Inverse heat conduction to model and optimise a geothermal field

The design of heat exchanger fields is a key phase to ensure the long-term sustainability of such renewable energy systems. This task has to be accomplished by modelling the relevant processes in the complex system made up of different exchangers, where the heat transfer must be considered within exchangers and outside exchangers. We propose a mathematical model for the study of the heat conduction into the soil as consequence of the presence of exchangers. Such a problem is formulated and solved with an analytical approach. On the basis of such analytical solution, we propose an optimisation procedure to compute the best position of the exchangers by minimising the adverse effects of neighbouring devices. Some numerical experiments are used to show the effectiveness of the proposed method also by taking into account a reference approximation procedure of the problem based on a finite difference method.


Introduction
Low-enthalpy geothermal applications provide sustainable, economic and safe solutions for the future energy supply. Their advantages are the very low environmental impact, the uninterrupted production during the summer and the winter, and the relatively simple installation that can be done in almost all geological settings.
The vertical exchangers are the most popular solution to exchange heat with the soil: a single or double U-shaped pipe is placed into a borehole [1], [2], and a carrier fluid is used to exchange heat with the ground material, which at suitable depth has constant temperature throughout the year. Their performances depend on a number of factors, such as exchanger depth, thermophysical properties of the ground material, thermal properties of the filling grout, fluid dynamics parameters of the fluid (velocity, viscosity, diffusivity, and so on). So optimising a geothermal exchanger means developing a fluid dynamics model that depends on those parameters and is able to consider their coupled effect [3]. However, a single geothermal exchanger is often not sufficient to satisfy the energy demand of a building, which requires the installation of multiple devices. A geothermal field consists in an array of devices tidily laid out in a given plot of land. For such systems, the complexity in the optimal design of a single device has to be jointly considered with the difficulty in the analysis of inter-device interference for optimising the geothermal field in a long time perspective. The modelling of such a complex system is a fundamental step in the design of a geothermal field, otherwise the energy extraction or injection for years may introduce a growing thermal anomaly in the ground, such as the local cooling or heating in the interior of the field [4]. Consequently, the thermal reservoir of the soil will vanish and the whole field will show a drastically reduced performance [5]. These thermal anomalies are partially compensated by the conductive heat flux through the boundaries of the field, but this is not very effective for systems with too much compact exchangers [4]. Thus, the design of an efficient geothermal field should avoid as much as possible thermal anomalies in the soil, even for the long-term heat extraction or injection. A key feature of the geothermal field (strictly connected to the soil temperature profile) is the above-mentioned arrangement of the devices, since the neighbouring exchangers may negatively interfere with each other [4]. In fact, a distance between devices of about 7-8 m is usually recommended [6]. However, such a prescription is rather empirical, since the complexity of the problem makes unreliable a general estimation of this distance for different geological settings and/or different systems [7]. A more precise design process has to take into account several system features as decision variables (such as the position of the exchangers, their geometrical features and the physical parameters of the ground), while its objective function should require the maximum possible produced energy. However, this quantity is very hard to be computed, and such an objective function would make the optimisation problem practically unsolvable. So, it is usually more effective a simplified objective function that can be computed directly from the soil temperature dynamics that is induced by the presence of the exchangers. The soil temperature dynamics is determined by the energy flux between exchangers and soil, which is influenced by a variety of factors [4], [8]- [10]: the number and the geometric characteristics of the devices, the arrangement of the devices, the lithotypes composing the stratigraphy of the soil from which depend both the heat diffusion rate and the presence of aquifers or rainwater infiltrations.
The rainwater infiltrations [11] and aquifer movements can influence the heat propagation in the soil. However, the estimation of soil moisture dynamics is not an easy task, it requires complex geological surveys and the evaluation of rainwater effects in a long-term perspective. Moreover, inaccurate estimations can introduce large bias in the design process with the consequent worsening of the results obtainable by using purely conduction models. Some studies have considered the advective transport under simplified conditions, for instance the well-known moving infinite line source model investigates the heat transfer only on a horizontal plane and neglects the axial effects due to the finite length of the exchanger [12]. However, the interesting study in [9] shows that the heat transport by soil moisture dynamics can be neglected if the Darcy velocity is lower than 10 −7 m/s. So the use of the conductive heat transfer model is a common practice to analyse the temperature of the soil around a geothermal exchanger, where groundwater dynamics is usually discarded, both in numerical computations [7], [13] and in analytical studies [14], [15].
Borehole heat exchanger (BHE) fields have been largely explored, since the time of the first lowtemperature geothermal applications in 1970s, so several theoretical and practical contributions have been already produced for BHE fields. The researches have analysed nonresidential groundsource heat pump installations [16], best practices for designing geothermal systems [17], and basic principles for the heating and cooling of buildings [18]. In particular, this last paper presents a strategy for the design of borehole fields known as Ashrae method, which is a fast procedure predicting the BHE field depth, based on the infinite cylindrical source solution (ICS) [19] and a penalty term to take into account the thermal interference between adjacent exchangers. The study [20] showed that previous methods fail in their task, and it proposed a new procedure for the accurate calculation of the penalty term that mimics the original Ashrae method. This improved procedure has also been tested by a sensitivity analysis [21] to show its reliability for an extended range of working conditions. Medium/large BHE fields with unbalanced seasonal loads are analysed in [22], where a new design method for BHE fields based on the superposition of the effects of different BHEs is proposed.
The methods described above provide not only models for the heat transfer in presence of exchangers but also strategies for the optimisation of their activity. In fact, once the model for the description of the heat exchange in the geothermal field has been formulated, it can be exploited in an optimisation procedure to support the design process. Usually, mathematical optimisation methods applied to borehole fields aim to strategically define the workloads of the devices [5] and/or to arrange the borehole positions [23], [24]. More precisely, optimisation techniques have been used in [23] to strategically operate and arrange in space borehole exchangers. On the other hand, an interesting fact is shown in [24] for homogeneous fields: an optimal load assignment yields the same result of an optimal BHE placement, and the combined optimisation approach produces only slightly better results. Thus, between the two control strategies, the placement optimisation has to be preferred to reduce the management effort of the BHE field. In general, the best position for a couple of exchangers in homogenous soils is the one that maximises the distance between them, in such a way that the mutual interference is minimised. When the array is made of more than two devices, the interference phenomenon among them becomes more complex and a regular lattice may be a suboptimal solution for the position of the devices, especially on heterogeneous soils having quite different diffusivities. So, such a geothermal system needs a detailed analysis.
In this paper, we consider a conductive model for the heat propagation into the soil and, on the basis of this model, we propose an optimisation method for the position of the exchangers that minimises the deviation of the soil temperature from the undisturbed temperature profile in a prescribed time interval. This model takes into account the mutual influence between the soil and the exchangers, since they constitute a coupled system. Thus, both the conductive heat transfer into the soil and the convective heat transfer into the exchangers are considered, as well as their joint action for the heat exchange inside a geothermal field. The key element in the proposed approach is the analytical solution of the heat conductive problem for a generic position of the exchangers. The derivation of this solution is based on standard arguments on the fundamental solution of the heat equation, and an inverse problem for estimating an ad-hoc source term in the conductive problem that mimics the effect of the exchangers. More precisely, this inverse problem allows a simple definition of the source term from the knowledge of the positions of the exchangers, providing in this way a powerful tool for the optimal arrangement of the exchangers. We finally perform a numerical experiment with the proposed method and the numerical results show good reliability and reasonable computational times.
In Section 2, the heat conductive problem for a geothermal field is formulated. In Section 3, the finite difference method is used to define a reference procedure for the approximate solution of this heat conductive problem. In Section 4, a forced heat conductive problem is introduced and the corresponding analytical solution is provided. In Section 5, an inverse source problem for the heat equation is defined and solved. In Section 6, an optimisation procedure for the placement of the borehole exchangers is described. Section 7 shows the results of a numerical experiment for the comparison of the analytical and numerical approaches, as well as some applications of the optimisation method. Finally, in Section 8, some remarks and further developments are provided.

