Semi-Analytical Method for 3D Transient Heat Transfer in Thermally Interacting Fields of Ground Heat Exchangers

: The study of heat transfer in ground heat exchangers (GHEs) considering the ﬂuid advection inside the pipes; the heat transfer between the ﬂuid and the ground through the grout material; and the thermal interaction between GHEs is a challenging task. The present paper presents a new semi-analytical method that takes into account the aforementioned effects to consider both the short-to long-term effects of GHEs. The heat transfer between the ﬂuid and grout was studied by a transient multipole expansion considering time-dependent ﬂuid temperatures and an advection model for the pipes obtained from an energy balance on the heat carrier ﬂuid. Thermal interactions were analyzed using an equivalent borehole method while penalizing the transient multipole expansion to include thermal interaction effects. Validation of the short-term predictions was performed by comparing the proposed model to experimental data found in the literature and to an FEA model. The proposed model was then compared with a FEA model in long-term simulations of a geothermal ﬁeld comprised of 24 GHEs for multi-annual simulation. The method resulted in substantial reduction in computational time while preserving good accuracy when compared with the FEA model.


Introduction
One of the renewable energy sources that has been growing steadily worldwide in the last decade is geothermal energy [1]. Ground coupled heat pump (GCHP) systems are among the most efficient geothermal systems used for heating and cooling buildings [2]. GCHP systems consist of one (or several) heat pump that supplies heating or cooling to a building through the ground via ground heat exchangers (GHEs). GHEs are a core part of the GCHP system and consist of boreholes containing one or several U-tubes through which a heat-carrier fluid flows, and it is backfilled with grout material. In each of the GCHP design, control and operation phases, it is necessary to understand the heat transfer process between not only one but several GHEs. For instance, meeting the energy demands of a building requires that the GHEs are adequately sized so that returning fluid temperatures remain within the operating bounds of the GCHP system. Fluid temperature variations can be predicted from a model of the heat transfer between the fluid and the soil given ground thermal loads, accounting for possible thermal interactions between GHEs located in close proximity. Thus, improvements in the accuracy of heat transfer models lead to improvements in the accuracy of the design of GHEs.
Heat transfer phenomena are studied on different time scales, from short-term to long-term effects. Short-term effects are due to the flow of the fluid inside the U-tube coupled with the heat conduction between the pipes and the grout. At this time scale that ranges from minutes to a few hours, the thermal capacity of the grout decreases the rate of heat propagation between the fluid and the soil, and not taking this phenomenon into account can lead to an overestimation of the fluid temperatures [3]. The short-term effects end when quasi-steady conditions are reached inside the GHE, i.e., when temperature gradients are constant under a constant ground load. Long-term effects, spanning days to hundreds of years, include thermal interactions between boreholes and axial heat transfer effects due to interaction with the surface and the ground below the GHEs. It is necessary to consider both the short-and long-term effects to properly design a field of GHEs and ensure proper operation of a GCHP system throughout its life cycle [4].
Long-term heat transfer outside GHEs is studied both analytically and numerically. The well-known classical analytical solutions correspond to the infinite line source (ILS) [5] and the cylindrical heat source (CHS) [6]. Analytical solutions can be superimposed in time and space to consider variable ground loads and thermal interactions between GHEs. An important limitation of these classical solutions is that they cannot represent axial effects. Eskilson [7] solved the tridimensional heat transfer problem using a finite difference method. GHEs were represented as hollow cylinders with uniform and equal-but timevarying-temperatures. The model was used to produce dimensionless unit-step response functions, called g-functions, that give the relation between the borehole wall temperature and the heat extraction rates. Like analytical solutions, g-functions can be superimposed on time to consider variable ground loads. The g-function is defined by: where T b is the borehole wall temperature, T 0 is the undisturbed ground temperature, k s is the thermal conductivity, Q is the total heat extraction rate into the borefield, L is the borehole length and N b is the number of boreholes. Eskilson [7] proposed the analytical finite line source (FLS) to estimate the g-function of a single borehole. The FLS was later superimposed on space to evaluate g-functions of fields of GHEs [8][9][10]. For larger fields, the FLS can lead to important errors in the g-functions due to a mismatch of the boundary conditions at the borehole wall between the FLS (uniform and equal heat extraction rates) and the g-functions (uniform and equal temperatures). These errors can be eliminated by discretizing the GHEs to consider the variation in heat extraction rates along the lengths of GHEs and their variations in time to reproduce the uniform temperature condition of Eskilson [11,12]. In reality, neither the heat extraction rates nor the borehole wall temperatures are uniform along the lengths of the boreholes, since heat transfer is driven by the temperature difference between the fluid and the borehole wall. A more proper boundary condition is that of an imposed inlet fluid temperature, considering the internal thermal resistances of the boreholes and the piping connections between them [13][14][15]. The g-functions-and thus the returning fluid temperatures-are dependent on the borehole thermal resistances and the fluid mass flow rate. Spatial and temporal superposition of the FLS solution are also used for network-based models that study, at every time step, the complete distribution of the borehole wall temperatures considering the hydraulic connections between GHEs [16][17][18]. The computational requirements of analytical models increase as the complexity of the boundary condition increases. In the case of the imposed inlet fluid temperature, a system of equations must be built and solved to evaluate heat extraction rate profiles and borehole wall temperature profiles along all GHEs in the field. Prieto and Cimmino [19] proposed a method to limit the sizes of the systems of equations by forming groups of similarly behaving boreholes and modeling each of these groups using a single equivalent borehole. This results in an efficient and robust method that allows the study of geothermal fields composed of thousands of GHEs.
Knowing that internal configurations have a direct impact on the thermal interactions among boreholes and the long-term responses of GHEs, some researchers have coupled short-terms models with long-term models to conduct simulations. Yavuzturk and Spitler [20] were the first to introduce an extension of the g-function considering short time steps while including thermal capacity inside the GHE but disregarding the thermal capacity in the fluid. Later, Xu and Spitler [21], based on the work from Yavuzturk and Spitler [20], developed a one-dimensional model using a finite volume method, where the fluid is considered as an annular volume within the pipe, along with an equivalent thermal capacity. Brussieux and Bernier [22] coupled a finite volume method for the interior of a GHE with the CHS solution for the exterior of a GHE, resulting in the definition of g * -functions for short-term effects. Claesson and Javed [8] coupled an equivalent geometry model for the interior of the GHE [23] and used the FLS solution to study thermal interactions between the GHEs. Li et al. [24] used the analytical infinite composite-medium line source (ICMLS) [25] to model the insides of GHEs and the g-functions to account for thermal interactions among GHEs. In all aforementioned cases, despite the coupling of the short-to long-term solutions, the methods do not consider the advection inside the pipes and consider that all the GHEs have the same borehole wall temperature, which is not the case, as demonstrated by Cimmino [13] and Monzó et al. [26]. Pasquier and Marcotte [27] developed a response model based on a thermal resistance and capacitances (TRC) method to generate the thermal response factor of a GHE and estimate the fluid temperatures. Laferrière et al. [28] considered the axially-discretized thermal resistances and capacitances (TRC) model of Bauer et al. [29] coupled to an advection model and to g-functions, using a single uniform temperature at the wall of a GHE that represents the entire bore field. Recently, Alaie et al. [30] developed a solution based on the TRC method proposed by Minaei and Maerefat [31] and the segment-to-segment response factors proposed by Cimmino and Bernier [12] to consider both short-to long-term responses. As opposed to the earlier methods, their method does not consider that borehole wall temperatures are the same for all GHEs and thus achieves better representations of the distribution of heat extraction rates and temperatures within the field of GHEs. Additionally, there exist mesh-dependent approaches, such as those of Al-Khoury and Bonnier [32] and Diersch et al. [33,34], which are able to impose many boundary conditions and use many arrangements of GHEs; however, these methods are computationally demanding when considering multiple boreholes and including short-to long-term effects.
Axially-discretized methods become computationally demanding when very large fields of GHEs are considered. Thus, there is a need for accurate and computationallyefficient models that consider short-term effects, including advection in the pipes, and are able to simulate not only a single GHE but thousands. Coupled short-and long-term models can also benefit from recent advancements in short-term models, such as the transient multipole methods. Rivero and Hermanns [35] proposed a solution by means of matched asymptotic expansion to evaluate the thermal resistances in a spectral domain. Prieto and Cimmino [36] developed a transient multipole expansion to describe the variations in fluid, pipe and borehole wall temperatures in the time domain. In both approaches, the fluid thermal capacity and the advection of the fluid through the pipes are not taken into account. This paper introduces a new methodology to study the transient heat transfer for both the interiors and exteriors of GHEs, considering fluid flow inside the pipes, fluidto-soil heat conduction via grout, thermal interactions among GHEs and axial effects. The equivalent borehole method of Prieto and Cimmino [19] is used to model the exteriors of the GHEs, including thermal interactions between GHEs. The transient multipole method of Prieto and Cimmino [36] is used to model the interiors of the GHEs, including temperature penalization to account for thermal interactions. A finite difference model of fluid advection inside the pipes is used to consider the axial distributions of fluid temperatures and heat extraction rates. The model was first compared against the sandbox test of Beier et al. [37] to validate the accuracy of the short-term response, and then against a multi-annual simulation of a thermally interacting field of 24 boreholes using a finite element analysis (FEA) solver based on Diersch et al. [33,34] to show that the method is able to include both short-term effects and thermal interactions between boreholes for larger periods of time. Figure 1 shows a borefield composed of N b = 2 vertical GHEs, i and j, of the same dimensions and separated by a horizontal distance r ij . All GHEs are connected in parallel and are fed a heat carrier fluid at the same inlet temperature and mass flow rate. Each GHE is composed of a single U-tube (although the presented model can be extended to any U-tube configuration). The geometrical parameters are the same for each GHE and include: the borehole length L, the buried depth D, the borehole radius r b , the pipe outer radius r p with wall thickness t h and the U-tube shank spacing e. The thermal properties are considered homogeneous, isotropic and constant. k b and k s are the thermal conductivities (k) of the grout (b) and soil (s), respectively. α b and α s are the thermal diffusivities (α) of the grout and soil, respectively. The horizontal conduction heat transfer inside the GHE is studied in a transient state using the transient multipole expansion method developed by Prieto and Cimmino [36]. This solution is extended to consider time-dependent fluid temperatures inside the U-tube by means of a coefficient updating scheme. The radial and axial heat conduction in the soil around GHEs, and the interactions between GHEs, are modeled with the FLS solution by considering thermal interactions between equivalent boreholes as proposed by Prieto and Cimmino [19]. The transient multipole method is coupled to the FLS solution to consider both of the short-and long-term effects around GHEs, and to a discretized finite difference model of the advection heat transfer process inside the U-tubes. The proposed model was compared to a finite element analysis (FEA) model, based on the methodology developed by Diersch et al. [33,34], with specific modifications to ensure solution stability, as will be discussed in Section 2.4. The short-term transient heat transfer problem in the local vicinity of a GHE j is modeled in 2D for a horizontal cross-section of the GHE, considering conduction in the grout (Ω b ) and the soil (Ω s ). The heat transfer problem relies on the transformation of coordinates based on geometrical characteristics of the GHE; thus, a point x = (ρ, φ) in the 2D-polar space could be expressed at any center O k with translated coordinates x k = (ρ k , φ k ), as shown in Figure 2. The combination of the convection from the fluid to the pipe inner surfaces and the conduction across the pipes is represented as Robin boundary conditions at the pipe outer surfaces ∂Ω k (at a radius ρ k = r k from a pipe axis). The thermal capacitance of the pipes is neglected. The temperature at a far-field boundary ∂Ω e (at r e → ∞) is considered constant and equal to the undisturbed ground temperature T 0 . For constant fluid temperatures T f k ,j inside N arbitrary positioned pipes and a uniform initial temperature T 0 in the grout and soil (i.e., in Ω b ∪ Ω s ), the transient heat transfer problem in radial coordinates (ρ, φ) is given by:

