Parameter identification algorithm for ground source heat pump systems

• This paper presents a new parameter identification (PI) algorithm for estimating effective and detailed thermal parameters of ground source heat pump systems using data obtained from the well-known thermal response test. The PI comprises an iterative scheme coupling a semi-analytical forward model to an inverse model. The forward model is formulated based on the spectral element method to simulate transient 3D heat flow in ground source heat pump (GSHP) systems, and the inverse model is formulated based on the interior-point optimization method to minimize the system objective function. Compared to existing interpretation tools for the thermal response test, the proposed PI algorithm has several advanced features, including: it can handle fluctuating heat pump power and inlet temperatures; interpret data obtained from multiple heat injection or extraction signals; produce accurate backcalculation for short and long duration experiments; and handle multilayer systems. The PI algorithm is tested against synthesized data, using a wide range of random noise, and versus an available laboratory experiment. The computational results show that the PI algorithm is accurate, stable and exhibiting relatively high convergence rate.


Introduction
The use of the ground source heat pump technology for heating and cooling of buildings is rapidly rising worldwide, and engineers are striving to improve its design. Efficient design of a GSHP system depends in part on the accuracy of the thermal parameters of the borehole heat exchanger and the soil mass [1], and in another on the optimization of the energetic, exergetic and economic performance of the system as a whole [2][3][4]. This paper focuses on the first part.
Thermal parameters of the borehole heat exchanger (BHE) components are usually known a priori, but thermal parameters of the soil mass and the borehole thermal resistance are not readily known and need to be determined. Laboratory and field experiments have been devised for estimating the soil thermal parameters. Thermal properties estimated from laboratory experiments are in principle more accurate.
However, the drawback in here is that only a small volume of the soil is tested, ignoring the inhomogeneity of the soil layers. Also, disturbances in the soil samples can lead to erroneous estimate [1]. Following this, the in-situ Thermal Response Test (TRT), first introduced by Mogensen [5], has become the tool for estimating the GSHP parameters [6]. The TRT became popular because in addition to determining the thermal properties of the soil without the need for taking samples to the lab, it determines the borehole thermal interaction coefficient, a property that is quite difficult to be quantified in the lab.
A typical TRT consists of a vertical borehole heat exchanger embedded in which a single U-tube, and equipped with thermometers to measure the fluid temperatures at the inlet and outlet of the U-tube. Modern TRT devises are equipped with fiber optics to measure the fluid temperature along the BHE [7]. Spitler and Gehlin [8] provided a historical review of the thermal response test, and Zhang et al. [9] and Witte [10] provided a comprehensive review on the state-of-the-art of TRT and its parameter identification algorithms.
An important element of the thermal response test is the interpretation procedure of the measured data. This experiment is relatively expensive but can only be useful if the utilized parameter identification algorithm is able to produce accurate estimate of the involved thermal parameters. Basically, parameter identification algorithms employ forward and inverse models to obtain effective or detailed parameters, depending on the rigor of the algorithm. The forward model is a mathematical expression which describes the physical system and predicts the response of the system for any given values of the model parameters. The inverse model is a mathematical optimization algorithm designed for inferring the values of the model parameters from measured and theoretical response quantities. The competence of a parameter identification algorithm depends on several factors, including: 1. The adequacy of the forward model to describe the physics of the problem, and its computational efficiency for utilization in iterative schemes for inverse calculations. 2. The consistency, stability, convergence rate and uniqueness of the inverse model for solving the objective functions.
Based on these factors, parameter identification algorithms utilized for TRT interpretation can be put in three categories: 1. Forward analytical -Inverse graphical interpretation (GI) algorithm (also denoted as curve fitting method). This algorithm employs the well-known infinite line source (ILS) model [11] for the forward calculation. The inverse calculation is conducted via curve fitting of the measured data to estimate the effective soil thermal conductivity and the borehole thermal resistance, as shown in Fig. 1. This method is widely used because of its simplicity; simple mathematical tools (such as MS Excel) can be used to interpret the measured data. Works fall in this category, among many others, are: Witte et al. [1], Busso et al. [12], Sanner et al. [13] and [14], Witte [15], Focaccia et al. [16], and Bujok et al. [17]. In Section 3 this algorithm is discussed in detail, and its performance is examined in Section 5.
The study shows that this algorithm has fundamental shortcomings, mainly in its non-uniqueness in determining two physical parameters (soil thermal conductivity and borehole thermal resistance) from one equation. 2. Forward analytical -Inverse optimization algorithms. This kind of algorithms employ analytical models such as the infinite line source (ILS), infinite cylindrical source (ICS) [11], thermal resistance and capacitance (TRC) [18], or alike, for the forward calculation. The inverse calculation is conducted using some sort of optimization techniques, as shown in Fig. 2. Works fall in this category, among others, are: Fujii et al [19], Li and Lai [20], Wagner et al. [21], Raymond and Lamarche [22], Choi and Ooka [23], Pasquier [24] and [25], Zhang et al. [26], Choi et al. [27], and Pasquier et al. [28].
In this category, the use of optimization techniques for minimizing The BHE pipe-out g The BHE grout s The soil film soil The soil mass the system objective functions overcomes some of the shortcomings of the previous category. However, most of the analytical forward models are not formulated to describe the detailed heat flow in the GSHP systems, making backcalculations based on such models limited in accuracy. 3. Forward numerical -Inverse optimization algorithms. This kind of algorithms employ numerical models based on the finite element method, finite difference method or finite volume method for the forward calculation. The inverse calculation is conducted using some sort of optimization techniques, as shown in Fig. 3. Works fall in this category, among others, are: Spitler et al. [29], Austin et al. [30], Wagner and Clauser [31], Signorelli et al. [32], and Marcotte and Pasquier [33]. In general, the parameter identification models of this category are more advanced, but the use of the numerical models for forward calculations make them computationally inefficient for utilization in iterative schemes for inverse problems.
In this paper, we adopt the second category (Forward analytical -Inverse optimization algorithm). A new parameter identification (PI) algorithm capable of estimating effective and detailed thermal parameters of GSHP systems is introduced. The PI comprises an iterative scheme, coupling a semi-analytical forward model to an inverse model.
The forward model is formulated based on the spectral element method [34] and the superposition principle to simulate transient 3D heat flow in multiple borehole heat exchangers embedded in multilayer systems. The spectral element method (SEM) solves linear partial differential equations for a homogeneous domain analytically, and for a layered domain it solves the governing equations by means of the finite element technique. It elegantly combines the exactness of the analytical methods to a great extent of generality of the numerical methods in describing the geometry and boundary conditions. It requires one element per layer, making the mesh size equivalent to the number of layers. These features make the proposed forward model accurate and computationally efficient, and thus suitable for utilization in iterative schemes for inverse problems.
The inverse model is formulated based on the interior-point method (IPM), a technique capable of optimizing multi-dimensional, sparse linear, quadratic or general nonlinear objective functions and constraints [35]. Optimization algorithms based on the IPM are usually implemented on the basis of the predictor-corrector technique, Cholesky decomposition and Newton's method. These features make the inverse model capable of identifying multiple variables problems with relatively high convergence rate.
Details of the PI algorithm is given in Section 4, followed by verifications and analyses in Section 5 and 6. But first, we elaborate on the challenges involved in interpreting the TRT measurements (Section 2) and discuss the current GI algorithm (Section 3).