Formulation of the heat conduction problem
We formulate the problem of the heat propagation into a geothermal field as an initialboundary value problem for the heat equation. We consider a three-dimensional bounded interval where N E is the number of exchangers placed in the field. The domain gives the space occupied by the soil surrounding the exchangers. All the exchangers have the same dimensions and the same workloads; they are approximated with parallelepipeds having square cross section with side length L E and depth H E . We consider the following heat conduction problem: is the set of the boundaries of all the exchangers, g is a function describing the initial soil temperature distribution that corresponds to the temperature profile of the undisturbed soil, α is the soil thermal diffusivity, T S is the temperature profile of the undisturbed soil before the system activation, which is supposed to remain constant in time, and T E is the temperature profile on the exchanger wall. We assume these compatibility conditions: g(x) = T S (x), x ∈ Γ S , and g(x) = T E (x, 0), x ∈ Γ E . The solution u of problem (2) describes how the temperature of the soil varies due to the presence of the exchangers. For simplicity, we analyse the heat transfer in the field up to a maximum depth H > H E ; moreover, we have used positive z to describe depth values, so the real geothermal field is modelled by a domain that is formally the symmetric with respect to the plane z = 0. This choice does not interfere with the description of the physical phenomenon.

Numerical solution of the heat problem
We briefly describe the approximation of the solution of problem (2) by a finite difference method. For simplicity, the space discretisation grid is built for the particular case of exchangers in a regular lattice, however different arrangements of the exchangers can be easily considered with the same procedure. Let N x , N y , N z be the partition size along the x, y and z directions, respectively. Let be the initial discretisation steps, which are adapted to fit the exchangers cross section. To fix idea, we consider the adaptation method on the x axis. The number N e x of grid points that are inside an exchanger is computed on the base of step h x , i.e., where [a] is the nearest integer to a ∈ R. To have exactly N e x grid points inside the exchanger, a new discretisation step is defined. Thus, the total number of grid points in the x direction may increase or decrease by few units with respect to the initial choice N x . For simplicity, we again denote the updated number of grid points by N x , N y , N z and the new discretisation steps by h x , h y , h z , respectively. The grid points are We denote the set of indices of all the inner points of the grid by I, i.e., the set of the indices of all the boundary points of the grid in Γ S by B S , i.e., and the set of the indices of all the boundary points of the grid in Γ E by B E , i.e., We note that the position of each exchanger is slightly adjusted along the x-y plane with respect to the discretisation grid, so, since the cross section of the exchangers contains an integer number of partitions, the definition of these sets of indices is particularly simple also for irregular lattices of exchangers. Moreover, there is one constraint the grid must satisfy: we prescribe that at least five inner points in I interpose between an exchanger and the nearest one, in order to have a reliable approximation of the laplacian operator. The discretisation of the time interval [0, t] consists in dividing it into N t partitions with step ∆t, that is and nodes t n = n∆t, n = 0, 1, . . . , N t .
We use this grid to compute the numerical solution U (x, t) of problem (2). We denote Applying the forward finite difference on the time derivative and the centred finite differences on the Laplacian term, we obtain the usual explicit Euler scheme: for (i, j, k) ∈ I, n = 0, . . . , N t − 1, where we can reorganise the terms to obtain for all (i, j, k) ∈ I, n = 0, . . . , N t −1. At the initial time, U 0 i,j,k = g i,j,k for all (i, j, k) ∈ I ∪B S ∪B E , where g i,j,k = g(x i , y j , z k ) for all (i, j, k) ∈ I ∪B S ∪B E . Finally, boundary conditions in problem (2) are used to obtain the solution U n+1 i,j,k for (i, j, k) ∈ B S ∪ B E . Therefore, we derived the complete solution U n i,j,k for all (i, j, k) ∈ I ∪ B S ∪ B E , n = 1, . . . , N t , for which the Courant-Friedrichs-Lewy stability condition [25] must be considered for the choice of the discretisation steps.