Mathematical Model
where α d is equal to the thermal diffusivity in the grout or soil (α b and α s ) for indices d = 1, 2, respectively; β k = 2πk b R k is the dimensionless fluid to outer pipe wall thermal resistance of pipe k; and R k is the fluid to outer pipe wall thermal resistance. T f k ,j is the fluid temperature at pipe k in a GHE j.
Following the method of Prieto and Cimmino [36], the transient heat transfer problem with non-homogeneous boundary conditions for T d,j is separated into two sub-problems: a transient heat transfer problem with homogeneous boundary conditions for T d h ,j , and a steady-state heat transfer problem with non-homogeneous conditions for T d ss ,j . The complete solution is then given by superposition, with T d,j = T d h ,j + T d ss ,j .
The homogeneous transient heat transfer problem is then: The non-homogeneous steady-state heat transfer problem is given by: The solution to the homogeneous transient heat transfer problem is obtained by means of separation of variables and expressed as a Fourier-Bessel expansion: where λ w d,j is the w-eigenvalue associated with the separation of variables sub-problem for each domain d satisfying the continuity condition λ 2,j = λ 1,j √ α b /α s for a GHE j. X d,j (x; λ w d,j ) are the quasi-orthogonal eigenfunctions related to the eigenvalue λ w d,j for each domain d in GHE j.
For the grout domain (d = 1), containing N arbitrarily positioned pipes, the quasiorthogonal eigenfunction is expressed as follows: and for the soil domain (d = 2): where J and H (1) are Bessel and Hankel functions of the first kind; M is the number of terms of the multipole expansion; and γ 0 l,j , γ k l.j and δ 0 l,j are coefficients introduced to match the boundary conditions in Equations (3b)-(3e). C m,j are the Fourier-Bessel coefficients satisfying the initial condition (Equation (3f)) and are given by: The solution to the non-homogeneous steady-state heat transfer problem defined in Equation (4) is treated as a multipole expansion for a T-complete basis. For the grout domain (d = 1), containing N arbitrarily positioned pipes, the solution is expressed as follows: (9) and for the soil domain (d = 2): where r min = min i∈{1,...,N} r i is the minimum external radius of the pipes; r max ≥ r b is an arbitrary radius and is equal to r b ; and h is the number of terms in the complete expansion. Similarly to multipole expansion for the quasi-orthogonal eigenfunctions, terms α m,j , β m,j , γ l m,j , γ m,j , δ l m,j and δ m,j are determined by matching the boundary conditions in Equations (4b)-(4e). According to Prieto and Cimmino [36], the values of α 0,j = T 0 and γ 0,j , α m,j and β m,j are set to 0. Here, the authors note that, for the present method, keeping all terms expressed in Equation (10) stabilizes the integration in Equation (8), where terms γ 0,j and α m,j , β m,j are small when matching the far-field boundary with r e = 100 m. For more details concerning the evaluation of eigenvalues and coefficients to match the boundary conditions in both problems, readers are referred to Prieto and Cimmino [36].

