{'fraction_non_alphanumeric': 0.03309560060902316, 'fraction_numerical': 0.009275582979405401, 'mean_word_length': 4.273851030110935, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Short-Term Behavior of a Geothermal Energy Storage: Numerical Applications

This paper is devoted to numerical simulations of the short-term behavior of the spatial temperature distribution in a geothermal energy storage. Such simulations are needed for the optimal control and management of residential heating systems equipped with an underground thermal storage. We apply numerical methods derived in our companion paper [15] in which we study the governing initial boundary value problem for a linear heat equation with convection. Further, we perform extensive numerical experiments in order to investigate properties of the spatio-temporal temperature distribution and of its aggregated characteristics.


Introduction
This paper is devoted to the computation of the spatial temperature distribution in a geothermal energy storage for short periods of time. We focus on underground thermal storages as depicted in Fig. 1.1 which can be found in heating systems of single buildings as well as of district heating systems. Such storages have gained 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 in new buildings and in renovations. Such facilities are used to mitigate and to manage temporal fluctuations of heat supply and demand and to move heat demand through time. It is well-known that 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 [8], Kitapbayev et al. [10]). 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.
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. [5] 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 multi-physics 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. [12] 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. [13] study the numerical simulation of aquifer thermal energy storages and focus on dynamic mesh optimisation for the finite element solution of the heat and flow equations. For an overview on thermal energy storages we refer to Dincer and Rosen [6] and for further contributions on the numerical simulation of such storages to [2,6,9,11,14,17]. This paper is based on our companion paper [15] where we give a detailed description of the mathematical model of an underground thermal energy storage and the derivation and theoretical justification of the numerical methods. The starting point is a 2D-model, see Fig. 1.2. A defined volume under or aside of a building is filled with soil and insulated to the surrounding ground. The storage is charged and discharged via pipe heat exchangers (PHXs) filled with some fluid (e.g. water). Thermal energy is stored by raising the temperature of the soil inside the storage. A special feature of the storage is its open architecture at the bottom. There is no insulation such that thermal energy can also flow into deeper layers as it can be seen in Fig. 1.2. This leads to 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.
A similar model has been already considered in Bähr et al. [4,3] where the authors focus on the numerical simulation of the long-term behavior over weeks and months of the spatial temperature distribution and the interaction between a geothermal storage and its surrounding domain. For the sake of simplicity the charging and discharging process using PHXs was not modeled in detail but described by a source term. In this work, we focus on the short-term behavior of the spatial temperature distribution. We believe that this is interesting 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. Contrary to [4,3] we include PHXs for a more realistic model of the storage's charging and discharging process. However, we do not consider the surrounding medium but reduce the computational domain to the storage depicted in Fig. 1.2 by a dashed black rectangle. Instead we set appropriate boundary conditions to mimic the interaction between storage and environment.
The temporal evolution of the spatial temperature distribution is governed by a linear heat equation with convection and appropriate boundary and interface conditions. A numerical solution of that PDE using finite difference schemes is sketched in Sec. 3. For more details we refer to our paper [15]. Management and control of a storage that is embedded into a residential system usually does not require the complete spatio-temporal temperature distribution but is based only on certain aggregated characteristics that can be computed in a post-processing step as explained in Sec. 4. Examples are the average temperatures in the storage medium, in the PHX fluid, at the outlet of the PHXs and at the storage's bottom boundary, respectively. From these quantities one can derive the amount of available thermal energy that can be stored in or extracted from the storage in a given short period of time.
In Sec. 5 we present results of extensive numerical experiments where we use simulations results for the temporal behavior of the spatial temperature distribution to determine how much energy can be stored in or taken from the storage within a given short period of time. Special focus is laid on the dependence of these quantities on the arrangement of the PHXs within the storage.
In another companion paper [16] we apply model reduction techniques known from control theory such as Lyapunov balanced truncation to derive low-dimensional approximations of the above mentioned aggregated characteristics. The latter is crucial if the costoptimal management of residential heating systems equipped with a geothermal storage is studied mathematically in terms of optimal control problems. It is well-known that most of the model reduction methods are developed for linear time-invariant (LTI) systems. However, the heat equation (2.1) which we derive in Sec. 2 contains a convection term that is driven by the velocity of the fluid in the PHXs. That velocity is time-dependent and typically piecewise constant during waiting, charging and discharging periods. Therefore, we are not in the framework of LTI systems and propose in Sec. 6 an LTI analogous model that mimics the most important features of the original non-LTI model of the geothermal storage.
The rest of the paper is organised as follows. In Sec. 2 we derive a linear heat equation with a convection term and appropriate boundary and interface conditions describing the dynamics of the spatial temperature distribution in the geothermal storage. Sec. 3 is devoted to the finite difference discretization of the heat equation. In Sec. 4 we introduce aggregated characteristics of the spatio-temporal temperature distribution and explain their numerical approximation. Sec. 5 presents numerical results of extensive case studies. We provide additional video material showing animations of the temporal evolution of the spatial temperature distribution in the geothermal storage. The videos are available at www.b-tu.de/owncloud/s/D68fmqXRcgbesKj . Finally, in Sec. 6 we derive an LTI analogous model of the geothermal storage and present some numerical results. An appendix provides a list of frequently used notations and some auxiliary results removed from the main text.