Analytical solution of the heat problem
Problem (2) and the corresponding numerical scheme (4) give a direct formulation of the heat diffusion in a geothermal field. This approximation approach provides also a quite efficient computational tool when the position of the exchangers is known. On the contrary, it is not so efficient to manage a geothermal field optimisation, because in this case problem (2) has to be solved several times with different positions of the exchangers. Thus we reformulate the heat diffusion in a geothermal field by a forced heat equation problem, whose solution can be computed analytically. Problem (2) and its numerical solution (4) are used in Section 7 as a benchmark for the comparison with the proposed model.
In Section 4.1, the diffusive problem is formulated by a forced heat equation, where the source term accounts for the contribution of the geothermal exchangers. In Section 4.2, we illustrate the procedure to obtain the analytical solution by means of the heat kernel, which can be easily used with the new formulation of the problem.

The analytical formulation
Let Ω = {x = (x, y, z) T ∈ R 3 : x, y ∈ R, 0 < z < H} ⊃ D be the domain for the space variables. The forced heat conductive problem is where f is the source term, g is the initial temperature distribution, T 0 , T H are the temperature values at the ground level and at depth H, respectively" which are supposed constant, α is the soil thermal diffusivity. We assume these compatibility conditions: The domain Ω corresponds to the three-dimensional slab of soil where the system of exchangers is placed, together with the exchangers themselves. So, problem (5) is strictly related with problem (2) and the action of the exchangers is described by f . The solution u of problem (5) formally describes the temperature of both the soil and the exchangers in the respective positions, but its main objective is the modelling of soil temperature among exchangers, since this is the main information used for the optimisation of the geothermal field.
Let u be the solution of problem (5), we write u as where is a suitable translation along the z direction such that v(x, t) satisfies the homogeneous boundary conditions, i.e., v(x, t) is the solution of the problem In the following, we describe an explicit formula for the solution of problem (8) and in turn of problem (5), as well as an approximation of the solution of problem (2), where the exchangers prescribe a proper source term f in problem (8).

The solution of the forced heat transfer problem
We illustrate the main steps for obtaining the analytical solution of problem (8). This procedure uses standard arguments on the theory of partial differential equations, so the reader can find all the details in [26].
The fundamental solution of the heat equation in R n is where || · || is the Euclidean norm; G represents the temperature at location x and time t resulting from an instantaneous point source of heat releasing a unit of thermal energy at location ξ and time τ . It solves the problem where δ(x − ξ) is the Dirac distribution on R n . Function (9) gives a powerful tool for obtaining solutions of forced heat diffusion problems [27], even in presence of inhomogeneous initial conditions and bounded domains. Let us consider the initial-value problem with inhomogeneous heat equation and an inhomogeneous initial condition. The solution is given by the following theorem: has solution Proof. The proof is based on the superposition principle, being the heat equation a linear partial differential equation. Thus, it is sufficient to add the solution of the associated homogeneous equation, i.e., with f = 0, satisfying the inhomogeneous initial condition and the solution of problem (11) equipped with initial condition g = 0, see [28] for the complete proof of this theorem.
Remark. This result for the unbounded spatial domain must be jointly considered with the case of bounded domain, since problem (8) is defined for (x, y) ∈ R 2 and z ∈ (0, H). Let us consider the homogeneous conductive problem on a rod of length H > 0, that is The solution of this problem can be computed in a way similar to problem (11), in fact, it has the following form where G is the fundamental solution of the heat equation in the domain (0, H) with homogeneous Dirichlet boundary conditions. From the separation of variables we have that In the following theorem we use the previous results to obtain the fundamental solution of the heat equation in the domain Ω of problem (8).

Theorem 2.
Let Ω = R 2 ×(0, H). Let G 1 , G 2 , G 3 be the fundamental functions of the heat equation in the space variables x, y, z, respectively; in particular, G 1 , G 2 are relative to problem (11) in R and G 3 to problem (13) in (0, H). Then G = G 1 G 2 G 3 is the fundamental function of the heat equation in Ω satisfying the boundary conditions of problem (8). Proof.
since G 1 , G 2 , G 3 verify the homogeneous heat equation. Also for the initial condition and in distributional sense, we have where ξ = (ξ, η, ζ) ∈ Ω. Finally, the boundary condition in x 0 is satisfied, in fact, since G 3 satisfies the null boundary condition in z = 0, and a similar argument holds for the boundary x H .
In our case, since the domain Ω of problem (8) is rectangular, Ω = R 2 ×(0, H), the fundamental solution G(x, t; ξ, τ ) of the heat equation in Ω is given by the product of the fundamental solutions related to proper one-dimensional problems. More precisely, from (9) and (15), we obtain Therefore, from Theorem 1 and relation (6), the complete solution of problem (5) can be expressed as follows where G(x, t; ξ, τ ) is given by (16).

Estimation of the source term
Formula (17) gives the solution of problem (5), where f is a given source function. In our case, this source term f is not known, but it has to mimic the presence of the exchangers. In this way, solution (17) also provides an approximation of the solution of problem (2) and it can be used in the optimisation of the geothermal field. In Section 5.1, we propose a method to obtain the source term as the solution of an inverse source problem. In Section 5.2, we illustrate a procedure, based on formula (17), for the explicit computation of the approximated solution of problem (2). In Section 5.3, we show a comparison among the exact solution of problem (2) and the solution of the corresponding source problem (5) in the case of a single exchanger, in order to supply a numerical evidence of the close equivalence between the two formulations.

The computation of the source term
The source term in problem (5) has to yield the thermal effect of the exchangers in the surrounding soil. More specifically, this source term has to produce a temperature distribution on the exchangers position that corresponds to the one on the exchangers wall, and a temperature distribution on the free soil that corresponds to the temperature profile of the undisturbed soil.
To this aim, we study two similar inverse source problems for the computation of the complete source. In fact, two types of sources are used: the exchanger source ε, with support on the exchanger E = {x = (x, y, z) T ∈ R 3 : x, y ∈ [−L E /2, L E /2], z ∈ [0, H]}, which has to produce the temperature surplus T E − T S on the exchanger; the free soil source σ, with support on the whole domain Ω, which has to produce the temperature profile of the free soil T S . Moreover, the contribution of each source ε, σ in the complete source f is: Note that the temperature of the exchanger depends on the surrounding soil temperature, and vice-versa, since there is a mutual influence that must be taken into account. We consider the following problems with source terms ε and σ: where the functions in the initial conditions are specifically chosen to provide a simplification in the analytical solution of the problems. From standard arguments on the heat equation, we have that the initial condition does not influence the solution for large values of t > 0. In particular, for our problems a characteristic value of t is 100 days (cold/warm season length) and the contribution of the initial term is negligible after 5 days. Thus, we can suppose that the solutions of problems (19), (20) are independent from the initial condition. The source terms ε(x, t), σ(x, t), x ∈ Ω, t ∈ (0, t) have to produce the following result for problems (19) and (20), respectively. The choice of the solutions e, u in Equations (22), (23), respectively, reveals that the soil source σ acts in any point of the domain Ω, exchangers included, while the exchanger source ε acts only in the points corresponding to the device position. We note that this choice for the source terms considers an exchanger placed at the domain center with respect to the x-y plane. Being the profiles T E , T S independent of x, y, we assume that the source terms ε, σ have the same symmetry, thus ε is valid for every exchanger in the field regardless of its position on the x-y plane. More precisely, the effect of an exchanger at (x, y) can be produced by ε(x − x, y − y, z, t); similar arguments hold for the source σ. This fundamental property is a simple consequence of the convolution form of G in (16), and allows us to avoid the computation of the sources at each point of the domain Ω, which would make the implementation of the proposed method computationally unaffordable.
We propose the following procedure to compute the solutions ε, σ of equations (22),(23), respectively. From formula (17) we have that the solution of problem (19) at x = 0, y = 0 (position of the exchanger) is where the space integrals in the x and y variables have been restricted to the horizontal cross section of exchanger E, assuming that E is the support of the source term ε. In a similar way the solution u of problem (20) can be translated by function g s and represented by where we assumed that the support of the source term σ isD = (−A, A) × (−B, B) × (0, H), that is the finite approximation of Ω. Formulas (24) and (25) whose coefficients e r , r = 1, 2, . . . , are given by From equation (24) and using the Fourier expansion (26), we obtain where Volterra integral equations (28) are considered for t ∈ [0, t] and are solved numerically by a simple quadrature method. We divide the time interval [0, t] into N t subintervals of length ∆t. So, the time variable t in discretised form becomes t n = n∆t, n = 1, . . . , N t , and the discretisation nodes of the time variable τ are τ p = (p − 0.5)∆t, p = 1, . . . , N t . From formula (28) by the midpoint quadrature rule, we obtain the following linear system for the source term coefficients: where R is the truncation index in the Fourier series, ε r (τ p ) are the unknowns, e r (t n ) are the known terms given by formula (27), and the matrix E r contains the approximation of the exponential time integral and the exact solutions of the space integrals with respect to x and y, i.e., where erf is the Gauss error function. An analogous procedure can be applied to s, thus integral equation (25) is discretised by the following linear system: where Systems (30), (31) are solved by means of LU factorisation. We note that an important property of the discretisation schemes (30), (31) is that they are decoupled with respect to the variable r, so the solution for each r = 1, . . . , R can be computed independently from the other ones.

The approximated solution for the BHE field
Once the source term has been estimated, formula (17) can be used to compute the solution of problem (5). We choose the initial temperature distribution g(x), x ∈D, equal to the function w(z) defined in (7); as already mentioned, this choice does not influence the temperature u for large values of t. For later convenience, we want to emphasise the dependence of u on the exchangers position, so we denote with p = (p 1 , p 2 , . . . , p 2N E ) ∈ R 2N E , the vector of the x, y coordinates of the exchangers centre. Thus, the solution of (5) becomes where x = (x, y, z) T ∈D, E l , N E , L E have been defined in Section 2, ε r , σ r have been defined in (29), (32), and We note that the source term f (ξ, τ ) in (33) actually does not depend on ξ, η but only on ζ, τ due to the simplifying assumption made in Section 5.1 about the independence of the source terms from the position of the exchangers. Moreover, the source has been split into ε (i.e., its Fourier coefficients ε r ) that accounts for the exchanger action and σ (i.e., its Fourier coefficients σ r ) that accounts for the soil action, according to Section 5.1. Consequently, we split the spatial integral into two parts: one over the exchangers support and the other over the soil support that is the wholeD. Another important consideration on the computation of u is that the source terms ε, σ neglect the interference among devices, which is instead considered in u by means of the superposition of functions ϕ. Formula (33) is approximated by a simple numerical scheme on D × [0, t]. We describe this approximation for the evaluation of u on a uniform grid; a similar approach can be used to evaluate u on different points. The time grids for t and τ are those defined in Section 5.1. The variable x is discretised as in Section 3, i.e., by the grid points (x i , y j , z k ) T = (−A + ih x , −B + jh y , kh z ) T , (i, j, k) ∈ I ∪ B S ∪ B E . Applying the midpoint rule to the time integral, the final form of the discretised solution is for all i = 0, . . . , N x , j = 0, . . . , N y , k = 0, . . . , N z , n = 1, . . . , N t , where t n−p = (n − p + 0.5)∆t, and the matrixes T , S are defined as follows,

Equivalence between heat conduction problems
The physical phenomena described by problem (2) and problem (5) are quite different. In fact, in problem (2) exchangers are outside the domain and their effect is considered by means of proper boundary conditions at the interface between the soil and the exchangers; whereas in problem (5) exchangers are part of the domain and their effect is reproduced as a heat source. Nevertheless, these two heat conduction problems can be equivalent; here, the equivalence of the two problems is based on the soil temperature generated by their solution. Now, we compare these two formulations in order to give a numerical evidence of their equivalence. For the sake of simplicity, this comparison deals with a single exchanger placed in a sufficiently large soil portion. This simple study is able to show that the effect of the exchanger on the soil temperature, which should be formulated as problem (2), is similar to the one of a proper heat source. In other words, the agreement of the results given by the two formulations allows the use of problem (5) to describe the heat transfer into a geothermal field, which has the advantage to make much more easier the solution computation especially with several exchangers.
In the case of a single exchanger, problem (2) can be stated in cylindrical coordinates, that is where the exchanger has circular cross section with radius a = L E /2. We solve problem (36) as done in Section 4, applying the translation T S (z) and the Green's function approach. The Green's function in the spatial domain outside a cylindrical hole is obtained by a multiplication theorem for cylindrical Green's functions similar to Theorem 2. In fact, if the problem has azimuthal symmetry, the multidimensional Green's function can be obtained by multiplying one-dimensional Green's functions [27], [29], that is where φ(β, r) = J 0 (βr/a) Y 0 (β) − Y 0 (βr/a) J 0 (β). Thus it can be shown that the solution of problem (36) is where e l , l = 1, 2, . . . , are the Fourier coefficients of the function T E − T S . See [29] for the general solution formula. The alternative formulation is given by problem (5) where the source term accounts for only one exchanger placed at the domain center. Thus, from formula (33), the solution v of this where x ∈D, p = (p 1 , p 2 ) = (0, 0), and ε r , σ r , ϕ have been defined in (29), (32), (34), respectively. Both solutions w and v are approximated onD × [0, t] by a numerical scheme similar to the one described at the end of Section 5.2. In addition, in formula (37), the domain of the integral over β is truncated at β = 5a/ √ α∆t and the resulting integral is approximated by the midpoint rule. We denote with w n i,j,k , v n i,j,k the discretised form of w, v, respectively, on the uniform grid used in Section 5.2. These results are reported in Section 7.

Optimisation of the exchangers position
We propose a method to compute the optimal placement of the exchangers in a geothermal field. In principle, this problem should maximise the energy exchanged with the soil but the corresponding objective function is too much complex for defining a practical design tool. However, a strictly related condition can be obtained by requiring a minimum deviation of the soil temperature with respect to the undisturbed temperature profile. The resulting optimisation process applied to the geothermal system is then made possible by using the relations discussed in Sections 4,5.
The objective function of the proposed optimisation problem is: where X = {x k = (x k , y k , z k ) ∈ D, k = 1, . . . , P } is the set of the evaluation points, P is the total number of evaluation points in D, T S is again the undisturbed soil temperature, u(x k , t; p) is the soil temperature at point x k and time t, defined in (33). We note that having required that the evaluation points lie in D means that (x k , y k ) = (p i , p i+1 ), for all k = 1, . . . , P , i = 1, 3, . . . , 2N E − 1, since we are interested in the evaluation of the temperature outside the exchangers. Moreover, the time interval where we search for the optimal arrangement can be in principle different from the time interval [0, t] where we need to analyse the temperature of the geothermal field. In particular, this optimisation process is considered for a time t =t with t ∈ (0, t). The optimisation problem consists in the minimisation of the objective function (39) constrained to C, i.e., min where the constraints are given by The solution of problem (40) is computed by using the steepest descent method, that is where ∇ is the gradient operator with respect to variables p, n ∈ N is the algorithm step, with p 0 the initial position of the exchangers, β is the parameter corresponding to the descent step length computed by a simple line search procedure. An active set strategy [30] is used to deal with the box constraints. We note that this is the simplest method among the class of gradient approaches, but it allows us to observe another important property of function F : its gradient can be calculated analytically avoiding the finite difference approximation; in fact, for i = 1, 3, . . . , 2N E − 1, we have: where ϕ has been defined in (34), while ε k has been defined in (29), and .
The derivative with respect to p i+1 can be computed by using a similar formula. The time integral in formula (42) is approximated by the midpoint rule. In addition to the analytical computation of the derivative, our method is also based on the choice of the evaluation points in set X that must be well separated from the position of the exchangers. A simple strategy for this choice is described in Section 7.3. Finally, the algorithm (41) stops when one of the following criteria is fulfilled: |p n+1 − p n | ∞ ≤ tol 1 or log 10 (F (p 0 ; X , t)) − log 10 (F (p n ; X , t)) ≥ tol 2 or n ≥ n, (44) where tol 1 , tol 2 ∈ R and n ∈ N are prescribed tolerances monitoring the distance between two successive solutions, the improvement in the objective function and the maximum number of allowed steps, respectively.

Results
We present some results obtained by numerical experiments to test the methods proposed in the previous sections, both for the analytical solution of the heat diffusion problem in a geothermal field and for the solution of the corresponding optimisation problem. The general context taken into account is a geothermal field operating in winter mode for about six months. As the geothermal plant is started, a conductive heat flux occurs from the warmer soil to the colder exchangers thus the soil temperature profile is modified by this heat loss, in turn the exchanger profile is modified by this heat gain. With the following simulations, we evaluate how the soil temperature field is affected by the presence of the exchangers. In Section 7.1, the physical and geometrical characterisation of the geothermal field is given. In Section 7.2, the simulation details are described. Finally, in Section 7.3, we show and discuss results obtained with the proposed methods.

The exchanger convective problem and data of the geothermal field
The temperature profile T E at the exchanger wall, i.e., z ∈ [0, H E ], has been obtained by exploiting the theory of completely developed fluid flows inside rectilinear pipes [31]. This hypothesis is reasonably valid for U-shaped exchangers, that constitute a simplified case but also one of the most common; other strategies should be considered for exchangers having different shapes. The convective thermal exchange between the fluid and the surrounding soil undergoes the first principle of Thermodynamics, that is where T is the mean temperature on the cross section of the pipe; t is a fixed time instant in [0, t]; x = (x, y, z) T is any point on the exchanger wall since there is azimuthal symmetry in the soil temperature; k is the thermal conductivity of the fluid; b is the pipe radius; ρ, c p are the density and the specific heat of the fluid, respectively; U is the mean velocity of the fluid flow and N u is the Nusselt number. Cauchy problem (45) holds for the fluid in the descending and ascending parts of the pipe; more precisely, we have: T = T d E , z = 0, T (t) = T in is the inlet fluid temperature, when considering the descending pipe; (H, t), when considering the ascending pipe. So, in this evaluation, we neglect the contribution provided by the small curve at the bottom of the exchanger. From standard arguments on ordinary differential equations, we have: where T f E = T d E (H E , t) is the temperature value at the bottom of the exchangers. On the other hand, the temperature profile of the undisturbed soil given by T S is where H is the maximum depth at which the soil temperature is expected to be influenced by seasonal variations. The soil temperature is assumed to linearly increase until depth H and then to have a constant value, which corresponds to the mean value in pelitic and pelitic-arenaceous lithotypes in the central Italy [32]. Note that the geothermal gradient is neglected since its effect is irrelevant for the depth of the geothermal fields under consideration.  [33], k = 4.1412 · 10 4 J/(day m K), ρ = 1.0411 · 10 3 kg/m 3 , c p = 3.6915 · 10 3 J/(kg K). The flow mean velocity has been estimated U = 3.6288 · 10 4 m/day, according to data in [33]. The evaluation of N u is obtained by N u = 0.012 Re 0.87 D − 280 P r 0.4 , where P r is the Prandtl number and Re D is the Reynolds number for pipe flows, see [3], [31] for a detailed discussion on empirical laws for N u. The soil thermal diffusivity α is 0.05 m 2 /day, unless otherwise specified, T 0 = 280 K, T H = 286 K, H = 15 m. The temperature in the geothermal field has been simulated for a time interval of 180 days, that is t = 180 days in problem (2), which in principle corresponds to the six-month period when the exchangers have winter operational mode. Finally, in the optimisation problem,t = 120 days has been used.

Settings of the numerical and analytical solutions
The numerical solution of problem (2) is computed by the finite difference scheme (4)  We note that the computation of the source terms and the analytical solution cannot be done consequently, in fact, problems (2) and (45) are coupled since the soil temperature and the exchanger temperature influence each other. Thus, at each time step we implement a fixed-point iteration on T E in this way: at the time step n + 1, from (46) the solution (T E ) n+1 k requires the knowledge of T d E and T u E thus of u n+1 i,j,k ; in turn, the solution u n+1 i,j,k in (35) requires the knowledge of (T E ) n+1 k . Such solutions are available at the previous time step n, thus we use u n i,j,k to first calculate the predictors (T E ) n+1,0 k and consequently u n+1,0 i,j,k , then we compute the correctors (T E ) n+1,ν k and u n+1,ν i,j,k for ν = 1, 2, . . . , until a stop criterion on T E in infinite norm is verified.

Numerical results
We firstly show the results of the comparison between the two formulations of the heat conduction in the geothermal field, see Section 5.3 for details. Figure 1 shows the temperature distribution at z = 20 m in the soil caused by the presence of a single exchanger. In particular, Figure 1(a) reports the temperature field (37) that is solution of problem (36), where the effect of the exchanger on the surrounding soil is classically described by a boundary condition on the interface between soil and exchanger. On the other hand, Figure 1(b) reports the temperature field (38) that is solution of problem (5), where the effect of the exchanger is described by means of a source term in the heat equation. The two temperature fields are almost the same. For a quantitative comparison of these two solutions, we calculate the relative error between the discretised form of w in (37) and v in (38) in norm 1 at each time step, that is where k = 40 is the index of the half-depth plane (that is at z = 20 m in these simulations) where the results have been shown by figures. Similar results are obtained for different choices of z due to the symmetry of the problem therefore have not been reported here. Figure 2 shows that the relative error is of the order of 10 −6 , so the two solutions are in good agreement. The results reported in Figures 1 and 2 gives a numerical evidence of the reliability of the proposed approach to solve analytically problem (2) by means of the reformulation (5). We expect that the two solutions are in good agreement also in the case of multiple exchangers. However, in this case, we cannot solve the heat conduction problem (2) by a simple formulation as the one given by problem (36); so, a numerical approach must be used to solve this problem, as the one adopted in (4). We show the results of the numerical experiment for the comparison between the finite difference solution of the boundary value problem (2) and the forced conduction problem (5). these figures are relative to the horizontal section of the field taken at half depth, that is z = 20 m. Figure 3 shows the temperature distribution at the initial time, where the soil has a uniform temperature that corresponds to the value of the undisturbed profile at half depth, i.e. 286 K, while the exchangers are colder (at 280 K that is the temperature of the incoming water) and correspond to the blue points. Figure 4 shows the temperature distribution at the final time step, i.e., after 180 days of operation of the geothermal system; in particular, Figure 4(a) shows the numerical solution computed by the finite difference scheme, i.e., formula (4), whereas Figure 4(b) shows the analytical solution computed by (35). In both these figures, cold areas arise in neighbourhoods of the exchangers. In particular, Figure 4(a) shows a more evident interference among exchangers but the areas in the cross points of the diagonals among four devices are almost undisturbed being at 285.9 K; in Figure 4(b), we can see a situation similar to Figure 4(a), with a slightly slower heat diffusion. The soil affected by interference in Figure 4(a) differs in temperature from the corresponding areas in Figure 4 qualitative point of view there is a good agreement between the two solutions. For a quantitative comparison of the two solutions, we consider the relative error (48) where w n i,j,k is replaced by u n i,j,k in (35) and v n i,j,k by U n i,j,k in (4). The diagram in Figure 5 shows this relative error for the various simulation time steps; more precisely, we can observe err 1 ≈ 2.2 · 10 −5 at the initial time and err Nt ≈ 4.2 · 10 −4 at the end. Different arrangements of the geothermal field show similar results, providing an experimental verification of the reliability of the analytical solution.
In the next experiment we consider the optimisation procedure described in Section 6. The uniform lattice, like the one in Figure 3, is taken as the initial guess for the iterative algorithm. The tolerances and the maximum number of steps in (44) have been chosen as follows: tol 1 = 0.1, tol 2 = 2, n = 150. The evaluation points in X are not fixed throughout the process but at each step n they are chosen at regular distance with respect to the exchangers position p n . Such strategy prevents unfortunate situations in which an exchanger moves to an evaluation point. Figure 6 shows the initial and final arrangements of the exchangers, the little blue circles are the starting positions and the red stars are the optimal positions obtained from problem (40). In Figure 6(b), the devices are moved away from each other but the symmetry of the initial configuration is maintained. The central devices are kept near their original position while the external ones are pushed towards the boundary to better exploit the useful space and allow a bigger distance among all the exchangers. Also the heat exchange between soil and devices in the optimal arrangement has been computed; Figure 7 reports these results. The temperature distributions computed by the finite difference method (Figure 7(a)) and the analytical approach (Figure 7(b)) are in good agreement. As in the previous experiment, the numerical solution develops a slightly bigger interference among exchangers, but the maximum difference between these two solutions in the soil portions affected by the influence of multiple exchangers is again no bigger than 0.1 degrees. Comparing the two solutions, an interesting feature is that in Figure 7(a) the colder areas around the boundary devices lose the circular symmetry and tend to develop mostly in the direction of the neighbouring exchangers, while the mitigating effect of the boundary condition slightly advances among the boundary devices; on the other hand, in Figure 7(b) all the devices develop around colder areas with circular symmetry since there is no boundary condition. Also for this case, the relative error (48) observed during the simulation ranges from 2.0 · 10 −5 to 3.8 · 10 −4 , see Figure 8. We note that the relative errors in the case of 16 exchangers  with regular arrangement ( Figure 5) and the case with optimal arrangement (Figure 8) are very similar, meaning that the accuracy of the analytical solution is independent of the position of the devices. We also note that the effectiveness of the optimal positioning is clearly visible by the soil temperature distribution. In fact, a brief comparison between Figures 4 and 7 shows that the thermal impact in the soil on the long-term conditions can be effectively reduced by the optimal positioning of the devices, even applying a constant heating load. In Figure 7, the thermal effect of the exchangers in the soil is mitigated as a consequence of the optimisation process, as can be seen from the bigger undisturbed soil portions among the exchangers. The advantage of our analytical approach is to make easier the solution computation of the heat problem, dividing it into the calculation of the direct solution by the Green's function and the source computation by an inverse problem technique, which is rather fast since it is independent of the position of the exchangers on the x-y plane. Moreover, this analytical approach makes feasible the proposed optimisation technique. We conclude the numerical experiment with two examples that emphasise the value of the physical-based approach used in the proposed optimisation procedure. Indeed, the previous example can be also solved by using a simple geometric optimisation approach, where the uniform distribution of the devices is computed with respect to the available space. On the contrary, a geometric approach cannot properly deal with a heterogeneous soil situation. So, in the following example, the field is assumed made up of two equivalent parts with extremely different values for thermal diffusion. The two soils are separated by the vertical median ofD; the soil in the left side has thermal diffusion 0.39 m 2 /day, and in the right side 0.022 m 2 /day; these values correspond to the mean diffusivity for dry gravel and rock salt [34], respectively. The stopping criteria in problem (40) are defined as before. Figure 9 shows the final configuration of 16 exchangers in the heterogeneous geothermal field. The initial configuration is the same of Figure 6(a). In Figure 9, the exchangers initially placed in the left side of the domain tend to maximise the distance among each other exploiting all the surrounding space. The two central devices increase their respective distance moving towards the exterior without going far from their original position. On the contrary, the exchangers initially placed in the right side of the domain keep their original position, except for the first and last exchanger of the third column, which move slightly away due to the influence of the two devices placed at the boundary points with x coordinate equal to 30. It is clear from Figure 9 that in each of the two parts with constant thermal diffusivity, the final configurations maintain the symmetry with respect to the axis y = 35, in accordance with the symmetry of the problem data; of course, other symmetry properties, like global symmetry, are no longer obtained in the final solution. Figure 10 shows the initial and final configurations of 30 exchangers into the heterogeneous geothermal field. The arrangement reached in Figure 10(b) is similar to that in Figure 9; in fact, the optimisation algorithm favours devices placed along the domain boundary over the central positions in the soil with higher diffusivity, and leave the arrangement quite unchanged in the other part. Another noteworthy feature is that the central column of exchangers occurs on the separation line between the two areas with largely different diffusivities, therefore these devices slightly shift to the right, where the thermal exchange is lower. In conclusion, these two tests on a heterogeneous field show that the proposed optimisation process is not a geometrical optimisation procedure, but it actually takes into account the physical characteristics of the problem. Finally, in Table 1 we show the computation time of the optimisation procedure for the examples considered in this section. All the simulations have been run in the machine HP ProLiant ML350 Gen9, equipped with Red Hat Enterprise Linux Workstation release 7.5 (Maipo) and 2 CPU Intel(R) Xeon(R) CPU E5-2620 v3 @2.40GHz, with 6 physical cores and 6 threads for each CPU. In particular, the simulation code has been written in Matlab and exploits one core. The computation time obviously increases as the number of exchangers increases. The less time-demanding procedure is that one for 16 exchangers in the homogeneous soil, where we have a shorter computation time of each iteration than for the heterogeneous soil case. This fact is mainly due to a simplified algorithm implementation in the homogeneous case. On the contrary, we can observe that the number of iterations for the heterogeneous case study is slightly lower than the one for the homogeneous case study. This is mainly due to the fact that in the homogeneous case the gradients in the steepest descent method are smaller than the ones in the  In Table 2, we also show a comparison between the computation time of a single iteration of the proposed optimisation procedure and the one of a similar optimisation procedure where the analytical method is replaced by the finite difference method for the heat diffusion computation. Each step of the proposed optimisation procedure consists in, as described in Section 6, the calculation of the solution (35) and its derivatives with respect to the exchangers position that are contained in formula (42). On the other hand, each step of the optimisation procedure based on the finite difference approach consists in the calculation of the solution U by using the time-marching scheme (4) and the finite difference approximation of its derivatives, that requires 4N E + 1 evaluation of the solution U , thus the computation of U in the whole space and time grid must be repeated 4N E + 1 times. We note that in this way we simply have a lower estimate of the computation time for the finite difference approach. In particular, Table 2 considers the configuration with homogeneous soil and 16 exchangers, operating for 120 days, where the computation is performed by using three grids with increasing fineness having smaller and smaller step size h x = h y along the x, y directions and a constant step size h z = 0.5 m along the z direction. Even the coarsest discretisation shows a big gap in the elapsed time for the two computational approaches, where the analytical method is about 76 times faster than the numerical approach. The reduction rate is increasingly higher as the discretisation steps reduce and for the finest grid the proposed optimisation procedure is about 680 times faster than the procedure based on the finite difference approach. We remind that such a gap in the computation time is also due to the fact that a reduction of the space step induces a reduction of the time step, according to the Courant-Friedrichs-Lewy condition. Instead, in the analytical method the usual time step ∆t = 0.8 days has been used, and it could be even increased to make the computation faster even keeping a comparable accuracy. In addition, the optimisation procedure based on the finite difference scheme reaches the optimal solution only if the space grid is sufficiently fine, that is h x , h y ≤ 0.25 m. Hence, Table 2 shows that the finite difference approach used within an optimisation procedure is definitely unfavourable in terms of computation cost.

Conclusions
The paper studies the geothermal fields of borehole heat exchangers for the exploitation of shallow geothermal energy resource. This study considers the coupling of the soil and exchangers thermal response, which consists in taking into account the mutual influence between the cooling of the soil and the corresponding warming of the exchanger (in winter operational mode). Two main contributions are obtained: the development of an analytical method to quantify the temperature evolution within the geothermal field, and the development of an optimisation procedure for the positions of the devices. These two results are strictly related, in fact, the heat exchange in the geothermal field is given in terms of the fundamental solution of the heat equation, allowing a direct connection between the exchangers position and the temperature field in the soil. The soil temperature obtained with the proposed method is compared with the result obtained by standard approximation methods, i.e., the Green's function approach in the case of a unique exchanger, and the finite difference method in the case of multiple exchangers. From this numerical experiment, we can conclude that the proposed method is reliable for describing the temperature distribution in a geothermal field. Moreover, the optimisation procedure is tested by using some reference cases, showing again a good reliability of the proposed procedure. It is also shown that the induced thermal anomalies are reduced by the optimisation process.
Despite the good results obtained with the proposed method, some assumptions made in the paper could be further refined. For instance, the undisturbed soil thermal profile could be adapted to be as close as possible to a real profile. In addition, the inverse problem for the estimation of the source term that mimics the exchangers presence should take into account the interference among sufficiently near devices. Finally, on-field experiments for the temperature measurements in a real geothermal field should be performed to obtain a practical validation of the proposed approach.