Time-Dependent Fluid Temperatures
Prieto and Cimmino [36] introduced a coefficient updating scheme to consider the variations in fluid temperatures as step-wise variations at each time step n for every t (n−1) < t ≤ t n (where t n = t (n−1) + ∆t). The initial temperatures in Equations (3f) and (8) are replaced by the temperature at the end of the latest time step as follows: Therefore, the Fourier-Bessel coefficients are calculated using this scheme at each time step: with where T n d ss ,j is the solution of the steady-state problem considering the fluid temperatures at the n-th time step, and T 0 1 = T 0 2 = T 0 in a GHE j. Note that when Equations (11) are substituted into Equation (8), the sums in Equations (11) disappear, since the integrals with terms As written, Equation (12) is computationally intensive, since at every time step the integrals expressed in Equation (13) must be calculated. However, using the linearity in the steady-state problem described by Equation (4), it is possible to decompose the temperature field as a weighted sum of the fluid temperatures at each pipe. The temperature field can be expressed as N subproblems, and subtracting T 0 from the global problem yields: where ψ dl,j are the weighting functions in the domain d for a subproblem l in a GHE j: where δ l k is the Kronecker delta function which takes the value 1 if l = k and 0 otherwise. Equations (14) and (15) constitute a decomposition on N subproblems representing a unitary pulse at each pipe for each ψ dl,j ; therefore, ψ dl,j is time independent. Solving ψ dl,j follows the same multipole expansion structure shown in Section 2.1.1.
Substituting Equation (14) into Equations (13) yields: Substituting Equations (16) into Equation (12), rearranging and simplifying terms result in: with The terms E ml,j are independent of the fluid temperatures and thereby only need to be evaluated once at the initialization of a simulation.

