Performance Evaluation of Borehole Heat Exchanger in Multilayered Subsurface

In layered subsurface, the soil around a vertical borehole heat exchanger (BHE) contains different geological layers. Non-uniformity and groundwater flow can affect the performance of BHE drastically. In this paper, through the field investigation of boreholes in Zhu Shan, Nanjing, China, a numerical model considering five strata is developed. Using thermal resistance and capacity models for inside borehole and a combination of a locally refined grid for discretizing and solving the soil mass governing equations, the numerical model is calculated and validated by field test data. The maximum temperature difference never exceeds 0.3 ◦C. The numerical model is also compared with the homogenous finite line source (FLS) model. Based on the numerical multilayered model, the axial temperature profile at different distances under different heating times are presented and explored. After 60 days heating at the distance of 0.2 m to heat injection borehole, the maximum temperature rise is 9.2 ◦C in unsaturated soil layer, but the temperature rise in aquifer layer and in fractured layer are only 7.6 ◦C and 6.7 ◦C, respectively. Furthermore, two modified numerical layered models, in which the groundwater flow in aquifer or fracture layer is negligible, are established to analyze how the different layered characteristics impact on performance of BHE. The results showed that ignoring the groundwater flow in aquifer layer made the outlet temperature 0.7 ◦C higher than that of the original numerical layered model.


Introduction
A large number of GSHP (ground source heat pump) systems have been applied in residential and commercial buildings worldwide.Vertical borehole heat exchangers (BHEs) with 100 m-150 m depth are the most important parts of GSHP system.The performance of GSHP systems is affected by soil stratigraphy, in which thermal conductivity, groundwater flow and initial temperature play an essential role [1].Detailed and accurate information of thermal behavior of subsoil layers crossed by perforation is a prerequisite for improving the ratio between the heat transfer optimization and cost of the installations [2].
To evaluate the performance of BHE, many analytical or numerical models have been developed.Eskilson [3] proposed the finite line source (FLS) model using the numerical finite-difference method.Zeng et al. [4] improved the FLS model by imposing a constant temperature at the ground surface.Lamarche and Beauchamp [5] developed alternative forms for FLS model with shorter computation time.Bandos, T.V. et al. [6] presented three-dimensional finite line-source (FLS) model for BHEs that considered the prevailing geothermal gradient and allowed arbitrary changes in ground surface found from 16 m to 40 m depth, and the main direction of the groundwater flow is from borehole 2 to borehole 3, as illustrated in Figure 1.
Sustainability 2017, 9, 356 3 of 16 are drilled.The radius of the hole, and inner and outer diameters of the U-tubes are 0.2 m, 0.026 m and 0.032 m, respectively.From the geologic exploration and method of salting tracer, the aquifer layer is found from 16 m to 40 m depth, and the main direction of the groundwater flow is from borehole 2 to borehole 3, as illustrated in Figure 1.The field investigations indicate that the geological profile of borehole 2# could mainly be divided into five layers, the range of different layers and descriptions are shown in Table 1.The thermal conductivities of the samples are measured in laboratory.Table 2 shows typical thermal parameters of soil-rock sample at different depths.The thermal response test (TRT) is also conducted under the flow velocity of 0.72 m/s in tube (Re = 20634) for 60 hours inside borehole 2#, the effective subsurface thermal conductivities (λTRT) is measured as 2.14 W/(m•K).The field investigations indicate that the geological profile of borehole 2# could mainly be divided into five layers, the range of different layers and descriptions are shown in Table 1.The thermal conductivities of the samples are measured in laboratory.Table 2 shows typical thermal parameters of soil-rock sample at different depths.The thermal response test (TRT) is also conducted under the flow velocity of 0.72 m/s in tube (Re = 20,634) for 60 hours inside borehole 2#, the effective subsurface thermal conductivities (λ TRT ) is measured as 2.14 W/(m•K).

Model Description
According to the experimental investigations, a numerical model including unsaturated soil layer, aquifer layer, impervious layer, fractured rock layer and dense layer is established (Figure 2).The groundwater flow exists in aquifer and rock fracture layer at depth of 16-40 m and 55-85 m, respectively.To simplify the calculation, groundwater is assumed to be along horizontal direction.Site specific geological parameters such as the subsurface heterogeneities, thickness of layer or groundwater flow are constrained by field observations.By using the method of cubic spline interpolation, the thermal properties of the subsurface varying with depth are obtained through experimental investigations in Table 2.