Dynamics of Spatial Temperature Distribution in a Geothermal Storage
The setting is based on our paper [15,Sec. 2]. For self-containedness and the convenience of the reader, we recall in this section the description of the model. The dynamics of the spatial temperature distribution in a geothermal storage can be described mathematically by a linear heat equation with convection term and appropriate boundary and interface conditions. We denote by Q the temperature in the geothermal storage depending on time as well as on the location in the storage.

2D-Model
We assume that the domain of the geothermal storage is a cuboid and consider a twodimensional rectangular cross-section. We denote by Q = Q(t, x, y) the temperature at time t ∈ [0, T ] in the point (x, y) ∈ D = (0, l x ) × (0, l y ) with l x , l y denoting the width and height of the storage. The domain D and its boundary ∂ D are depicted in Fig. 2.1. D is divided into three parts. The first is D M and 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 D F , it represents the PHXs filled with a fluid (water) with constant material parameters ρ F , κ F and c F p . The fluid moves with timedependent velocity v 0 (t) along the PHX. For the sake of simplicity we restrict ourselves 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 D J between D M and D F . That interface is split into upper and lower interfaces D J and D 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 in (2.4) and (2.5) below. The above can be summarized in the following  2. Fluid velocity is piecewise constant, i.e. v 0 (t) = v 0 > 0, pump on, 0, pump off.
3. Perfect contact at the interface between fluid and medium.
Remark 2.2 Results obtained for our 2D-model, where D represents the rectangular crosssection of a box-shaped storage can be extended to the 3D-case if we assume that the 3D storage domain is a cuboid of depth l z with a homogeneous temperature distribution in zdirection. A PHX in the 2D-model then represents a horizontal snake-shaped PHX densely filling a small layer of the storage.
Heat equation. The temperature Q = Q(t, x, y) in the external storage is governed by the linear heat equation with convection term where ∇ = ∂ ∂ x , ∂ ∂ y 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.
, v y (x, y)) denotes the velocity vector with (v x , v y ) being the normalized directional vector of the flow. According to Assumption 2.1 the material parameters ρ, κ, c p depend on the position (x, y) and take the values ρ M , κ M , c M p for points in D M (medium) and ρ F , κ F , c F p in D 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 (2.1) can be written as ∂ y 2 is the Laplace operator and a = a(x, y) is the thermal diffusivity which is piecewise constant with values a † = κ † ρ † c † p with † = M for (x, y) ∈ D M and † = F for (x, y) ∈ D F , respectively. The initial condition Q(0, x, y) = Q 0 (x, y) is given by the initial temperature distribution Q 0 of the storage.

Remark 2.3
In real-world geothermal storages PHXs are often designed in a snake form located in the storage domain at multiple horizontal layers. There may be only a single inlet and a single outlet. We will mimic that design by a computationally more tractable design characterized by multiple horizontal straight PHXs as it is sketched in Fig. 1.2. This allows to control the PHXs in different layers separately. For a topology with single inlet and outlet snake-shaped PHXs the outlet of a straight PHX in one layer can be connected with the inlet of the straight PHX in the next layer.