Thermal Interactions between Ground Heat Exchangers
A temperature penalization approach is proposed to extend the transient multipole model to 3D, including axial heat transfer and thermal interactions between boreholes. Following the method of Cimmino and Bernier [12], each GHE is divided into n s segments (totaling (n s + 1) layers), and each segment is modeled as a line heat source with a uniform heat extraction rate positioned at the axis of the GHE segment. Figure 3 shows a uniform discretization for two boreholes i and j separated with a distance r ij . Each segment u has a length L u = L/n s . The finite line source (FLS) solution can then be superimposed to evaluate a penalty on the temperatures predicted by the transient multipole method. Additionally, the equivalent borehole method of Prieto and Cimmino [19] is used to simplify the thermal interaction problem and obtain a computationally-efficient simulation model of the bore field by only requiring the simulation of a limited number G of equivalent boreholes instead of the entire field of N b boreholes.

Spatial and Temporal Superposition of the Finite Line Source Solution
The temperature drop at the borehole wall of a segment u of a borehole i at time t n can be evaluated from the spatial and temporal superposition of the analytical FLS solution: where ∆T n b,i,u is the temperature drop at the wall of segment u of borehole i at time t k ; T n b,i,u is the temperature at the wall of segment u of borehole i at time t k ; Q p j,v is the heat extraction rate per unit length of segment v of borehole j from time t p−1 to t p ; and h ij,uv is the segment-to-segment thermal response factor for the borehole wall temperature change over segment u of borehole i caused by heat extraction from segment v of borehole j.
The segment-to-segment response factor is given by the FLS solution: where r ij is the radial distance between boreholes i and j (with r ii = r b ), erf(x) is the error function and erfint(x) is the integral of the error function.
Recently, Prieto and Cimmino [19] proposed a methodology to account for thermal interactions among boreholes as interactions between equivalent boreholes sharing similar borehole wall temperatures and heat extraction rates. The field is divided into G equivalent boreholes using a hierarchical aglomerative algorithm. Each equivalent borehole I represents a group G I composed of N b,I boreholes. The temperature drop at the borehole wall of a segment u of an equivalent borehole I at time t n is then: where ∆T n b,I,u is the temperature drop at the wall of segment u of equivalent borehole I at time t k ; Q p J ,v is the heat extraction rate per unit length of segment v of equivalent borehole J from time t p−1 to t p ; andh IJ ,uv is the equivalent segment-to-segment response factor for the borehole wall temperature change over segment u of equivalent borehole I caused by heat extraction from segment v of equivalent borehole J .
The equivalent segment-to-segment response factor is given by the average segmentto-segment response factors for the borehole wall temperature change along segments u of all boreholes in group G I due to the heat extraction at all segments v of all boreholes in group G J :h To reduce the computation time involved in the Equation (21), a load aggregation scheme is employed and presented in Appendix A. This scheme is already implemented in the pygfunction library [38], which was used here both for the aggregation of loads (Equation (21)) and for the calculation of segment-to-segment response factors (Equation (22)).

Penalization of the Transient Multipole Solution
A penalization technique is proposed to couple the transient multipole method presented in Section 2.1 with the equivalent borehole method. The transient multipole method can be used to evaluate the borehole wall heat extraction rate on ground layers, and the equivalent borehole method can be used to evaluate the borehole wall temperature change over segments (see Figure 3). Thus, an interpolation scheme is introduced to evaluate the heat extraction rate over a segment u from the transient multipole solution of its two delimiting layers, s and s + 1 (with u = s): where Q n I,u is the heat extraction rate per unit borehole length of a segment u of equivalent borehole I. Expressions for the temperature derivatives at the borehole wall are given in Appendix B.
The borehole wall temperature drop at a layer s is obtained from the borehole wall temperature drops at the two adjoining segments u − 1 and u (with u = s): where C n j,I,s is the j-th Fourier-Bessel coefficient at time step t n for a layer s in an equivalent borehole I, and λ Although the complete solution in Equation (28) looks computationally intensive due to eigenvalues, eigenfunctions, integrals and weighting functions, in practical applications, all boreholes possess a uniform geometry along their length and across the field. These terms thus only need to be evaluated once at the initialization of the simulation, except for the Fourier-Bessel coefficients that depend on the layer and the time step. Another interesting thing in Equation (28) is that even though T 0 , k s and α s are considered constant, the undisturbed ground temperature could change at each equivalent borehole, enabling the inclusion of thermal gradients and layered solutions when the soil thermal properties are not isotropic.

