Analysis of Relaxation Time of Temperature in Thermal Response Test for Design of Borehole Size

: A new practical method for thermal response test (TRT) is proposed herein to estimate the groundwater velocity and e ﬀ ective thermal conductivity of geological zones. The relaxation time of temperature (RTT) is applied to determine the depths of the zones. The RTT is the moment when the temperature in the borehole recovers to a certain level compared with that when the heating is stopped. The heat exchange rates of the zones are calculated from the vertical temperature proﬁle measured by the optical-ﬁber distributed temperature sensors located in the supply and return sides of a U-tube. Finally, the temperature increments at the end time of the TRT are calculated according to the groundwater velocities and the e ﬀ ective thermal conductivity using the moving line source theory applied to the calculated heat exchange rates. These results are compared with the average temperature increment data measured from each zone, and the best-ﬁtting value yields the groundwater velocities for each zone. Results show that the groundwater velocities for each zone are 2750, 58, and 0 m / y, whereas the e ﬀ ective thermal conductivities are 2.4, 2.4, and 2.1 W / (m · K), respectively. The proposed methodology is evaluated by comparing it with the realistic long-term operation data of a ground source heat pump (GSHP) system in Kazuno City, Japan. The temperature error between the calculated results and measured data is 6.4% for two years. Therefore, the proposed methodology is e ﬀ ective for estimating the long-term performance analysis of GSHP systems.


Introduction
Ground source heat pump (GSHP) systems utilizing the ground as a heat source or heat sink have been recognized as being high performance and environmentally friendly [1][2][3]. The studies on the GSHP system with various types of heat exchanger have studied how to supply heating and cooling to houses and the buildings [4,5]. Among them, borehole heat exchangers (BHEs), composed of vertical closed-loops, are preferred because of their high efficiency and minimal installation space requirement in expensive neighborhoods [6]. The reason for the high efficiency of the vertical loop type is that they can use the stable heat source from the relatively deeper ground, compared with the horizontal loop type which uses the heat source of the ground surface affected by the weather condition. The design of BHEs primarily depends on the thermal properties of soil as the design parameter. The main The methods to determine the thermal properties in a multilayer using the optical fiber distributed temperature sensing (DTS) technique have been proposed [24][25][26][27][28][29][30]. Fujii et al. [27] utilized optical fiber DTS to measure the temperature-time profile of circulating fluids in U-tubes during a TRT. The measured profiles were history matched with the cylindrical source function [8]. McDaniel et al. [30] provided a detailed description of variability in the subsurface heat transfer using the optical fiber DTS with laboratory thermophysical measurements. Sakata et al. [24] proposed a multilayer-concept TRT using optical fiber DTS and determined the stepwise ground thermal conductivity based on the depth of each sublayer. Their method was simple and practical use in TRT analysis without numerical interpretations. However, their method using conventional TRT analysis could not incorporate the effects of groundwater. It is important to understand the effective thermal conductivity and groundwater flow speed of each geological layer, as they will be used as the design parameters of the GSHP system. Therefore, a practical method for the thermal response test (TRT) is proposed herein to estimate the groundwater velocity and effective thermal conductivity of each geological layer. For simplicity, the layers were categorized into vertical zones within the depth of a borehole. Theoretical methodologies and experimental measurements are required to estimate the groundwater flow velocity and effective thermal conductivity.
First, the temperatures in the U-tube were measured according to the depth at an interval of 0.5 m from the optical fiber DTS. Furthermore, the temperature was monitored during both the heat injection and the stoppage of heat supply. Subsequently, a relaxation time of temperature (RTT) was applied to determine the depths of the zones. The RTT is the moment when the temperature in the borehole recovers to a certain level compared with that when the heating is stopped, including the effect of groundwater flow. Furthermore, the heat exchange rate of the zones was needed, which can be calculated from the vertical temperature profile of the circulating fluid during heat supply. Finally, the temperature increments of the circulating fluid were calculated according to the groundwater velocities using the MLS theory based on the calculated heat exchange rates. These results were then compared with the measured data from each zone and the best-fitting value yielded the groundwater velocities. The proposed methodology was evaluated by comparing it to the realistic long-term operation data of a GSHP system in Kazuno City, Japan. Figure 1 shows the ground plan and geological column section of the test site. The test site was located in Kazuno City (40 • 19 N and 140 • 78 E), Japan. The soils of the test site were investigated from the geological column section (Figure 1b); they were primarily composed of gravel, gravelly sand, and sandy gravel. The effective thermal conductivity of the soils should be 1-3 W/m·k) based on the geometric column section [31][32][33][34]. The thickness-weighted average value of the effective thermal conductivity was estimated to be 2.4 W/(m·k) when the soils were saturated in water. The air-conditioned space using the GSHP system was 143 m 2 in the three-floor building. The GSHP system was composed of an inverter-driven heat pump unit, and its maximum/minimum capacities for heating and cooling were 31.5/7.4 and 28.0/7.0 kW, respectively. Four BHEs were connected to the heat pump unit. The depth of each borehole was 100 m, and the borehole diameter was 144 mm. The inside and outside diameters of the U-tube made of the high-density poly ethene (HDPE) were 32 and 25 mm, respectively. A double U-tube without spacers was installed in each borehole. The borehole backfill material was dry silica sand, and it was inserted in the borehole from the ground surface. The shank spacing between the pipe center and the borehole center was estimated to be 0.045 m. The thermal conductivity of the pipe was 0.45 W/(m·K). The water table level in the test site was −7 m from the ground surface.  Figure 2 shows the schematics of the TRT. To determine the effect of groundwater flow, one long-term continuous TRT was performed from 12 to 29 January 2017. During the TRT, the temperature increase and the temperature recovery were monitored; the heat injection period was 198 h from 12 to 20 January and the relaxation period (after stopping the heating) was 202 h from 20 to 29 January. Two optical fiber-distributed temperature sensors were inserted in the supply and return sides of the U-tube. These sensors measured the temperature of the circulating fluid from the inlet/outlet side to the bottom side of the U-tube at intervals of 0.5 m. The measured temperatures in the U-tube were used to calculate the heat exchange rate of the heat injection period. Furthermore, they provided the temperature behavior during both the heat supply and the stoppage of the heat supply. The average initial temperature of the borehole was 12.0 °C before the TRT was conducted in a No. 3 BHE, as shown in Figure 1.   Figure 2 shows the schematics of the TRT. To determine the effect of groundwater flow, one long-term continuous TRT was performed from 12 to 29 January 2017. During the TRT, the temperature increase and the temperature recovery were monitored; the heat injection period was 198 h from 12 to 20 January and the relaxation period (after stopping the heating) was 202 h from 20 to 29 January. Two optical fiber-distributed temperature sensors were inserted in the supply and return sides of the U-tube. These sensors measured the temperature of the circulating fluid from the inlet/outlet side to the bottom side of the U-tube at intervals of 0.5 m. The measured temperatures in the U-tube were used to calculate the heat exchange rate of the heat injection period. Furthermore, they provided the temperature behavior during both the heat supply and the stoppage of the heat supply. The average initial temperature of the borehole was 12.0 • C before the TRT was conducted in a No. 3 BHE, as shown in Figure 1.  Figure 2 shows the schematics of the TRT. To determine the effect of groundwater flow, one long-term continuous TRT was performed from 12 to 29 January 2017. During the TRT, the temperature increase and the temperature recovery were monitored; the heat injection period was 198 h from 12 to 20 January and the relaxation period (after stopping the heating) was 202 h from 20 to 29 January. Two optical fiber-distributed temperature sensors were inserted in the supply and return sides of the U-tube. These sensors measured the temperature of the circulating fluid from the inlet/outlet side to the bottom side of the U-tube at intervals of 0.5 m. The measured temperatures in the U-tube were used to calculate the heat exchange rate of the heat injection period. Furthermore, they provided the temperature behavior during both the heat supply and the stoppage of the heat supply. The average initial temperature of the borehole was 12.0 °C before the TRT was conducted in a No. 3 BHE, as shown in Figure 1.