TRT interpretation challenges
Despite the bulk of research work in this field, there are yet many challenges that affect the accuracy of the TRT interpretation procedures, including: A) Heat pump power fluctuation Heat pumps used in GSHP systems usually produce fluctuating thermal loads, leading eventually to fluctuations in the measured inlet and outlet temperatures. As it will be shown in Section 3 and the subsequent sections, the GI algorithm is directly related to the heat power, and if this power is not constant, the temperature can be misrepresented, giving rise to inaccurate curve fitting. Two solutions have been employed for this problem. One is modifying the hardware by using electric resistance heaters for heating the circulating fluid [13,36]; and another is modifying the interpretation method by using numerical forward models which are usually capable of handling fluctuated heat fluxes. Witte and van Gelder [37], for instance, employed TRNSYS simulator to estimate the soil thermal parameters using a time varying heat flux injection. Raymond et al. [38] and Choi and Ooka [23] tackled this problem by using a temporal superposition technique combining the ILS model to an optimization algorithm.

D) TRT duration
In practice, the TRT is conducted for 72 h in an effort to let the system reaches the steady state and guarantee better interpretation [37]. This duration is particularly crucial for the GI algorithm, as will be explained in Section 3 and verified in Section 5.1. However, performing the experiment for 72 h can be expensive. Spitler et al. [29] have shown that the cost of performing TRT increases significantly with the increase of the TRT duration. As a consequence, recently, several studies are initiated to investigate the possibility of shortening the TRT duration. Bujok et al. [17] performed several TRT in a field containing 16 boreholes to analyze the effect of shortening the TRT durations on the accuracy of the thermal parameter identification. Their study shows that shortening the test to 24 h can bring an acceptable amount of inaccuracy. Pasquier [25] introduced an interpretation technique that can make use of the measured temperatures during the first three hours to produce an effective thermal conductivity accurate within 10% of a reference value.