Advection Inside the Pipes
The model described in Sections 2.1.1 and 2.2 allows the evaluation of ground temperatures based on known fluid temperature profiles T f k ,I (z, t) at each pipe k at every equivalent borehole I. This model needs to be coupled to a model of advection inside the pipes to complete the simulation model and predict fluid and ground temperatures based on ground loads.
Assuming quasi-steady plug flow inside the pipes and neglecting heat conduction in the direction of flow (i.e., only considering heat convection between the fluid and the inner pipe wall), energy balances on a layer at downward-flowing pipe p and an upward-flowing pipe q are given by: where z is the axial position along the pipe (positive in the downward direction); ρ f and c p f are the fluid density and specific heat at constant pressure, respectively; andṁ n I is the average mass flow rate in a U-tube containing pipes p and q for an equivalent borehole I at the particular time step t n .
The fluid to outer pipe thermal resistance, R p , is given by the sum of the film resistance at the inner pipe wall and the conduction resistance across the pipe wall: where Nu is the Nusselt number, Re is the Reynolds number, Pr is the Prandtl number and k p is the pipe thermal conductivity. An implicit finite difference scheme is used to solve the problem. The GHE is discretized into layers, according to Figure 3. The discretized energy balance on a control volume delimited by two consecutive ground layers is given by: The number of unknowns for the fluid temperatures per equivalent borehole, T n k,I,s , at each time step t n in Equation (31), is equal to N(n s + 1). A continuity condition is imposed at the bottom of the equivalent borehole: where pipes p = q but belong to the same U-tube.
The system of equations is completed by imposing the ground loadsQ n , for the complete borefield. In this case, boreholes are assumed to be connected in parallel: where T n in and T n out are the inlet and outlet temperatures at time step n in the borefield anḋ m n T = N bṁ n b is the total mass flow rate, withṁ n b the mass flow rate per borehole. As the borefield is parallel-connected, the inlet temperature is the same for each equivalent borehole I: T n in = T n f p ,I,1 The outlet temperature of the borefield, T n out , is obtained by mixing of the fluid at the outlet of all equivalent boreholes: The problem is then complete, and it is possible to estimate the fluid temperatures at each equivalent borehole while imposing the heat extracted or rejected from the ground. An iterative scheme is used to solve the simulation model at each time step n:

1.
Calculate fluid temperatures, T n f p ,I,s , using the pipe temperatures T n−1 1,I,s r p at the previous time step (n − 1) as a first guess in Equation (31) with their boundary conditions and ground loads at time step n; 2.
Update the coefficients in Equation (27)  Compute the new fluid temperatures as in step (1) with the new pipe temperature calculated in step (4); 6.
Repeat from (2) if the difference of T n f p ,I,s between the previous and current iterations is greater than 10 −4 • C.
The method converges quickly, usually after two or three iterations. In this case, a linear system is solved with a guessed or updated value for T n 1,I,s r q ; thus, the unknowns in the system of equations for estimating the fluid temperatures are equal to G N(n s + 1). For the limiting case when the ground load is 0 W, the mass flow is set to 0 kg/s. Therefore, the solution of Equations (29) does not require solving a system of equations, since the term containing the derivative with respect to z is eliminated.

Finite Element Modeling
A reference finite element analysis (FEA) model was developed to validate the proposed method. Diersch et al. [33,34] developed a complete FEA model that includes the transient effects inside and outside GHEs and thermal interactions between them. This is an extension of the model previously proposed by Al-Khoury and Bonnier [32]. The same mathematical formalism is used in the present paper with some modifications.
The method was created on double-continuum finite element (i.e., one continuum media for the soil and one for the GHEs) space based on Galerkin discretization with 6-node prismatic elements for the 3D model. Each GHE is simulated as a line element surrounded by six nodes that virtually represent the borehole wall, located at a distance greater than r b , as presented in Diersch et al. [34]. Figure 4 shows the discretizations used for the two simulations used to validate the proposed semi-analytical methodology. First, the discretization shown in Figure 4a corresponds to a single borehole and was used for the comparison against experimental data. Then, Figure 4b shows the discretization for a borefield composed of 24 GHEs that was used in multi-annual simulations.  For each line element, short-term effects were modeled using the TRC model developed by Bauer et al. [29,39]. The correction noted by Bauer et al. [39] and presented in FEFLOW white papers (Section 1.5.5) [40] was used when negative internal thermal resistances were present. Depending on the number of U-tubes, each node over the line element presents two pairs of temperatures: two for the inlet temperature and the grout temperature where the inlet pipe is located, and two for the outlet temperature and the grout temperature where the outlet pipe is located. The TRC model is coupled to the FE-discretized equations considering a steady-state problem for a given borehole wall temperature, as solved by Eskilson and Claesson [41], to predict the fluid temperatures inside the pipes. This procedure is the same as that described by Diersch et al. [33].
A Picard iteration was employed to solve the system of equations and ensure the temperatures converge. Diersch et al. [34] suggest that the GHE domain should be separated from the whole domain and only later included to influence the soil by means of Schur complement operation. However, in the present FEA model, the Picard iteration achieved convergence after one or two iterations without the need to separate domains. Additionally, due to the initialization of the model, a few iterations with a small time step (i.e., ∆t/10 for a fixed time step ∆t) are required to stabilize the heat transport between the fluid and grout. After this stabilization, larger time steps (∆t) can be used.