Model Description
According to the experimental investigations, a numerical model including unsaturated soil layer, aquifer layer, impervious layer, fractured rock layer and dense layer is established (Figure 2).The groundwater flow exists in aquifer and rock fracture layer at depth of 16-40 m and 55-85 m, respectively.To simplify the calculation, groundwater is assumed to be along horizontal direction.Site specific geological parameters such as the subsurface heterogeneities, thickness of layer or groundwater flow are constrained by field observations.By using the method of cubic spline interpolation, the thermal properties of the subsurface varying with depth are obtained through experimental investigations in Table 2.The numerical system is decomposed into two subdomains, one representing a soil mass, and another representing borehole in Figure 3. Thermal resistance and capacity models (TRCMs) [35] is used to reduce the number of nodes and computation time.The thermal short-circuiting resistance between the upward and downward tubes given by Yong Li et al. [36] is also introduced to obtain a more accurate heat transfer model inside borehole.
In Figure 3, to save the computational amount, the 3-D (X-Y-Z) model is simplified as 2-D (Z-X) transient numerical model in rectangular coordinate system.The axial heat transfer along the vertical z-axis is carefully taken into account, the groundwater in aquifer or fracture rock layer is assumed to flow along the horizontal direction.The numerical system is decomposed into two subdomains, one representing a soil mass, and another representing borehole in Figure 3. Thermal resistance and capacity models (TRCMs) [35] is used to reduce the number of nodes and computation time.The thermal short-circuiting resistance between the upward and downward tubes given by Yong Li et al. [36] is also introduced to obtain a more accurate heat transfer model inside borehole.
In Figure 3, to save the computational amount, the 3-D (X-Y-Z) model is simplified as 2-D (Z-X) transient numerical model in rectangular coordinate system.The axial heat transfer along the vertical z-axis is carefully taken into account, the groundwater in aquifer or fracture rock layer is assumed to flow along the horizontal direction.

Governing Equations
In this section, the equations for heat transfer inside and soil in layered subsurface will be presented.
In Figure 3, the heat transfer governing equations for circulated fluid inside downward and upward tubes node are expressed as: Downward tube: ( ) ( ) Upward tube: where Tf1 and Tf2 are fluid temperature downward and upward tubes, respectively.Tb is the borehole temperature.uf is flow velocity in the tube.λf and (ρc)f are the thermal conductivity and volume thermal capacity of the fluid, respectively.Af is the cross-sectional area of tube.Rfb is the unit length borehole resistance including pipe resistance, which is calculated as: Rb is pipe-to-borehole thermal resistance given by Sharqawy, M.H. et al. [37]: R12 is the thermal short-circuiting resistance [36] between upward and downward tubes:

Governing Equations
In this section, the equations for heat transfer inside and soil in layered subsurface will be presented.
In Figure 3, the heat transfer governing equations for circulated fluid inside downward and upward tubes node are expressed as: Downward tube: Upward tube: where T f1 and T f2 are fluid temperature downward and upward tubes, respectively.T b is the borehole temperature.u f is flow velocity in the tube.λ f and (ρc) f are the thermal conductivity and volume thermal capacity of the fluid, respectively.A f is the cross-sectional area of tube.R fb is the unit length borehole resistance including pipe resistance, which is calculated as: R b is pipe-to-borehole thermal resistance given by Sharqawy, M.H. et al. [37]: R 12 is the thermal short-circuiting resistance [36] between upward and downward tubes: The heat transfer in soil is greatly influenced by the layered system.According to the experimental investigations, the numerical layered subsurface will be conducted with five geological layers.Through the method of cubic spline interpolation, λ s (z) and (ρc) s (z), obtained through the experimental investigations shown in Table 2, are converted into a function of depth z.
At shallow depth, the medium is unsaturated soil.The heat transfer model is expressed as: where T 1 (z, x, t) is the temperature distribution in unsaturated soil layer.λ 1 and (ρc) 1 are the thermal conductivity and volume thermal capacity of unsaturated soil, respectively.In aquifer layer, the advective heat transfer should be taken into account.The heat transport is given as: (ρc) where T 2 (z, x, t) is the temperature distribution in aquifer layer.λ por and (ρc) por are the bulk thermal conductivity and volume thermal capacity of medium, respectively.λ w , (ρc) w , and u wr are thermal conductivity, volume thermal capacity and groundwater flow velocity in horizontal direction, respectively.φ 2 is porosity of the medium, and its value is set as 0.25.
In the impervious layer, the subsurface with low permeability is hard and compact rock, and the heat transfer is mainly conductive model; the equation is written as: where T 3 (z, x, t) is the temperature distribution in impervious layer.λ s3 and (ρc) s3 are the thermal conductivity and volume thermal capacity of medium in this layer, respectively.In fractured rock layer, the moving groundwater flow exists in the fracture.S.E.A. Gehlin et al. [38] indicated that groundwater flow in fractures, even at relatively low specific flow rates, may cause significantly enhanced heat transfer for BHE.The rock zone is assumed to be completely impermeable: the groundwater only flows through the fracture zone horizontally following Darcy's law shown in Figure 4.
Sustainability 2017, 9, 356 6 of 16 The heat transfer in soil is greatly influenced by the layered system.According to the experimental investigations, the numerical layered subsurface will be conducted with five geological layers.Through the method of cubic spline interpolation, λs(z) and (ρc)s(z), obtained through the experimental investigations shown in Table 2, are converted into a function of depth z.
At shallow depth, the medium is unsaturated soil.The heat transfer model is expressed as: where T1(z, x, t) is the temperature distribution in unsaturated soil layer.λ1 and (ρc)1 are the thermal conductivity and volume thermal capacity of unsaturated soil, respectively.In aquifer layer, the advective heat transfer should be taken into account.The heat transport is given as: where T2 (z, x, t) is the temperature distribution in aquifer layer.λpor and (ρc)por are the bulk thermal conductivity and volume thermal capacity of medium, respectively.λw, (ρc)w, and uwr are thermal conductivity, volume thermal capacity and groundwater flow velocity in horizontal direction, respectively.φ2 is porosity of the medium, and its value is set as 0.25.
In the impervious layer, the subsurface with low permeability is hard and compact rock, and the heat transfer is mainly conductive model; the equation is written as: where T3(z, x, t) is the temperature distribution in impervious layer.λs3 and (ρc)s3 are the thermal conductivity and volume thermal capacity of medium in this layer, respectively.In fractured rock layer, the moving groundwater flow exists in the fracture.S.E.A. Gehlin et al. [38] indicated that groundwater flow in fractures, even at relatively low specific flow rates, may cause significantly enhanced heat transfer for BHE.The rock zone is assumed to be completely impermeable: the groundwater only flows through the fracture zone horizontally following Darcy's law shown in Figure 4.In rock zone, the heat transfer model is written as: In rock zone, the heat transfer model is written as: Sustainability 2017, 9, 356 7 of 16 where T 4 (z, x, t) is the temperature distribution in rock zone.λ s4 , (ρc) s4 varying with depth are the thermal conductivity and volume thermal capacity of rock, respectively.In fracture zone, the heat transfer between the fracture water and surrounding rock is in the form of convective scheme.With the aid of instantaneous local thermal equilibrium presented by Cheng et al. [39], and Ghiassemi et al. [40], the heat transfer in fracture with water moving is expressed as: where T x is the temperature distribution of fracture.λ por-x , (ρc) por-x , u xx are the bulk thermal conductivity, volume thermal capacity and fracture water flow velocity of fracture, respectively.
λ por-x and (ρc) por-x are calculated by: (ρc φ 4 is the porosity in fracture, its value is set as 0.65.In the dense layer, the rock becomes more hard and compact, the fracture aperture becomes smaller and smaller, the groundwater in fracture keeps static, and the advective heat transfer form can be ignored; thus, the heat transfer model is expressed as Equation ( 12) in rock and Equation (13) in the fracture, respectively.
In rock: In fracture with static water: where T 5 (z, r, t) and T zx are the temperature distribution of rock and fracture, respectively.λ s5 and (ρc) s5 vary with depth are the thermal conductivity and volume thermal capacity of rock, respectively.λ por-zx and (ρc) por-zx are the bulk thermal conductivity and volume thermal capacity of fracture, respectively.

Initial and Boundary Conditions
Considering the surface temperature variations and geothermal gradient, the initial ground temperature is calculated by: where z is ground depth, t is the running time, T sur is the average atmospheric temperature, and A s is the amplitude observed during the 24 h cycles that constitute the measured temperature wave time.t 24h is equal to 24 h.When t = 0, T f1 (z, t = 0) = T f2 (z, t = 0) = T j (z, x, t = 0) = T g (z, t = 0) (where j = 1, 2, . . ., 5 indicates the different layers).
At the interface between the subsurface and the borehole wall, heat exchange is expressed as following mixed Neumann's boundary condition: where T b is the borehole wall temperature at depth of z.K eq is equivalent convective heat transfer coefficient inside the borehole and is obtained by: T f (z, t) is the arithmetic average of the inlet and outlet fluid temperature in cross-section, it is given as the mean fluid temperature: where z = 0, the model is subjected to the boundary conditions: where T air is atmospheric temperature, which is measured every minute.h air is the convective heat transfer coefficient between subsurface and air.Inside borehole, where z = 0, T fi and T fo are the inlet and outlet temperature of the borehole, respectively.They are measured by the data logging system.T fi is chosen as the input parameters for numerical model, and the outlet temperature T fo would be computed to compare with the field tested data. where In the x-direction, x ∞ = 200•r b is chosen as farfield boundary condition: In vertical direction, according to Eskilson [3], the distance traveled by a heat front from the heat source after time t can be estimated as ∆z = 3 √ at max .The maximum time t max is assumed to be 1000 days, where z ∞ = H + 3 √ at max , Moreover, at the interface between two layers, two boundary conditions must be satisfied.The continuity of temperature: and the continuity of the heat flux: The amount of heat flow ( Q) injected into the borehole is calculated as: Equations ( 1)- (25), which govern heat transfer of BHE in layered subsurface, are solved numerically using finite volume method, with an upwind difference scheme.A combination of a locally refined grid and a multigrid with hierarchal tree data structure was used for solving the soil mass governing equations to make the model computationally efficient.The model domain is vertically discretized in accordance with the geological layers.The mesh is refined around the boundary between layers.A combination of the TDMA and Gauss-Seidel method in an iterative solver are used to calculate outlet temperature, and the outlet temperature T fo at time k + 1 is computed each iteration until the temperature difference between two iterations becomes less than 10 −3 .

Numerical Verification
In this section, the input parameters including the geometrical parameters and thermal parameters of borehole are elaborated in accordance with field test.The numerical model for BHE in layered subsurface is validated by using the tested data, and the FLS model for the homogeneous profile is also implemented to compare with the numerical model.

Input Parameters of the Numerical Model
Thermal properties of borehole, pipe, grouting material and heat carrier fluid are obtained from the product specifications.Table 3 shows the geometrical and thermal parameters of borehole.The flow velocity in tube is obtained from the field test by magnetic flowmeter.In Figure 3, the geological profile in soil is subdivided into five different bedded layers.The parameters of the five different bedded layers are showed in Table 1.By the method of cubic spline interpolation, λ s (z) and (ρc) s (z) are converted into the function of depth z according to the experimental data given in Table 2.The groundwater flow velocity is given as 2 × 10 −7 m/s and 4 × 10 −5 m/s in the aquifer and fracture layer, respectively.
According to the investigations of rock-soil sample at the depth of 55 m-100 m, there are about five fractures every meter.The widths of fractures are mainly 0.0015 m and 0.001 m at the depth of 55-85 m and 86-100 m, respectively.

Numerical Verification
Under the field test, the heating water was circulated through borehole 2# for 60 days and followed by a 65-day self-regeneration period.The inlet and outlet temperature were recorded every minute, therefore, the heat injection into borehole 2# is calculated.
With aid of the related parameters in Tables 1-3, according to Equations ( 1)-( 25), the simulated outlet temperature is calculated and compared to the tested outlet temperature.The results are presented in Figure 5.
It is obvious that the simulated outlet temperature conforms with the experimental outlet temperature very well.The deviation of the outlet temperature never exceeds 0.3 • C. The deviations of the numerical model are 0.25 • C and 0.28 • C on the 30th day and 75th day, respectively.

Comparison with FLS Model
The modified FLS model shown in Equation ( 26) developed by S. Koohi-Fayegh et al. [41] is also chosen to conduct comparisons with the numerical model.
where ΔT(x, y, z) is the temperature rise at time t on the position of (x, y, z).r is the distance to the borehole.H is the borehole length.qi−qi−1 is the incremental load between two successive hours.λs and as are the ground homogeneous thermal conductivity and thermal diffusivity, respectively.The effective subsurface thermal conductivities (λTRT) 2.14 W/(m•K) from TRT and volumetric heat capacity 3.52 × 10 6 J/(m 3 •K) are used in the FLS model.The comparison between FLS model and the numerical layered model is shown in Figure 5.In Figure 5, outlet temperature of FLS is 0.8 ℃ higher than the tested outlet temperature on the 10th day and the temperature different extends to 1.9 °C on the 60th day.During the self-regeneration period, the outlet temperature of FLS is also 1.2 °C higher than the tested outlet temperature on the 80th day.The FLS model, ignoring the heterogeneities of subsurface thermal properties, the groundwater flow and the rock fracture flow, leads to underestimate the performance of BHE.5.2 Temperature profile Based on the numerical MLM model, after 60 days heating injection period, the axial temperature rise distribution at distance of 0.2 m and 3 m to borehole are obtained in Figures 6 and  7, respectively.The refining temperature response distribution of the fracture zone is given as well.The curved temperature profile clearly indicates that heat transfer is not uniform in axial direction.
In Figure 6, at distance of 0.2 m to borehole, the temperature rise is 9.2 °C at the depth of 10 m in unsaturated soil layer.This is due to the smaller soil thermal diffusion in unsaturated soil layer.In aquifer layer, the temperature rise at the depth of 30 m is 7.6 °C, which is 1.2 °C and 0.3 °C lower than that at depth of 15 m and 45 m, respectively.This indicates that the groundwater flow makes the heat spread further and avoids heat accumulation.In fractured rock layer with moving water, the temperature rise in rock and fracture is different.In Figure 6a, the temperature difference between fracture and rock is 0.3 °C at the depth of 65-66 m, while, in Figure 6b, the temperature

Comparison with FLS Model
The modified FLS model shown in Equation ( 26) developed by S. Koohi-Fayegh et al. [41] is also chosen to conduct comparisons with the numerical model. ) where ∆T(x, y, z) is the temperature rise at time t on the position of (x, y, z).r is the distance to the borehole.H is the borehole length.q i −q i − 1 is the incremental load between two successive hours.λ s and a s are the ground homogeneous thermal conductivity and thermal diffusivity, respectively.The effective subsurface thermal conductivities (λ TRT ) 2.14 W/(m•K) from TRT and volumetric heat capacity 3.52 × 10 6 J/(m 3 •K) are used in the FLS model.The comparison between FLS model and the numerical layered model is shown in Figure 5.In Figure 5, outlet temperature of FLS is 0.8 °C higher than the tested outlet temperature on the 10th day and the temperature different extends to 1.9 • C on the 60th day.During the self-regeneration period, the outlet temperature of FLS is also 1.2 • C higher than the tested outlet temperature on the 80th day.The FLS model, ignoring the heterogeneities of subsurface thermal properties, the groundwater flow and the rock fracture flow, leads to underestimate the performance of BHE.5.2 Temperature profile Based on the numerical MLM model, after 60 days heating injection period, the axial temperature rise distribution at distance of 0.2 m and 3 m to borehole are obtained in Figures 6 and 7, respectively.The refining temperature response distribution of the fracture zone is given as well.The curved temperature profile clearly indicates that heat transfer is not uniform in axial direction.
In Figure 6, at distance of 0.2 m to borehole, the temperature rise is 9.2 • C at the depth of 10 m in unsaturated soil layer.This is due to the smaller soil thermal diffusion in unsaturated soil layer.In aquifer layer, the temperature rise at the depth of 30 m is 7.6 • C, which is 1.2 • C and 0.3 • C lower than that at depth of 15 m and 45 m, respectively.This indicates that the groundwater flow makes the heat spread further and avoids heat accumulation.In fractured rock layer with moving water, the temperature rise in rock and fracture is different.In Figure 6a, the temperature difference between fracture and rock is 0.3 • C at the depth of 65-66 m, while, in Figure 6b, the temperature difference is only 0.06 • C at the depth of 95-96 m in dense layer, in which the water in the fracture is static.This demonstrates that the fractured moving water should not be neglected when discussing the heat transfer of BHE in layerd system.Figure 7 presents the axial temperature rise distribution at the distance of 3 m to borehole.It shows that the maximum temperature response is 1.7 • C at the depth of 30 m in aquifer layer, which is 5.8 • C lower than that at the distance of 0.2 m to the borehole in Figure 6.The temperature response in unsaturated soil layer, impervious layer and dense layers are almost the same as 0.6 • C, which is 1.1 • C lower than that in aquifer layer.It indicates that the heat accumulates downwards along the groundwater flow.Combined with Figure 6, it is found that the heat accumulates around the borehole in unsaturated soil layer, impervious layer and dense layer.However, the heat is dissipated further by groundwater flow in aquifer.The groundwater in the aquifer leads to less heat accumulation around the heat injected borehole and improves the performance of BHE.The thermal interaction between two neighboring boreholes would become complex because of the layered subsurface and groundwater flow.Choosing the distance between the boreholes must take the multilayered system into account.
The temperature rise of fracture with moving water is about 0.2 • C higher than that in the rock region in Figure 7a, however, the temperature difference between the fracture and rock is only 0.04 • C in dense layer in Figure 7b.It indicates that the moving water, even fractured moving water would reduce the heat accumulation around the borehole.
In order to evaluate how the heating time affects the performance of BHE, Figure 8 presents the ground temperature rise profile at the distance of 1 m under different running times.It can be found that the temperature rise increases rapidly with time in the unsaturated soil layer, which accords with the temperature variation pattern of heat conduction.In aquifer layer, the temperature rise is larger than other layers at the beginning, but the temperature response rises slowly and tends to be stable.After the 365 days, the temperature rise at the depth of 30 m is 4.7 • C, which are 1.2 • C and 1.0 • C lower than that at the depth of 50 m and 90 m, respectively.In fractured rock layer with moving water, the temperature of the fracture is about 0.4 • C higher than that in the rock region after 30 days, but with running time increasing, the rock temperature is higher than the temperature of the fracture.The temperature difference is about 0.8 • C when the running time reaches to 365 days, it indicates that the moving water in fracture could also improve the heat transfer of BHE.However, in the dense layer with unmoving water, the temperature difference between rock and fracture becomes smaller with running time and tends to be zero.
Combined with the Figures 6-8, it is found that the axial temperature rise distribution is different in different layers.This is due to the heterogeneous thermal properties at different layers as well as the groundwater.The multilayered subsurface and groundwater flow leads to heat transfer of BHE become complex in axial direction, when designing the boreholes field, the multilayered subsurface system must be given careful consideration.
difference is only 0.06 °C at the depth of 95-96 m in dense layer, in which the water in the fracture is static.This demonstrates that the fractured moving water should not be neglected when discussing the heat transfer of BHE in layerd system.
Figure 7 presents the axial temperature rise distribution at the distance of 3 m to borehole.It shows that the maximum temperature response is 1.7 °C at the depth of 30 m in aquifer layer, which is 5.8 °C lower than that at the distance of 0.2 m to the borehole in Figure 6.The temperature response in unsaturated soil layer, impervious layer and dense layers are almost the same as 0.6 °C, which is 1.1 °C lower than that in aquifer layer.It indicates that the heat accumulates downwards along the groundwater flow.Combined with Figure 6, it is found that the heat accumulates around the borehole in unsaturated soil layer, impervious layer and dense layer.However, the heat is dissipated further by groundwater flow in aquifer.The groundwater in the aquifer leads to less heat accumulation around the heat injected borehole and improves the performance of BHE.The thermal interaction between two neighboring boreholes would become complex because of the layered subsurface and groundwater flow.Choosing the distance between the boreholes must take the multilayered system into account.
The temperature rise of fracture with moving water is about 0.2 °C higher than that in the rock region in Figure 7a, however, the temperature difference between the fracture and rock is only 0.04 °C in dense layer in Figure 7b.It indicates that the moving water, even fractured moving water would reduce the heat accumulation around the borehole.
In order to evaluate how the heating time affects the performance of BHE, Figure 8 presents the ground temperature rise profile at the distance of 1 m under different running times.It can be found that the temperature rise increases rapidly with time in the unsaturated soil layer, which accords with the temperature variation pattern of heat conduction.In aquifer layer, the temperature rise is larger than other layers at the beginning, but the temperature response rises slowly and tends to be stable.After the 365 days, the temperature rise at the depth of 30 m is 4.7 °C, which are 1.2 °C and 1.0 °C lower than that at the depth of 50 m and 90 m, respectively.In fractured rock layer with moving water, the temperature of the fracture is about 0.4 °C higher than that in the rock region after 30 days, but with running time increasing, the rock temperature is higher than the temperature of the fracture.The temperature difference is about 0.8 °C when the running time reaches to 365 days, it indicates that the moving water in fracture could also improve the heat transfer of BHE.However, in the dense layer with unmoving water, the temperature difference between rock and fracture becomes smaller with running time and tends to be zero.
Combined with the Figures 6-8, it is found that the axial temperature rise distribution is different in different layers.This is due to the heterogeneous thermal properties at different layers as well as the groundwater.The multilayered subsurface and groundwater flow leads to heat transfer of BHE become complex in axial direction, when designing the boreholes field, the multilayered subsurface system must be given careful consideration.