Boundary and Interface Conditions
For the description of the boundary conditions we decompose the boundary ∂ D into several subsets as depicted in Fig. 2.1 representing the insulation on the top and the side, the open bottom, 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 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 -Robin condition describing heat transfer on the bottom where λ G > 0 denotes the heat transfer coefficient and Q G (t) the underground temperature.
-Dirichlet condition at the inlet if the pump is on (v 0 (t) > 0), i.e. the fluid arrives at the storage with a given temperature Q 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 k f ∂ Q ∂ n and a convective heat flux given by v 0 (t)ρ F c F p Q. 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 D J between fluid and soil leading to a coupling condition Here, Q F , Q 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, i.e., we have

Discretization of the Heat Equation
We now sketch the discretization of the heat equation (2.1) together with the boundary and interface conditions given in (2.2) through (2.5). For details we refer to our paper [15,Sec. 3 and 4]. We proceed in two steps. In the first step we apply semi-discretization in space and approximate only spatial derivatives by their respective finite differences. This approach is also known as 'method of lines' and leads to a high-dimensional system of ODEs for the temperatures at the grid points. In the second step also time is discretized resulting in an implicit finite difference scheme.

Semi-Discretization of the Heat Equation
The spatial domain depicted in Fig. 2.1 is discretized by the means of a mesh with grid points (x i , y j ) as in Fig. 3 Here, N x and N y denote the number of grid points while h x = l x /N x and h y = l y /N y are the step sizes in x and y-direction, respectively. We denote by 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. 3. Each interface between medium and fluid contains grid points.
We approximate the spatial derivatives in the heat equation, the boundary and interface conditions by finite differences as in [15, Subsec. 3.1-3.3] where we apply upwind techniques for the convection terms. The result is the system of ODEs (3.1) (given below) for a vector function Y : [0, T ] → R n collecting the semi-discrete approximations Q i j (t) of the temperature Q(t, x i , y j ) in the "inner" grid points, i.e., all grid points except those on the boundary ∂ D and the interface D J . For a model with n P PHXs the dimension of Y is n = (N x − 1)(N y − 2n P − 1), see [15].
Using the above notation the semi-discretized heat equation (2.1) together with the given initial, boundary and interface conditions reads as with the initial condition Y (0) = y 0 where the vector y 0 ∈ R n contains the initial temperatures Q(0, ·, ·) at the corresponding grid points. The system matrix A results from the spatial discretization of the convection and diffusion term in the heat equation (2.1) together with the Robin and linear heat flux boundary conditions. It has the tridiagonal structure As a result of the discretization of the Dirichlet condition at the inlet boundary and the Robin condition at the bottom boundary, we get the function g : [0, T ] → R 2 called input function and the n × 2 input matrix B called input matrix. The entries of the input matrix B lr , l = 1, . . . , n, r = 1, 2, are derived in [15,Subsec 3.4] and are given by The entries for other l are zero. Here, K denotes the mapping (i, j) → l = K(i, j) of pairs of indices of grid point (x i , y j ) ∈ D to the single index l ∈ {1, . . . , n} of the corresponding entry in the vector Y . The input function reads as Recall that Q I is the inlet temperature of the PHX during pumping and Q G is the underground temperature.

Full Discretization
After discretizing the heat equation (2.1) w.r.t. spatial variables we will now also discretize the temporal derivative and derive a family of implicit finite difference schemes. We introduce the notation N τ for the number of grid points in the t-direction, τ = T /N τ the time step and t k = kτ, k = 0, ..., , v 0 (t) at time t = t k . Further, we denote by Y k = (Y k 1 , . . . ,Y k n ) the discrete-time approximation of the vector function Y (t) at time t = t k . Discretizing the temporal derivative in (3.1) with the forward difference gives Substituting (3.4) into (3.1) and replacing the r.h.s. of (3.1) by a convex combination of the values at time t k and t k+1 with the weight θ ∈ [0, 1] gives the following general for which we provide in our paper [15,Sec. 4] a detailed stability analysis. For our numerical experiments in Sec. 5 we use an explicit scheme which is obtained for θ = 0 and given by the recursion as with the initial value Y 0 = Y (0) and the notation I n is the n × n identity matrix. 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. In [15,Theorem 4.2], we show that the above explicit scheme is stable if the time step τ satisfies the condition