Validation
The proposed method and the reference FEA model were first validated by comparison to the experimental results of Beier et al. [37] for a single borehole. Then, the proposed method was validated against the FEA model for a multi-year simulation of a borefield comprised of 24 boreholes. Inlet and outlet fluid temperatures and borehole wall temperatures were compared in both cases. The accuracy of the method was evaluated by computing the root mean square error (RMSE) and temperature differences between models and an experiment.
For the proposed method, the discretization along the GHE was done using 12 segments based on the recommendation of Cimmino and Bernier [12]. Following the recommendations of Prieto and Cimmino [36], the order of the transient multipole expansion was set to M = 7 and h = 15, corresponding to the homogeneous transient heat equation and non-homogeneous steady-state heat equation, respectively, and the first 100 eigenvalues were used for each simulation. All calculations were done on a personal computer with 16 GB of RAM and a 6-core processor (12 threads) running at normal speed of 2.60 GHz and maximum speed of 4.80 GHz.

Sandbox Test
Beier et al. [37] designed a sandbox test widely used in the context of short-term response validation of a GHE. The experiment consists of an average heat injection of 1051 W while maintaining an average mass flow rate of 0.197 kg/s in an 18 m long single U-tube GHE for 3106 min (approximately 52 h). The dimensions and parameters used for the experiment are shown in Table 1. Since in the experiment the GHE included an aluminum housing that affects heat transfer to the soil, Pasquier and Marcotte [27] use a modified thermal conductivity of the grout that preserves the borehole thermal resistance reported in the experiment. Therefore, the thermal conductivity calculated by Pasquier and Marcotte [27] was used in the proposed method.  The time step used in the simulation corresponds to the resolution time reported in the experiment and was equal to 60 s. Regarding the FEA model, 7260 6-node prismatic elements were used, equivalent to 4144 nodes in total. A square domain was proposed to consider the soil part with 6 m per side to avoid that the boundary conditions affect the thermal response of the GHE. The length of the FEA along the GHE axis was equal to 18.5 m. The ground load and mass flow rate reported in Beier et al. [37] were variable in time. Instead of using the average of ground load and mass flow rate, time varying values were used. Figure 5 shows the ground load and mass flow rate over the duration of the experiment.   The comparison between the proposed model and the experimental data is shown in Figure 6a, resulting in RMSEs for inlet and outlet fluid temperatures and for the average borehole temperature equal to 0.134, 0.131 and 0.105 • C, respectively. The maximum absolute errors for inlet and outlet fluid temperatures occurred at the 76th minute of the simulation and were 0.556 • C and 0.519 • C, respectively. In contrast, the maximum absolute error for the average borehole wall temperature occurred at the 196th minute and was equal to 0.253 • C. The fact that the errors did not occur at the same time for inlet and outlet fluid temperatures and the borehole wall is explained by the fact that the heat transfer process inside the borehole reached a quasi-steady state (or constant heat flux) at a time in the order of r 2 b /α s . Figure 6b shows the comparison between the FEA model and the proposed method computing inlet and outlet fluid temperature and average borehole wall temperature, similar to what was shown for comparison with the experimental data. For this case, the RMSEs for the inlet and outlet fluid temperatures and average borehole wall temperature were equal to 0.196, 0.196 and 0.183 • C, respectively. The maximum absolute errors for each temperature, T in , T out and T b , occurred at the end of the simulations with values equal to 0.317, 0.317 and 0.270 • C. Figure 7 shows a comparison of the fluid temperature profiles along the pipes between the proposed model and the FEA model at four different times, corresponding to 1, 10, 100 and 1000 min. Arrows indicate fluid flow direction. For the inlet pipe, the maximum absolute errors for 1, 10, 100 and 1000 min were equal to 0.003, 0.047, 0.013 and 0.125 • C. Similarly, for the outlet pipe, the maximum absolute errors were equal to 0.002, 0.048, 0.013 and 0.125 • C, corresponding to 1, 10, 100 and 1000 min of the simulation. The maximum absolute errors for inlet and outlet pipes were expected to be approximately equal, since the energy balance is preserved by the ground load in the discretization of the advection equation. The calculation time for Beier's experiment by means of the proposed method was equal to 214 s, including 210 s for the initialization and 4 s for the solution of the system of equations and coefficient updating scheme.

Multi-Annual Simulation of a Geothermal Field
Two simulations were used to evaluate the accuracy of the proposed method for long-term and short-term temperature predictions. A first multi-annual simulation was performed for a period of 10 years using a 1 h time step to evaluate the long-term accuracy of the method, followed by a period of 7 days using a time step of 1 min to evaluate the short-term accuracy of the method. The borefield was comprised of 24 boreholes in a rectangular configuration. Every GHE received the same inlet temperature and the same mass flow rate. The parameters used in the simulation are shown in Table 2.

