Effect of Borehole Material on Analytical Solutions of the Heat Transfer Model of Ground Heat Exchangers Considering Groundwater Flow

Abstract: Groundwater flow is one of the most important factors for the design of a ground heat exchanger (GHEX) since the thermal environment of the ground around the buried GHEX is significantly affected by heat convection due to the groundwater flow. Several preceding studies have been conducted to develop analytical solutions to the heat transfer model of GHEX with consideration of groundwater flow. One of these solutions is the combined heat transfer model of conduction and convection. However, the developed combined analytical models are inapplicable to all of the configurations of ordinary GHEXs because these solutions assume that the inner part of the borehole is thermally inert or consists of the same material as that of the surrounding ground. In this paper, the applicability of the combined solid cylindrical heat source model, which is the most suitable to energy piles until now, was evaluated by performing a series of numerical analyses. In the numerical analysis, the inner part of the borehole was modeled as two different materials (i.e., permeable ground formation and impermeable fill such as concrete) to evaluate applicability of the analytical solution along with different diameter-length (D/L) ratios of borehole. In a small value of the D/L ratio, the analytical solution to the combined heat transfer model is in good agreement with the result of numerical analysis. On the other hand, when increasing the D/L ratio, the analytical solution significantly overestimates the effect of groundwater flow on the heat transfer of GHEXs because the analytical solution disregards the existence of the impermeable region in the borehole. Consequently, such tendency is more critical in the GHEX with a large D/L ratio such as large-diameter energy piles.