E) Nonhomogeneous soil mass
Most TRT parameter identification algorithms assume that the soil mass consists of a homogeneous layer, which can be described by effective soil thermal conductivity and borehole thermal resistance. Obviously, this assumption is not realistic as the borehole can go as deep as 200 m into the ground, which most probably consists of layers with varying soil types, geometry and thermal properties [37,8]. An attempt to backcalculate the thermal properties of a multilayer system necessities two main requirements: 1) The TRT needs to include temperature measurements, not only at the inlet and outlet of the borehole, but also along its depth. Currently, this is gradually becoming a common practice, especially by the introduction of what is known as the Distributed Thermal Response Test (DTRT), where measurements are made along the BHE using optical fiber [19,7,39].
2) The forward model needs to be able to simulate multilayer systems.
Raymond and Lamarche [22] presented a noteworthy attempt to estimate the thermal parameters of a multilayer soil mass system. They utilized the MLU simulator (an analytical tool for hydraulic pumping test in layered domains) for the forward calculation. The inverse calculation is conducted using the Levenberg-Marquardt algorithm [40]. They backcalculated thermal parameters for three layers and showed that the results were, for most cases, within 20% accuracy.
As indicated above, these challenges have been addressed in the literature, but most of the models can only deal with one or two challenges at a time. In this paper we address these challenges, together with the possibility to identify parameters other than the conventional thermal conductivity and borehole thermal resistance.

Current graphical interpretation (GI) algorithm
The GI algorithm estimates the effective soil thermal conductivity and borehole thermal resistance based on the infinite line source model (ILS) and the least square curve fitting method. For a sufficiently long time, Carslaw and Jaegar [41] solved the ILS partial differential equation for a constant heat flow rate per unit length, q, giving: is the temperature of the medium (soil mass in this case) at radial distance r from the heat line source, (W/(m.K)) is the medium thermal conductivity, (m 2 /s) is the thermal diffusivity and = 0.5772 is Euler's constant. At the borehole boundary surface, = r r b , the temperature is [42], and [31]: are the inlet and outlet temperatures, and R b is the borehole thermal resistance.
Substituting Eq. (1) into Eq. (2), and rearranging, gives In a compact format, Eq. (4) can be written as which represents a line in a semi-log scale with a its slope and b its T -intercept. Relating a to the coefficient of the first term on the righthand side of Eq. (4), and knowing the heat pump q, the soil thermal conductivity can readily be determined. Then, relating b to the second and third terms on the right-hand side of Eq. (4), and known from the previous step together with all other known parameters, R b can be determined. Apparently, the GI algorithm is easy to implement, making it attractive for daily engineering practice. However, it suffers from fundamental shortcomings, including: 1. This method is undetermined. It entails determining two physical parameters from one equation in two steps: not by solving two equations simultaneously. An error in determining would inevitably lead to error in R b . Basically, the parameter identification procedure requires either determinate systems (number of equations is equal to the number of unknowns) or over-determinate systems (number of equations is more than the number of unknowns). 2. The time when the curve fitting begins can have a crucial effect on the values of a and b, and hence on and R b . 3. The solution given in Eq. (1), provided by Carslaw and Jaegar [41], is applicable for "sufficiently" long time, necessitating conducting the experiment for several days to be valid.