Parameter Value
Borehole buried depth D = 4 m Borehole length L = 100 m Borehole radius r b = 0.075 m U-tube pipe outer radius r p = 0.0211 m U-tube pipe thickness t h = 0.00406 m U-tube shank spacing e = 0.05 m Ground thermal conductivity k s = 2.4 W/(mK) Ground thermal diffusivity α s = 1.2 × 10 −6 m 2 /s Grout thermal conductivity k b = 0.81 W/(mK) Grout thermal diffusivity α b = 2.1 × 10 −7 m 2 /s Fluid thermal conductivity k f = 0.48 W/(mK) Fluid thermal capacity (ρc p ) f = 4.02 × 10 6 J/(m 3 K) U-tube pipe thermal conductivity k p = 0.42 W/(mK) Undisturbed ground temperature T 0 = 12.5 • C Total mass flow rateṁ = 6 kg/s Number of segments n s = 12 The borefield was divided into equivalent boreholes that shared similar borehole wall temperatures and heat extraction rates using the methodology proposed by Prieto and Cimmino [19]. Figure 8 shows the borefield studied, in this case a rectangular configuration where each consecutive borehole was 6 m equidistant both horizontally and vertically, and the equivalent boreholes were identified. Equivalent boreholes denoted by the numbers I, I I, I I I each represent 12, 8 and 4 GHEs, respectively. The FEA model was composed of 123,828 6-node prismatic elements with 65,934 nodes in a ground domain, 300 m per side and an axial length equal to 140 m.  Figure 9 shows the complete simulation for the first 10 years of heat extraction and heat rejection with the ground (Q), including inlet and outlet (T in and T out ) temperature and the average borehole wall temperature for each equivalent borehole (T b,I , T b,I I , T b,I I I ) calculated with the proposed method. Figure 9a shows the ground loads (Q > 0 corresponds to heat extraction with the ground, negative otherwise) for the 10-year period, generated using the synthetic ground load profile of Bernier et al. [42]. The maximum heat extraction and rejection were equal to 93, 502 and −76, 704 W, respectively. The average net heat transfer for all 10 years was equal to 3761 W, which means an imbalance towards heat extraction. The impact of this imbalance is shown in Figure 9b,c, exhibiting decreasing behavior in both fluid temperatures and borehole wall temperature at the end of the simulation. The maximum and minimum inlet fluid temperatures were 19.98 and 0.95 • C, respectively, and the maximum and minimum outlet fluid temperatures were 17.29 and 4.25 • C, respectively. At the end of the simulation, the fluid had a temperature of 9.17 • C for both T in and T out , since there was no heat exchange at this time. This represents a difference between the temperature at the beginning of the simulation and the end of the simulation equal to 3.33 • C. For the average borehole wall temperature for each equivalent borehole, the maximum and minimum temperatures in the 10-year period were, respectively, 14.82 and 7.28 • C for equivalent borehole I; 14.86 and 7.37 • C for equivalent borehole I I; and 14.89 and 7.45 • C for equivalent borehole I I I. At the end of the simulation, the equivalent boreholes (I, I I, I I I) presented average borehole wall temperatures equal to 9. 36, 9.44 and 9.52 • C, respectively. Although the variation in the borehole wall temperature between the equivalent boreholes was small, it existed and shows the need to simulate in detail a field composed of several GHEs for larger periods of time. It is also worth mentioning that if the ground loads maintain their imbalance for a longer time, the average borehole wall and fluid temperatures will tend to decrease until they reach steady-state.   Table 3. Lastly, the RMSEs for the inlet and outlet fluid temperatures were 0.13 and 0.14 • C over the 10th simulation year, respectively. Therefore, the inlet and outlet temperatures show good agreement when compared with the FEA model.   The behavior between the soil and the GHE can be studied through the behavior of the borehole wall temperature. For this purpose, according to the groups identified in Figure 8, the average borehole wall temperature for each of the equivalent boreholes was calculated in the FEA model to compare with the proposed methodology. Figure 11 shows the average borehole wall temperature for each of the equivalent boreholes for the proposed method and the FEA model, and the error between each of them for the 10th simulation year. For Figure 11a-c, the average borehole wall temperature for each equivalent borehole was estimated by means of the proposed method; for Figure 11d The results corresponding to the time steps of the minimum and maximum differences and at the end of simulation are presented in Table 4. The largest error (in magnitude) was reached when the minimum temperature was reached. The RMSEs between FEA and proposed method were 0.10, 0.09 and 0.10 • C, for each equivalent borehole, I, I I and I I I, respectively. The errors show that the proposed model is able to predict with reasonable accuracy the thermal response between GHE and the soil.    The calculation time for this simulation was 310 s: 210 s dedicated to the initialization procedures (i.e., integrals, eigenvalues, coefficients, etc.) and 100 s to the solution of the systems of equations and iterations for the period of 10 years with a time step equal to 1 h.

Short-Term Effects after 10 Years
For the validation of the proposed method with respect to the short-term effects, an additional period of 7 days was simulated following the 10 years of simulation. The time step during this period was set to 1 min. The FEA model was initialized using the final results of the 10-year simulation as initial conditions. Figure 12 shows the inlet and outlet fluid temperatures and the borehole wall temperatures for each of the equivalent boreholes for the FEA model and the proposed method, and a comparison of the results. Figure 12a,b shows the temperatures obtained using the proposed method, and Figure 12c,d shows the same temperatures obtained with the FEA model. Figure 12e,f shows the differences between the FEA and proposed method. Table 5 presents the inlet and outlet fluid temperatures at the time of the maximum and minimum errors and at the end of the simulation, and the differences between the FEA and the proposed model. The RMSEs between the two models were 0.23 and 0.24 • C for the inlet and outlet fluid temperatures, respectively. This shows that the proposed method is able to accurately estimate the fluid temperatures when the short-term effects are considered. Table 6 compares the average borehole wall temperatures at each equivalent borehole at the times of the maximum and minimum errors and at the end of the simulation, and the differences between the FEA and the proposed model. The RMSE for each equivalent borehole, I, I I and I I I, was 0.16, 0.16 and 0.15 • C, respectively. As previously noted, these results show good agreement when compared with the FEA model.     −0.14 −0.14 −0.13 Figure 13 shows the fluid temperature profiles and the borehole wall temperature profiles for each of the equivalent boreholes along their lengths estimated with the proposed method and the FEA for the time step where the largest difference (in magnitude) between the two methods was obtained, which in this case corresponded to 7029 min. The arrows indicate the fluid flow direction for each pipe. The inlet temperatures for each equivalent borehole for the FEA and proposed method were equal to −0.05 and 0.4 • C, respectively. This amounts to the error previously computed as −0.46 • C between FEA and the proposed method. The outlet temperatures for the proposed method were equal to 4.15 • C for I, 4.19 • C for I I and 4.23 • C for I I I. On the other hand, for the FEA the temperatures were equal to 3.69 • C for I, 3.73 • C for I I and 3.76 • C for I I I. Figure 13b shows the equivalent borehole wall temperature profiles along the lengths of the GHEs for the FEA and proposed model. According to results in Table 6, the difference between FEA and proposed method for each equivalent borehole was preserved, and were similar to the differences in the average temperature, being −0.24, −0.24 and −0.23 • C boreholes I, I I and I I I, respectively. A small difference was observed at all equivalent boreholes, top and bottom, equal to −0.22 • C on both parts.