Aggregated Characteristics
The numerical methods introduced in Sec. 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 a few aggregated characteristics of the temperature distribution which can be computed via post-processing. 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 given subdomain of the storage which are related to the amount of stored energy in that domain. Let B ⊂ D be a generic subset of the 2D computational domain. We denote by |B| = B dxdy the area of B. Then W B (t) = l z B ρc p Q(t, x, y)dxdy represents the thermal energy contained in the 3D spatial domain For B = D † , † = M, F, we can use that the material parameters on D † equal the constants ρ = ρ † , c p = c † p . Thus, for the corresponding gain of thermal energy we obtain denotes the average temperature in the medium ( † = M) and the fluid ( † = F), respectively. We denote by Q S the average temperature in the whole storage. It can be obtained from Q M and Q F by Further, the total gain in the storage denoted by G S is obtained by

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 C ⊂ ∂ D be a generic curve on the boundary, then we denote by |C| = C ds the curve length. The rate at which the energy is injected or withdrawn via the PHX is given by is the average temperature at the outlet boundary. Here, we have used that in our model we have horizontal PHXs such that |∂ D I | = |∂ D O | and a uniformly distributed inlet temperature at the inlet boundary ∂ D I . Note that the fluid moves at time t with velocity v 0 (t) and arrives at the inlet with temperature Q I (t) while it leaves at the outlet with the average temperature Q 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 is the average temperature at the bottom boundary. Note that the second equation in the first line follows from the Robin boundary condition. The quantity describes the amount of heat transferred via the bottom boundary of the storage.

Energy Balance
In our model we assume perfect thermal insulation at all boundaries except the inlet, outlet and the bottom boundary. At the outlet we impose a homogeneous Neumann condition describing zero diffusive heat transfer. At the inlet we also have a zero diffusive heat transfer under the reasonable assumption that the temperature in the supply pipe is constant and equals Q I (t), thus the normal derivative ∂ Q ∂ n is zero. This implies that gains and losses of thermal energy in the storage are caused either by injections or withdrawals via the PHXs or by heat transfer via the open bottom boundary. Thus, we can decompose the total gain G S to obtain the following energy balance

Numerical Computation of Aggregated Characteristics
In this subsection we consider the approximate computation of aggregated characteristics introduced in the previous subsections by using finite difference approximations of the temperature Q = Q(t, x, y). The approximations are given in terms of the entries of the vector function Y (t) satisfying the system of ODEs (3.1) and containing the semi-discrete finite difference approximations of the temperature in the inner grid points of the computational domain D. Recall that the temperatures on boundary and interface grid points can be determined by linear combinations from the entries of Y (t). The extension to approximations based of the solution of the fully discretized PDE (3.5) is straightforward using the relation Let us start with the average temperatures Q M and Q F , where the temperature Q(t, x, y) is averaged over unions of disjoint rectangular subsets of the computational domain D. Assume that B ⊂ D is a generic rectangular subset with corners defined by the grid points We assume further that the domain B contains at least one layer of horizontal and vertical inner grid points, respectively. Thus we require i + 2 ≤ i and j + 2 ≤ j. We denote by x, y)dxdy the average temperature in B. Rewriting the double integral as two iterated single integrals and applying trapezoidal rule to the single integrals the average temperature Q B can be approximated by (for details see Appendix B.1) where N B = {(i, j) : i = i, . . . , i, j = j, . . . , j} and the coefficients d i j of the above quadrature formula are given by (boundary points, except corners) 1 4 , for i = i, i, j = j, j (corner points).
and we can rewrite approximation (4.4) as with an 1×n−matrix D and an 1×n−matrix D, whose entries are given for l = 1, . . . , n, l = 1, . . . , n by respectively. Finally, substituting Y = CY into (4.6) yields a representation of the average temperature Q B as a linear combination of entries of the vector Y which reads as Based on the above representation we can derive similar approximations for the average temperatures Q M and Q F in the medium and the fluid, respectively. Our model assumptions imply that for a storage with n P PHXs the domain D F splits into n P disjoint rectangular subsets D F j , j = 1, . . . , n P (PHXs), whereas D M consists of n P + 1 of such subsets between the PHXs and the top and bottom boundary of D which we denote by D M j , j = 0, . . . , n P . Then we can apply (4.4) to derive the approximation An approximation of the form Q M ≈ C M Y can be obtained analogously. Further, from Eq. (4.1) the approximation for the average temperature in the whole storage can be derived as In Appendix B.2 we derive approximations Q O ≈ C O Y and Q B ≈ C B Y for the average temperatures at the outlet and the bottom boundary, respectively. Here, the line integrals in the definitions (4.2) and (4.3) of these two characteristics are approximated by trapezoidal rule.

Numerical Results
In this section we present results of numerical experiments based on the finite difference discretization (3.5) of the heat equation (2.1). We determine the spatio-temporal temperature distribution in the storage. Further, we study the impact of the PHX topology and vary the number and arrangement of the PHXs. In Subsecs. 5.2, 5.3 and 5.4 we present results for a storage with one, two and three PHXs, respectively. For these experiments we also compute and compare certain aggregated characteristics which are introduced in Sec. 4 and computed via post processing of the temperature distribution.
Note that we provide additional video material showing animations of the temporal evolution of the spatial temperature distribution for which in the following we can present snapshots only. The videos are available at www.b-tu.de/owncloud/s/D68fmqXRcgbesKj .

Experimental Settings
The model and discretization parameters are given in Table 5.1. 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 Q 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 Q I (t) = Q 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 Q I (t) = Q I D = 5°C leading to a cooling of the storage. At the outlet we impose a vanishing diffusive heat flux, i.e. 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 reducing 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 One Horizontal Straight PHX
In this experiment we run simulations with one horizontal PHX located at different vertical positions p between the bottom (p = 0 cm) and the top (p = l y = 100 cm) of the storage. We compare the spatial temperature distributions as well as aggregated characteristics such as the average temperature in the storage Q S (t), the average outlet temperature Q O (t), and the gain or loss of energy G S (0, T ) in the storage during a period of T = 36 hours. Charging is realized by sending fluid through the PHX for 36 hours. It arrives at the inlet with constant temperature Q I C (t) = 40°C. We start with an initial temperature Q(0, x, y) = 10°C, uniformly distributed in the storage. In the experiment with discharging we start with an uniformly distributed initial temperature 35°C. For 36 hours the storage is cooled by the moving fluid arriving at the storage inlet with constant temperature Q I D (t) = 5°C.   domains is only slightly changed. Due to the direction of the moving fluid from left to right, warming and cooling in the left part of the storage is slightly stronger than in the right part. A closer inspection of the results shows that except in the experiment with the PHX close to the bottom boundary (p = 10 cm), after 36 hours of charging the temperatures in the vicinity of that boundary are below the underground temperature Q G = 15°C. Thus in addition to the injection of heat via the PHX we also have an inflow of thermal energy from the warmer underground into the storage. This results in a "boundary layer" which is slightly warmer than in the inner storage region. The reverse effect can be observed during discharging where close to the bottom boundary the temperature is always above Q G = 15°C. This induces a heat flux from the storage to the colder underground which contributes together with the extraction of heat via the PHX to the total loss of thermal energy in the storage. In Fig. 5.3 we plot in the upper panels the gain G S (respectively loss −G S ) of thermal energy during 36 hours of charging (respectively discharging) against time for vertical positions p = 10, 20, . . . , 90 cm. The lower panel shows these quantities at the end of the 36 hour charging and discharging period, depending on the vertical PHX position p. In the  Further, for all p we observe a decaying slope of the curves in the upper plots. This can be explained by the "thermal saturation" in the vicinity of the PHX and the slow diffusive propagation of the heat to the more distant regions of the storage. It shows that (dis)charging the storage becomes less efficient after longer periods of operation. Injecting (extracting) a certain amount of energy takes longer and needs more electricity consumed by the pumps. This effect suggests to interrupt (dis)charging and include waiting periods in which the heat (cold) in the vicinity of the PHXs can propagate to other regions of the storage. The impact of such waiting periods will be studied in more detail in Subsec. 5.3. we observe remarkable deviations. Here charging and discharging is considerably slower and gains and losses of thermal energy are smaller. For a PHX close to the top this can be explained by the saturation of the storage domain in the vicinity of the PHX. During charging (discharging) the boundary and its insulation prevent the propagation of heat into (from) the inner storage regions. On the bottom boundary that effect is combined with heat transfer to the underground. During charging a part of the injected heat is lost to the underground while during discharging the vicinity of the PHX is also cooled by the colder underground. Thus as expected, for an efficient operation of the storage the PHX should be located in the central region of the storage.

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 cm and 90 cm. Recall that placing a single PHX at p = 50 cm showed quite good performance in the last subsection. First we study the spatial temperature distribution and some aggregated characteristics during (dis)charging for T = 36 hours. Then we introduce waiting periods allowing the injected heat (cold) to spread within the storage.    = 90 cm). As in the experiment with only one PHX 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. Thus, 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 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.

Charging and Discharging Without Waiting Periods
Next we will have a look at aggregated characteristics. In Fig. 5.5 the average temperatures in the storage Q S are plotted against time for distances of the PHXs d = 10, 20, . . . , 90 cm. Fig. 5.6 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 figures reveal that apart from the first 4 hours 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 period 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 ] into three subintervals of length 8, 12, 16 hours. In each subinterval (dis)charging is followed by a waiting period of the same length as it can be seen in Fig. 5.7 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 Q S and at the outlet Q O , respectively, during charging and discharging. We compare a storage with two PHXs of distance d = 40 cm and a storage with more close-by PHXs d = 10 cm and two PHXs at distance d = 90 cm. Recall that in the previous subsection we have seen that d = 40 cm allows for 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 Q 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 Q S . In the two lower panels of Fig. 5.7 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 distance d = 10 cm and PHXs at distance 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. In this example we add a third PHX to the storage architecture and study two different PHX arrangements. We proceed with the experimental design including the same waiting periods considered in the previous subsection but now we "glue" together the two periods of charging and discharging each of length 36 h. The result is a total period of length T = 72 h starting with a storage at temperature Q(0, x, y) = 10°C. Within the the first 36 hours the storage is charged by the moving fluid arriving at the PHX inlet with temperature Q I C (t) = 40°C. In the second 36 hours it is discharged using the inlet temperature Q I D (t) = 5°C. The charging, waiting and discharging periods can be seen in Fig. 5.9. Contrary to the above experiments, discharging now starts not with a temperature 35°C but with a nonuniformly temperature distribution which is obtained after 36 h of charging (and waiting). In this more realistic setting, temperatures typically are higher in the vicinity of the PHXs and lower in other regions. Fig. 5.8 shows snapshots of the spatial temperature distribution during the last charging period (at t = 27h), during the subsequent waiting period (at t = 35h) and during of the last discharging period (at t = 63h), respectively. We compare two storage architectures with three PHXs. In the first, the PHXs are located symmetrically w.r.t. the vertical mid level. For the second, the central PHX was moved upwards such that we get a non-symmetric arrangement with two quite close-by PHXs in the upper region. The snapshots show a strong saturation between the two upper PHXs of the non-symmetric PHX arrangement while for symmetric PHXs the temperature distribution is much more uniform, in particular during the waiting period as it can be seen in the middle panel for time t = 35 h.

Storage With Three Horizontal Straight PHXs
In Fig. 5.9 we present aggregated characteristics which are plotted over time and observe similar patterns as in the experiment with a two PHX storage considered in the previous subsection. During the waiting periods after charging the average outlet and PHX temperatures decay at a faster rate for symmetric PHXs than for non-symmetric PHXs. Vice versa they increase faster in waiting periods after discharging. This is a consequence of the stronger saturation for non-symmetric PHXs which prohibits a faster cooling (warming) of the PHX during waiting. For symmetric PHXs the average storage temperature during charging increases faster and during discharging decreases faster than for non-symmetric PHXs. This explains the similar patterns for the gain of stored energy which are plotted in the right panel. It shows that the storage with symmetric PHXs (dis)charges faster than the storage with non-symmetric PHXs.