Proposed PI algorithm
The proposed PI algorithm consists of a forward model and an inverse model, described hereafter.

Forward model
The background theory of the forward model has been thoroughly presented in a string of publications, mainly Al-Khoury [43], BniLam and Al-Khoury [44][45][46] and BniLam et al. [47]. Here, we present the governing heat equations, the associated thermal resistance between GSHP components and their correspondence to the commonly used formulation for the borehole thermal resistance.
The governing heat equations of a ground source heat pump consisting of a single U-tube borehole heat exchanger embedded in a soil mass, Fig. 4, can be described as Pipe-in Pipe-out Grout Soil mass The soil film (Eq. (9)) is introduced in this model to couple the borehole heat exchanger to the soil mass. The essence of this coupling has been discussed in BniLam and Al-Khoury [45] and BniLam et al. [47].
The initial condition is The boundary condition at the inlet of pipe-in might be any of two types: a Neumann boundary condition: (12) or a Dirichlet boundary condition: (13) where q in is the prescribed heat flux and T in is the prescribed inlet temperature, that might have any arbitrary distribution in time.
Eqs. (6)-(9) are coupled via local thermal interaction terms describing heat flow at the contact surfaces between neighboring BHE components. At the boundary between pipe-in and grout, the heat flow is described in terms of the thermal interaction coefficient, b ig , expressed as where dS ig is the surface area of pipe-in in contact with the grout, and R ig is the corresponding thermal resistance, expressed analytically as in which r pi and r po are the inner and outer radius of pipe-in, respectively; L is the length of the spectral element; p is the thermal conductivity of pipe-in material; and = h D Nu is the convective heat transfer coefficient, where D is the inner diameter of the pipe, Nu and are the Nusselt number and thermal conductivity of the circulating fluid. A similar formulation is used for b og : pipe out -grout thermal interaction.
The thermal interaction coefficient for the soil film -soil mass is expressed as where dS s is the surface area of the soil film in contact with the soil mass, and R ss is the corresponding thermal resistance, expressed analytically as in which r f is the radius of the soil film and r b is the borehole radius. The thermal interaction coefficient for grout -soil film is expressed as where dS gs is the surface area of grout in contact with the soil film, and R gs is the corresponding thermal resistance. Unlike R ig , R og and R ss , R gs does not have an analytical expression due to the acentric positions of pipe-in and pipe-out inside the grout. This coefficient is usually approximated using an equivalent centric pipe with an equivalent crosssectional area of pipe-in and pipe-out [43]. Here, R gs is backcalculated using the inverse model.
The above local thermal interaction terms are unique to this model, Eqs. (6)- (9). Different models can have different local thermal interaction arrangements and formulations. Nevertheless, regardless of the local interaction terms, models that are utilized to interpret the thermal response test (TRT) are required to back calculate the borehole thermal resistance, R b . The usual practice is to lump the local thermal interaction coefficients into an equivalent coefficient using some sort of arrangements, such as the Delta-or Y -configuration analogy to Ohm's law [43]. The governing BHE heat equations (Eqs. (6)-(9)) are a typical Y -configuration, where pipe-in and pipe-out are interacting with each other via the grout, which is in contact with the soil, Fig. 5a. Fig. 5b shows the thermal resistance configuration, where R ig and R og are in parallel to each other and in series with R gs and R ss . The equivalent thermal resistance for this configuration is given in Fig. 5c that can be expressed as where R ig , R og and R ss are determined analytically using Eqs. (15) and (17), respectively, and R gs is backcalculated.