Introduction
The ground source heat pump (GSHP) system is one of the most effective heating and cooling systems because it utilizes the constant subsurface temperature regardless of season to improve the thermal performance of heat pump and reduce energy consumption [1].Compared with conventional heating, ventilating, and air conditioning (HVAC) systems, the GSHP system has been favorably evaluated due to its safety, economic feasibility, and effectiveness [2].In the GSHP system, heat exchange pipes are installed underground where heat exchange occurs between a working fluid and the subsurface ground.By extracting or releasing thermal energy from or into the ground, the GSHP system increases the temperature level of the working fluid [3].The buried heat exchange pipes, called the ground heat exchanger (GHEX), are then connected via a heat pump for heating and cooling the building [4].Therefore, the optimum design of the GHEX is crucial because the thermal performance of GHEX dominates the overall performance of the GSHP system.
Usually, the GHEX uses geothermal energy indirectly with the aid of a heat carrier that circulates in a closed system installed in a horizontal [5,6] or vertical configuration [7].The closed-loop vertical GHEX is the most popular heat exchanger in practice, and is commonly equipped with a single U-type or a double U-type heat exchange pipe in a borehole of 100 ~300 mm diameter and 50 ~200 m drilling depth.However, the conventional closed-loop vertical GHEX can involve a relatively high construction cost due to the requirement of additional borehole drilling and extra construction area [8][9][10].Thus, a novel type of GHEX, the energy pile, has recently attracted attention.The energy pile is constructed by adopting an existing pile foundation in order to reduce the initial construction cost, such as the drilling cost, and the requirement of additional space for construction [11][12][13].The energy pile has relatively shorter length (i.e., 20~50 m) and larger diameter (i.e., 350~2000 mm) than those of the conventional closed-loop vertical GHEX because it is designed for a structural purpose.
Among a variety of design methods for GHEXs, the analytical model of heat transfer (i.e., heat conduction) around a buried heat source has been widely used in studying field performance of the GHEX [14] due to its clarity and convenience for computation.The finite line heat source model [15] is a typical model which is available to simulate the thermal behavior around the line-shape of a heat source.Especially, the thermal properties of the surrounding medium such as thermal conductivity and thermal diffusivity can be obtained with simple computation using the line heat source model.As a modified analytical heat transfer model, the solid cylinder heat source model [16] is applicable to an energy pile that has a cylinder-shaped heat source with a large diameter and short length, assuming that the homogeneous thermo-physical properties of the internal heat source medium are the same as the surrounding medium.Additionally, Cui et al. [17] proposed the ring-coil heat source model that is suitable to an energy pile with coil-shape heat exchange pipes.The ring-coil heat source model assumes that the coil-shaped heat sources are located at the same interval in the borehole, and a uniform heat flux is supplied.
However, the current homogeneous models cannot consider the thermal resistance of the borehole, which means these models calculate the transient thermal responses of the surrounding medium of GHEXs by applying heat flux to the borehole wall, not to the working fluid or heat exchange pipes.In other words, these models assume that the heat transfer process inside a borehole is approximately steady-state.Therefore, short-term thermal behaviors estimated by the homogeneous models marginally deviates from the actual situation, and this phenomenon becomes more critical for a large-diameter energy pile.Lamarche [18] discussed limitations of the homogeneous models assuming that the mean fluid temperature is identical to the temperature in the borehole wall.In addition, Park et al. [19] evaluated the validity of the homogeneous models for a large-diameter cast-in-place energy pile by comparing the analytical solutions with field measurements of the short-term thermal response.
Recently, composite analytical models, considering the heat capacity inside the borehole, have been introduced as an alternative solution to the homogeneous models for the more accurate simulation of the short-term heat transfer process in GHEXs [20].Lei et al. [21] developed an advanced analytical solution considering a periodic cylinder heat source in composite media.However, these analytical models for heat conduction cannot take into account the existence of groundwater flow, even though the groundwater flow is one of the most important factors for designing ground heat exchangers (GHEXs) since the thermal environment of the ground around the buried GHEX is significantly affected by heat convection due to groundwater flow.When the influence of groundwater flow is obvious, the advection and conduction of water in a porous medium should be simultaneously considered along with the heat conduction between solid particles.Several studies have been conducted to enhance the previous analytical models for simulating the thermal behavior of GHEX with consideration of groundwater flow.Molina-Giraldo et al. [22] developed the moving finite line heat source model, and Zhang et al. [23] introduced the combined heat transfer model of conduction and convection toward the Energies 2016, 9, 318 3 of 19 line heat source, cylindrical heat source, and ring-coil heat source model.Especially, Zhang et al. [24] proposed the combined solid cylinder heat source model with consideration of groundwater flow by introducing the effective heat transfer coefficient to take into account the conduction of water and solid particles in a porous medium, and the moving heat source equation to implement advection exerted by groundwater flow.In addition, Tye-Gingras and Gosselin [25] proposed the generic ground response functions for GHEXs in the presence of groundwater flow by developing the combined analytical solutions.However, the developed combined analytical heat transfer models are inapplicable to all the configurations of ordinary GHEXs such as a large-diameter energy pile because these models assume the inner part of the borehole is thermally inert or consists of the same material as that of the ground.
In this paper, the applicability of the combined solid cylinder heat source model to a large-diameter energy pile was evaluated by performing numerical analysis on two different assumptions about the inner part of the heat source.First, the inner part of the heat source is assumed to be the same material as that of the surrounding ground formation, which is identical to the composite analytical model.Second, the inner part of the heat source is modeled as an impermeable borehole to evaluate applicability of the analytical solution.A series of numerical analyses was carried out with three different velocities of groundwater flow (0, 2.5 ˆ10 ´6 and 5.0 ˆ10 ´6 m/s) and two different diameter-length (D/L) ratios of heat source (0.025 and 0.00075).