Different Multi-Layers Characteristics
To obtain how the different multilayered characteristics affect the performance of BHE, the numerical simulations are implemented for other two modified cases: Case (1) ignoring the groundwater flow in aquifer layer; and Case (2) assuming no fracture water in fractured rock layer.The thermal parameters in others layers are same as the above-mentioned parameters in Tables 1-3.The axial temperature response distribution of the three numerical models at the distance of 1 m after 60 days heating time is shown in Figure 9, and the outlet temperature of the models is given in Figure 10.
In Figure 9, the temperature rise of the Case (1) model at the depth of 25 m is about 0.5 °C lower than that of the original numerical multilayered model (MLM) in aquifer layer, but the temperature rise in the impervious layer, fractured rock layer and dense layer is about 0.3 °C higher than that of the original MLM model.The temperature is slightly higher under the Case (1) model.In Figure 10, during the beginning 10 days, the outlet temperature of the Case (1) model is almost consistent with the original MLM model, but the outlet temperature difference becomes larger with the heating time.On the 50th day, the outlet temperature of Case (1) model is 0.7 °C higher than that of the original MLM model.

Different Multi-Layers Characteristics
To obtain how the different multilayered characteristics affect the performance of BHE, the numerical simulations are implemented for other two modified cases: Case (1) ignoring the groundwater flow in aquifer layer; and Case (2) assuming no fracture water in fractured rock layer.The thermal parameters in others layers are same as the above-mentioned parameters in Tables 1-3.The axial temperature response distribution of the three numerical models at the distance of 1 m after 60 days heating time is shown in Figure 9, and the outlet temperature of the models is given in Figure 10.
In Figure 9, the temperature rise of the Case (1) model at the depth of 25 m is about 0.5 °C lower than that of the original numerical multilayered model (MLM) in aquifer layer, but the temperature rise in the impervious layer, fractured rock layer and dense layer is about 0.3 °C higher than that of the original MLM model.The temperature is slightly higher under the Case (1) model.In Figure 10, during the beginning 10 days, the outlet temperature of the Case (1) model is almost consistent with the original MLM model, but the outlet temperature difference becomes larger with the heating time.On the 50th day, the outlet temperature of Case (1) model is 0.7 °C higher than that of the original MLM model.

