Modeling and simulation of the input–output behavior of a geothermal energy storage

This paper investigates mathematical modeling and numerical methods for simulations of the input–output behavior of a geothermal energy storage. Such simulations are needed for the optimal control and management of residential heating systems equipped with an underground thermal storage. There, a given volume under or aside of a building is filled with soil and insulated to the surrounding ground. The thermal energy is stored by raising the temperature of the soil inside the storage. It is charged and discharged via pipe heat exchangers filled with a moving fluid. Simulations of geothermal energy storages aim to determine how much energy can be stored in or taken from the storage within a given period of time. The latter depends on the dynamics of the spatial temperature distribution in the storage, which is governed by a linear heat equation with convection and appropriate boundary and interface conditions. We consider semidiscretization and full discretization of that PDE using finite difference schemes and study associated stability problems. Finally, we present numerical results.


INTRODUCTION
Thermal storage facilities help to mitigate and to manage temporal fluctuations of heat supply and demand for heating and cooling systems of single buildings as well as for district heating systems.They allow heat to be stored in the form of thermal energy and be used hours, days, weeks or months later.This is attractive for space heating, domestic or process hot water production, or generating electricity.Note that thermal energy may also be stored in the way of cold.Thermal storages can significantly increase both the flexibility and the performance of district energy systems and enhancing the integration of intermittent renewable energy sources into thermal networks (see Guelpa and Verda [1], Kitapbayev et al. [2]).Since heat production is still mainly based on burning fossil fuels (gas, oil, coal), these are important contributions for the reduction of carbon emissions and an increasing energy independence of societies.
For an overview on thermal energy storages, we refer to Dincer and Rosen [3].Zalba et al. [4] provide a review for the history of thermal energy storages with solid-liquid phase change and focused in three aspects: materials, heat transfer, and applications.An overview of the European thermal energy storage potential is presented in Arce et al. [5].The authors show that thermal energy storages make an important contribution to the reduction of CO 2 -emissions.In Soltani et al. [6], the authors provide a comprehensive review on the evolution of geothermal energy production from its beginnings to the present time by reporting production data from individual countries and collective data of worldwide production.
The efficient operation of geothermal storages requires a thorough design and planning because of the considerable investment cost.For that purpose, mathematical models and numerical simulations are widely used.We refer to Dahash et al. [7] and the references therein.In that paper, the authors investigate large-scale seasonal thermal energy storages allowing for buffering intermittent renewable heat production in district heating systems.Numerical simulations are based on a multiphysics model of the thermal energy storage, which was calibrated to measured data for a pit thermal energy storage in Dronninglund (Denmark).Another contribution is Major et al. [8], which considers heat storage capabilities of deep sedimentary reservoirs.The governing heat and flow equations are solved using finite element methods.Further, Regnier et al. [9] study the numerical simulation of aquifer thermal energy storages and focus on dynamic mesh optimization for the finite element solution of the heat and flow equations.For further contributions to the numerical simulation of such storages, we refer to [3,6,[10][11][12][13].
In this paper, we focus on geothermal storages as depicted in Figure 1.Such storages gain more and more importance and are quite attractive for residential heating systems since construction and maintenance are relatively inexpensive.Furthermore, they can be integrated both into new buildings and in renovations.We will work with a 2D model of a geothermal thermal energy storage (see Figure 2), where a defined volume under or aside of a building is filled with soil and insulated to the surrounding ground.Thermal energy is stored by raising the temperature of the soil inside the storage.It is charged and discharged via pipe heat exchangers (PHX) filled with some fluid such as water.These PHXs can be connected to a short-term storage such as a water tank or directly to a solar collector and (heat) pumps move the fluid carrying the thermal energy.A special feature of the storage in this work is that it is not insulated at the bottom such that thermal energy can also flow into deeper layers as it can be seen in Figure 2.This can be considered as a natural extension of the storage capacity since that heat can to some extent be retrieved if the storage is sufficiently discharged (cooled) and a heat flux back to storage is induced.Of course, there are unavoidable diffusive losses to the environment but due to the open architecture, the geothermal storage can benefit from higher temperatures in deeper layers of the ground and serve as a production unit similar to a downhole heat exchanger.Note that in many regions in Europe the temperature in a depth of only 10 m is almost constant around 10 • C over the year.
Geothermal storages enable an extremely efficient operation of heating and cooling systems in buildings.Further, they can be used to mitigate peaks in the electricity grid by converting electrical into heat energy (power to heat).Pooling several geothermal storages within the framework of a virtual power plant gives the necessary capacity, which allows to participate in the balancing energy market.
This paper extends and complements the results in Bähr et al. [14,15], where the authors focus on the numerical simulation of the long-term behavior of the spatial temperature distribution in a geothermal storage over weeks and months and the interaction between a geothermal storage and its surrounding domain.For simplicity, charging and discharging was described by a simple source term but not by PHXs.
In the present work, we focus on the computation of the of the spatial-temporal temperature distribution and the associated input-output behavior of the geothermal storage on smaller time scales as in [14,15].This is needed for storages embedded into residential heating systems and the study of the storage's response to charging and discharging operations on time scales from a few minutes to a few days.We extend the setting in [14,15] and include PHXs for a more realistic model of the storage's charging and discharging process.However, for the sake of simplicity, we do not consider the surrounding medium but reduce the computational domain to the storage depicted in Figure 2 by a black rectangle.Instead, we set appropriate boundary conditions to mimic the interaction between storage and environment.
For the management and control of a storage, which is embedded into a residential heating system, one needs to know the amount of available thermal energy that can be stored in or extracted from the storage in a given period of time.Such questions can only be answered if one knows the spatial temperature distribution, in particular around the PHXs.Charging and discharging is not efficient or even impossible if there are only small differences between the temperatures inside and in the vicinity of the PHXs.Long periods of (dis)charging may lead to saturation in the vicinity of the PHXs.As a consequence, (dis)charging is no longer efficient and should be stopped since propagation of heat to regions away from the PHXs takes time.
The dynamics of the spatial temperature distribution is governed by a linear heat equation with convection and appropriate boundary and interface conditions.We solve that PDE using finite difference schemes; see Duffy [16].For the convection terms, we apply upwind techniques.In a first step, we study the semidiscretization with respect to spatial variables leading to a system of linear ODEs.In a second step, we consider full space-time discretization and derive implicit finite-difference schemes.The current paper provides the following theoretical contributions.First, we prove that the chosen semidiscretization ensures a system of linear ODEs with a stable system matrix.Second, we provide a detailed stability analysis for the implicit finite-difference schemes of the fully discretized PDE and establish a stability condition.
We perform numerical experiments, where simulation results for the temporal behavior of the spatial temperature distribution are used to determine how much energy can be stored in or taken from the storage within a given period of time.Special focus is laid on the dependence of these quantities on the arrangement of the PHXs within the storage.We refer to [17] in which we apply model reduction techniques known from control theory such as balanced truncation to derive low-dimensional approximations of aggregated characteristics of the temporal behavior of the spatial temperature distribution.The latter is crucial if the geothermal storage is embedded into a residential heating system and the cost-optimal management of such systems is studied mathematically in terms of optimal control problems.
The rest of the paper is organized as follows.In Section 2, we describe the dynamics of the spatial temperature distribution in the geothermal storage which is governed by a linear heat equation with a convection term and appropriate boundary and interface conditions.In Section 3, we present the semidiscretization with respect to spatial variables of the initial boundary value problem for that heat equation.For the resulting system of linear ODEs, we show that the system matrix is stable.The full space-time discretization is studied in Section 4, where we derive implicit finite-difference schemes and provide the associated stability analysis.In Section 5, we introduce aggregated characteristics of the spatio-temporal temperature distribution and explain their numerical approximation.Finally, Section 6 presents numerical results of some case studies.An appendix provides a list of frequently used notations, auxiliary results from matrix analysis, proofs, and some technical details that were removed from the main text.

DYNAMICS OF SPATIAL TEMPERATURE DISTRIBUTION IN A GEOTHERMAL STORAGE
In this section, we describe the dynamics of the spatial temperature distribution in a geothermal storage mathematically by a linear heat equation with convection term and appropriate boundary and interface conditions.

Two-dimensional model
We assume that the domain of the geothermal storage is a cuboid and consider a two-dimensional rectangular cross-section.We denote by T = T(t, x, ) the temperature at time t ∈ [0, t E ] in the point (x, ) ∈  = (0, l x ) × (0, l  ) with l x , l  denoting the width and height of the storage.The domain  and its boundary  are depicted in Figure 3.  is divided into three parts.The first is  M , which is filled with a homogeneous medium (soil) characterized by material parameters  M ,  M , and c M p denoting mass density, thermal conductivity, and specific heat capacity, respectively.The second is  F ; it represents the PHXs filled with a fluid (water) with constant material parameters  F ,  F , and c F p .The fluid moves with time-dependent velocity v 0 (t) along the PHX.For the sake of simplicity, we restrict to the case, often observed in applications, where the pumps moving the fluid are either on or off.Thus, the velocity v 0 (t) is piecewise constant taking values v 0 > 0 and zero, only.Finally, the third part is the interface  J between  M and  F .That interface is split into upper and lower interfaces  J and  J , respectively.Observe that we neglect modeling the wall of the PHX and suppose perfect contact between the PHX and the soil.Details are given below in (7) and (8).We summarize as follows: Assumption 1.
1. Material parameters of the medium  M ,  M , c M p in the domain  M and of the fluid  F ,  F , c F p in the domain  F are constants.

Heat equation
The temperature T = T(t, x, ) in the external storage is governed by the linear heat equation with convection term where ∇ = denotes the gradient operator.The first term on the right-hand side describes diffusion while the second represents convection of the moving fluid in the PHXs.Further, v = v(t, x, ) = v 0 (t)(v x (x, ), v  (x, )) ⊤ denotes the velocity vector with (v x , v  ) ⊤ being the normalized directional vector of the flow.According to Assumption 1, the material parameters , , c p depend on the position (x, ) and take the values  M ,  M , c M p for points in  M (medium) and  F ,  F , c F p in  F (fluid).Note that there are no sources or sinks inside the storage and therefore the above heat equation appears without forcing term.Based on this assumption, the heat equation ( 1) can be written as where Δ =  2 x 2 +  2  2 is the Laplace operator and a = a(x, ) is the thermal diffusivity which is piecewise constant with values a † =  †  † c † p with † = M for (x, ) ∈  M and † = F for (x, ) ∈  F , respectively.The initial condition T(0, x, ) = T 0 (x, ) is given by the initial temperature distribution T 0 of the storage.

Boundary and interface conditions
For the description of the boundary conditions, we decompose the boundary  into several subsets as depicted in Figure 3 representing the insulation on the top and the side, the open bottom, and the inlet and outlet of the PHXs.Further, we have to specify conditions at the interface between PHXs and soil.The inlet, outlet, and the interface conditions model the heating and cooling of the storage via PHXs.We distinguish between the two regimes "pump on" and "pump off."For for simplicity, we assume perfect insulation at inlet and outlet if the pump is off.This leads to the following boundary conditions.

• Homogeneous Neumann condition describing perfect insulation on the top and the side
where ×{l  } and  denotes the outer-pointing normal vector.
• Robin condition describing heat transfer on the bottom with  B = [0, l x ] × {0}, where  G > 0 denotes the heat transfer coefficient and T G (t) the underground temperature.
For more interpretation, we refer to Remark 2. • Dirichlet condition at the inlet if the pump is on (v 0 (t) > 0), that is, the fluid arrives at the storage with a given temperature T I (t).If pump is off (v 0 (t) = 0), we set a homogeneous Neumann condition describing perfect insulation.
• "Do Nothing" condition at the outlet in the following sense.If the pump is on (v 0 (t) > 0), then the total heat flux directed outwards can be decomposed into a diffusive heat flux given by  F T  and a convective heat flux given by v 0 (t) F c F p T. Since in real-world applications the latter is much larger than the first, we neglect the diffusive heat flux.This leads to a homogeneous Neumann condition If the pump is off, then we assume (as already for the inlet) perfect insulation, which is also described by the above condition.
• Smooth heat flux at interface  J between fluid and soil leading to a coupling condition Here, T F , T M denote the temperature of the fluid inside the PHX and of the soil outside the PHX, respectively.Moreover, we assume that the contact between the PHX and the medium is perfect which leads to a smooth transition of a temperature; that is, we have Remark 1.If the contact between the PHX and the medium is not perfect (e.g., in case of contact resistance), then the transition of the temperature at the interface  J will not be smooth, that is, T F ≠ T M .This leads to a temperature jump between the PHX and the medium.That phenomenon occurs in the heat transfer between the medium and an insulation as shown in [14,15].
Remark 2. Imposing the Robin condition (4) at the bottom boundary aims to mimic the thermal behavior at the bottom boundary.A more realistic description requires embedding the storage domain  into a larger computational domain including the surrounding regions as in Figure 2.This allows for warming and cooling in the vicinity of the storage resulting from the outflow and inflow of the storage heat.Contrary to that, condition (4) assumes an exogenously given underground temperature T G independent of the temperature in the storage.
The heat transfer coefficient  G describes the resistance to the heat flux at the boundary.For the limiting case  G → 0 we get a homogeneous Neumann condition, that is, perfect insulation, while in the limit for  G → ∞ condition ( 4) is the Dirichlet condition T = T G (t).The underground temperature in general shows seasonal fluctuations that can be described by where K G 1 is the intensity of the fluctuation, K G 2 is the average ground temperature, t 0 a time or phase shift, and t a the number of time units per year.Since our focus is on the short-term behavior, we assume in the sequel that the underground temperature is constant over time, that is, K G 1 = 0.

SEMIDISCRETIZATION OF THE HEAT EQUATION
This and the next section are devoted to the finite difference discretization of the heat equation ( 2) with the boundary and interface conditions (3) through (8).We proceed in two steps.In the first step, we apply semidiscretization in space and approximate only spatial derivatives by their respective finite differences.This approach leads to a high-dimensional system of ODEs with a stable system matrix for the temperatures at the grid points.The latter will be used as starting point for model reduction in our paper [17].In the second step (see Section 4), also time is discretized resulting in a family of implicit finite difference schemes for which we perform a stability analysis.

Semidiscretization of the heat equation
We now apply the finite difference method (see Duffy [16]) combined with upwind techniques for the convection terms for semidiscretization of the heat equation (2).
Let N x and N  be the number of grid points and h x = l x ∕N x and h  = l  ∕N  the step sizes in x-direction and -direction, respectively.The spatial domain is discretized by means of a mesh with grid points (x i ,   ) as in Figure 4 where We denote by T i (t) ≃ T(t, x i ,   ) the semidiscrete approximation of the temperature and by v 0 (t) which we identify with the corresponding sets of grid points.We denote by  S =  F ∪  M the set of grid points in the inner domain  S =  F ∪  M .In addition, we decompose the set of grid points on the interface  J =  J ∪  J between the fluid and medium into Here,  J and  J denote the lower and upper interface, respectively; see Figure 3. Further, we decompose the set   of grid points on the boundary domain  according to the decomposition of  given in Figure 3 into The spatial derivatives in the PDE (1) are approximated by linear combinations of values of T at the grid points (x i ,   ) in  S at time t.We use central second-order finite difference for the diffusion term: For the convection term, we use the upwind discretization to get We have to point out that the above upwind approximations of the convection terms need to be applied only to the set of grid points  F in the fluid domain  F , since there is no convection outside the fluid and we can set v x i = v  i = 0.For the sake of simplification and tractability of our analysis, we restrict ourselves to the following assumption on the arrangement of PHXs and impose conditions on the location of grid points along the PHXs.

Assumption 2.
1.There are n P ∈ N straight horizontal PHXs, the fluid moves in positive x-direction.2. The interior of PHXs contains grid points.

Each interface between medium and fluid contains grid points.
Then for grid points in the domain  S , the semidiscrete scheme is given by For grid points (i, ) ∈  F in the "fluid" domain  F Assumption 2 implies that v x i = 1 while v  i = 0 and the above coefficients are given by In the "medium" domain  M , the convection terms disappear and the coefficients of the scheme (9) become time-independent and are given for and . Note that for the neighboring grid points to the interfaces, we have to slightly modify the above scheme (9) due to the extra contribution from the interfaces; see Equations ( 16) and ( 17) below.

Semidiscretization of the boundary conditions
In this paragraph, we consider the discretization of boundary conditions.We start with the homogeneous Neumann conditions ( 3) and ( 6) for the top, left, right, and the outlet boundary, where the normal vector  is equal to (0, 1) ⊤ , (−1, 0) ⊤ , (1, 0) ⊤ , and (1, 0) ⊤ , respectively.Using first-order differences for the normal derivative, we obtain for all t ∈ [0, Next, we discretize the Robin condition (4) at the bottom boundary  B .We have  = (0, −1) ⊤ such that for all grid points (i, 0) ∈   B , we have for all t ∈ [0, t E ] On the inlet boundary  I , we have according to (5) a Dirichlet boundary condition during pumping and a Neumann condition if the pump is off.Then for all grid points (0, ) The relations ( 12) through ( 14) represent linear algebraic equations, which allow to express the grid values T i (t) in the boundary grid points (i, ) ∈   in terms of the corresponding values in the neighboring points in the interior of the domain and the input data to the boundary conditions.Thus, in the finite difference scheme, these values T i (t) can be removed from the set of unknowns.

Semidiscretization of interface condition
Now we consider grid points on the interface  J between fluid and medium, which are by Assumption 1 straight lines in x-direction.That interface can be decomposed as  J =  J ∪  J , with  J and  J representing the lower and upper interface, respectively; see Figure 5.We define the outer normal by  = (0, 1) ⊤ on the upper interface and by  = (0, −1) ⊤ for lower interface.Note that we have n P PHXs and each PHX has two interfaces.Then, we have in total 2n P interface subdomains.
For a grid point (x i ,   ) on the interface  J , the perfect contact condition ( 8) implies that at a given time t the temperature of the fluid T F (t, x i ,   ) is equal to the temperature T M (t, x i ,   ) of the medium at that point.As usual, T i (t) denotes the semidiscrete approximation of that temperature.Then discretization of the interface condition (7) leads to We obtain the following coupling between the grid values in an interface grid point (i, ) ∈  J and its neighbors in vertical direction a time t ∈ [0, t E ], The above relations show that the grid values T i (t) in the interface grid points (i, ) ∈  J can be expressed as linear combinations of the grid values in the two vertical neighboring points in the fluid and medium.Thus, in the finite difference scheme these values T i (t) can be removed from the set of unknowns.Now, let (i, ) ∈  J be an interface point on the lower interface.Then substituting the above expressions for T i (t) into the finite differences scheme (9) applied to the lower neighbor (i,  − 1) ∈  M in the medium leads to whereas for the upper neighbor (i,  + 1) ∈  F in the fluid, it holds with Similar expressions can be derived for points (i, ) ∈  J on the upper interface.

Matrix form of the semidiscrete scheme
We are now in a position to establish a semidiscretized version of the heat equation ( 2) in terms of a system of ODEs by summarizing relations ( 9), (16), and (17).To this end, we recall that the temperature at the boundary grid points can be obtained by the linear algebraic equations ( 12) through ( 14) derived from the boundary conditions.Further, the values at the interface points are obtained by the interpolation formulas in (15) derived from the perfect contact condition.Thus, we can exclude these grid points from the subsequent considerations where we collect the semidiscrete approximations of the temperature T(t, x i ,   ) at the remaining points of the grid in the vector function The enumeration of the entries of Y is such that we start with the first inner grid point (1, 1) next to the lower left corner of the domain.Then we number grid points consecutively in vertical direction where we exclude the 2n P points of the interfaces of the n P PHXs such that we have q = N  − 2n P − 1 points in each "column" of the grid.Thus, Y (i−1)q+1 corresponds to grid point (i, 1) for i = 1, … , N x − 1, and the last entry Y n to the inner grid point (N x − 1, N  − 1) next to the domain's upper right corner.The dimension of Y is n = (N x − 1)q = (N x − 1)(N  − 2n P − 1).The enumeration described above can be expressed formally by a mapping  ∶  S → {1, … , n} with (i, )  → l = (i, ) which maps pairs of indices (i, ) of grid point (x i ,   ) ∈  to the single index l of the corresponding entry in the vector Y .
Using the above notations we can rewrite relations ( 9), ( 16) and ( 17) as the following system of ODEs for the vector function Y representing the semidiscretized heat equation ( 2) together with the given boundary and interface conditions.
with the initial condition Y (0) =  0 .Here, the vector  0 ∈ R n contains the initial temperatures at the corresponding grid points with  0l = T 0 (x i ,   ) where l = (i, ), l = 1 … , n.The system matrix A results from the spatial discretization of the convection and diffusion term in the heat equation ( 2) together with the Robin and linear heat flux boundary conditions.It has tridiagonal structure consisting of (N x − 1) × (N x − 1) block matrices of dimension q given by The inner block matrices A M , i = 2, … N x − 2 of dimension q have tridiagonal structure and are sketched for the case of one PHX in the top row of Table 1.The matrix entries  F ,  F are given in (10),  M ,  M in (11),  M I ,  M I in (16), and  F I ,  F I in (17).The first and last diagonal entry reads as , respectively.They are obtained if the discretized top and bottom boundary conditions (12) and ( 13) are substituted into (9).The first and the last diagonal block matrices A L and A R ∈ R q×q of the matrix A given in (19) are also tridiagonal and sketched in the bottom row of Table 1.Its entries result from the discretization of boundary conditions at the left and right boundary.The entries in the first and last row are related to the inner grid points next to the four corners of the domain and obtained by substituting homogeneous Neumann condition (12) and Robin condition (13) into (9).A detailed derivation can be found in Appendix B of [18] and leads to Recall that  M ,  M ,  M are given by (11),  F± ,  F ,  F in (10).Further,  M and  F are given in (15) and the off-diagonal coefficients  M I ,  F I are given in ( 16) and ( 17).The block matrices on the subdiagonals D ± ∈ R q×q , i = 1, … , N x − 1, are diagonal matrices of the form where  M is given in (11) and  F± in (10).Here, we denote by | the location of the interfaces where we only sketched the case of one PHX.The n × 2 input matrix B is a result of the discretization of the inlet and Robin boundary conditions, its entries B lr , l = 1, … , n, r = 1, 2, are derived in Appendix B in [18] and are given by The entries for other l are zero.The input function g ∶ [0, t E ] → R 2 is defined by Recall that T I is the inlet temperature of the PHX during pumping and T G is the underground temperature.

Stability of matrix A
The finite difference semidiscretization of the heat equation ( 2) given by the system of ODEs ( 18) is expected to preserve the dissipativity of the PDE.This property is related to the stability of the system matrix A = A(t) in the sense that all eigenvalues of A lie in left open complex half plane.That property will play a crucial role for model reduction techniques for (18) based on balanced truncation which we study in [17].The next theorem confirms the expectations on the stability of A.
Theorem 1 (Stability of Matrix A).Under Assumption 1 on the model and Assumption 2 on the discretization, the matrix A = A(t) given in (19) is stable for all t ∈ [0, T], that is, all eigenvalues (A) of A lie in left open complex half plane.
Proof.Lemma 2 below shows by using Gershgorin's circle theorem (see Lemma 7 in Appendix B), that the eigenvalues are either located in left open complex half plane or zero.Further, Lemma 4 below shows that A(t) is nonsingular for all t ∈ [0, T] and thus excludes the case (A) = 0. Thus, for all eigenvalues it holds (A) ∈ C − and A is stable.□ We now state and prove the above mentioned two Lemmas 2 and 4 for which we introduce for a generic matrix Note that the maximum norm of M is given by ||M|| ∞ = max i S i (M).The quantities R i (M) appear as radii of Gershgorin's circles of M and the J i (M) are used to describe diagonal dominance of M. For the convenience of the reader, we provide in Appendix B Gershgorin's circle theorem and the notions and properties of diagonal dominant and strongly connected matrices known from matrix theory.We recall that the time-dependence of A(t) is a result of the discretization of convection terms in the heat equation ( 2).The latter depend on the time-dependent velocity v 0 (t) and with some abuse of notation we can write A = A(t) = A(v 0 (t)) and B = B(t) = B(v 0 (t)).Recall that we assume in Assumption 1 that v 0 (t) is piecewise constant with v 0 (t) = v 0 during charging and discharging when the pump is on whereas v 0 (t) = 0 if the pump is off.Therefore, the matrices A, B share this property.They take only the two values A P = A(v 0 ), B P = B(v 0 ) during pumping and A N = A(0), B N = B(0) if the pump is off.Thus, for studying properties of A(t) on [0, T] it is sufficient to look at the properties of A P and A N .
A closer look to the entries of the block matrices forming the system matrix A, that is, to A M , A L , A R given in Table 1 and to D ± given in (20), shows that for the diagonal entries and the row characteristics R i , J i given in (23) one has to distinguish 14 different cases.Instead of n rows it is sufficient to consider only 14 representative rows whose indices we denote by i l , l = 1, … , 14. Table C1 provides a list of diagonal entries A i l i l and the row characteristics R i l (A), J i l (A), S i l (A) in terms of the model and discretization parameters.

Lemma 1. The matrix A = A(t) is weakly diagonal dominant for all t ∈ [0, T].
Proof.Inspecting the quantities J i l (A) and Table C1 it can be seen that it holds J i l (A) ≥ 0, hence by Definition 2 in Appendix B the matrix is diagonal dominant.□ Note that A is weakly but not strictly diagonal dominant since not all of its rows are strictly diagonal dominant.

Lemma 2. The Gershgorin circles of the matrix A
Here, C − denotes the set of complex numbers with negative real part.
Proof.Let us examine Gershgorin's circles of A for the 14 different representative rows denoted by with centers C i l = A i l i l and the radii R i l (A), l = 1, … , 14, given in Table C1.Since all entries of A are real, the centers C i l = A i l i l < 0 of the discs are on the negative real axis.Lemma 1 shows that A is diagonal dominant, that is, Proof.Let (p, q) be a pair of distinct integers with p, q ∈ {1, … , n}.Then we can choose the sequence of distinct integers k 1 , k 2 , … , k m , such that m = |p − q| + 1 and k  = p +  − 1 (for p < q) and k  = p −  + 1 (for p > q).It holds A k  k +1 ≠ 0 since these entries are located on the upper and lower subdiagonal of A for which we have where  S is the set of grid points in the fluid and medium  F ∪  M and  J N the set of neighboring grid points to the interface.Since  F∕M given in (10) and (11) and given in ( 16) and ( 17) are positive, we have A k  k +1 ≠ 0,  = 1, 2, … , m.Thus, the matrix A is strongly connected.□ Lemma 4. The matrix A = A(t) is non-singular for all t ∈ [0, T].
Proof.From Lemma 1 and 3, it is known that A(t) is weakly diagonal dominant and strongly connected for all t ∈ [0, T].Table C1 shows that there exist strictly diagonal dominant rows.Hence, Better's corollary (see Lemma 9 in Appendix B) implies that A(t) is nonsingular.□

FULL DISCRETIZATION
After discretizing, the heat equation ( 2) w.r.t.spatial variables we will now also discretize the temporal derivative and derive an explicit and a family of implicit finite difference schemes for which we perform a stability analysis.

𝜃-Implicit finite difference scheme
We introduce the notation N  for the number of grid points in t-direction,  = t E ∕N  the time step, and Recall that Y contains the temperatures T = T(t, x, ) at the points of the grid excluding points on the boundary and interface.Discretizing the temporal derivative in (18) with the forward difference gives Substituting ( 24) into (18) and replacing the r.h.s. of ( 18) by a linear combination of the values at time t k and t k+1 with the weight  ∈ [0, 1] gives the following general -implicit finite difference scheme from which we derive for k = 0, … , N  − 1 the recursion with the initial value Y 0 = Y (0) and the notation I n for the n × n identity matrix.The above general -implicit scheme leads for  = 1∕2 and 1 to special cases which are known in the literature as Crank-Nicolson scheme and backward Euler or fully implicit scheme, respectively.The limiting case  = 0 is not an implicit but a fully explicit scheme also known as forward Euler scheme.

Stability of the finite difference scheme
In this subsection, we investigate the stability of the finite difference scheme (25) in the maximum norm and present in Theorem 2 below a stability condition to the time discretization.The use of such stability results is twofold.First, it ensures "robustness" w.r.t.round-off errors of the problems's input data, which are the initial condition and the inlet and underground temperature, in the sense that we can run the recursion for an arbitrarily long time without a total loss of accuracy.Second, stability of the scheme is a key ingredient in any analysis of convergence of the exact solution of the finite difference scheme to the exact solution of the given initial boundary value problem for the PDE for an infinite refinement of space and time discretization.
Note that a complete convergence analysis is beyond the scope of this paper.In particular, we do not investigate consistency issues.Consistency roughly says that the finite differences scheme approximates correctly the PDE.The proof of consistency is straightforward and based on Taylor series expansions.We refer to the Lax-Richtmyer equivalence theorem; see Sanz-Serna and Palencia [19], Theorem 2.5.3 in Thomas [20], stating that a consistent finite difference scheme for a well-posed linear initial boundary value problem is convergent if and only if it is stable.Hence, for a consistent scheme, convergence is synonymous with stability.
Our stability result is given in terms of maximum norms, which are defined for a vector and for a square matrix Definition 1 (Stability of difference scheme in the maximum norm).The finite difference scheme ( 25) is stable in the maximum norm if there exist constants C 0 , C g > 0 such that Theorem 2 (Stability of -implicit scheme).Under Assumption 1 on the model and Assumption 2 on the discretization, it holds 1.For  ∈ [0, 1), the semi-implicit finite difference scheme (25) is stable if the time step  satisfies the condition 2. For  = 1, the fully implicit finite difference scheme ( 25) is unconditionally stable, that is, stable for any  > 0.
The constants C 0 , C g in ( 26) can be chosen as The proof is based on the following two lemmas.

Lemma 5. For the maximum norm of the matrix A = A(t), it holds
Proof.A(t) is piecewise constant taking only the two values A P and A N .From the last column of Table C1 showing the 14 different row sums S i (A P∕N ) of the two matrices, it can be easily seen that yielding the estimate in the lemma.□ Lemma 6.Under Assumption 1 on the model and Assumption 2 on the discretization, it holds for all k = 0, … , N  − 1 and  ∈ [0, 1] that 1. the matrices G k+1 given in (25) are invertible and ||(G k+1 ) −1 || ∞ ≤ 1 with equality for  = 0; 2. the matrices H k given in (25) satisfy ), where  is given in (27); 3. the vectors F k given in (25) The proof of Lemma 6 can be found in Appendix D. It is based on the auxiliary results from matrix theory given in Appendix B and a closer inspection of the entries of the matrices G k , H k using the results from Table C1.
Proof of Theorem 2. From the invertibility of G k (see assertion 1 of Lemma 6) and the iteration of the recursion (25) we obtain for k = 1, … ,   the explicit representation where we define 0 ∏ =1 (•) = I n .Taking the maximum norm on both sides and applying the triangular and Cauchy-Schwarz inequality gives Substituting the estimates for where we used k ≤ N  = T.According to the second assertion of Lemma 6 the above estimate holds for all  > 0 if  = 1 and for  ≤ 1

AGGREGATED CHARACTERISTICS
The numerical methods introduced in Section 3 allow the approximate computation of the spatio-temporal temperature distribution in the geothermal storage.In many applications, it is not necessary to know the complete information about that distribution.An example is the management and control of a storage, which is embedded into a residential heating system.Here, it is sufficient to know only the response of a few aggregated characteristics of the temperature distribution to charging and discharging operations.These quantities can be computed via a post-processing procedure.In this section, we introduce some of these aggregated characteristics and describe their approximate computation based on the solution vector Y of the finite difference scheme.

Aggregated characteristics related to the amount of stored energy
We start with aggregated characteristics given by the average temperature in some subdomain of the storage, which are related to the amount of stored energy in that domain.
Let  ⊂  be a generic subset of the 2D computational domain.We denote by || = ∫∫  dxd the area of .Then W  (t) = l z ∫∫  c p T(t, x, )dxd represents the thermal energy contained in the 3D spatial domain is the gain of thermal energy during the period [t 0 , t 1 ].While positive values correspond to warming of , negative values indicate cooling and −G  (t 0 , t 1 ) represents the magnitude of the loss of thermal energy.For  =  † , † = M, F, we can use that the material parameters on  † equal the constants  =  † , c p = c † p .Thus, for the corresponding gain of thermal energy, we obtain ) , where denotes the average temperature in the medium ( † = M) and the fluid ( † = F), respectively.We denote by T S the average temperature in the whole storage.It can be obtained from T M and T F by Further, the total gain in the storage denoted by G S is obtained by G S = G S (t 0 , t 1 ) = G M (t 0 , t 1 ) + G F (t 0 , t 1 ).

Aggregated characteristics related to the heat flux at the boundary
Now we consider the convective heat flux at the inlet and outlet boundary and the heat transfer at the bottom boundary.
Let  ⊂  be a generic curve on the boundary, then we denote by || = ∫  ds the curve length.The rate at which the energy is injected or withdrawn via the PHX is given by ] , where is the average temperature at the outlet boundary.Here, we have used that in our model we have horizontal PHXs such that | I | = | O | and a uniformly distributed inlet temperature at the inlet boundary  I .Note that the fluid moves at time t with velocity v 0 (t) and arrives at the inlet with temperature T I (t) while it leaves at the outlet with the average temperature T O (t).For a given interval of time [t 0 , t 1 ], the quantity describes the amount of heat injected (G P > 0) to or withdrawn (G P < 0) from the storage due to convection of the fluid.Next, we look at the diffusive heat transfer via the bottom boundary and define the rate ) , where is the average temperature at the bottom boundary.

Numerical computation of aggregated characteristics
Approximations of the aggregated characteristics introduced above can be derived by the using finite difference approximations of the temperature T = T(t, x, ).They can be given in terms of the entries of the vector function Y (t) satisfying the system of ODEs (18) and containing the semidiscrete finite difference approximations of the temperature in the inner grid points of the computational domain .Recall that the temperatures on boundary and interface grid points can be determined by linear combinations of the entries of Y (t).
Let T † (t) with † = M, F, S, O, B be one of the average temperatures introduced above.Then the numerical approximation of the defining single and double integrals by quadrature rules leads to approximations by linear combinations of the entries of Y of the form where C † is some 1 × n-matrix.For the details we refer to subsection 4.2 and Appendix B in [21].The extension to approximations based of the solution of the fully discretized PDE ( 25) is straightforward using the relation

NUMERICAL RESULTS
In this section we present results of numerical experiments based on the finite difference discretization (25) of the heat Equation ( 2).We work with explicit schemes ( = 0).The advantage of an explicit scheme is that it avoids the time-consuming solution of systems of linear equations but one has to satisfy stronger conditions on the time step  to ensure stability of the scheme, see Theorem 2. We determine the spatio-temporal temperature distribution in the storage and study the impact of the heat exchanger pipe topology.In Section 6.2, we present results for a storage with two PHXs.
For these experiments we also compute and compare certain aggregated characteristics which are introduced in Section 5 and computed via post processing of the temperature distribution.Note that Subsections 5.2 and 5.4 in [21] contain additional results for a storage with one and three PHXs.

Experimental settings
The model and discretization parameters are given in Table 2.The storage is charged and discharged via PHXs filled with a moving fluid and thermal energy is stored by raising the temperature of the storage medium.We recall the open architecture of the storage which is only insulated at the top and the side but not at the bottom.This leads to an additional heat transfer to the underground for which we assume a constant temperature of T G (t) = 15 • C. In the simulations the fluid is assumed to be water while the storage medium is dry soil.During charging a pump moves the fluid with constant velocity v 0 arriving with constant temperature T I (t) = T I C = 40 • C at the inlet.If this temperature is higher than in the vicinity of the PHX, then a heat flux into the storage medium is induced.During discharging the inlet temperature is T I (t) = T I D = 5 • C leading to a cooling of the storage.At the outlet, we impose a vanishing diffusive heat flux, that is, during pumping, there is only a convective heat flux.We also consider waiting periods where the pump is off.This helps to mitigate saturation effects in the vicinity of the PHXs, which reduce the injection and extraction efficiency.During that waiting periods, the injected heat (cold) can propagate to other regions of the storage.Since pumps are off, we have only diffusive propagation of heat in the storage and the transfer over the bottom boundary.

Storage with two horizontal straight PHXs
In this experiment, we run the simulations with two horizontal PHXs located symmetrically to the vertical mid level of p = 50 cm and separated by a distance d varying between 10 and 90 cm.Note that in Subsection 5.2 in [21] we observed that placing only a single PHX at p = 50 cm shows a quite good performance.First, we study the spatial temperature distribution and some aggregated characteristics during (dis)charging for t E = 36 h.Then we introduce waiting periods allowing the injected heat (cold) to spread within the storage.It can be seen that warming and cooling in the left part of the storage is slightly stronger than in the right part.It mainly takes places in a vicinity of the PHX whereas after 36 h temperatures in more distant regions are only slightly changed.The spatial temperature distributions differ considerably for the three arrangements of two PHXs.For a small distance (d = 10 cm), we observe a strong saturation at a temperature level close to the inlet temperature in the small region between the PHXs while the region at the top is almost at the initial temperature and the region at the bottom is only slightly warmed (cooled) by the underground.For the PHXs at distance d = 90 cm, we observe an extreme saturation in the small layer between the upper PHX and the top boundary while the lower PHX is also warming (cooling) the underground.Next, we will have a look at aggregated characteristics.In the left panel of Figure 7, the average temperature in the storage T S during charging is plotted against time for distances of the PHXs d = 10,20, … , 90 cm.The right panel presents the gain G S and loss −G S of thermal energy in the storage at the end of the charging and discharging period, respectively.The results reveal that apart from the first 4 h, there is a strong impact of the PHX distance.The most efficient mode of operation is obtained for the PHXs distance of d = 40 cm.Here, the gain (loss) of thermal energy during charging (discharging) is at maximum.These quantities strongly decay for smaller and larger distances because of the saturation effect which becomes stronger if PHXs are arranged closer to each other or closer to the top and bottom boundary of the storage.

Charging and discharging with waiting periods
The above experiments have shown how saturation effects can be mitigated by choosing an appropriate vertical distance of the two PHXs.This option is only available in the design of the storage architecture and not during the operation of an already existing storage.Therefore, we now want to examine another option, which is the interruption of (dis)charging cycles allowing the heat injected to (extracted from) the vicinity of the PHXs to propagate to the other storage regions.The idea is that after a sufficiently long waiting period the saturation in the vicinity of the PHX is considerably reduced such that (dis)charging can resumed with higher efficiency.Although the introduction of such waiting periods will increase the time needed to inject (extract) a given amount of thermal energy, it reduces the saturation effect and helps to save operational costs for electricity used for running the pumps.
In our experiments, we divide the time interval [0, t E ] into three subintervals of length 8, 12, and 16 h.In each subinterval, (dis)charging is followed by a waiting period of the same length as it can be seen in Figure 8 where charging, waiting, and discharging periods are represented by red, green, and blue background color.The top panels show the average temperatures in the storage T S and at the outlet T O , respectively, during charging and discharging.We compare storages with two PHXs of distance d = 10,40, 90 cm.Recall that in the previous subsection we have seen that d = 40 cm allows for a much more efficient operation than for d = 10,90 cm.As expected, during the waiting periods, the average temperatures at the outlet and in the PHX decay after charging and rise after discharging.This is due to the diffusion of heat in the storage, in particular the heat flux induced by the different temperatures inside and outside the PHX.During waiting, the average temperature in the storage T S is almost constant since injection or extraction of heat is stopped.However, the heat transfer to and from the underground at the bottom boundary continues also during waiting but the waiting periods are too short to produce a visible change of T S .In the two lower panels of Figure 8, we compare the storage operation with and without waiting periods.We plot the gain G S (loss −G S ) of thermal energy in the storage during charging (discharging) over time.Note that for operation with waiting (dis)charging takes place only 50% of the time.However, for the "optimal" PHX distance d = 40 cm, the resulting gain (loss) reaches more than 80% of the values for uninterrupted operation.For the less efficient cases of PHXs at distances d = 10 cm and d = 90 cm that cause strong saturation effects the differences are smaller and the gaps are quickly reduced to almost zero after resuming (dis)charging.

CONCLUSION
We have investigated mathematical modeling and numerical simulations for the dynamics of spatial temperature distribution in a geothermal energy storage and the associated input-output behavior.The underlying initial boundary value problem for the heat equation with a convection term has been discretized using finite difference schemes.In a first step, we studied the semidiscretization with respect to spatial variables.For the resulting system of linear ODEs, we proved that the system matrix is stable.In a second step, the full space-time discretization has been considered.Here, we derived explicit and implicit finite-difference schemes and investigated associated stability problems.Finally, some numerical experiments showing the dependence of the charging and discharging efficiency on the topology and arrangement of two PHXs have been presented.Our numerical experiments show how these simulations can support the design and operation of a geothermal storage.Examples are the dependence of the charging and discharging efficiency on the topology and arrangement of PHXs and on the length of charging and discharging periods.
Based on the findings of this paper, we study in [17] model reduction techniques to derive low-dimensional approximations of aggregated characteristics of the temperature distribution describing the input-output behavior of the storage.The latter is crucial if the geothermal storage is embedded into a residential heating system and the cost-optimal management of such systems is studied mathematically in terms of optimal control problems as in [22].
Matrices that are weakly but not strictly diagonal dominant can be singular.A criterion for nonsingularity is based on the following property of a matrix and the subsequent lemma.That property was introduced in Definition 6.2.7 in Horn and Johnson [25] and termed "property SC."In the literature it is also known as "strongly connected."Definition 3 (Strongly connected matrix).A matrix M ∈ C n×n is called strongly connected (or of property SC) if for each pair of distinct integers p, q ∈ {1, … , n} there is a sequence of distinct integers k 1 = p, k 2 , … , k m = q such that each entry M k 1 k 2 , M k 2 k 3 , … ,M k m−1 k m is nonzero.
For strongly connected matrices, Corollary 6.2.9 in Horn and Johnson [25] gives the following criterion for nonsingularity.
Lemma 9 (Better's corollary).Suppose that the matrix M ∈ C n×n is strongly connected, weakly diagonally dominant and there exists one strictly diagonal dominant row.Then M is nonsingular.

APPENDIX C: ROW CHARACTERISTICS OF MATRIX A
We give in Table C1 for the system matrix A the centers C i l = A i l i l and radii R i l (A) of the Gershgorin circles, the differences J i l (A) as well as the row sums S i l (A) for the 14 representative rows of A. )  (25), it follows that for k = 0, … , N  − 1