MLS
The MLS model [8] was used to analyze the TRT data in the area with the groundwater flow. The alternative flow speed, U, was introduced from the relationship between the Darcy velocity, v, and the moving medium of the line source theory [35] as shown in Equation (1). The time-dependent temperature increase in the polar coordinates (r, ϕ) can be obtained using the heat flux and alternative flow speed, as shown in Equation (2): The average temperature increase of the circulating fluid (∆T f ) was calculated using the borehole thermal resistance and the temperature increase (∆T) when r was assigned to the borehole radius (r bh ) in Equation (3). Meanwhile, ∆T f without the effect of U, can be approximately expressed as a linear equation with a logarithmic time when αt r 2 > 5 is satisfied in Equation (4). Subsequently, λ eff can be estimated as shown in Equation (5): where Ei is the exponential integral function, and k is slope of the temperature. The temporal superposition principle can provide the temperature response of the assumed borehole by considering the time-dependent heat flux [36][37][38][39]. When this method is utilized with the MLS model, the average temperature changes of the circulating fluid resembling the TRT temperature results are reproduced, which vary by in response to heat injection. The MLS model applied to the temporal superposition can express the average fluid temperature in accordance with the time-varying heat flux shown in Equation (6): The TRT was carried out at the same test site and borehole of the previous research [23]. R bh in Equation (3) was applied to the results of R bh in the previous research [23]. This R bh was calculated by the consideration of the groundwater flow in the borehole backfilled with the permeable material. The values of the R bh were reduced when the rapid groundwater passed through the inside borehole.