Inverse model
The procedure for parameter identification entails minimizing the objective function of a system subjected to predefined constraints, such that: where f x ( ) is the objective function that needs to be minimized; is the vector of the unknown parameters that need to be estimated, such as the soil thermal conductivity and the borehole thermal resistance; and represent the constraints of x.
For a ground source heat pump, the objective function can be expressed as the Euclidean distance (norm 2) between the measured and computed temperatures. Using the spectral analysis, the Euclidean distance is described in the frequency domain, as can be at any frequency, n , at any depth, z, and for any component of the GSHP system, including the soil mass. However, in typical TRT, only the temperatures at the inlet of pipe-in and outlet of pipe-out are measured. The summation over n in Eq. (21) indicates that the objective function can contain any number of frequencies. This constitutes an important advantage of the spectral analysis compared to the time domain. In the time domain the whole signal must be considered and the whole system must be calculated in every iteration, while in the spectral analysis the system can be calculated for only few frequencies. This stems from the fundamental property of the frequency domain which describes the characteristic of the system to respond to a range of frequencies regardless of the time of occurrences. In the time domain, the sequence of occurrences determines the response of the physical system, specifically for that event. The verification and numerical examples given in later sections are conducted using 20 frequencies for each reading set. For more information regarding the frequency response of a signal, the reader is referred to Doyle [34].
Solving the objective function, Eq. (21), requires a minimization algorithm capable of solving complex, multidimensional, nonlinear equations. The interior-point optimization method [48], which can handle large and sparse problems, is employed for this purpose. It satisfies bounds at all iterations and can recover from NAN (Not a Number) or Inf (Infinity) results. This algorithm has been implemented and optimized in the MATLAB nonlinear programming solver called fmincon [49].
The proposed PI algorithm can be summarized as

Performance of PI algorithm
The performance of the proposed PI algorithm is examined by means of numerical experiments using synthetic data generated by the forward model to represent measurements of a typical TRT experiment. Synthetic data is usually smooth, but to make it comparable to measured data, random noise with different contamination levels is applied. The advantage of using synthetic data is that the input parameters are known a priori that enables quantifiable verification of the inverse model. The PI algorithm is tested in terms of its accuracy; stability, which is a measure of the ability of an algorithm to converge even in the presence of contaminated data; and convergence rate, which is a measure of the speed at which a convergence sequence reaches the desired accuracy.
Typical TRT heat pumps exhibit fluctuations in their power signal, with standard deviation (std) ranging between 2% and 10% [50] and [1]. To mimic this, we contaminate the synthetic heat pump power using an additive white Gaussian noise (AWGN), such that = + P P t ( ,std) con sy (22) in which P con is the contaminated heat pump signal, P sy is the synthetic heat pump power and is a zero-mean random AWGN with standard deviation, std. Here, the synthetic heat pump power is assumed 5000 W.
To generate a wide range of randomness, 150 random AWGN with three levels of std are applied: The heat pump is assumed to be connected to a 100 m BHE, constituting a single U-tube embedded in a soil mass. Details of the BHE parameters are given in Table 1.
Two numerical TRT experiments have been conducted: 1) for a BHE embedded in a half space; and 2) for a BHE embedded in a three-layer system. Parameter identifications have been performed on the contaminated data using the GI method and the proposed PI algorithm. The performance of these two algorithms is compared based on: 1. accuracy in predicting the parameters; 2. noise effect; 3. duration of experiment; and 4. number of backcalculated parameters.