Different Multi-Layers Characteristics
To obtain how the different multilayered characteristics affect the performance of BHE, the numerical simulations are implemented for other two modified cases: Case (1) ignoring the groundwater flow in aquifer layer; and Case (2) assuming no fracture water in fractured rock layer.The thermal parameters in others layers are same as the above-mentioned parameters in Tables 1-3.The axial temperature response distribution of the three numerical models at the distance of 1 m after 60 days heating time is shown in Figure 9, and the outlet temperature of the models is given in Figure 10.
In Figure 9, the temperature rise of the Case (1) model at the depth of 25 m is about 0.5 • C lower than that of the original numerical multilayered model (MLM) in aquifer layer, but the temperature rise in the impervious layer, fractured rock layer and dense layer is about 0.3 • C higher than that of the original MLM model.The temperature is slightly higher under the Case (1) model.In Figure 10, during the beginning 10 days, the outlet temperature of the Case (1) model is almost consistent with the original MLM model, but the outlet temperature difference becomes larger with the heating time.On the 50th day, the outlet temperature of Case (1) model is 0.7 • C higher than that of the original MLM model.In Figure 9, it is obvious that the axial temperature would be higher when ignoring the influence of fracture water under the Case (2) model: the temperature response at depth of 70 m is 0.8 °C higher compared with the original numerical model.In Figure 10, the outlet temperature of Case (2) model is 0.9 °C higher than that of the original numerical model on the 60th day, which shows that ignoring the heat transfer of moving fracture water would underestimate the performance of BHE.