Analogous Linear Time-Invariant System
This section is motivated by our companion paper [16] in which we aim to approximate the dynamics of certain aggregated characteristics for the infinite dimensional spatial distribution of the temperature Q = Q(t, x, y) describing the storage's input-output behavior by a low-dimensional system of ODEs. Recall that the dynamics of the spatial distribution of Q is governed by the heat equation (2.1). We applied semi-discretization to that PDE and obtained the finite-dimensional approximation (3.1) which reads as dY (t) dt = A(t)Y (t)+B(t)g(t) and constitutes a high-dimensional system of ODEs for the vector function Y containing the temperatures in the grid points. In [16] that system of ODEs is the starting point for the application of model reduction techniques to find a suitable lowdimensional system of ODEs from which the aggregated characteristics can be obtained with a reasonable degree of accuracy.
Eq. (3.1) represents a system of n linear non-autonomous ODEs. Since some of the coefficients in the matrices A, B resulting from the discretization of convection terms in the heat equation (2.1) depend on the velocity v 0 (t), it follows that A, B are time-dependent. Thus, (3.1) does not constitute a linear time-invariant (LTI) system. The latter is a crucial assumption for most of model reduction methods such as the Lyapunov balanced truncation technique that is considered in our paper [16]. We circumvent this problem by replacing the model for the geothermal storage by a so-called analogous model which is LTI.
The key idea for the construction of such an analogue is based on the observation that under the assumption of this paper our "original model" is already piecewise LTI. This is due to our assumption that the fluid velocity is constant v 0 during (dis)charging when the pump is on, and zero during waiting when the pump is off. This leads to the following approximation of the original by an analogous model which is performed in two steps.
Approximation Step 1. For the analogous model we assume that contrary to the original model the fluid is also moving with constant velocity v 0 during pump-off periods. During these waiting periods in the original model the fluid is at rest and only subject to the diffusive propagation of heat. In order to mimic that behavior of the resting fluid by a moving fluid we assume that the temperature Q I at the PHX 's inlet is equal to the average temperature of the fluid in the PHX Q F . From a physical point of view we will preserve the average temperature of the fluid but a potential temperature gradient along the PHX is not preserved and replaced by an almost flat temperature distribution. It can be expected that the error induced by this "mixing" of the fluid temperature in the PHX is small after sufficiently long (dis)charging periods leading to saturation with an almost constant temperature along the PHX.
In the mathematical description by an initial boundary value problem for the heat equation (2.1), the above approximation leads to a modified boundary condition at the inlet. During waiting the homogeneous Neumann boundary condition in (2.3) is replaced by a non-local coupling condition such that the inlet boundary condition reads as The above condition is termed 'non-local' since the inlet temperature is not only specified by a condition to the local temperature distribution at the inlet boundary ∂ D I but it depends on the whole spatial temperature distribution in the fluid domain D F . Semi-discretization of the above boundary condition using approximation (4.8) of the average fluid temperature Q F = C F Y formally leads to a modification of the input term g(t) of the system of ODEs (3.1) given in (3.3). That input term now reads as Further, the non-zero entries B l1 of the input matrix B given in (3.2) are modified. They are now no longer time-dependent but given by the constant B l1 = a F h 2 x + v 0 h x which was already used during pump-on periods.
Approximation Step 2. From (6.1) it can be seen that the input term g during pumping depends on the state vector Y via C F Y and can no longer considered as exogenous. Formally, the term C F Y has to be included in AY which would lead to an additional contribution to the system matrix A given by B •1 C F where B •1 denotes the first column of B. Thus, the system matrix again would be time-dependent and the system not LTI. In order to obtain an LTI system we therefore perform a second approximation step and treat Q F as an exogenously given quantity (such as Q I C , Q I D , Q G ). This leads to a tractable approach for model reduction by the Lyapunov balanced truncation technique applied in [16]. The latter generates low-dimensional systems depending only on the system matrices A, B but not on the input term g. Further, from an algorithmic or implementation point of view this is not a problem since given the solution Y of (3.1) at time t, the average fluid temperature Q F (t) can be computed as a linear combination of the entries of Y (t).
Numerical Results. In Figs. 6.1 and 6.2 we present some numerical results where we compare the spatio-temporal temperature distribution and its aggregated characteristics of the original and the associated analogous model. These results are based on the experimental design in the Subsec. 5.4 for a storage architecture with three symmetric PHXs and waiting periods. Fig. 6.1 compares snapshots of the spatial temperature distribution in the storage for the original and analogous model. One snapshot is taken during charging and the other at the end of the last waiting period after preceding discharging periods. At first glance there are no visible differences. A look at the aggregated characteristics in Fig. 6.2 shows negligible approximation errors for the average temperature in the storage Q S and the fluid Q F . However, the approximation of the average outlet temperature Q O suffers slightly from the replacement of a resting fluid by a moving fluid during the waiting period. The resulting "mixing of the temperature profile" inside the PHX adjusts the outlet to the average in the PHX. This can be seen in the right panel where the relative error for the outlet temperature dominates the errors for the two other average temperatures in the storage and the fluid. The experiment indicates that apart from some noticeable approximation errors in the PHX during waiting periods, in particular at the outlet, the other deviations are negligible. Finally, it can be nicely seen that during the (dis)charging periods the errors decrease and vanish almost completely, i.e., in the long run there is no accumulation of errors.
Remark 6.1 The poor precision of the outlet temperature approximation by the analogous model during waiting periods is of no relevance for the management and operation of the geothermal storage within a residential heating system. Here, the outlet temperature is required only during charging and discharging but not during the waiting periods. The interesting quantity for which a good approximation precision is required is the average temperature in the storage and this is provided by the analogous model.