Conclusions
This paper presented a novel approach to the study of heat transfer in a field of thermally interacting GHEs where the heat capacity effects of the fluid, grout and soil are included. In the proposed method, the heat transfer inside the GHE is studied by the transient multipole expansion and coupled with the outside of the GHE by penalizing the interior solution with the equivalent borehole method to include the interactions between the GHEs. The advection inside the pipes is modeled by the energy balance along the pipes and discretized using a finite difference scheme. The method drastically reduces the dimensions of the system of equations due to the use of the equivalent borehole method to limit the number of modeled boreholes.
The main contributions of this paper are listed below: • The method is able to be used in efficient short-term simulations by extending the transient multipole expansion [36] with an efficient coefficient updating scheme to compute fluid and GHE temperatures. • It is no longer necessary to separate the solutions for short-and long-term effects as done by some researchers [8,24,28], since the heat extraction rates are calculated by the extension of the transient multipole expansion; therefore, the transition between time scales is done naturally by the approaching exact solution. Additionally, no condition is assumed at the borehole wall. • Calculation times show that it is possible to perform computationally efficient simulations, even when both short-to long-term effects are included. The majority of this calculation time is due to the initialization (210 s on every simulation). The initialization time for the multipole method is independent of the time step and the number of boreholes. The initialization time for the equivalent borehole method is only weakly dependent on the number of boreholes and the number of time steps. The simulations for Beier's experiment, the 10-year borefield simulation and the 7-day borefield simulation were done in 4 s, 100 s and 16 s, respectively.
The method was first validated against experimental data provided by Beier et al. [37] for a single GHE and against an FEA simulation of the same experiment. The maximum errors for the inlet fluid temperature, outlet fluid temperature and average borehole wall temperature were 0.556, 0.519 and 0.253 • C, respectively, when compared to the experimental data; and 0.317, 0.317 and 0.270 • C, respectively, when compared to the FEA method. The validation showed that the proposed model compares well with the FEA method, and that both methods compare favorably with the experiment. Differences between the proposed method and the FEA method were expected, since the two methods use different models for the interiors of the GHEs. A 10-year simulation of a field of 24 GHEs was then presented to evaluate the long-term accuracy of the proposed method when compared to the FEA method. A 10 year period was simulated with a 1 h time step, followed by a 7-day period with a time step of 1 min to evaluate the short-term accuracy of the proposed method in the context of long-term simulations. For the 10-year period, the last year presented maximum and minimum differences with the FEA method of 0.18 and −0.39 • C, respectively, for both the inlet and outlet fluid temperatures. For the 7-day period, the maximum and minimum differences with the FEA method were −0.24 and −0.46 • C, respectively, for both the inlet and outlet fluid temperatures.
The method is not restricted to applications where the GHEs are connected in parallel, but could also be extended to consider series connections between the boreholes by modifying the model in Section 2.3. The method could also be extended to include stratified solutions for anisotropic soils, adapting, for example, the solution by Abdelaziz et al. [43]. Additionally, by specific modifications of the FLS solution, the method could include effects such as groundwater flow by making use of proposed solutions such as the moving infinite line source (MILS) by Diao et al. [44], the moving finite line source (MFLS) by Molina-Giraldo et al. [45] or the stratified MFLS by Luo et al. [46]. This will, however, require future work to adapt the equivalent borehole method to consider groundwater flow during the clustering process.
For p ≥ 2, aggregated loads at a time t n are calculated according to the aggregated loads on the previous time step t n−1 : The first aggregated load (p = 1) has its value equal to the average load since the last time step. Here, loads are constant during a time step:

Appendix B. Additional Expressions for the Transient Multipole Method
In addition to the expressions presented in Prieto and Cimmino [36], the proposed method requires the evaluation of the average temperatures at the pipe walls ( T 1,I r q ) and the average temperature at the borehole wall ( T 1,I r b ), along with their derivatives (e.g., . Their evaluation necessitates expressions for averages and derivatives of the eigenfunctions (X 1 ) for the homogeneous problem, and the average temperatures (T 1 ss ) and derivatives for the steady-state non-homogeneous problem.
For the homogeneous problem, the average and derivatives of the eigenfunction X 1 at the pipe boundary (∂Ω l ) are given by: For the steady-state non-homogeneous problem, the average and derivatives of the multipole expansion at the pipe boundary (∂Ω l ) are given by: