A New Approach for Characterizing Pile Heat Exchangers Using Thermal Response Tests

Abstract

Pile heat exchangers offer a cost effective route to implementation of ground-source heat pump systems for many large commercial buildings compared with traditional boreholes. Such projects typically use thermal response tests to determine the key input parameters for system design, namely soil thermal conductivity and heat exchanger thermal resistance. However, this brings challenges for pile heat exchanger based systems, where in situ thermal response tests are known to be less reliable due to the large thermal capacity of the pile. This paper presents a new “black box” resistance capacitive model for applications to pile thermal response tests. The approach is tested against case study data and shown to perform well. Additional test duration savings are shown to be possible if a novel combination of borehole and pile thermal response tests is applied together to determine design parameters.

1. Introduction

Ground-Source Heat Pump (GSHP) systems can decrease the emission of greenhouse gases resulting from heating, cooling and hot water provision. Energy geostructures, which combine structural and thermal function, are an opportunity to reduce the installation costs of ground heat exchangers (GHE), which are the in-ground component of the GSHP system. Energy geostructures have been installed for several decades [1]. This includes the equipping of foundation piles [2,3], embedded retaining walls [4], or tunnels [5] with plastic heat transfer pipes so they can act as a GHE. Of the types of energy geostructure, pile heat exchangers (PHE), e.g., Figure 1, are the most common. This is partly because most are likely to be associated with an overlying building and hence have ready users of the heat they can supply, but also because their typically cylindrical shape makes them superficially similar in geometry to the more common borehole heat exchanger (BHE). Thermal analysis methods for BHE are well developed and hence are typically drawn on for application to PHE.

However, the radius of a PHE, which can exceed 50 cm, is much greater than the radius of a borehole heat exchanger (BHE) (typically 8 to 10 cm). PHE are also usually much shorter than BHE (typically 10 to 20 m against 100 m to 200 m). Therefore many approaches to the analysis of BHE, which assume the presence of a very long and thin heat source, may not be applicable to PHE. Consequently, more appropriate methods for thermal modeling such as fully discretized models (e.g., finite element analysis) to investigate the pipe arrangements and thermal performance [6,7] or semi-analytical models for thermal dynamic simulation of the whole energy system with hourly time steps [8] are often applied.

Outside of research, analysis methods for PHE are typically used for two main purposes. First, forward simulation is required for long-term design. In this scenario long duration (decades) analyses, often with hourly fluctuating thermal demand are required, and the PHE is often located beneath the overlying building. The thermal properties and the temperature limits are input, and the output is the available thermal power. In the second case, reverse simulation is carried out for thermal response test (TRT) interpretation. In this scenario, short duration (days to weeks) analysis is carried out for in situ characterization and the power applied and temperatures measured are input, while the thermal properties are the output. Typically, the PHE is open to the air because the overlying building is yet to be constructed.

During a TRT, a heat-carrying liquid (usually water) circulates in the pipes of the ground heat exchanger with a constant flow volume Q_v [m³·s⁻¹] while an electric heater supplies a constant power P [W] to the fluid. The entry T_in and exit T_out temperatures of the pile are recorded for the duration of the TRT [9]. The “classical” interpretation of a TRT is based on a purely resistive thermal model developed for BHE [9], i.e., a model overlooking the thermal inertia of the GHE. It assumes that the evolution of the temperature of the heat-carrying liquid T_f is described by the sum of the resistive component and transient infinite line source (ILS) step response, which takes the borehole to be a line emitting power with a constant rate:

$$T_f=\frac{T_{in}+T_{out}}{2}=T_0+p{R}_b+\frac{p}{4\pi {\lambda }_m}{E}_1\left(\frac{1}{4t^{*}}\right)$$

(1)

In Equation (1), T₀ is the undisturbed ground temperature, p is the linear power (W·m⁻¹) defined as the ratio between the power P and the depth of the exchanger H, $R_b$ is the borehole thermal resistance, ${\lambda }_m$ is the effective thermal conductivity of the ground, $E_1\left(x\right)={\int }_x^{\infty }\frac{e^{-u}}{u}du$, t* is the normalized time or Fourier number:

$$t^{*}=\frac{{\lambda }_m}{{\left(\rho C_p\right)}_m{r}_b{}^2}t$$

(2)

(ρC_p)_m is the volumetric calorific capacity of the ground [J·K⁻¹·m⁻³] and r_b is the borehole radius. (ρC_p)_m is estimated according to the lithology, for example by using the standard SIA-384/6 [10]. The “classical” interpretation of the TRT allows one to back-calculate:

• The non-disturbed initial temperature of the ground T₀ (°C);
• The thermal conductivity of the ground λ_m [W·K⁻¹·m⁻¹];
• The thermal resistance of the borehole R_b [K·m·W⁻¹].

The first item is obtained based on an inspection of the initial fluid temperature before heating has commenced, while for 2 and 3 above, an approximation of Equation (1) is generally used, with $\gamma$ the Euler constant ($\gamma \approx$ 0.5773):

$$T_f\approx T_0+p{R}_b+\frac{p}{4\pi {\lambda }_m}\left[\ln \left(4t^{*}\right)-\gamma \right]$$

(3)

Equation (3) is valid as soon as a stationary thermal regime in the borehole is reached, i.e., t* > t_min*, which leads to the exclusion of the temperatures measured before t_min of the interpretation. λ_m and R_b can be identified graphically from the slope and intercept of the curve $T_f=f\left(\ln \left(t^{*}\right)\right)$. An alternative is to minimize the root mean square error (RMSE) $\epsilon \left(\left\{X\right\}\right)$ (°C) between the measured temperature T_fl and that calculated by the model T_fl,mod:

$$\epsilon \left(\left\{X\right\}\right)={\left[\frac{1}{t_{max}-t_{min}}{{\displaystyle \int }}_{t_{min}}^{t_{max}}{\left(T_{f,exp}\left(t\right)-T_{f,mod}\left(\left\{X\right\},t\right)\right)}^{2}dt\right]}^{\frac{1}{2}}$$

(4)

where $\left\{X\right\}$ contains the parameters to back-calculate.

Initially developed for BHE, the model described in Equations (1) and (3) do not account for the thermal inertia of the backfilling material. Past work shows that this purely resistive model (Figure 2a) is not suitable for pile heat exchangers of large diameters [11,12], since once λ_m and R_b have been fitted, the temperature change at small time scales (e.g., 1–10 h) is overestimated by several °C, leading to an underestimation of the transfer capacity of the PHE. Moreover, the duration required to thermally load the backfill material t_min* is barely compatible with the operational constraints of a construction site. Indeed, considering the typical criteria t_min* = 5 [13] and assuming a soil with thermal properties ${\lambda }_m$ = 1.4 W·K⁻¹·m⁻¹ and (ρC_p)_m = 2.2 MJ·K⁻¹·m⁻³, this leads to t_min = 14 h for r_b = 8 cm (typical for a BHE) while t_min reaches 196 h (c.a. 8 days) for r_b = 30 cm (typical for a PHE).

Despite these limitations, many authors report application of TRT to PHE to determine the piles and ground properties [11,12]. In these studies, the measured fluid temperature is typically interpreted with the ILS approach, which leads to the exclusion of a significant part of the temperature measurements before t_min* reaches 5. Alternative approaches have been reported in [11,14,15,16]. Loveridge et al. [16] used pre-defined pile step responses, which lead to values of ${\lambda }_m$ and R_b similar to that determined by the ILS, while improving the forecast of the temperature changes at the early stage of the TRT. However, the pile step responses had to be chosen from a small number of representative numerically pre-calculated cases for different PHE geometries. When used for TRT back-calculation, selection from these cases is made before the analysis ([15,17]). It is also not possible to look for the characterization of the heat transfer inside the PHE, since the nature of this is an assumption of the case chosen. This means the approach is not appropriate in all cases and may contain errors related to choice of the closest geometry. Zarrella et al. [14] used the computational Capacity Resistance Model (CaRM) [18] to account for the grouting material inertia. CaRM takes a resistive-capacitive approach and has been developed for specific geometries including for single-U, double-U and coaxial heat exchangers. It has the key limitations that the pipe arrangements in the borehole are predefined (e.g., “pipes close to the borehole wall”) and resistances describing the inner thermal transfer must be pre-calculated with a distinct method, such as finite element code. Once the inner parameters have been tuned, CaRM was in very good agreement with a TRT performed on a BHE of 140 mm diameter. Alberdi et al. compared several analytical approaches, as well as a finite element model representing an horizontal PHE cross section [16]. They showed that models that considered both the short length of the PHE and took account of its internal capacity gave the best results when compared with a benchmark 3D numerical simulation.

However, the main approaches, which account for the pile thermal capacity, have limitations as described above, making them applicable only in certain circumstances, or applied with greater error when generalized. Consequently, there remains a need for a flexible model, which can be specific to a given pile geometry, accounting for concrete thermal inertia and is easy to use and can be embedded in a numerical procedure of back-calculation for TRT interpretation.

In this paper we present a simple model of a PHE with thermal inertia, which deals with PHE of large diameters (r_b $\gg$ 10 cm) and is specific to a given pile geometry (Section 2). The model is tested by back analysis of two experimental TRT data sets, where the PHE parameters are back-calculated and compared with the classical ILS interpretation with Equations (1) and (3) (Section 3). We show that the temperature rise at small time scales (i.e., t ≈ 1 h) can be better accounted for with the new simple approach. Additionally, the expected errors are comparable or less compared with standard borehole TRT interpretation. However, the time to convergence of the new model output means that the soil thermal conductivity cannot be determined from a short TRT on a pile. However, we go on to demonstrate that a novel combination of BHE and PHE TRTs interpreted using the new model can lead to a reduction in pile testing time. On this basis we estimate the minimum TRT duration for PHE (Section 4) potentially saving time on practical operations.

2. Methods

Section 2.1 describes the field test data used in the study and Section 2.2 describes the new resistive capacity model for PHE TRT interpretation. Analysis of the data using the new model and the classical approach introduced in Section 1 is then set out in Section 3.

2.1. Experimental Data

Two TRT, summarized in

Table 1. Characteristics of the TRT.

	Set B	Set C
Depth of pile H [m]	18.3	31.0
Radius of pile rb [m]	0.225	0.300
Geothermal equipment	Double-U (tested as single-U)	Double-U
External diameter of pipes [cm]	3.00	2.50
Thickness of pipes [cm]	0.29	0.23
Distance between two tubes diametrically opposed pipes [m]	0.157	0.425
Initial temperature of the ground T0 [°C]	24.97	14.23
Power applied P [kW]	2.27	1.69
Linear power pf = P/H [W·m−1]	123.7	54.6
Volume flow in the pile [m3·h−1]	2.46	1.15
Duration of the heating [h]	103.4	354.1

, and respectively referred to as set B and set C, are analyzed in this paper. Set B was undertaken on a pile of 22.5 cm radius at Richmond, Texas, USA, mostly located in saturated sand [15,19]. Set C was a TRT undertaken on a pile of 30.0 cm radius located in London clay [12]. The inner details of the piles are given in Table 1 and illustrated in Figure 3. Note that for set B, the power was shut down at t = 103.4 h, and restarted 4 h later until t = 140.2 h. The classical interpretation was carried out for t = 103.4 h. For both TRT, the heat-carrier fluid was water, whose properties are assumed to have a density ρ_fl = 1000 kg·m⁻³ and a heat capacity C_p,fl = 4180 kJ·kg⁻¹·K⁻¹. The temperatures, flow-rate and power were monitored with a time step of 5 min for set B and 1 min for set C. Data of set C were averaged with a time step of 5 min for analysis purposes.

The thermal inertia of the concrete means that the evolution of the fluid temperature goes away from purely linear behavior depending on the logarithm of time foreseen by the ILS model, i.e., Equation (1) (cf. Figure 4). This is even more noticeable for set C, the pile of greater diameter. The temperature measured in set B also exhibits slight periodic variations of 24 h [c.a. 0.5 °C], suggesting the fluid temperature is influenced by the atmosphere. One explanation may be that the pipes connecting the PHE to the TRT module may not have been sufficiently insulated. Note that in what follows, (ρC_p)_m = 2.4 MJ·K⁻¹·m⁻³ is assumed for both soils based on typical values [10].

2.2. Resistive-Capacitive Model for a Pile Heat Exchanger

A model has been developed to take into account the thermal inertia of the backfilling material. The heat transfer inside the concrete relies on a linear capacitance $\mathcal{C}$ [J·K⁻¹·m⁻¹] located between two resistances R₂ and R₃ [K·m·W⁻¹] (Figure 2b). The transfer outside the pile considers a hollow cylinder, and relies on the “classical” infinite cylindrical source (ICS) model [20].

The resistive-capacitive model (RC) in Figure 2b is qualified as semi-analytical since the transfer inside the pile is treated numerically whilst the transfer in the ground is treated analytically. The capacitance is simply estimated by:

$$\mathcal{C}=\pi {\left(\rho C_p\right)}_c{r}_b{}^2$$

(5)

(ρC_p)_c is the volumetric calorific capacity of the concrete [J·K⁻¹·m⁻³], estimated as 2.11 MJ·K⁻¹·m³, a value compatible with preceding studies [21,22]. A parameter x represents the position of the capacitive node (Figure 2b) within the borehole total resistance R_b:

$$x=\frac{R_2}{R_b}; R_b=R_2+R_3$$

(6)

The definition of the fluid mean temperature T_f remains unchanged:

$$T_f=\frac{T_{in}+T_{out}}{2}$$

(7)

A power balance on the fluid leads to:

$$p_f=\frac{P}{H}=\frac{{\left(\rho C_p\right)}_f{Q}_v\left(T_{out}-T_{in}\right)}{H}=\frac{C_p{}_f\dot{m}\left(T_{out}-T_{in}\right)}{H}$$

(8)

where p_f [W·m⁻¹] is the ratio between the power P [W] applied by the TRT test module and the pile depth H [m], $C_p{}_f$ and ${\rho }_f$ the fluid mass-specific heat capacity [J·K⁻¹·kg⁻³] and density [kg·m⁻³], respectively, $Q_v$ [m³·s⁻¹] and $\dot{m}$ [kg·s⁻¹] the volume and mass flow-rate, respectively.

While applied to the concrete:

$$p_f=\frac{T_f-T_{\mathcal{C}}}{R_2}$$

(9)

(ρC_p)_f is the fluid capacity [J·K⁻¹·m⁻³]. A power balance on the central node at temperature $T_{\mathcal{C}}$ leads to:

$$\mathcal{C}\frac{d{T}_{\mathcal{C}}}{dt}=p_f-p_b=\frac{T_f-T_{\mathcal{C}}}{R_2}-\frac{T_{\mathcal{C}}-T_b}{R_3}$$

(10)

In what follows, the equations are discretized with a time step Δt. Δt = 15 min has been chosen to remain in line with the experimental measurements.

An analytical solution to the heat equation, the infinite cylindrical source (ICS) $G\left(t^{*}\right)$ [20], allows the calculation of the temperature at the borehole wall (i.e., at the ground/pile interface):

$$G\left(t^{*}\right)=\frac{1}{{\pi }^2}{\int }_0^{\infty }\frac{e^{-{\beta }^2{t}^{*}}-1}{J_1{}^2\left(\beta \right)+Y_1{}^2\left(\beta \right)}\left(J_0\left(\beta \right)Y_1\left(\beta \right)-J_1\left(\beta \right)Y_0\left(\beta \right)\right)\frac{d\beta }{{\beta }^2}$$

(11)

Under steady-state conditions (constant heat flux), the borehole wall temperature is related to the step response G(t*):

$$T_b=T_0+\frac{p_b}{{\lambda }_m}G\left(t^{*}\right)$$

(12)

During a TRT, the thermal power provided by the test module is constant. However, this power first heats up the concrete before being progressively transferred to the ground. Therefore the steady-state condition cannot be assumed straightaway. The superposition principle is introduced to estimate T_bⁿ at time step n [23]:

$$T_b{}^n=T_0+\{\begin{array}{c}{p}_b{}^1{G}^1 \mathrm{if} n=1\\ \frac{1}{{\lambda }_m}\left(p_b{}^1{G}^n+{\displaystyle \sum _{l=1}^{n-1}}\left(p_b{}^{l+1}-p_b{}^l\right)G^{n-l} \right) \mathrm{if} n>1\end{array}$$

(13)

In Equation (13), the superscripts refer to the time steps, e.g., $T_b{}^n\approx T_b\left(n\Delta t\right)$. p_b is the linear power received at the borehole wall:

$$p_b{}^n=\frac{T_{\mathcal{C}}{}^n-T_p{}^n}{R_3}$$

(14)

The Equations (13) and (14) combine as follows:

$$\left(1+\frac{G^1}{R_3}\right)T_P{}^n=T_0+\frac{p_b{}^{’ n}}{{\lambda }_m}+\frac{T_{\mathcal{C}}{}^n{G}^1}{R_3}$$

(15)

With:

$$p_b{}^{’ n}=\{\begin{array}{c}0 \mathrm{if} n=1\\ p_b{}^1\left(G^2-G^1\right) \mathrm{if} n=2 \\ p_b{}^1{G}^n+{\displaystyle {\displaystyle \sum }_{l=1}^{n-2}}\left(p_b{}^{l+1}-p_b{}^l\right)G^{n-l}-p_b{}^{n-1}{G}^1 \mathrm{if} n>2\end{array}$$

(16)

The Equations (7)–(10) and (15) recombine in a matrix system as follows:

$$\begin{gathered} \left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& {\mathcal{C}}_2& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\frac{d}{dt}\left\{\begin{array}{c}{T}_f\\ T_c\\ T_p\\ T_{in}\\ T_{out}\end{array}\right\}+\left[\begin{array}{ccccc}\frac{1}{R_2}& -\frac{1}{R_2}& 0& 0& 0\\ -\frac{1}{R_2}& \frac{1}{R_2}+\frac{1}{R_3}& -\frac{1}{R_3}& 0& 0\\ 0& 0& -\frac{G^1}{R_3}& \frac{G^1}{R_3}+{\lambda }_m& 0\\ 2{\lambda }_m& 0& 0& -{\lambda }_m& -{\lambda }_m\\ 0& 0& 0& \frac{\dot{m}{C}_{p,f}}{H}& \frac{\dot{m}{C}_{p,f}}{H}\end{array}\right]\left\{\begin{array}{c}{T}_f\\ T_c\\ T_p\\ T_{in}\\ T_{out}\end{array}\right\} \\ =\left\{\begin{array}{c}{p}_f\\ 0\\ {\lambda }_m{T}_0+p_b{}^{\text{'}}\\ 0\\ p_f\left(t\right)\end{array}\right\} \end{gathered}$$

(17)

At each time step n, the vector of the temperatures at the following time step {T}ⁿ⁺¹ is determined by an implicit Euler scheme, i.e., with the following approximation:

$$\frac{\mathrm{d}\left\{T\right\}}{\mathrm{d}t}\approx \frac{\left\{T^{n+1}\right\}-\left\{T^n\right\}}{\Delta t}$$

(18)

With:

$$\left\{T\right\}={\left\{\begin{array}{ccccc}{T}_f& T_c& T_p& T_{in}& T_{out}\end{array}\right\}}^T$$

(19)

leading to the following system solved in Matlab^®:

$$\left(\frac{1}{\Delta t}\left[\mathcal{C}\right]+\left[{\mathsf{\Lambda }}^n\right]\right)\left\{T^{n+1}\right\}=\frac{1}{\Delta t}\left[\mathcal{C}\right]\left\{T^n\right\}+\left\{{𝒫}^{n+1}\right\}$$

(20)

$\left[\mathcal{C}\right]$ is a capacitance matrix [J·m⁻¹·K⁻¹], $\left[{\mathsf{\Lambda }}^n\right]$ a conductance matrix [W·m⁻¹·K⁻¹] and $\left\{{𝒫}^{n+1}\right\}$ a vector [W·m⁻¹] defined as:

$$\left[{\mathsf{\Lambda }}^n\right]=\left[\begin{array}{ccccc}\frac{1}{R_2}& -\frac{1}{R_2}& 0& 0& 0\\ -\frac{1}{R_2}& \frac{1}{R_2}+\frac{1}{R_3}& -\frac{1}{R_3}& 0& 0\\ 0& 0& -\frac{G^1}{R_3}& \frac{G^1}{R_3}+{\lambda }_m& 0\\ 2{\lambda }_m& 0& 0& -{\lambda }_m& -{\lambda }_m\\ 0& 0& 0& \frac{{\dot{m}}^n{C}_{p,f}}{H}& \frac{{\dot{m}}^n{C}_{p,f}}{H}\end{array}\right]$$

(21)

$$\left[\mathcal{C}\right]=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& {\mathcal{C}}_2& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$$

(22)

$$\left\{{𝒫}^n\right\}=\left\{\begin{array}{c}{p}_f{}^n\\ 0\\ {\lambda }_m{T}_0+p_b{}^{’ n}\\ 0\\ p_f{}^n\end{array}\right\}$$

(23)

The typical execution time of the numerical model is a fraction of a second.

Note that the model considers that all temperatures are constant along the pile and neglects the axial heat transfer. However, since the yearly variations of temperatures affect the ground up to a few meters, additional research is still required on how these changes, which result in non-uniform initial ground temperature T₀ measurements, could affect the TRT result.

2.3. Analysis Approach

The initial temperature was determined from the early TRT phases, before heat was injected to the ground. This circulation phase lasted 9 min for set B and 60 min for set C. Averaging the measured temperature yields T₀ = 24.97 °C for set B and 14.23 °C for set C. For set B the value was confirmed to be appropriate with reference to three temperature sensors in a nearby borehole. Less than 0.5 °C variation in ground temperature was observed over 18 m depth.

In this paper, the RMSE (Equation (4)) was minimized with the Matlab^® software. $\left\{X\right\}$ was determined with the fmincon function, using the SQP (sequential quadratic programming) solver. The active-set solver was tested, but appeared to stop before converging to a local minimum of $\epsilon \left(\left\{X\right\}\right)$ in some cases, providing unreliable results. A second stage has the benefit of using the fitnlm and coefCI functions to fit nonlinear regression models and to determine the 95% confidence intervals (CI) on the fitted parameters. The fitnlm function was run with the results of fmincon as an initial point.

The first step was to interpret the TRT with the classical method described in the introduction. Setting the criterion t_min* = 5, the minimum duration t_min was 60 h for set B and 214 h for set C whilst estimating, respectively, λ_m at 2.8 W·K·m⁻¹ and 1.4 W·K·m⁻¹. Indeed, as t_min depends on λ_m which is sought by the interpretation (see Equation (2)), t_min was determined by a manual trial and error process, increasing or decreasing t_min, then minimizing $\epsilon \left(\left\{{\lambda }_m,R_b\right\}\right)$, until the condition on t_min* was reached. The investigation was also undertaken with t_min* = 3. Regarding set C, previous studies [12] reported λ_m between 1.35 and 1.45 W·K·m⁻¹ depending on the values chosen for t_min* (from 5 to 7) and t_max (from 250 h to 350 h) so this new analysis is consistent with previous work.

Having established t_min, the root mean square error between the field data and the modeled data was minimized over all time steps from t_min until t_max, where t_max was the last data time step used in the analysis and was increased in stages until the end of the test period. In this way, the effect of the data window could be investigated and for an ideal data set convergence of derived parameters (R_b, λ_m) would be seen with time.

Secondly, to capture the transient phase at the beginning of the TRT, the resistive-capacitive model presented in Section 2.2 was used to back-calculate the ground thermal conductivity λ_m and the borehole resistances R_b and x parameter by minimizing the misfit Equation (4), i.e., {X} = {λ_m, R_b, x} with the fitnlm function. Note that the first optimization with fmincon was run with {X} = {λ_m, R₂, R₃}. Without the limitations of the stationary assumption for the PHE in the classical method, t_min was fixed to 1 h for both datasets, to capture the transient heat transfer in the pile.

3. Results

3.1. Classical Interpretation

As t_max is increased the root mean square error, thermal conductivity and thermal resistance should converge as the model error is reduced as the PHE thermal capacity is overcome with time. Set B (Figure 5a,c,e) presents a different behavior since the values of λ_m and R_b do not converge when the duration of the test increases. There are two factors, which may contribute to explaining this result. First, the cyclic daily perturbation of the measured temperature between the pile and the near surface interfere with the interpretation, resulting in a 24 h cyclic pattern in both λ_m and R_b. It has been observed elsewhere that for TRTs performed on piles, effects like this can lead to increasing errors since any perturbations in temperature due to the environment have a proportionally greater impact later in the tests when the rate of absolute temperature change is lower [24].

Secondly, it has been reported [24] after t* = 10 the pile thermal response will be affected by end effects leading to systematic overestimation of thermal conductivity. For set B, t* was 9.9 at the end of test, so this effect could also impact the late time accuracy.

However, the variation of λ_m (between 2.7 and 3.2 W·K⁻¹·m⁻¹, that is 15%) remains acceptable in as much as the uncertainty in the interpretation of a TRT is of the order of 10% [23]. The misfit soars for t_max = 105 h since it encompasses temperatures measured after the power shutdown. The daily perturbation of temperature results in narrower intervals of confidence.

For set C, λ_m and R_b converge if the duration of the TRT t_max is above approximately 220 h for t_min* = 3 (cf. b, d, f). λ_m and R_b stabilize, respectively, at 1.44 W·K⁻¹·m⁻¹ and at 0.136 K·m·W⁻¹. The interpretation is similar with t_min* = 5 as soon as t_max = 260 h: using t_min* = 3 allows the reduction of the duration of the TRT by only 10%. Beyond this threshold at ≈260 h, λ_m and R_b vary by less than 3% if t_max increases: the thermal transfer at the borehole wall is stationary, which validates the resistive model. t_min was estimated to ≈220 h for λ_m ≈ 1.4 W·K⁻¹·m⁻¹, which means that about 84% (≈(260 − 220)/260) of the TRT duration only served to thermally charge the concrete—and only 40 h was actually used for the interpretation.

The dependence of λ_m with t_max has been reported in [24] for a TRT exhibiting daily perturbations performed on a 25 m deep BHE, as for set B. Set C does not exhibit such a behavior. One explanation may be that set C PHE is 70% deeper than set B PHE, while λ_m is half, resulting in a shallower affected thermal zone. A more realistic explanation is that set C measurements do not need to be corrected since specific care has been taken to insulate the pipes. For set B, the perturbation of the fluid temperature, estimated as a deviation from the fitted Jacob approximation of the ILS (Equation (3)), can reach c.a. 0.4 °C, while it remains below 0.1 °C for set A.

The engineer sizing a ground-source heat pump is concerned with the ability of the models representing every system component to forecast the evolution of the fluid temperature at small, medium and large time scales. The capability of the purely resistive model coupled to the ILS (Equation (1)) or the ICS (Equation (11)) has been assessed for both experimental datasets (Figure 6). The TRNSYS Type 557, which is based on the duct storage model (DST) for ground heat exchangers [25], has been included in the comparison as well. The DST uses superposition of three elements: a heat balance on the fluid, a local resistive process close to an individual pipe (so-called “duct”) rather than the borehole, and a global process.

All the models use the same values of ground conductivity λ_m and borehole resistance R_b, obtained by the “classical” interpretation with t_min* = 5 for both sets, t_max = 103 h for set B and t_max = 350 h for set C. Equation (3) is also represented as it is used for the interpretation, though it may lead to a significant underestimation of the temperature at small time scales due to its logarithmic behavior. After 1 h of thermal load the ILS model (Equation (1)) overestimates the temperature of the fluid T_fl by ≈5.0 °C for set B and by ≈3.3 °C for set C. For both datasets, neglecting the thermal inertia of the concrete leads to an underestimation of the pile heat transfer capability, and therefore possibly a cost-ineffective geothermal equipment of many more piles than required. Indeed, the heat pump could exchange a larger amount of thermal energy with the underground than a purely resistive model would forecast.

When parameterizing the duct storage model, the resistance R_b was taken to be equivalent to the “fluid to ground resistance” and for simplicity no account was taken of any additional pipe-to-pipe interactions. The DST is remarkably able to reproduce the temperature evolution. This is, to some extent, surprising since DST is a purely resistive model in terms of internal heat transfer and it does not allow defining the pipe locations in the borehole or pile [25]. However, DST has been primarily designed for borehole thermal energy storage. Consequently, the ground heat exchangers are assumed to be located within a cylinder at the nodes of a hexagonal grid, which will not be realistic for most PHE fields. Meanwhile, thanks to the superimposition principle [21], analytic solutions are far more flexible to handle arbitrary borehole locations. A framework can even be built to select which piles to equip in a field [8].

3.2. Interpretation with the New Resistive-Capacitive (RC) Model

Figure 7 shows the back-analyzed thermal conductivity for the two datasets using the RC model coupled to the infinite cylindrical source. For set B, the conductivity oscillates in the range 3.1–3.3 W·K⁻¹·m⁻¹ as long as t < 103.4 h, so before the power shutdown (Figure 7c). λ_m exhibits smaller variations than the classical interpretation, suggesting that the new approach is better capable of dealing with daily perturbations in the measured temperature since it uses temporal superposition. However, though the RC model has the ability to deal with time-varying temperature changes, it fails to properly estimate the ground conductivity when encompassing data after the power shutdown, λ_m then soaring to the range 3.7–3.9 W·K⁻¹·m⁻¹.

For set C, as soon as t_max > 250 h, the conductivity obtained by the inversion of the RC model converges towards λ_m = 1.43 W·K⁻¹·m⁻¹, a value almost identical to that obtained by the normal interpretation (Figure 7d). The convergence of R_b and x beyond 250 h can also be observed (Figure 7f). Hence, while the RC model appears to allow successful interpretation of set B, which was not possible with the classical approach, the RC interpretation does not allow the reduction of the duration of the TRT compared to the resistive model in this case.

3.2.1. Error Analysis

An in-depth analysis of the error on the back-calculated parameters has been performed by Witte for a conventional TRT on a BHE, given that analytical expressions of λ_m and R_b are known [23]. For comparison, the estimated error $\delta X_i$ in the present study on the three back-calculated parameters $X_i$ is determined as follows:

$$\delta X_i=\sqrt{{\displaystyle \sum }_j^{}{\left(\frac{\partial X_i}{\partial y_j}\mathsf{\Delta }{y}_j\right)}^2}$$

(24)

where $y_j$ account for the input parameters and $\mathsf{\Delta }{y}_j$ for the related error on this parameter. The partial derivatives $\frac{\partial X_i}{\partial y_j}$ are computed numerically by changing the input value by a small amount ($\pm$1%), then minimizing the RMSE to get the updated value of $X_i$, and computing the relative change of $X_i$. The measurement error has been computed for a heat injection duration of 100 h for set B and 350 h for set C. The values of $\mathsf{\Delta }{y}_j$ are given in

Table 2. Considered values for the error on the input parameters.

Input Parameter	Error Values
	Set B	Set C
Ground heat capacity [MJ·K−1·m−3]	0.4	0.3
Concrete heat capacity [MJ·K−1·m−3]	0.2	0.2
Pile diameter [m]	0.025	0.03
Height of the equipped pile [m]	0.05	0.05
Initial ground temperature [K]	0.3	0.1
Fluid heat capacity [J·K−1·kg−1]	1.0	9.5
Power [W]	43.7	5.0
Flow rate [m3·s−1]	3.44 × 10−5	3.29 × 10−6
Inlet temperature [K]	0.3	0.1
Outlet temperature [K]	0.3	0.1

, with the following comments:

• According to SIA standards, the volume-specific heat capacities of wet clay and wet sand are, respectively, in the range 2.0–2.8 MJ·K⁻¹·m⁻³ and 2.2–2.8 MJ·K⁻¹·m⁻³ [10]. We considered the error on the ground capacity as the half of these intervals, i.e., 0.3 MJ·K⁻¹·m⁻³ for set B and 0.4 MJ·K⁻¹·m⁻³ for set C. Typical ranges for concrete heat capacity could not be found in the literature, but given the values reported in previous studies [26,27], an error of 0.2 MJ·K⁻¹·m⁻³ was chosen.
• The error for the pile diameter and height were determined according to the UK specification for construction tolerances [28]. In this respect, it should be noted that the dimensions of a constructed pile should not be less than the specified dimensions. A tolerance on these dimensions of up to the lesser of 50 mm or 5% is permissible.
• For set B, the test was performed with reference to the ASHRAE standard [29]. This states that the accuracy of temperature measurement must be less than 0.3 °C, for power measurements less than 2% and for flow rate measurements less than 5%. These are conservative values, since the test may have been performed with more accurate instruments.
• For set C, the client specification had tighter accuracy requirements, which can be reasonably applied. They would be error for temperature of 0.1 °C, flow measurement of 0.01 m·s⁻¹ and power to 5W.

For both datasets, the most significant sources of error are related to the heat capacities of the ground and the concrete (refer to Figure 8 which represent the distribution of $\frac{\partial X_i}{\partial y_j}\mathsf{\Delta }{y}_j$ terms). As expected, temperatures and power play a more significant role for set B, where the standards are less restrictive. Errors on the pile diameter affect the inner parameters R_b and x significantly. The error on the ground thermal conductivity is 0.06 W·K⁻¹·m⁻¹ for set C (about 4%) and 0.20 W·K⁻¹·m⁻¹ for set B (about 6%) reflecting the more restrictive standards for set C. These values are well within the expected bounds for thermal response testing for boreholes, giving confidence to the new approach. For comparison, for a classical TRT on a BHE, Witte reported a typical error of 5% on the ground thermal conductivity [23].

3.2.2. Pile Thermal Capacity

The interest of the RC model on the resistive model is in the transitory phase. Here we arbitrarily define this phase as t* < 2.5, which results in 27 h for set B and 107 h for set C (Equation (2)). Once their inner parameters have been fitted, the classical and RC methods are compared in Figure 6. Given the similar results between the two methods in terms of conductivity prediction, for comparative purposes, the same input parameters are used in Figure 6. Note that both methods also use the normalized duration t_max* of experimental data, so t_max = 103 h for set B (t_max* ≈ 9.9) and t_max = 350 h for set C (t_max* ≈ 8.5). After 1 h of thermal load, the RC model overestimates the temperature of the fluid T_fl by only ≈1.0 °C for set B and ≈1.3 °C for set C. Beyond the transitory phase (t* > 2.5), the two approaches accord well with the measured temperature and both correctly forecast the measured temperature, especially when the t* > 5 criteria is reached. Further, the RC model allows an understanding of the dynamics of heat transfer within the pile (Figure 9). The pile is half-loaded at t* ≈ 0.4 for set B and t* ≈ 0.2 for set C.

4. Discussion and Recommendations

Whatever the interpretation method, the duration t_max for set C must be of the order of 250 h, i.e., about 10 days. This duration is not necessarily compatible with the management of a construction site. Therefore it is clear that the ground thermal conductivity may be obtained from a much shorter TRT carried out on a special purpose BHE as recommended by [17]. Yet this proposed approach retains a drawback, namely that the pile itself has not been characterized to allow appropriate forward simulation for the design process. No information has been obtained about its thermal resistance. While a closed form analytical expression has now been developed to obtain pile thermal resistance [30] this still requires information about concrete thermal conductivity and assumes accuracy of pipe placements compared with design.

Therefore we propose that by using the RC model the combination of a short BHE TRT and an additional short TRT on one of the constructed piles may provide the best solution. One TRT is performed on a BHE having a small radius (e.g., r_b = 8 cm) to determine λ_m within a few days, while simultaneously a TRT is performed on the PHE to obtain the inner parameters for the RC model. This undertaking assumes that the ground is relatively homogeneous and lateral variations of composition and ground properties can be overlooked. This approach is illustrated with respect to datasets B and C. The best-fit values of λ_m are assumed to have been obtained independently from a BHE TRT. In this case values of 3.24 W/mK and 1.43 W/mK are used, respectively, for sets B and C, based on the results in Section 3. The RC model is then optimized for the pile TRT data with only two rather than three unknowns. The fitted values of x and R_b start to stabilize around t_max ≈ 30–40 h for set B at x ≈ 0.56–0.58 and R_b = 0.080 K·m·W⁻¹ and t_max ≈ 100 h for set C at x ≈ 0.76–0.78 and R_b = 0.121–0.124 K·m·W⁻¹ (Figure 10). For both datasets, the difference, with the value of R_b determined in Section 3, does not exceed 4%.

It is desirable to determine the minimum test duration t_max beyond which the determined values of R_b and x are reliable. Here we use as a criteria the ability of the RC model to forecast the evolution of the fluid temperature at the very end of the data record, from t_min’ to from t_max’. A second RMSE indicator ε’ is build considering t_min’ = 80 h and t_max’ = 100 h for set B, and t_min’ = 300 h and t_max’ = 350 h for set C. For set C, ε’ reaches a minimum around t_max* = 100 h (t_max* = 2.4), which suggests this duration is sufficient to estimate the temperature evolution from 300 to 350 h, as confirmed by the comparison to the recorded temperature (Figure 11).

5. Conclusions

A resistive-capacitive PHE model has been developed to capture the transient energy transfers within a pile. This model has been shown to perform well and in particular:

• Numerical back-calculation of the model parameters on two thermal response tests yield similar values of ground conductivity and thermal resistance as the well-established infinite line source model.
• Inclusion of temporal superposition with the model allows reliable results to be obtained even when tests are affected by ambient air interference.
• The RC model better represents the transient phase of pile warm-up in the early part of the test (approximately up to a Fourier number t* = 1 to 2).
• The errors associated with the calculation of thermal conductivity are all less than 10% and well within expected ranges for boreholes thermal response tests interpreted with the classic infinite line source.
• Standard back-calculation using the RC model does not allow to significantly reduce the TRT duration below t* = 5.
• However, if the thermal conductivity can be obtained by another means, the time for the RC model to converge is much reduced meaning that pile resistance can be obtained from a pile TRT in a duration corresponding to a Fourier number t* ≈ 2 to 2.5.

Given these characteristics of the RC model the following novel approach for pile characterization is recommended:

• Use a borehole at the same site and of the same length as the piles to carry out a BHE TRT to determine the effective soil thermal conductivity using the classical approach.
• Carry out a short duration pile TRT according to Fourier number t* ≈ 2 to 2.5, or around 100 h for the cases demonstrated in this paper.
• Interpret the pile TRT using the RC model to determine both the pile thermal resistance and the inner resistances of the RC model, which can then be used in forward simulation for design purposes.