TABLE C1 Diagonal entries A
where we have used that B k = B(k) takes only the values B P and B N .□

FIGURE 1 FIGURE 2
FIGURE 1 Geothermal storage: in the new building, under a building (left) and in the renovation, aside of the building (right), see www.ezeit-ingenieure.eu,www.geo-ec.de.[Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 3
FIGURE 3 Two-dimensional model of the geothermal storage: decomposition of the domain  and the boundary .[Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 5
FIGURE 5 Interfaces between the fluid and soil.[Colour figure can be viewed at wileyonlinelibrary.com]

Lemma 3 .
of the Gershgorin circles never exceed |C i l | and it holds D i l ⊂ C − ∪ {0}.□ The matrix A = A(t) is strongly connected for all t ∈ [0, T].

Figure 6
Figure 6 shows for three different distances d of the two PHXs the spatial distribution of the temperature in the storage after 36 h of charging (left) and discharging (right).In the top panels the PHXs are very close (d = 10 cm).The panels in

FIGURE 6 FIGURE 7
FIGURE 6 Spatial distribution of the temperature in the storage with two horizontal PHXs of vertical distance d after 36 h of charging (left) and discharging (right).Top: d = 10 cm.Middle: d = 40 cm.Bottom: d = 90 cm.[Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 8
FIGURE 8 Charging and discharging during 36 h with several waiting periods for a storage with two horizontal PHXs at distance d = 10, d = 40, and d = 90 cm.Top: Aggregated characteristics T S and T O .Bottom: Gain G S / loss −G S of stored energy.Left: Charging.Right: Discharging.[Colour figure can be viewed at wileyonlinelibrary.com] of Gershgorin circles), radii of Gershgorin circles R i l , differences J i l and row sums S i l of the matrix A, l =

TABLE 1
Sketch of block matrices for the case of one PHX.Top: inner block matrices A M .Bottom: matrices A L and A R .

TABLE 2
Model and discretization parameters.