Estimation of Groundwater Velocity and Effective Thermal Conductivity in the Multilayer
RTT (t r ) is defined as the moment when the temperature in the borehole recovers to a certain level, where ∆T f , at a time when the heat supply is stopped (t ), is reduced by the ratio of the temperature decrement (ω): ∆T f (t r ) = ω∆T f (t ). Here, ω ranges between 0 and 1 (0 < ω < 1), whereas t r primarily depends on t , λ eff , and v, which affects the heat transfer intensity around the borehole Equation (8). Figure 3 illustrates the relationship between the dimensionless temperature variation (∆T * f (t)) and Energies 2020, 13, 3297 6 of 20 dimensionless RTT (t r /t ) according to the groundwater velocity when ω is 0.5. The expression for ∆T * f (t) is shown in Equation (8), based on Equation (6). When the groundwater flow speed was high, the RTT decreased: Energies 2020, 13, x FOR PEER REVIEW 6 of 20 borehole Equation (8). Figure 3 illustrates the relationship between the dimensionless temperature variation (∆ f * ( )) and dimensionless RTT ( / ′) according to the groundwater velocity when is 0.5. The expression for ∆ f * ( ) is shown in Equation (8), based on Equation (6). When the groundwater flow speed was high, the RTT decreased: The layers in this research were decided by the measured position of the optical-fiber DTS, and the length of a layer was 0.5 m ( Figure 2). For the simplicity of a complex geological structure, the layers were categorized into vertical zones within the depth of a borehole according to the effects of the groundwater flow. The vertical distribution of the RTT ( , ) was calculated from the measured recovery temperature of each layer. To classify the zones, the boundaries of the RTT ( , ) were determined by the groundwater velocity. The layers were assigned the zone of th when RTT ( , ) was shorter than RTT ( , ). Here, is the number of zones. The layers with the rapid groundwater flow were sequentially classified from starting with the first Zone. The boundary of RTT ( , ) was applied until the sum of the length of zones equaled the depth of the borehole.
The heat exchange rate was calculated by the temperature changes of the circulating fluid in the layer using Equation (9). These temperature changes in the layer were observed from the optical fiber DTS during the heat injection period. The heat exchange rate of each zone was calculated using the average heat exchange rate of the grouped layers: Finally, eff ang of each zone were estimated thought the comparison between the calculated results (∆ f,cal ̅̅̅̅̅̅ ( ′ , eff, , , Z N )) and the measured temperature increment from the optical fiber DTS in each zone (∆ f ̅̅̅̅̅̅ ( ′ )). Figure 4 shows a flow chart of the calculation performed in this method. The layers in this research were decided by the measured position of the optical-fiber DTS, and the length of a layer was 0.5 m ( Figure 2). For the simplicity of a complex geological structure, the layers were categorized into vertical zones within the depth of a borehole according to the effects of the groundwater flow. The vertical distribution of the RTT (t r,L i ) was calculated from the measured recovery temperature of each layer. To classify the zones, the boundaries of the RTT (t r,b N ) were determined by the groundwater velocity. The layers were assigned the zone of Nth when RTT (t r,L i ) was shorter than RTT (t r,b N ). Here, N is the number of zones. The layers with the rapid groundwater flow were sequentially classified from starting with the first Zone. The boundary of RTT (t r,b N ) was applied until the sum of the length of zones equaled the depth of the borehole.
The heat exchange rate was calculated by the temperature changes of the circulating fluid in the layer using Equation (9). These temperature changes in the layer were observed from the optical fiber DTS during the heat injection period. The heat exchange rate of each zone was calculated using the average heat exchange rate of the grouped layers: Finally, λ eff ang v of each zone were estimated thought the comparison between the calculated results (∆T f,cal t , λ eff,N , v N , q Z N ) and the measured temperature increment from the optical fiber DTS in each zone (∆T fzN (t )). Figure 4 shows a flow chart of the calculation performed in this method.  Figure 5 shows the measurement data of the TRT. The temperature difference between the inlet, in and outlet, out was almost constant, i.e., 4.6 °C, as shown in Figure 5a. The flow rate and heat injection were approximately 20 L/min and 6.06 kW during the heat injection period, respectively, as shown in Figure 3b. From the conventional TRT analysis, the temperature gradient ′ was 0.51 °C/ln(t), and the apparent effective thermal conductivity was 9.5 W/(m · K) for 12-60 h [40].   Figure 5 shows the measurement data of the TRT. The temperature difference between the inlet, T in and outlet, T out was almost constant, i.e., 4.6 • C, as shown in Figure 5a. The flow rate and heat injection were approximately 20 L/min and 6.06 kW during the heat injection period, respectively, as shown in Figure 3b. From the conventional TRT analysis, the temperature gradient k was 0.51 • C/ln(t), and the apparent effective thermal conductivity was 9.5 W/(m·K) for 12-60 h [40].  Figure 5 shows the measurement data of the TRT. The temperature difference between the inlet, in and outlet, out was almost constant, i.e., 4.6 °C, as shown in Figure 5a. The flow rate and heat injection were approximately 20 L/min and 6.06 kW during the heat injection period, respectively, as shown in Figure 3b. From the conventional TRT analysis, the temperature gradient ′ was 0.51 °C/ln(t), and the apparent effective thermal conductivity was 9.5 W/(m · K) for 12-60 h [40].    Figure 6 indicates the dimensionless average temperature variation of the circulating fluid. The temperature increase during the heat injection was the measured results of the Pt-100 sensors and the optical fiber-distributed temperature sensors at the inlet and outlet pipes. The recovery temperature during the relaxation period was halted and was the average temperature of the optical fiber distributed temperature sensors in all the layers. ∆T * f was calculated from Equations (7) and (8) when λ eff was 3 W/(m·k) and v was 0 m/y. In addition, the other parameters were obtained based on the same TRT conditions. These results indicated the effect of groundwater flow on the temperature behaviors of the circulating fluid.

Determination of Zone Depths
Energies 2020, 13, x FOR PEER REVIEW 8 of 20 Figure 6 indicates the dimensionless average temperature variation of the circulating fluid. The temperature increase during the heat injection was the measured results of the Pt-100 sensors and the optical fiber-distributed temperature sensors at the inlet and outlet pipes. The recovery temperature during the relaxation period was halted and was the average temperature of the optical fiber distributed temperature sensors in all the layers. ∆ f * was calculated from Equations (7) and (8) when eff was 3 W/(m • k) and was 0 m/y. In addition, the other parameters were obtained based on the same TRT conditions. These results indicated the effect of groundwater flow on the temperature behaviors of the circulating fluid. The ground was categorized into three vertical zones: Zone 1 comprised the layers with strong groundwater flow effects; Zone 2 comprised the layers with intermediate groundwater flow effects; and Zone 3 comprised the layers with the weak groundwater flow effects. These zones were distinguished by the boundary of the RTT ( ,b ). Two ,b were discovered from the variation in according to . Figure 7 shows the RTT according to the groundwater velocity when was 0.5. Here, eff was considered from 1.5 to 3 W/(m • K) at an interval of 0.5, which is the expected range at the test site. ' was 198 h, which is the same TRT duration as during the heat supply. The ,b1 for classifying Zones 1 and 2 in this study was determined when the difference in the variation of was less than 10 −4 %. The ,b1 was 2.1 h when the groundwater velocity was 200 m/y. Meanwhile, the ,b2 for classifying Zones 2 and 3 was decided from the variation in . In the zone with the weak effect of groundwater flow, the variations in were small. Studies [13,19] have indicated that it was difficult to determine a less than 30 m/y through the TRT. In this study, ,b2 was determined when the variations in were less than 10%. In addition, the highest value of eff in the test site (3 W/(m • k)) was applied to determine ,b2 . Using this ,b2 , Zone 3 exhibited a eff lower than 3 W/(m • k). Consequently, the ,b2 was 10.7 h.  The ground was categorized into three vertical zones: Zone 1 comprised the layers with strong groundwater flow effects; Zone 2 comprised the layers with intermediate groundwater flow effects; and Zone 3 comprised the layers with the weak groundwater flow effects. These zones were distinguished by the boundary of the RTT (t r,b ). Two t r,b were discovered from the variation in t r according to v. Figure 7 shows the RTT according to the groundwater velocity when ω was 0.5. Here, λ eff was considered from 1.5 to 3 W/(m· K) at an interval of 0.5, which is the expected range at the test site. t was 198 h, which is the same TRT duration as during the heat supply. The t r,b1 for classifying Zones 1 and 2 in this study was determined when the difference in the variation of t r was less than 10 −4 %. The t r,b1 was 2.1 h when the groundwater velocity was 200 m/y. Meanwhile, the t r,b2 for classifying Zones 2 and 3 was decided from the variation in t r . In the zone with the weak effect of groundwater flow, the variations in t r were small. Studies [13,19] have indicated that it was difficult to determine a v less than 30 m/y through the TRT. In this study, t r,b2 was determined when the variations in t r were less than 10%. In addition, the highest value of λ eff in the test site (3 W/(m·k)) was applied to determine t r,b2 . Using this t r,b2 , Zone 3 exhibited a λ eff lower than 3 W/(m·k). Consequently, the t r,b2 was 10.7 h. Figure 8 shows the vertical distribution of t r , which was calculated from Equations (7) and (8) based on the temperature recovery in the borehole of each layer. First, Zone 1 with significant effects of groundwater flow was decided when t r,L i was shorter than t r,b1 . Subsequently, Zone 3 with weak effects of groundwater flow was determined when t r,L i was higher than t r,b2 . Finally, the undecided layers were allocated to the intermediate zone (Zone 2). Hence, the depths of Zones 1, 2 and 3 were 40, 28 and 32 m, respectively. Table 1 lists the conditions of the zones.

Determination of Zone Depths
Energies 2020, 13, 3297 9 of 20 the ,b2 for classifying Zones 2 and 3 was decided from the variation in . In the zone with the weak effect of groundwater flow, the variations in were small. Studies [13,19] have indicated that it was difficult to determine a less than 30 m/y through the TRT. In this study, ,b2 was determined when the variations in were less than 10%. In addition, the highest value of eff in the test site (3 W/(m • k)) was applied to determine ,b2 . Using this ,b2 , Zone 3 exhibited a eff lower than 3 W/(m • k). Consequently, the ,b2 was 10.7 h.  Energies 2020, 13, x FOR PEER REVIEW 9 of 20 Figure 8 shows the vertical distribution of , which was calculated from Equations (7) and (8) based on the temperature recovery in the borehole of each layer. First, Zone 1 with significant effects of groundwater flow was decided when , was shorter than ,b1 . Subsequently, Zone 3 with weak effects of groundwater flow was determined when , was higher than ,b2 . Finally, the undecided layers were allocated to the intermediate zone (Zone 2). Hence, the depths of Zones 1, 2 and 3 were 40, 28 and 32 m, respectively. Table 1 lists the conditions of the zones.  Based on the standard procedure of TRTs in Japan [40], the data measured from 12 to 60 h were used to analyze the TRT results. Figure 9 shows the vertical temperature distribution during the heat injection period and the heat transfer rate of each zone. The heat exchange rate was calculated from Equation (9). As shown in Figure 9b, the heat exchange change rate in Zone 1 decreased gradually. This was caused by the temperature difference between the ground and the circulating fluid. The ground temperature was stabilized by the effect of the groundwater flow, whereas the temperature of the circulating fluid increased at the beginning of the heat supply. Meanwhile, the temperature increase rate of the circulating fluid decreased gradually as the TRT progressed, causing the temperature difference between the ground and the circulating fluid to stabilize. Hence, the average heat exchange rates of Zones 1, 2 and 3 were calculated to be 92.8, 31.9, and 27.6, respectively.
Based on the standard procedure of TRTs in Japan [40], the data measured from 12 to 60 h were used to analyze the TRT results. Figure 9 shows the vertical temperature distribution during the heat injection period and the heat transfer rate of each zone. The heat exchange rate was calculated from Equation (9). As shown in Figure 9b, the heat exchange change rate in Zone 1 decreased gradually. This was caused by the temperature difference between the ground and the circulating fluid. The ground temperature was stabilized by the effect of the groundwater flow, whereas the temperature of the circulating fluid increased at the beginning of the heat supply. Meanwhile, the temperature increase rate of the circulating fluid decreased gradually as the TRT progressed, causing the temperature difference between the ground and the circulating fluid to stabilize. Hence, the average heat exchange rates of Zones 1, 2 and 3 were calculated to be 92.8, 31.9, and 27.6, respectively.  Figure 10 shows the temperature increment according to the groundwater velocity and effective thermal conductivity in each zone. The blue line represents the ∆ f ( ′ ) calculated based on the heat exchanger rates in each zone. Here, the effective thermal conductivity in the advection-influenced zones was considered as the thickness-weighted average value estimated from the geological column section. In advection-influenced conditions, ∆ f is not significantly affected by the effective thermal conductivity but depends primarily on the effect of the groundwater flow ( Figure 7). Meanwhile, in the zone with a weak effect of groundwater flow, the effective thermal conductivity primarily affected ∆ f . In this study, eff of 1-3 W/(m • k)) was applied to the calculation of ∆ f (Figure 10c). The red line represents the measured ∆ f ( ′ ). It is the average temperature of the circulating fluid in each zone and was obtained from the optical fiber DTS, and their values were 6.38 °C, 6.35 °C, and 6.32 °C for Zones 1, 2 and 3, respectively. The groundwater velocity was determined at a point where the red line intersected with the blue lines. This intersection of two lines implies that the ∆ f,cal ( ′ ) applied to the specified groundwater velocity with the other parameters was the ∆ f,data ( ′ ) in Figure 10a,b. In Figure 10c, the eff was estimated similarly. These estimated parameters in each zone can represent the soil parameters. The groundwater velocity obtained was 2750 m/y, 58 m/y, and 0 m/y, whereas the thermal conductivity was 2.4 W/(m • k), 2.4 W/(m • k), and 2.1 W/(m • k) for Zones 1, 2 and 3, respectively.  Figure 10 shows the temperature increment according to the groundwater velocity and effective thermal conductivity in each zone. The blue line represents the ∆T f (t ) calculated based on the heat exchanger rates in each zone. Here, the effective thermal conductivity in the advection-influenced zones was considered as the thickness-weighted average value estimated from the geological column section. In advection-influenced conditions, ∆T f is not significantly affected by the effective thermal conductivity but depends primarily on the effect of the groundwater flow ( Figure 7). Meanwhile, in the zone with a weak effect of groundwater flow, the effective thermal conductivity primarily affected ∆T f . In this study, λ eff of 1-3 W/(m·k)) was applied to the calculation of ∆T f (Figure 10c). The red line represents the measured ∆T f (t ). It is the average temperature of the circulating fluid in each zone and was obtained from the optical fiber DTS, and their values were 6.38 • C, 6.35 • C, and 6.32 • C for Zones 1, 2 and 3, respectively. The groundwater velocity was determined at a point where the red line intersected with the blue lines. This intersection of two lines implies that the ∆T f,cal (t ) applied to the specified groundwater velocity with the other parameters was the ∆T f,data (t ) in Figure 10a,b. In Figure 10c, the λ eff was estimated similarly. These estimated parameters in each zone can represent the soil parameters. The groundwater velocity obtained was 2750 m/y, 58 m/y, and 0 m/y, whereas the thermal conductivity was 2.4 W/(m·k), 2.4 W/(m·k), and 2.1 W/(m·k) for Zones 1, 2 and 3, respectively.

Comparison between Calculated Results and TRT Data
To validate the proposed method, the calculated temperatures were compared with the TRT data. These results were calculated by using the MLS model with the estimated parameters in Equations (10)- (12). Figure 11 shows the temperature change of the calculated results and the TRT data. The calculated results matched well with the TRT data for 198. The comparison results validated the estimated design parameters:

Comparison between Calculated Results and TRT Data
To validate the proposed method, the calculated temperatures were compared with the TRT data. These results were calculated by using the MLS model with the estimated parameters in Equations (10)- (12). Figure 11 shows the temperature change of the calculated results and the TRT data. The calculated results matched well with the TRT data for 198. The comparison results validated the estimated design parameters: Energies 2020, 13, 3297 12 of 20 Energies 2020, 13, x FOR PEER REVIEW 12 of 20 Figure 11. Comparison between the calculated results and the TRT data.

Analysis Results of Thermal Parameters According to TRT End Time
The TRT is generally operated within a limited time owing to the required construction time and cost. For this reason, the TRT data for 60 h [40], 48 h [21], and 96 h [41] was used to estimate the eff and after the determination of the zone. The estimated parameters from the proposed method were validated by comparing the TRT data for 200 h. Figure 12 shows the temperature change according to the design parameters determined by the end time of the TRT, which was referred from previous studies: 60 h [40] (Case 1), 48 h [21] (Case 2), and 96 h [41] (Case 3). Table 2 indicates the results for each case. The temperature increments were calculated based on the heat exchange rates of each sublayer. These heat exchange rates were determined by the measured temperature, which might contain measurement errors or be affected by the surrounding environment. Nevertheless, the heat exchange rate converged gradually as the TRT progressed. Based on the estimated thermal parameter at each end time, the temperature results of Cases 1 and 3 were similar but differed slightly from those of Case 2, compared with the measurement data from the DTS. For this reason, the heat exchange rate converged over time. Although a longer TRT time can provide more accurate design parameters, increased construction time and effort are required. Because the optimal TRT time depends on the test location, future studies will be performed to determine the optimal TRT time based on the TRT results of other sites.

Analysis Results of Thermal Parameters According to TRT End Time
The TRT is generally operated within a limited time owing to the required construction time and cost. For this reason, the TRT data for 60 h [40], 48 h [21], and 96 h [41] was used to estimate the λ eff and v after the determination of the zone. The estimated parameters from the proposed method were validated by comparing the TRT data for 200 h. Figure 12 shows the temperature change according to the design parameters determined by the end time of the TRT, which was referred from previous studies: 60 h [40] (Case 1), 48 h [21] (Case 2), and 96 h [41] (Case 3). Table 2 indicates the results for each case. The temperature increments were calculated based on the heat exchange rates of each sublayer. These heat exchange rates were determined by the measured temperature, which might contain measurement errors or be affected by the surrounding environment. Nevertheless, the heat exchange rate converged gradually as the TRT progressed. Based on the estimated thermal parameter at each end time, the temperature results of Cases 1 and 3 were similar but differed slightly from those of Case 2, compared with the measurement data from the DTS. For this reason, the heat exchange rate converged over time. Although a longer TRT time can provide more accurate design parameters, increased construction time and effort are required. Because the optimal TRT time depends on the test location, future studies will be performed to determine the optimal TRT time based on the TRT results of other sites.
Energies 2020, 13, x FOR PEER REVIEW 12 of 20 Figure 11. Comparison between the calculated results and the TRT data.

Analysis Results of Thermal Parameters According to TRT End Time
The TRT is generally operated within a limited time owing to the required construction time and cost. For this reason, the TRT data for 60 h [40], 48 h [21], and 96 h [41] was used to estimate the eff and after the determination of the zone. The estimated parameters from the proposed method were validated by comparing the TRT data for 200 h. Figure 12 shows the temperature change according to the design parameters determined by the end time of the TRT, which was referred from previous studies: 60 h [40] (Case 1), 48 h [21] (Case 2), and 96 h [41] (Case 3). Table 2 indicates the results for each case. The temperature increments were calculated based on the heat exchange rates of each sublayer. These heat exchange rates were determined by the measured temperature, which might contain measurement errors or be affected by the surrounding environment. Nevertheless, the heat exchange rate converged gradually as the TRT progressed. Based on the estimated thermal parameter at each end time, the temperature results of Cases 1 and 3 were similar but differed slightly from those of Case 2, compared with the measurement data from the DTS. For this reason, the heat exchange rate converged over time. Although a longer TRT time can provide more accurate design parameters, increased construction time and effort are required. Because the optimal TRT time depends on the test location, future studies will be performed to determine the optimal TRT time based on the TRT results of other sites.

Long-Term Simulation Results According to Time-Variable Building Load
The estimated design parameters were validated through a comparison between the measured data and the calculated results for the long-term period. The measured data above were those of the circulating fluid when the GSHP system was operated with four BHEs (the GSHP system and the target building are described in Section 2). The measured period was from 1 November 2017 to 22 October 2019. The temperature change in the circulating fluid is calculated by using the GSHP simulation tool applied to the estimated parameters and the time-variable building loads as the input parameters. The GSHP simulation was "Ground Club," developed by Hokkaido University [42]. The simulation tool can incorporate the effects of groundwater flow as well as the multiple BHEs in multilayers [43][44][45]. Figure 13 shows the heating and cooling load of the target building. These building loads were calculated using the electric power consumption of the heat pump and the heat extraction and injection of the BHEs based on Equations (13) and (14). The annual heating and cooling loads were 106.6 and 20.1 GJ from 2017 to 2018, respectively, and 108.7 and 18.1 GJ from 2018 to 2019, respectively. The maximum, minimum, and average outdoor temperatures of the test site (Figure 13b) were 35.6 • C, −18.4 • C, and 9.6 • C, respectively: Energies 2020, 13, x FOR PEER REVIEW 13 of 20 (m/y) 0 0 0

Long-Term Simulation Results According to Time-Variable Building Load
The estimated design parameters were validated through a comparison between the measured data and the calculated results for the long-term period. The measured data above were those of the circulating fluid when the GSHP system was operated with four BHEs (the GSHP system and the target building are described in Section 2). The measured period was from 1 November 2017 to 22 October 2019. The temperature change in the circulating fluid is calculated by using the GSHP simulation tool applied to the estimated parameters and the time-variable building loads as the input parameters. The GSHP simulation was "Ground Club," developed by Hokkaido University [42]. The simulation tool can incorporate the effects of groundwater flow as well as the multiple BHEs in multilayers [43][44][45]. Figure 13 shows the heating and cooling load of the target building. These building loads were calculated using the electric power consumption of the heat pump and the heat extraction and injection of the BHEs based on Equations (13) and (14). The annual heating and cooling loads were 106.6 and 20.1 GJ from 2017 to 2018, respectively, and 108.7 and 18.1 GJ from 2018 to 2019, respectively. The maximum, minimum, and average outdoor temperatures of the test site ( Figure  13b Figure 14 shows the average temperatures of the circulating fluid. These temperatures were those measured from the inlet and outlet of the heat pump and the calculated results of the simulation tool. The annual temperature error between the calculated results and the measured data was 6.4% for two years. These comparison results demonstrated the appropriateness of the estimated design parameters for designing the GSHP systems. Table 3 presents the uncertainties of the measured and evaluated parameter in the experimental estimation of the heat transfer rate.
Energies 2020, 13, x FOR PEER REVIEW 14 of 20 Figure 14 shows the average temperatures of the circulating fluid. These temperatures were those measured from the inlet and outlet of the heat pump and the calculated results of the simulation tool. The annual temperature error between the calculated results and the measured data was 6.4% for two years. These comparison results demonstrated the appropriateness of the estimated design parameters for designing the GSHP systems. Table 3 presents the uncertainties of the measured and evaluated parameter in the experimental estimation of the heat transfer rate.

Estimation of the Circulating Fluid According to the Borehole Size
The proposed method has the advantage of the consideration for the optimal BHE size in the areas where the groundwater flows rapidly in the specific layers. The conventional TRT analysis methods or the previous research [23] can provide only the weight-average values of the effective thermal conductivity and the groundwater velocity with respect to the depth of the BHE where the TRT was carried out. If the BHE size was changed, the estimated values from the TRT could not use it for the design of the BHE. However, the proposed methodology could provide the effective thermal conductivities and the groundwater velocities in multi-layer. In particular, it is possible to effectively design a BHE size by grasping the thermal properties of the zones and the groundwater velocities. The following paragraph presents the comparison results between the conventional TRT analysis method and the proposed method. Table 4 shows the estimated parameter values according to the TRT analysis methods.

Estimation of the Circulating Fluid According to the Borehole Size
The proposed method has the advantage of the consideration for the optimal BHE size in the areas where the groundwater flows rapidly in the specific layers. The conventional TRT analysis methods or the previous research [23] can provide only the weight-average values of the effective thermal conductivity and the groundwater velocity with respect to the depth of the BHE where the TRT was carried out. If the BHE size was changed, the estimated values from the TRT could not use it for the design of the BHE. However, the proposed methodology could provide the effective thermal conductivities and the groundwater velocities in multi-layer. In particular, it is possible to effectively design a BHE size by grasping the thermal properties of the zones and the groundwater velocities. The following paragraph presents the comparison results between the conventional TRT analysis method and the proposed method. Table 4 shows the estimated parameter values according to the TRT analysis methods. Figure 15 shows that the average temperature of the circulating fluid during the heating and cooling seasons according to the borehole size. T f in the condition of Case 1 decreased by −3.14 • C from T 0 during the heating season. Meanwhile, T f in the condition of Case 2 decreased by −4.42 • C from T 0 during the heating season. This difference shows the importance considering the groundwater flow effects in the TRT analysis method. If the BHEs are designed based on the conventional TRT analysis result (Case 2), more BHEs are required to supply the appropriate temperature of the fluid entering the heat pump to operate the GSHP system. Meanwhile, T f in the condition of Case 1 decreased by −4.61 • C from T 0 despite the BHEs of 100 m being reduced to that of 40 m. This result was obtained because T f was calculated using the finite cylinder source model and the parameters of Zone 1 in the condition of Case 1 was considered. Zone 1 in the condition of Case 1 was located from the ground surface to −40 m below ground and the v of Zone 1 was 2750 m/y. By contrast, ∆T in the condition of Case 2 decreased by −9.36 • C when the length of the BHEs was reduced by 40% from 100 m. Table 4. Estimated parameter values according to the TRT analysis methods.
Case    Figure 16 shows the average temperature variations of the circulating fluid. These temperatures were calculated for two years based on the conditions of Cases 1-1 and 2-1. The heating and cooling loads were of the same condition as that of the target building described in Section 5.3. The temperature variation in the condition of Case 2-1 exhibited larger fluctuations of −2.7 °C to 1.9 °C compared with that in the condition of Case 1-1.  Figure 16 shows the average temperature variations of the circulating fluid. These temperatures were calculated for two years based on the conditions of Cases 1-1 and 2-1. The heating and cooling loads were of the same condition as that of the target building described in Section 5.3. The temperature variation in the condition of Case 2-1 exhibited larger fluctuations of −2.7 • C to 1.9 • C compared with that in the condition of Case 1-1. Figure 17 shows the average temperature variations of the circulating fluid. The condition of Case 3-1 was obtained from a previous study [23], which provided the values of λ eff and v in a single layer of a 100 m BHE at the same test site. The result shows that the temperature error between Cases 1-1 and 3-1 was 0.3%. Both methods were regarded as suitable for the design of BHEs in areas with rapid groundwater flows. Meanwhile, Figure 18 shows the average temperature variations of the circulating fluid when the borehole size was 40 m. The temperature variation in the condition of Case 3-3 exhibited larger fluctuations of −3.0 • C to 1.8 • C compared with that in the condition of Case 1-3. The seasonal temperature error in the heating season was 9.2%. The proposed method obtained layers with strong groundwater flow effects and provided the temperature variation prediction of the circulating fluid according to the BHE size. It was effective in designing the appropriate BHE size.  Figure 17 shows the average temperature variations of the circulating fluid. The condition of Case 3-1 was obtained from a previous study [23], which provided the values of eff and in a single layer of a 100 m BHE at the same test site. The result shows that the temperature error between Cases 1-1 and 3-1 was 0.3%. Both methods were regarded as suitable for the design of BHEs in areas with rapid groundwater flows. Meanwhile, Figure 18 shows the average temperature variations of the circulating fluid when the borehole size was 40 m. The temperature variation in the condition of Case 3-3 exhibited larger fluctuations of −3.0 °C to 1.8 °C compared with that in the condition of Case 1-3. The seasonal temperature error in the heating season was 9.2%. The proposed method obtained layers with strong groundwater flow effects and provided the temperature variation prediction of the circulating fluid according to the BHE size. It was effective in designing the appropriate BHE size.    Figure 17 shows the average temperature variations of the circulating fluid. The condition of Case 3-1 was obtained from a previous study [23], which provided the values of eff and in a single layer of a 100 m BHE at the same test site. The result shows that the temperature error between Cases 1-1 and 3-1 was 0.3%. Both methods were regarded as suitable for the design of BHEs in areas with rapid groundwater flows. Meanwhile, Figure 18 shows the average temperature variations of the circulating fluid when the borehole size was 40 m. The temperature variation in the condition of Case 3-3 exhibited larger fluctuations of −3.0 °C to 1.8 °C compared with that in the condition of Case 1-3. The seasonal temperature error in the heating season was 9.2%. The proposed method obtained layers with strong groundwater flow effects and provided the temperature variation prediction of the circulating fluid according to the BHE size. It was effective in designing the appropriate BHE size.    Figure 17 shows the average temperature variations of the circulating fluid. The condition of Case 3-1 was obtained from a previous study [23], which provided the values of eff and in a single layer of a 100 m BHE at the same test site. The result shows that the temperature error between Cases 1-1 and 3-1 was 0.3%. Both methods were regarded as suitable for the design of BHEs in areas with rapid groundwater flows. Meanwhile, Figure 18 shows the average temperature variations of the circulating fluid when the borehole size was 40 m. The temperature variation in the condition of Case 3-3 exhibited larger fluctuations of −3.0 °C to 1.8 °C compared with that in the condition of Case 1-3. The seasonal temperature error in the heating season was 9.2%. The proposed method obtained layers with strong groundwater flow effects and provided the temperature variation prediction of the circulating fluid according to the BHE size. It was effective in designing the appropriate BHE size.

Conclusions
A TRT analytical method was proposed herein to estimate the groundwater velocity and the effective thermal conductivity of geological zones. Furthermore, the RTT was applied to determine the zone depths by considering the temperature recovery during the relaxation period (Sections 3.2 and 4.2). Temperature increments at the TRT end time were calculated according to the groundwater velocities using the MLS model (Sections 3.2 and 4.3). These results were compared with the average temperature increments measured from each zone and its best-fitting value yielded the groundwater velocities.
The groundwater velocity was estimated to be 2750, 58, and 0 m/y, whereas the effective thermal conductivity was evaluated to be 2.4, 2.4, and 2.1 W/(m·k) for Zones 1, 2 and 3, respectively (Section 4.3).
These estimated velocities and other parameters were applied to the simulation tool to calculate the temperatures of the circulating fluid. The calculated temperatures were validated by comparing them to the measurement data of the circulating fluid in the target building. The annual temperature error between the calculated results and the measured data was 6.4% (Section 5.3). In addition, the seasonal average temperatures of the circulating fluid were estimated according to the BHE size by applying them to the parameters estimated using the TRT analysis methods (Section 6). It was discovered that the proposed method was effective in obtaining layers that were significantly affected by the groundwater flow and demonstrated the appropriate BHE size. The long-term performance can be predicted well using the groundwater velocity and the effective thermal conductivity for each zone obtained from the proposed method as the calculating condition.