Conclusion
We have investigated the numerical simulation of the short-term behavior of the spatial temperature distribution in a geothermal energy storage. The underlying initial boundary  value problem for the heat equation with a convection term has been discretised using finite difference schemes. In a large number of numerical experiments we have shown 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 heat exchanger PHXs and on the length of charging, discharging and waiting periods.
Based on the findings of this paper, our companion paper [16] studies 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. Q(x, y j ) + 1 2 Q(x, y j ) dx.
Approximating the inner integrals again by trapezoidal rule we get Q(x i , y j ) + 1 2 Q(x i , y j ) , j = j, . . . , j.
Substituting into the above expression for J yields J ≈ h x h y 1 4 Q(x i , y j ) + Q(x i , y j ) + Q(x i , y j ) + Q(x i , y j ) Q(x i , y j ) = h x h y 1 4 Q i j (t) + Q i j (t) + Q i j (t) + Q i j (t) Since the area of the rectangle B is given by (i − i)( j − j)h x h y the average temperature Q B can be approximated by with the coefficients µ i j given in (4.5).

B.2 Numerical Approximation of Q O and Q B
Now we consider the average temperatures Q O and Q B where the temperature Q(t, x, y) is averaged over one-dimensional curves on the boundary ∂ D. Assume that C ⊂ ∂ D is a generic curve on one of the four outer boundaries. For the ease of exposition we restrict C to be a line between the grid points (x i , y 0 ) and (x i , y 0 ) on the bottom boundary, where 0 ≤ i, i + 2 ≤ i ≤ N x . We denote by Q C = Q C (t) = 1 |C| C Q(t, x, y) ds the average temperature in C. Applying trapezoidal rule to the line integral we obtain (suppressing the time variable t) C Q(x, y) ds = Since the length of the curve C is given by (i−i)h x the average temperature Q C can be approximated by where N C = {(i, j) : i = i, . . . , i, j = 0} and the coefficients µ i j of the above quadrature formula are given by 1, for i < i < i, j = 0, (inner points) 1 2 , for i = i, i, (end points).
Using the same notation and approach as above we can rewrite approximation (B.1) as where the matrices D and D are defined as in (4.7) with N 0 B and N 0 B replaced by N 0 C and N C 0 , respectively. Note that in our finite difference scheme the grid values of boundary points are not contained in Y . Thus, we have N 0 C = ∅ and D = 0 1×n . Finally, substituting Y = CY into (B.2) yields a representation of the average temperature Q C as a linear combination of entries of the vector Y which reads as