Conclusions
When boreholes are drilled in the layered subsurface, the traditional homogenous FLS model may lead to underestimate the performance of BHE and larger length of the borehole.
The site stratigraphy and the thermal properties of soil-rock samples at different depths are investigated through laboratory measurements.A numerical multilayered model considering five layers including unsaturated soil layer, aquifer layer, impervious layer, fractured rock layer and dense layer is established.The numerical system is decomposed into soil mass and borehole.Thermal resistance and capacity models (TRCMs) and thermal short-circuiting effect between the upward and downward tubes is introduced to obtain a more accurate and lower computation heat transfer model for inside borehole.A multi-grid with hierarchal tree data structure is used for solving the layered soil mass governing equations to make the model computationally efficient.
The numerical model and the modified FLS model are compared with the test outlet temperature.It shows that the maximum temperature different between the FLS and tested outlet temperature is 1.9 °C, but the maximum deviations between the numerical model and measured  In Figure 9, it is obvious that the axial temperature would be higher when ignoring the influence of fracture water under the Case (2) model: the temperature response at depth of 70 m is 0.8 °C higher compared with the original numerical model.In Figure 10, the outlet temperature of Case (2) model is 0.9 °C higher than that of the original numerical model on the 60th day, which shows that ignoring the heat transfer of moving fracture water would underestimate the performance of BHE.

Conclusions
When boreholes are drilled in the layered subsurface, the traditional homogenous FLS model may lead to underestimate the performance of BHE and larger length of the borehole.
The site stratigraphy and the thermal properties of soil-rock samples at different depths are investigated through laboratory measurements.A numerical multilayered model considering five layers including unsaturated soil layer, aquifer layer, impervious layer, fractured rock layer and dense layer is established.The numerical system is decomposed into soil mass and borehole.Thermal resistance and capacity models (TRCMs) and thermal short-circuiting effect between the upward and downward tubes is introduced to obtain a more accurate and lower computation heat transfer model for inside borehole.A multi-grid with hierarchal tree data structure is used for solving the layered soil mass governing equations to make the model computationally efficient.
The numerical model and the modified FLS model are compared with the test outlet temperature.It shows that the maximum temperature different between the FLS and tested outlet temperature is 1.9 °C, but the maximum deviations between the numerical model and measured In Figure 9, it is obvious that the axial temperature would be higher when ignoring the influence of fracture water under the Case (2) model: the temperature response at depth of 70 m is 0.8 • C higher compared with the original numerical model.In Figure 10, the outlet temperature of Case (2) model is 0.9 • C higher than that of the original numerical model on the 60th day, which shows that ignoring the heat transfer of moving fracture water would underestimate the performance of BHE.

Conclusions
When boreholes are drilled in the layered subsurface, the traditional homogenous FLS model may lead to underestimate the performance of BHE and larger length of the borehole.
The site stratigraphy and the thermal properties of soil-rock samples at different depths are investigated through laboratory measurements.A numerical multilayered model considering five layers including unsaturated soil layer, aquifer layer, impervious layer, fractured rock layer and dense layer is established.The numerical system is decomposed into soil mass and borehole.Thermal resistance and capacity models (TRCMs) and thermal short-circuiting effect between the upward and downward tubes is introduced to obtain a more accurate and lower computation heat transfer model for inside borehole.A multi-grid with hierarchal tree data structure is used for solving the layered soil mass governing equations to make the model computationally efficient.
The numerical model and the modified FLS model are compared with the test outlet temperature.It shows that the maximum temperature different between the FLS and tested outlet temperature is 1.9 • C, but the maximum deviations between the numerical model and measured outlet temperature are 0.25 • C under heating period and 0.28 • C under heat recovery period, respectively.This indicates the reliability of the developed numerical multilayered model.
The axial temperature profiles are obtained based on the numerical multilayered model.After 60 days heating at the distance of 0.2 m to heat injection borehole, the maximum temperature rise is 9.2 • C at the depth of 10 m in unsaturated soil layer, but the temperature rise at the depth of 30 m in aquifer layer and 78 m in fractured layer are only 7.6 • C and 6.7 • C, respectively.This is due to soil thermal diffusion difference and groundwater flow.At the distance of 3 m to heat injection borehole, the temperature response at depth of 30 m in aquifer is 1.8 • C, which is 5.8 • C lower than that at the distance of 0.2 m.This finding indicates heat transfer efficiency of the BHEs is obviously higher within aquifer layers.
The temperature rises rapidly with time in unsaturated soil layer.This is due to the lower thermal diffusivity.However, the temperature response in aquifer layer rises slowly and gradually tends to be stable with time.In the rock fracture layer, with the running time increasing, the rock temperature is higher than the fracture temperature with moving water.The temperature difference is about 0.8 • C after 365 days.However, in dense layer with unmoving water, the temperature difference between rock and fracture becomes smaller with running time and tends to be zero.
Two modified numerical MLM models are established to analyze how the different layered characteristics influence on the performance of BHE.The results show that ignoring the groundwater flow in aquifer layer makes the outlet temperature 0.7 • C higher than that of the original MLM model.At depth of 70 m, ignoring the fracture moving water would lead to the soil temperature rise at radial distance of 1 m 0.8

Figure 1 .
Figure 1.Boreholes locations and experiential system flow diagram.

Figure 1 .
Figure 1.Boreholes locations and experiential system flow diagram.

Figure 3 .
Figure 3.The numerical model description and simplification for multilayered borehole.

Figure 3 .
Figure 3.The numerical model description and simplification for multilayered borehole.
Figure for simplified model

Figure 4 .
Figure 4.The simplified heat transfer model of fracture rock.

Figure 4 .
Figure 4.The simplified heat transfer model of fracture rock.

Figure 5 .
Figure 5. Compared test outlet temperature with multilayered model and FLS model.

Figure 5 .
Figure 5. Compared test outlet temperature with multilayered model and FLS model.

Figure 6 .
Figure 6.The temperature rise distribution at the radial distance of 0.2 m to borehole after 60 days.Figure 6.The temperature rise distribution at the radial distance of 0.2 m to borehole after 60 days.

Figure 6 .
Figure 6.The temperature rise distribution at the radial distance of 0.2 m to borehole after 60 days.Figure 6.The temperature rise distribution at the radial distance of 0.2 m to borehole after 60 days.

Figure 7 .
Figure 7. Temperature rise distribution at the radial distance of 3 m to borehole after 60 days.

Figure 8 .
Figure 8. Temperature rise distribution of different running time at the distance of 1 m to borehole.

Figure 7 . 16 Figure 7 .
Figure 7. Temperature rise distribution at the radial distance of 3 m to borehole after 60 days.

Figure 8 .
Figure 8. Temperature rise distribution of different running time at the distance of 1 m to borehole.

Figure 8 .
Figure 8. Temperature rise distribution of different running time at the distance of 1 m to borehole.

Figure 9 .
Figure 9.The temperature response distribution at the radial distance of 1 m under different multilayered models.

Figure 10 .
Figure 10.The outlet temperature variation under different multilayered models.

Figure 9 . 16 Figure 9 .
Figure 9.The temperature response distribution at the radial distance of 1 m under different multilayered models.

Figure 10 .
Figure 10.The outlet temperature variation under different multilayered models.

Figure 10 .
Figure 10.The outlet temperature variation under different multilayered models.

Table 1 .
The depth and lithology for five bedded layers.

Table 2 .
The typical thermal parameters of soil-rock sample at different depth.

Table 2
lists only smaller part of the thermal parameters of soil-rock sample.

Table 1 .
The depth and lithology for five bedded layers.

Table 2 .
The typical thermal parameters of soil-rock sample at different depth.

Table 2
lists only smaller part of the thermal parameters of soil-rock sample.

Table 3 .
Geometrical and thermal parameters of borehole.
• C higher than that of the original MLM model.