Analytical Model for Pure Conduction
The concept behind fundamental analytical models used for heat transfer originates from estimating temperature response from an instantaneous point heat source in an infinite medium.Based on this concept, Carslaw and Jaeger [26] developed the Green function as shown in Equation (1): where a is the thermal diffusivity (m 2 ¨s´1 ), τ is the time (s) and `x1 , y 1 , z 1 ˘is the location of an instantaneous point heat source.
With the Green function, any shape of heat sources can be considered by integrating the point heat source in the spatial domain.If the point heat sources are accumulated to form a line-shaped heat source while neglecting its radial dimension, the temperature response (i.e., temperature increment, θ " T ´T0 where T is the temperature at a certain time (K) and T 0 is the initial temperature (K)) induced by transient heat conduction in a homogeneous medium at the time, τ, can be calculated by integrating the Green function with respect to time and spatial domain, as shown in Equation (2): where Ei pzq " ş z 8 e u u du is the exponential integral function, q l is the heating rate (W/m) generated constantly in the line-shaped heat source, ρ is the density (kg/m 3 ), c is the specific heat (J/kgK), k is the thermal conductivity (W/mK), and r is the radial coordinate from the central axis of the heat source (m).
Accordingly, the solid cylinder heat source model and ring-coil heat source model can be obtained.The solid cylinder heat source model improves the existing infinite line heat source model by implementing numerous line heat sources as a cylindrical shape while the thermo-physical properties Energies 2016, 9, 318 4 of 19 inside the heat source medium are regarded as a homogeneous medium identical to the outside heat source medium [16]: where r 0 is the radius of cylindrical heat source (m) and h 1 and h 2 are the top and bottom locations of heat source in the z-direction, respectively.The ring-coil heat source model accumulates the point heat sources to form coil-shaped heat sources, which are located at the same interval, and a uniform heat flux is supplied from the initial time [17]: θ " where b is the coil pitch (m) and m is the number of coil-shaped heat sources.All of the aforementioned models are developed based on the following critical assumptions: (1) Both the inside and outside of the heat source circumference are considered as one homogeneous medium, and the thermo-physical properties are constant with temperature changes.(2) The homogeneous medium has a uniform initial temperature, T 0 , and the heating rate per unit length of heat source, q l (W/m), is constant from the initial state, τ " 0. (3) The heat source is buried in the homogeneous medium in accordance with z-axis.The thickness, mass and heat capacity of the heat source are neglected.(4) The temperature boundary of the medium (i.e., ground surface) is maintained as constant initial temperature, T 0 .

Analytical Model for Heat Transfer with Groundwater Flow
Groundwater flow seeping through buried heat sources is able to induce a combined heat transfer mechanism: heat conduction of the solid and water in porous medium, and advection due to the moving water.In order to take into account the effect of groundwater flow in the analytical heat transfer model, Zhang et al. [24] developed the combined solid cylinder heat source model by introducing two additional concepts: the effective heat transfer coefficient, U, (to consider the conduction of water and solid particles), and the moving heat source equation (for simulating the advective effect caused by groundwater flow).
To simulate the thermal advection using the analytical model, either the heat source should move through a fixed medium or the heat production of the homogeneous medium should move under the fixed heat source [27].In the combined analytical model for heat transfer, the heat source is assumed to be motionless while the medium moves only in the x direction with a constant flow velocity, u (m/s).The energy conservation equation in the fixed coordinate px, y, zq can be expressed as Equation ( 5) [22]: Energies 2016, 9, 318 5 of 19 The moving heat source equation can be obtained by solving Equation ( 5) with the Green function.The moving heat source equation is shown in Equation ( 6) [23]: Meanwhile, the governing equation of heat transfer in the saturated porous medium when the Darcy velocity is uniform with respect to the x direction in the entire domain is expressed as Equation ( 7): where U " uρ w c w {ρ c, which refers to the effective heat transfer coefficient.The combined heat transfer equation can then be formulated by substituting U for u in Equation ( 5).In the same way of the pure conduction model, Zhang et al. [24] obtained the combined solid cylinder heat source model by integrating the moving heat transfer function with time and space as shown in Equation ( 8): θ " In addition to the assumptions for the pure conduction model, two additional assumptions are introduced for the combined analytical heat transfer model as follows.
(1) The groundwater seeps through at a constant velocity, u, exclusively in the x-direction in the medium (2) The inner part of the heat source is the same as the surrounding medium, which means the velocity of groundwater flow through the inner part of the borehole with same velocity of the surrounding medium.
In the combined analytical heat transfer model, the second assumption (i.e., the inner part of the heat source is assumed to be the same as the surrounding medium) is the primary reason for the overestimation of the groundwater flow effect on the heat transfer of GHEXs in the analytical solution, because the borehole of GHEX is actually filled with less permeable materials such as grout or concrete than the surrounding ground formation.In this study, the applicability of the combined solid cylinder heat source model was evaluated by comparing it with the numerical analysis result considering the permeability and configuration of boreholes.

Overview of Numerical Model
A series of numerical analyses was performed using a commercial 3-D finite-element analysis program, COMSOL Multiphysics.COMSOL Multiphysics is a useful tool for solving Partial Differential Equations (PDEs) of multi-physical phenomena including fluid flow and heat transfer.In order to analyze the heat transfer in the saturated porous medium, the mass, momentum and energy conservation should be satisfied [28].The heat transfer equation in the saturated porous medium, Equation (7), was considered as a governing equation in the numerical model.The solution of Equation (7) for the energy conservation with respect to the heat conduction and convection is formulated to yield the temperature distribution in the time and space domain.The backward different formula (BDF) was adopted to deal with the time domain, while the finite element formulation was utilized to deal with the space domain.Subsequently, the PARDISO (Parallel sparse direct solver) method [29], which is the solver for large sparse symmetric and non-symmetric linear systems of equations, was used.
In this paper, the numerical analysis was performed under two different assumptions about how to deal with the inner medium of the heat source.First, the inner part of the heat source is regarded as the same material as that of the surrounding ground formation, which is the same assumption as that for the combined analytical heat transfer model.Second, the inner part of the heat source is modeled as an impermeable borehole.In addition, different diameter-length (D/L) ratios of heat source and varied groundwater velocities were considered in each case of numerical analysis.With respect to the configuration of heat source, two different D/L ratios were considered: the cylindrical shape of heat source with a diameter of 750 mm and depth of 30 m (i.e., the D/L ratio of 0.025 typically representing large diameter energy piles), and the thread-like shape of heat source with a diameter of 150 mm and depth of 200 m (i.e., the D/L ratio of 0.00075 that is a common configuration of closed-loop vertical GHEXs).For the groundwater velocity condition, three different cases were applied with consideration of the common groundwater velocity in the Korea peninsula, 0, 2.5 ˆ10 ´6 and 5.0 ˆ10 ´6 m/s.The predetermined conditions of each numerical analysis are summarized in Table 1.The surrounding ground formation was modeled as a saturated porous medium and the groundwater flow was applied only in the horizontal direction with a constant velocity.The change in the thermal behavior of the ground formation was simulated more than one month until the ground temperature converged to a certain value.The no-flux boundary condition was applied to the exterior boundaries of the ground formation domain with a diameter of 20 m, which is large enough to avoid the boundary effect on the temperature change inside domain.In other words, the temperature change in the ground formation cannot reach the outside boundary during the period of numerical simulation.The no-flux boundary condition was also imposed on the ground surface to comply with the boundary condition adopted in the finite analytical model that is derived by the image method of heat flux [16].In addition, heat flux was applied to the entire surface of the borehole wall (i.e., heat source wall) to reproduce the same condition as that of the analytical model with the constant heating power of 7 kW.The heating power was arbitrarily determined because the main purpose of this paper is to observe the effect of groundwater flow and borehole permeability on the thermal behavior of the ground rather than to observe the heating power effect itself.Then, the temperature change in the ground formation was numerically calculated at the predetermined location of 0 m, 0.15 m, 0.30 m, and 0.50 m in the horizontal direction from the heat source to compare with the analytical solutions.The thermal properties of the ground formation and borehole material applied in the numerical analysis are summarized in Table 2, and an overview of the numerical model is illustrated in Figure 1.purpose of this paper is to observe the effect of groundwater flow and borehole permeability on the thermal behavior of the ground rather than to observe the heating power effect itself.Then, the temperature change in the ground formation was numerically calculated at the predetermined location of 0 m, 0.15 m, 0.30 m, and 0.50 m in the horizontal direction from the heat source to compare with the analytical solutions.The thermal properties of the ground formation and borehole material applied in the numerical analysis are summarized in Table 2, and an overview of the numerical model is illustrated in Figure 1.

Effect of Impermeable Borehole on Ground Thermal Behavior
The results of numerical analysis without considering the groundwater effect are shown in Figures 2-5.Figures 2 and 3 show the temperature contours obtained from the numerical model after one month operation with two different D/L ratios of heat source, respectively.Figures 4 and 5 compare the numerical analysis results with the analytical solutions (i.e., temperature variation at the location of 0 m, 0.15 m, 0.30 m and 0.50 m from the heat source wall) for two different D/L ratios of the heat source.In comparing the results, the time domain was expressed as dimensionless time, Fo " ατ{r b 2 .Without groundwater flow, only heat conduction occurs, which means that the combined analytical solution becomes the pure conduction solution.The comparison indicates that the developed numerical model can simulate pure conduction showing a good agreement with the analytical solution in both D/L ratios when the inner part of the heat source (i.e., borehole) is assumed to be the same as the surrounding ground formation.However, when the inner part of the heat source is modeled as impermeable material, the analytical solution shows a marginal deviation from the numerical analysis result, especially for the D/L ratio of 0.025, because the effect of different thermo-physical properties is enhanced with an increase in borehole volume.This result coincides with the result by Park et al. [19], where the validity of homogeneous models on short-term thermal response for a large-diameter cast-in-place energy pile was evaluated by comparing analytical solutions with field test results.

Effect of Impermeable Borehole on Ground Thermal Behavior
The results of numerical analysis without considering the groundwater effect are shown in Figures 2-5.Figures 2 and 3 show the temperature contours obtained from the numerical model after one month operation with two different D/L ratios of heat source, respectively.Figures 4 and 5 compare the numerical analysis results with the analytical solutions (i.e., temperature variation at the location of 0 m, 0.15 m, 0.30 m and 0.50 m from the heat source wall) for two different D/L ratios of the heat source.In comparing the results, the time domain was expressed as dimensionless time, Fo = ατ/  2 .Without groundwater flow, only heat conduction occurs, which means that the combined analytical solution becomes the pure conduction solution.The comparison indicates that the developed numerical model can simulate pure conduction showing a good agreement with the analytical solution in both D/L ratios when the inner part of the heat source (i.e., borehole) is assumed to be the same as the surrounding ground formation.However, when the inner part of the heat source is modeled as impermeable material, the analytical solution shows a marginal deviation from the numerical analysis result, especially for the D/L ratio of 0.025, because the effect of different thermophysical properties is enhanced with an increase in borehole volume.This result coincides with the result by Park et al. [19], where the validity of homogeneous models on short-term thermal response for a large-diameter cast-in-place energy pile was evaluated by comparing analytical solutions with field test results.
(a) (b)      The results of numerical analysis with the groundwater velocity of 2.5 × 10 −6 m/s are shown in Figures 6-9 and those with the groundwater velocity of 5.0 × 10 −6 m/s are shown in Figures 10-13.Figures 6, 7, 10 and 11 show the temperature contours obtained from the numerical model, and Figures 8, 9, 12 and 13 present a comparison of the analytical solution and numerical analysis result.In the comparison, the temperature variation at horizontal distance of 0, 0.15, 0.30 and 0.50 m from heat source wall is compared to each other.For the thermal response estimated at the front side of the heat source, however, only the thermal behaviors of ground at the borehole wall and a horizontal distance of 0.15 m were compared with the analytical solutions because the temperature change of the ground at the distances of 0.30 and 0.50 m from the front side of the heat source is not significant enough to compare with the analytical solutions.The results of numerical analysis with the groundwater velocity of 2.5 ˆ10 ´6 m/s are shown in Figures 6-9 and those with the groundwater velocity of 5.0 ˆ10 ´6 m/s are shown in Figures 10-13.Figures 6, 7, 10 and 11 show the temperature contours obtained from the numerical model, and Figures 8, 9, 12 and 13 present a comparison of the analytical solution and numerical analysis result.In the comparison, the temperature variation at horizontal distance of 0, 0.15, 0.30 and 0.50 m from heat source wall is compared to each other.For the thermal response estimated at the front side of the heat source, however, only the thermal behaviors of ground at the borehole wall and a horizontal distance of 0.15 m were compared with the analytical solutions because the temperature change of the ground at the distances of 0.30 and 0.50 m from the front side of the heat source is not significant enough to compare with the analytical solutions.
In the comparison, the temperature variation at horizontal distance of 0, 0.15, 0.30 and 0.50 m from heat source wall is compared to each other.For the thermal response estimated at the front side of the heat source, however, only the thermal behaviors of ground at the borehole wall and a horizontal distance of 0.15 m were compared with the analytical solutions because the temperature change of the ground at the distances of 0.30 and 0.50 m from the front side of the heat source is not significant enough to compare with the analytical solutions.The groundwater firstly seeps into the front side of the heat source and then later passes through the back side of the heat source as shown in Figure 1, which means the temperature responses at the front and back side of the heat source show different tendencies.At the front side of the heat source, the groundwater flow significantly reduces the temperature by eliminating the heat.This tendency is enhanced with an increase in groundwater velocity.On the other hand, the groundwater flow causes a larger temperature increment at the back side of the heat source than at the front side by accumulating the heat taken up from the front side.In addition, the underground thermal environment is varied with the different inner material of the heat source.When the inner part of the heat source is modeled as impermeable material, the thermal convection due to groundwater flow is reduced compared to the borehole modeled as the same material as that of the ground formation.Therefore, the temperature at the back side of the heat source more rapidly increases when the borehole is modeled as impermeable material because the groundwater bypasses the impermeable borehole.Figure 14 summarizes the temperature responses at the heat source wall with consideration of different groundwater velocities and borehole materials.When the inner part of borehole is modeled as the same as the surrounding ground formation, the presence of groundwater flow causes a steady-state, where the temperatures converge at a certain value, because thermal energy is removed by advection.In case of higher groundwater velocity, the temperature converges in the earlier stage.At the front side of the heat source, the discrepancy between the permeable and impermeable borehole conditions occurs only at a steady-state condition.Especially, as for the highest groundwater velocity (i.e., u = 5.0 ˆ10 ´6 m/s), the analytical solution underestimates the temperature variation by 10% for the D/L ratio of 0.00075 and 25% for the D/L ratio of 0.025, which indicates the effect of groundwater flow on the temperature response of GHEX becomes marked with increasing D/L ratio.
On the other hand, the back side of the heat source is more directly affected by the impermeability of borehole because the groundwater bypasses it.Therefore, the impermeable borehole condition leads to a more rapid temperature increase at the back side of the heat source than the permeable borehole condition.The root mean square error (RMSE) of the numerical analysis results between the permeable and impermeable borehole conditions at the back side of the heat source are shown in Figure 15 corresponding to the groundwater velocity and the location of temperature measurement.The maximum discrepancy is estimated near the heat source wall to be 12.33 ˝C at the D/L ratio of 0.025 and 0.59 ˝C at the D/L ratio of 0.00075 when u = 5.0 ˆ10 ´6 m/s.advection.In case of higher groundwater velocity, the temperature converges in the earlier stage.At the front side of the heat source, the discrepancy between the permeable and impermeable borehole conditions occurs only at a steady-state condition.Especially, as for the highest groundwater velocity (i.e., u = 5.0 × 10 −6 m/s), the analytical solution underestimates the temperature variation by 10% for the D/L ratio of 0.00075 and 25% for the D/L ratio of 0.025, which indicates the effect of groundwater flow on the temperature response of GHEX becomes marked with increasing D/L ratio.On the other hand, the back side of the heat source is more directly affected by the impermeability of borehole because the groundwater bypasses it.Therefore, the impermeable borehole condition leads to a more rapid temperature increase at the back side of the heat source than the permeable borehole condition.The root mean square error (RMSE) of the numerical analysis results between the permeable and impermeable borehole conditions at the back side of the heat source are shown in Figure 15 corresponding to the groundwater velocity and the location of temperature measurement.The maximum discrepancy is estimated near the heat source wall to be 12.33 °C at the D/L ratio of 0.025 and 0.59 °C at the D/L ratio of 0.00075 when u = 5.0 × 10 −6 m/s.Consequently, the discrepancy between the analytical solution and the numerical analysis result in case of the impermeable borehole becomes significant with increasing groundwater velocity and D/L ratio (refer to Figures 9 and 13), especially at the back side of the heat source.On the other hand, in the case of D/L ratio of 0.00075, the discrepancy between the analytical solution and the numerical analysis result is relatively smaller than that of D/L ratio of 0.025 (refer to Figures 8 and 12).
In this paper, the developed numerical model successfully reflects the effect of ground water and permeability of borehole.Therefore, the numerical model can overcome the limitation of the analytical solution proposed by Zhang et al. [24] that unavoidably assumes the inner part of heat Consequently, the discrepancy between the analytical solution and the numerical analysis result in case of the impermeable borehole becomes significant with increasing groundwater velocity and D/L ratio (refer to Figures 9 and 13), especially at the back side of the heat source.On the other hand, in the case of D/L ratio of 0.00075, the discrepancy between the analytical solution and the numerical analysis result is relatively smaller than that of D/L ratio of 0.025 (refer to Figures 8 and 12).
In this paper, the developed numerical model successfully reflects the effect of ground water and permeability of borehole.Therefore, the numerical model can overcome the limitation of the analytical solution proposed by Zhang et al. [24] that unavoidably assumes the inner part of heat source occupied by the same material as that of the surrounding ground formation.As a result, the critical assumption adopted in the combined analytical solution leads to the following outcomes: (1) The combined analytical solution largely deviates from the actual underground thermal behavior around buried GHEXs with a large diameter and relatively small length (i.e., large value of D/L ratio).( 2) The accuracy of the combined analytical solution is significantly reduced at the back side of the heat source, which is more affected by the impermeable borehole than the front side of the heat source.(3) As the groundwater velocity increases, the error caused by the critical assumption in the analytical solution increases.
Even though the analytical model for heat transfer is widely used in practice for designing GHEXs due to its clarity and convenience of computation, the analytical solution may significantly overestimate the effect of groundwater flow on the heat transfer of GHEXs (i.e., underestimate the temperature increment), especially in large D/L ratios of borehole, because the analytical model disregards the existence of the impermeable zone in the borehole.Consequently, such tendency will be more critical in the case of large-diameter energy piles at a high groundwater velocity.Because large-diameter energy piles have been adopted in many public buildings recently [30,31], the development of a novel design method applicable to the large-diameter energy pile will become an important issue in the GSHP system field.

Conclusions
The combined analytical heat transfer model for heat transfer proposed by Zhang et al. [24] is widely used in practice for designing GHEXs due to its clarity and convenience of computation.However, this analytical model has a critical limitation because it assumes that the inner part of the borehole consists of the same material as that of the surrounding ground formation, while an actual borehole of GHEX is filled with less permeable materials such as grout or concrete than the ground formation.In this paper, applicability of the combined solid cylindrical heat source model was evaluated by performing a series of numerical analyses considering the effect of groundwater flow and permeability of boreholes.The key conclusions of this paper are summarized as follows: (1) A series of numerical analyses was performed with different borehole materials (i.e., same material as that of the surrounding ground formation and impermeable material), groundwater velocities (i.e., 0 m/s, 2.5 ˆ10 ´6 m/s, and 5.0 ˆ10 ´6 m/s) and D/L ratios of heat source (i.e., 0.00075 and 0.025).The developed numerical model successfully reflects the effect of ground water and permeability of borehole.(2) In the case of no groundwater flow, the developed numerical model can simulate pure conduction showing a good agreement with the analytical solution in both D/L ratios when the inner part of the heat source is assumed to be the same as the surrounding ground formation.However, when the inner part of the heat source is modeled as impermeable material, the analytical solution slightly deviates from the numerical analysis result, especially for the D/L ratio of 0.025, because the effect of different thermo-physical properties is enhanced with an increase in borehole volume.(3) The groundwater flow significantly reduces the temperature at the front side of the heat source by eliminating the heat.On the other hand, the groundwater flow causes a larger temperature increment at the back side of the heat source than at the front side by accumulating the heat taken up from the front side.However, if the borehole is impermeable, the thermal convection due to groundwater flow is significantly reduced compared to the permeable borehole, which leads to a rapid increase in temperature at the back side of the heat source.(4) The accuracy of the combined solid cylindrical heat source model is reduced at the back side of the heat source, which is more affected by the existence of impermeable borehole than the front side of the heat source.Moreover, as the groundwater velocity and the D/L ratio of heat source increase, the error caused by the critical assumption in the combined analytical solution also increases.

Figure 1 .
Figure 1.Schematic overview of the numerical model.

Figure 1 .
Figure 1.Schematic overview of the numerical model.

Figure 8 .
Figure 8.Comparison of analytical solution and numerical analysis results (u = 2.5 ˆ10 ´6 m/s, D/L ratio = 0.00075).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 8 .Figure 9 .
Figure 8.Comparison of analytical solution and numerical analysis results (u = 2.5 × 10 −6 m/s, D/L ratio = 0.00075).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 9 .Figure 9 .Figure 10 .Figure 11 .
Figure 9.Comparison of analytical solution and numerical analysis results (u = 2.5 ˆ10 ´6 m/s, D/L ratio = 0.025).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 12 .
Figure 12.Comparison of analytical solution and numerical analysis results (u = 5.0 × 10 −6 m/s, D/L ratio = 0.00075).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 12 .
Figure 12.Comparison of analytical solution and numerical analysis results (u = 5.0 ˆ10 ´6 m/s, D/L ratio = 0.00075).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 13 .
Figure 13.Comparison of analytical solution and numerical analysis results (u = 5.0 × 10 −6 m/s, D/L ratio = 0.025).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 13 .
Figure 13.Comparison of analytical solution and numerical analysis results (u = 5.0 ˆ10 ´6 m/s, D/L ratio = 0.025).Temperature increment at (a) wall of front side; (b) 0.15 m of front side; (c) wall of back side; (d) 0.15 m of back side; (e) 0.30 m of back side; (f) 0.50 m of back side.

Figure 14 .
Figure 14.Temperature responses at heat source wall with different groundwater velocity and different permeability of borehole.Temperature increment at (a) front side of heat source (D/L ratio = 0.00075); (b) back side of heat source (D/L ratio = 0.00075); (c) front side of heat source (D/L ratio = 0.025); (d) back side of heat source (D/L ratio = 0.025).

Figure 14 .
Figure 14.Temperature responses at heat source wall with different groundwater velocity and different permeability of borehole.Temperature increment at (a) front side of heat source (D/L ratio = 0.00075); (b) back side of heat source (D/L ratio = 0.00075); (c) front side of heat source (D/L ratio = 0.025); (d) back side of heat source (D/L ratio = 0.025).

Figure 15 .
Figure 15.RMSE of numerical analysis between permeable and impermeable borehole condition at the back side of heat source.

Figure 15 .
Figure 15.RMSE of numerical analysis between permeable and impermeable borehole condition at the back side of heat source.

Table 1 .
Conditions of each numerical analysis case.

Table 2 .
Properties of materials for numerical analysis.

Table 2 .
Properties of materials for numerical analysis.