Numerical TRT in a half space
In this set up, the BHE is assumed to be embedded in a half space constituting a homogeneous soil layer with thermal conductivity = 2 s W/(m.K), and volumetric heat capacity = c 2.6 s s MJ/(m 3 .K). 150 numerical TRT experiments were conducted for 72 h. They were divided into three sets, each representing 50 numerical experiments subjected to randomly contaminated power signal (Eq. (22)) with a specific standard deviation (Eq. (23)). Fig. 6 shows an example of one of the three sets. It shows randomly contaminated power with 300 W standard deviation (Fig. 6a), and the associated computed inlet and outlet temperatures (Fig. 6b). This figure evidently demonstrates that the use of AWGN makes the synthetic data comparable to measured data.
The borehole thermal resistance is calculated based on the forward model results, via [22]: where T is the average fluid temperature (Eq. (3)), and T b the borehole temperature (T s in Eq. (8)). Figure (6c) is a plot of R b versus time, which clearly shows that, upon the end of the transient period, R b becomes nearly constant at 0.2 m.K/W.

Backcalculation based on GI algorithm
The GI algorithm is utilized to backcalculate the effective soil thermal conductivity, s , and the borehole thermal resistance, R b , for the 150 numerical experiments. Fig. 7 shows an example of one of the numerical experiments, generated from processing the computed results of Fig. 6. It displays the average fluid temperature T versus time, along with the best fitted curves. Four fitted curves, depending on the start time, were conducted: after 1 h; after 6 h; after 15 h; and after 30 h. The fitted curves are shown in normal scale plots and in semi-log scale plots (inside the blue boxes). Obviously, the fitted curve which starts after 1 h is the least accurate due to the high transient gradient at the beginning of the experiment. m.K/W) the figure qualitatively indicates that the fitted curve which starts after 30 h and = std 100W (2% std) gave the best estimate. The boxplots is utilized as a means for descriptive statistics to observe the estimation errors of the best fitted curve case based on five quantities: minimum, first quartile, median, third quartile, and maximum. Fig. 9 shows boxplots of the 30 h fitting case; and Table 2    These observations indicate that the GI algorithm depends strongly on the noise level.
It is worth noting that the error in the backcalculated borehole thermal resistance has the same statistical description as that for the

N. BniLam and R. Al-Khoury
Applied Energy 264 (2020) 114712 8 thermal conductivity (see Fig. 9a and b). This proves the argument which was raised in Section 3. As the estimation of R b using the GI algorithm is dependent on s , the error in estimating s is reflected on R b .

Backcalculation based on PI algorithm
Here, the proposed PI algorithm was utilized to backcalculate s and R b from the 150 numerical experiments. Fig. 8 shows the backcalculated results (black lines). The estimated s is 1.99 W/(m.K) (input value is 2 W/(m.K)), and R b is 0.207 m.K/W (input value is 0.2 m.K/W). The figure also shows that the PI algorithm is stable such that it demonstrates minor oscillations with noise. Fig. 10 shows the convergence rate of one of the experiments with = std 300W. It shows that both parameters were identified in 25 iterations. This relatively high convergence rate can be attributed to several reasons, including the use of the spectral element method for the forward calculation, and the use of the interior-point optimization method and Euclidean norm for the inverse calculation. 20 frequencies were utilized and on average it took 5 s in a normal laptop for the whole backcalculation procedure.

Numerical TRT in a layered system
In this subsection, a numerical TRT constituting a BHE embedded in a three-layer soil mass is conducted. Fig. 11 shows the geometry and material properties of the soil layers. The BHE physical and thermal properties are as those given in the previous example. Also, 150 numerical experiments with heat pump power and noise levels as for the previous example were conducted. The borehole thermal resistance is calculated using Eq. (24).
The PI algorithm is examined in two ways: 1) effective parameter identification, where the thermal parameters are estimated assuming that the soil mass consists of a single layer; and 2) detailed parameter identification, where the thermal parameters of all layers are estimated assuming that sufficient temperature readings along the pipe are available.

Effective parameter identification
The soil mass is assumed half space, characterized by effective (average) thermal parameters. The average thermal conductivity of the soil mass, as given in Fig. 11, can readily be calculated to give 2.16 W/ (m.K), and the average borehole thermal resistance to give 0.206 m.K/ W. If the layer thicknesses are considered, calculation of the weighted averages would give 2.3 W/(m.K) and 0.203 m.K/W, respectively. Fig. 12 shows the backcalculated effective parameters using both the proposed PI algorithm and the GI algorithm. The backcalculation due to the PI algorithm gave 2.188 W/(m.K) for the thermal conductivity, and 0.187 W/(m.K) for the borehole thermal resistance; which are close to both averages of the input values. Additionally, the figure reveals that the PI algorithm has minor oscillations for all noise levels. On the other hand, the GI algorithm gave different results depending on the start time of the curve fitting and the noise level. In all cases it gave underestimated values.
We furthermore examined the performance of the PI algorithm for different experiment durations: 12, 24, 48 and 72 h. Table 3 presents the estimated values of the effective s and R b , backcalculated from the PI algorithm and GI algorithm for different curve fitting periods (from 1, 6, 15 and 30 h). This independency on the duration of the experiment  is attributed to the capability of the forward model to simulate the highly transient period at the beginning of the experiment, and to the choice of the minimization procedure. On the other hand, the GI algorithm gave varying estimates with somewhat better performance for the longer experiment duration and shorter curve fitting. This explains why the TRT is conducted in practice for 72 h.

Detailed parameter identification
Detailed parameter identification entails estimating the thermal parameters of all involved layers and/or parameters including the volumetric heat capacities, c s s . As mentioned earlier, the basic principle of the parameter identification is that the system of equations must be either determinate or over-determinate. This implies that the number of measured set of readings must be equal or more than to the number of backcalculated parameters. To estimate the thermal parameters of the three-layer system of Fig. 11, we assume having four measurement sets along the BHE. One set is the standard inlet and outlet temperatures, and the other three measurement sets are taken in the U-tube, on the middle of each layer. Nine parameters are backcalculated: s and R b for each layer together with their associated volumetric heat capacities, c s s . The TRT operation time is assumed 12 h. Table 4 shows the estimated parameters compared to the exact (input values). It reveals that the estimated thermal conductivities for the three layers are nearly exact, the borehole thermal resistances have a maximum 5.5% error, and the volumetric heat capacities have a maximum 4.6% error.

Model verification
Witte et al. [1] and [37] had conducted several in-situ thermal response tests to backcalculate the thermal conductivities of a multilayer soil mass. The borehole heat exchanger is 30 m in depth and 0.25 m in diameter, filled with soil materials of the surroundings, and embedded in which a U-tube, 0.025 m in diameter. The working fluid in the U-tube consists of water with 17% ethylene glycol solution. Detailed description of the geometry and materials can be found in Witte et al. [1]. Fig. 13 and Table 5 (Column TRT) give an overview of the geometry and show that the soil mass consists of 8 layers of different soil types and thicknesses.
Witte et al. [1] also provided the range of thermal conductivities of each soil type, as given in literature, and laboratory results for nine soil samples taken from the site. Witte and van Gelder [37] conducted nonconventional experiments comprising three successive heat pump power pulses: two for heat injection and one for heat extraction. The power pulses lasted 96 h, followed by a switching off period until 169 h. The heat pump power together with the temperatures at the inlet and outlet of the U-tube were recorded for 169 h every 4 min. Figure 14a shows the heat pump power versus time.
Here, we make use of Witte et al. [37] TRT measurement to backcalculate the thermal conductivities of the involved soil layers with two varying level of complexities: 1. Effective parameter identification of the average soil thermal conductivity, assuming that all soil layers can be grouped into one. The backcalculation results are given in   involves noise. The backcalculation results are given in Table 5 (Column Parameter Identification). It shows that the backcalculated thermal conductivities are rather close to the laboratory results. The weighted average of the backcalculated values of all layers is 2.16 W/(m.K).
Figure (14b) shows the measured inlet and outlet temperatures, together with the theoretical outlet temperature calculated based on the backcalculated thermal conductivities. The theoretical outlet temperatures were calculated in two ways: based on the effective thermal conductivity (2.01 W/(m.K)), and based on the detailed thermal conductivities, as given in Table 5 (Column Parameter Identification). The figure obviously shows a very good matching between the measured temperatures and the computed from both ways. This gives the impression that backcalculating the effective thermal conductivity is sufficient to produce accurate results; avoiding thus the need for detailed backcalculations. This argument has been highlighted in literature, see for example Lee [51]. However, matching accurately the fluid outlet temperature does not necessarily reflect the accuracy in calculating the temperature distribution in the soil mass. Figure (14c) shows the associated temperature profiles in the soil mass in the z direction, 10cm away from the BHE. The temperature profiles are computed at the end of heat injection to the soil and the end of heat extraction. The figure clearly shows that calculations based on the detailed parameter identification exhibit different temperature distributions in the soil layers, depending on their thermal conductivities. This implies that relying on the effective thermal parameters might be useful for individual GSHP systems, but for regional GSHP systems where the thermal interaction between adjacent boreholes is crucial, knowing the detailed thermal properties of the soil layers is unavoidable.

Conclusions
The thermal response test (TRT) is an effective tool for quantifying the thermal parameters of the ground that form the basis for designing shallow geothermal systems. Despite the bulk of interpretation algorithms to interpret the TRT results, their accuracy and validity are limited to the duration of the test and the noise level in the measured data. This paper introduces a new parameter identification algorithm    Table 5 TRT field and laboratory data versus backcalculated.
e is the mean of two readings at 4 and 6 m. **Normal average of all the values.
that overcomes most of current limitations. It constitutes a semi-analytical forward model based on the spectral element method and an inverse model based on the interior-point optimization. Important features of the PI algorithm are: 1. It can handle fluctuating heat pump power. 2. It can interpret data obtained from multiple heat extraction or injection signals for heating or cooling modes. 3. It can produce accurate backcalculation for short and long duration experiments. 4. It can handle multilayer systems. 5. It can identify multiple parameters. 6. It is accurate and computationally efficient, stable and has a high convergence rate. Additionally, the numerical experiment on the multilayer system and the verification example demonstrated that the proposed PI algorithm has remarkable uniqueness even for multiple parameter identification. These features (accuracy, computational efficiency, stability, high convergence rate and uniqueness) constitute the essence of any parameter identification technique. They are manifested in the proposed algorithm due to the novel coupling between the spectral element method and the interior-point optimization method together with the Euclidean norm.
In practice, these features can help in: 1. Reducing the TRT cost by applying short duration experiments. The proposed forward model is accurate at both, the highly transient period of TRT and its steady state. This makes the backcalculation possible for almost any duration. 2. Facilitating better GSHP design. Knowing the detailed soil mass parameters yields accurate calculation of the temperature distribution in the soil mass, making the PI particularly suitable for regional GSHP systems design. Combining this parameter identification technique to a heat pump optimization tool such as that of Jin and Spitler [52], and to energetic, exergetic and economic performance optimization tools, such as those of Ozgener and Hepbasli [2], Conti [3] and Verrax [4], would lead to optimal GSHP designs.
In spite of the extensive testing of the proposed PI algorithm against noisy data and experimental measurements, it yet requires a study examining its uniqueness in handling accurate laboratory and field measurements constituting different soil types, heat pulse durations and noise levels.