Grid-Coupled Geothermal and Decentralised Heat Supply Systems in a Holistic Open-Source Simulation Model for 5GDHC Networks

Abstract

In order to reach climate protection goals at national or international levels, new forms of combined heating and cooling networks with ultra-low network temperatures (5GDHC) are viable alternatives to conventional heating networks. This paper presents a simulation library for 5GDHC networks as sustainable shared energy systems, developed in the object-oriented simulation framework OpenModelica. It comprises sub-models for residential buildings acting as prosumers in the network, with additional roof-mounted thermal systems, dynamic thermo-hydraulic representations of distribution pipes and storage, time-series-based sources for heating and cooling, and weather conditions adjustable to user-specified locations. A detailed insight into an in-house development of a sub-model for horizontal ground heat collectors is given. This sub-model is directly coupled with thermo-hydraulic network simulations. The simulation results of energy balances and energetic efficiencies for an example district are described. Findings from this study show that decentralised roof-mounted solar thermal systems coupled to the network can contribute 21% to the total source heat provided in the network while annual thermal gains from the distribution pipes add up to more than 18% within the described settings. The presented simulation library can support conceptual and advanced planning phases for renewable heating and cooling supply structures based on environmental sources.

1. Introduction

The climate protection goals of the German Federal Government increasingly put the focus on technological advancements in the building heating sector. The necessity for reductions in both useful and primary energy is evident from statistics in many years. In 2021, the share of heating applications in the final energy demand of residential buildings was as much as 90% in Germany [1]. Air-to-water heat pumps (AWHPs) constitute over 92% of the total amount of heat pumps for heating installations in Germany in 2023 [2]. One negative implication of a vast roll-out of AWHPs, promoted by monetary funding for investing building owners, is the higher electrical power demand of AWHPs compared to heat pumps exploiting more steady and warmer environmental sources in peak load times during winter. On the other hand, ground-coupled heat pump systems require higher financial and organisational expenses for accessing the heat source. A district-wide shared utilisation of regenerative heat sources and the shift of simultaneous heating and cooling demands make the concept of 5th-generation heating and cooling (5GDHC) networks attractive. They have the potential to increase the overall energetic efficiency in comparison to individual building-level solutions and support a development towards more flexible and sustainable energy supply systems. In the face of rising temperatures due to climate change, an increase in cooling demands in residential buildings further promotes the use of this excess heat source in shared energy systems.

1.1. 5GDHC Network Systematisation

In general, low-temperature networks with system temperatures below 60 °C are the most recent and favoured type of heating network in terms of primary supply from renewable heat sources. Heating and cooling networks with decentralised heat pump installations are seen as a central component for the decarbonisation of the heating sector for room heating and domestic hot water (DHW) production, with a special focus on new buildings [3]. Regarding the systematic classification of low-temperature networks there are multiple approaches found in the literature [3,4,5]. A widely used arrangement of different development stages in district heating and cooling networks is provided by Lund et al. [6], who categorise heating networks with system temperatures below 70 °C as 4th-generation networks. According to this description, most contemporary heating networks fall into the class of 3rd-generation heating networks [7]. Whereas 5GDHC networks are operated close to the temperature of the environment or the soil regime in which the distribution pipes are buried. Apart from the fact of having more suitable conditions for the feed-in of thermal energy from renewable sources at very low temperature levels than 4th-generation networks, the advantages of 5GDHC networks also encompass the shift of simultaneous heating and cooling demands between prosumers connected to the grid, as is outlined in existing research work [3,5,8]. This ability for bidirectional energy flow paired with undirected medium flow in the network is characteristic of 5GDHC networks [3]. These 5GDHC networks can be designed as active or passive networks or hybrid forms of both. Active networks feature central pumping units to ensure proper circulation of the medium whereas passive networks rely on circulation pumps contained in decentralised substations or connected buildings, respectively [9]. An example of the realisation of a hybrid network can be found in [10].

1.2. Current Scope of Modelling and Simulating 5GDHC Networks

Heissler [11] uses a co-simulation for a common simulation of a thermo-hydraulic representation of the grid and building models in a 5GDHC network. In the study, demand scenarios and network topologies as well as heat demand per trail length are varied. The building models are developed in TRNSYS [12] and incorporate buffer storage on the primary and the secondary sides of the heat pump, hot water storage to supply internal heating loops, and an additional solar loop for roof-mounted units. Bidirectional energy flows and undirected medium flows are handled by two pump models with reverse flow direction. Sub-models for the pumps and distribution pipes are modelled in Dymola [13].

Abugabbara [14] describes a simulation model for 5GDHC networks designed in Dymola [13] and based on the Modelica.Fluid library [15]. The work also provides a case study for an example network in Sweden with both heating and cooling demands. The focus of the study is put on the effect of shared simultaneous heating and cooling demands and the benefits compared to conventional dedicated heating networks or cooling networks. Sub-models for the buildings in this simulation model are equipped with substations able to operate in heating mode, in active cooling mode by use of reversible heat pumps, or in passive cooling mode. Bidirectional energy flows and undirected medium flows are handled so that substations may draw from the cold line when cooling demands are dominant. A balancing unit is included in the presented simulation model to maintain the temperatures of the warm line and the cold line in the network within set limits by feed-in or extraction of energy.

In [16], MATLAB and Simulink are used as software environments for a simulation tool for 5GDHC networks. The aim of the authors is to provide a flexible software tool for an early stage of network planning and an analysis of energy flows in the network. Building models are represented by simple water-to-water heat pumps as substations with an efficiency mapping according to the specified sink-side feed temperatures and the network temperature. As described in [14], a balancing unit serves as a generic source for heating and cooling in the network.

In [10], the simulation and subsequent construction of a multi-source 5GDHC network is described. Whereas horizontal ground heat collectors (GHCs) and ground-coupled ice storage are modelled in the software environment DELPHIN (version DELPHIN5), the network itself and the building models are designed in an in-house development environment called Sim-Vicus. All sub-models are connected in a complete-system simulation via functional mockup interfaces. The direct coupling of the soil model in the sub-models of the GHC and the ice storage with the network calculation are seen as a promising advancement in the field of simulation of 5GDHC networks. This aspect also stands out as a main innovation of the present work.

Blacha et al. [17] present a simulation model and a case study for low-temperature networks with simultaneous heating and cooling demands similar to the application in [14]. The building models include heat pumps and decentralised distribution pumps, which provide the mass flow to the substations. Furthermore, active cooling modes can be simulated with the help of chiller units. Demand profiles for heating and cooling can be transferred to the building models as time series. A central balancing unit connected to the network keeps the network temperatures in the warm and cold lines within defined ranges. Further references for simulation models capable of handling the complex interdependencies of components in 5GDHC networks can be found in the literature review in [18].

With the focus on the Modelica language for integrated simulations at building or district level, there are several existing libraries worth mentioning. A comprehensive review of these and other libraries can be found in [19]. Among the libraries developed within the Annex60 project [20] which has further been known under its new name IBPSA since 2017 are the following:

AixLib is a development by the RWTH, Aachen. The library incorporates a vast collection of sub-models at low and high levels of detail for building performance analysis simulations. The building models at a low level of detail are composed of serial and parallel connections of thermal capacities and resistances. The aforementioned plug-flow pipe model is included in the library as well as models for all kinds of thermal and electrical components in energy systems on the building and district scales.

The Buildings library, developed by Lawrence Berkeley National Laboratory contains detailed sub-models for the thermal characterisation of buildings in interaction with the environment. Apart from large collections of material characteristics in the Buildings.Fluid database, the library also contains extensive models for heating, ventilation andair conditioning systems.

IDEAS, as a free library developed by KU Leuven, puts the emphasis on simultaneous and coupled simulations of electrical and thermal systems on a smaller building or component scale. However, the library also includes a model for fields of borehole heat exchangers, which are commonly used as heat sources and storage on the scale of a district heating and cooling network. This model uses interfaces for media transport from the FluidHeatFlow library, part of the Modelica standard library.

1.3. Paper Aims

This paper aims to present the developed open-source library components in OpenModelica which can be used for the design of holistic simulation models for 5GDHC networks. The library so far comprises sub-models for regenerative heat sources at both district grid scale and individual building scale, sub-models for the thermo-hydraulic representation of distribution pipes, and sub-models for decentralised substations and prosumers. Thus, several effects crucial to the dimensioning of heat sources and overall performance analysis can be investigated. An in-house development of a sub-model for a GHC that can be connected to the network enables analysis of the interaction between this ground-coupled geothermal heat source and additional central or decentralised heat sources in 5GDHC networks. Aspects like the degree of regeneration of this geothermal heat source by other heat sources distributed in the network and the impact of roof-mounted solar thermal systems on the overall efficiency of the network operation can be addressed. Sufficient dimensioning and regeneration methods are particularly important to maintain sustainable agricultural use of land above the heat source. A low-threshold for access to the simulation of 5GDHC networks is ensured as the input parameters to building models, piping models, and heat source models are kept concise and relate to commonly used methodologies. Special focus is put on the GHC model as only a few numerical calculation approaches are known from the literature and no general and standardised approach is available. This paper not only provides a detailed insight into the modelling of the GHC, it also gives a description of the straight coupling of the model with a connected 5GDHC network topology and attached prosumers. Another focus is put on two case studies performed: The first case study highlights the behaviour of the developed GHC model in a district simulation under varied boundary conditions from climatic constraints, soil conditions, and reduced GHC area. The impact of free cooling from residential buildings on the performance of a GHC is analysed as well. In a second case study, a more comprehensive district simulation comprising the GHC as a primary heat source is presented. Here, the effect of additional, network-coupled heat sources on the efficiency of decentralised heat pumps and on the overall system performance is examined.

The originality of this paper is pointed out in the following: In the scope of the object-oriented Modelica language, no other work is known to the authors that provides a direct coupling of detailed GHC models to thermo-hydraulic network representations. The comprehensive survey of modelling approaches which are combined in the developed GHC model serves as a guidance and reference work in this field. Finally, there are only a few studies that combine a detailed geothermal heat source model, heating and cooling distribution networks, and prosumer models equipped with decentralised heat sources in a single open-source simulation framework. The results from two case studies further highlight the comparability to other established research findings and outline the possibilities and the soundness of the simulation library.

The presented work can contribute to the propagation of the still-young technology of 5GDHC networks. The understanding of interconnections between different heat sources in a shared energy network as well as between heating and cooling demands on a district scale are important to the planning process of centralised heat supply systems. This also applies to the economical implications of the chosen technologies like ground-coupled geothermal heat sources and storage. A promotion of this form of sustainable heating and cooling supply system as an alternative to conventional heating systems helps to transfer climate protection policy into long-term practical realisation. What is more, a broader use of open-source simulation frameworks is supported by the outcomes of this research work. The provided collection of used modelling approaches for the representation of the main physical components in 5GDHC networks can be a great assistance for other researchers when developing or adopting similar models. The results from the case studies are expected to underline the rising need for highly individualised and simulation-assisted planning procedures. This is of special interest for heating networks of recent generations due to the importance of the temporal resolution of heating and cooling demand and supply. In addition, for operators of 5GDHC networks, more controllable energy consumption patterns, e.g., by sharing among prosumers, can assist an adequate dimensioning of connections to superior electrical grids and heat sources within the network. These aspects are seen by the authors as highly beneficial to the current technical and societal challenges arising with the need for legal guidelines for municipal planning of heat supply like in Germany.

2. Materials and Methods

The following section presents an overview of the developed holistic simulation environment for 5GDHC networks in OpenModelica v.1.20.0-64bit. It is further used to present several sub-models of the library. Focus is put onto crucial components like the piping models and in-house developments of heat sources, especially on the horizontal ground heat collector.

OpenModelica is used as the simulation framework for this work for several reasons. As stated in [18], this framework is widely used for integrated building and network simulations. For this reason, a vast collection of validated sub-models and basic components, e.g., for modelling heat transfer problems, exists. Furthermore, it comes as an open-source distribution, which makes it possible to share developments with a growing community. In current times of increasing need for simulation-assisted heat planning procedures at district- and community-levels, this is seen as an advantageous argument for this framework.

2.1. Basic Modelling Structures and Libraries

The structure of the presented simulation library is aligned with the division into subsystems found in conventional heating networks or energy supply networks in general. It features sub-models for supply systems of heat and cold, thermo-hydraulic sub-models for distribution pipes, and comprehensive sub-models for decentralised prosumers attached to the network. The prosumer models themselves can be equipped with different decentralised components of heat suppliers such as roof-mounted solar thermal systems.

The interfaces of the components are modelled with the use of the Modelica standard library FluidHeatFlow [21]. This library uses constant material properties, assumes incompressibility of the fluids, and enables the user to explicitly and easily define new media. For the application of 5GDHC simulations at a network scale, the simplification of constant material properties is seen as a reasonable trade-off due to comparatively small ranges of network temperatures.

2.2. Sub-Models for Heat and Energy Sources

The following sections deal with the sub-models for heat sources in 5GDHC networks included in the simulation library. The list of described models is not exhaustive though, and only the most relevant models in the context of the case studies described later on are presented.

2.2.1. Horizontal Ground Heat Collector (GHC)

As one ground-coupled heat source and seasonal buffer storage, a model of a horizontal ground heat collector is incorporated in the simulation library. Due to a lack of standardised or validated models available for implementation in Modelica, an in-house development is used.

GHCs are usually installed at depths of around 0.8 m to 2 m and make active use of the latent heat stored in the water content of the soil. The GHC model is based on the thermal representation of the soil regime surrounding the collector pipes and the heat exchange between the circulating medium in the pipes and the soil. Heat transport due to convection and irradiance at the ground surface as well as heat conduction within the soil regime are modelled. Moisture transport itself is not modelled explicitly in the model at hand, but it is used as a time-series input to the sub-models of the soil cells, impacting material properties and both the sensible and latent heat transfer.

Discretisation and Computational Domain

For numerical calculations, the soil regime is commonly discretised in the horizontal x-direction and vertical z-direction with varying resolution, from fine granularity around the collector pipes and near the ground surface to a more coarse segmentation at greater depths; compare [22,23,24]. The aims and outputs of this model are the temperature distribution within the soil regime and the energy balance of the operational medium passing through the collector pipes. In order to reduce the effort for numerical calculations, the computational domain is limited by symmetry borders in the x-direction. In this way, both geometrical and thermal symmetries can be exploited. The computational domain can be reduced in size, and hence in the number of discretised soil cells. The boundaries of the computational domain in the x-direction represent adiabatic planes. The GHC is modelled as if the distance between the parallel collector pipes is constant throughout the entire GHC area. With this approach, the representation of collector pipes within the soil regime does not take into account the different temperature levels of the circulating medium in neighbouring pipes stemming from counterflow arrangements. Hence, complex pipe routing such as bifilar windings, as shown in Figure 1, are simplified and reduced to a layout of quasi-indefinite parallel collector pipes (symmetrical in the x-direction and denoted as the second level of abstraction in Figure 1).

Otherwise, an explicit study of the boundary effects of single modules with bifilar windings that are connected in series or in parallel, or 3-dimensional representations of the soil regime to account for asymmetries, would become necessary. In Figure 1, distribution pipes can be seen on either side of the collector pipe windings. In order to adequately consider the additional pipe length of these distribution pipes for pressure losses and contributions to heat transfer with the surrounding soil, a fraction of the distribution pipes is added to the total pipe length. This grows proportionally with the specified length of the collector pipes. For reasons of simplification, their length-specific heat flow is assumed to be identical to the value calculated explicitly for the collector pipes (see following paragraphs).

A GHC can be modelled with 2-dimensional or 1-dimensional discretisation of the soil regime dependent on the user input. In case a 1-dimensional discretisation is desired, the computational domain of the soil regime extends to 1 m in the x-direction. Consequently, energy balances are calculated as average values within this extent. The collector pipe model adjoins one soil zone above and one soil zone below.

In the case of a 2-dimensional discretisation, the computational domain in the x-direction extends from the centre of a collector pipe to half the distance to the neighbouring pipe to the right. Hence, only one half of the collector pipe is modelled and takes part in the energy exchange with the soil regime. The collector pipe borders soil zones above, to the right, and below.

In both approaches, the maximum depth in the z-direction can be set by user inputs. It is advisable to choose a minimal depth limitation of 10 times the installation depth of the collector pipes, e.g., 20m–25m. At this depth limitation, a fixed temperature is set as the boundary condition to the model. This is commonly chosen as a constant value of 9 °C to –11 °C; see [24]. A juxtaposition of both discretisation approaches is shown on the right in Figure 2.

The resolution of the soil discretisation can be controlled by user inputs to the model. For this purpose, the soil regime is split into 15 separate zones in the z-direction, which can each be equipped with an independent discretisation resolution and an individual total extent. This results in a varying number of soil cells in each zone based on the user inputs. In the case of a 2-dimensional discretisation, the soil regime can further be split into three zones in the x-direction, also with adjustable discretisation. However, the model features a dynamic and automatic discretisation control depending on the specified installation depth of the collector pipes: In the vicinity of the collector pipe and near the ground surface, the (relatively) highest heat fluxes and temperature gradients are expected. Therefore, the discretisation resolution is refined in these regions. Ramming [22] and Hirsch [23] point out that the discretisation resolution close to the collector pipe is essential to an accurate representation of the ice-phase formation within the soil texture. Ramming [22] uses a resolution for soil cells of half a pipe diameter $d_{pipe}$ in the near field of the collector pipe. Hirsch [23] provides a study on discretisation resolution and subsequently uses a resolution of $\Delta x=\Delta z={\displaystyle \frac{d_{pipe}}{6}}$.

As a compromise between accuracy and computational effort, a default value of $\Delta x=\Delta z={\displaystyle \frac{d_{pipe}\pi }{4n}}$ is used in the GHC model, where $n^2$ constitutes the number of quadratic soil cells with the same area of the cross-section of the pipe. Near the ground surface, the discretisation resolution is chosen as 5 cm by default. With increasing distance below the installation depth of the collector, the resolution increases up to $\Delta z=1m$ near the bottom of the computational domain.

In Figure 3, the thermal connection between neighbouring soil cells is depicted in a region around the collector pipe in the 2-dimensional model. Thermal connections between cells within the same zone and between soil zones are built using heat port models from the standard library Modelica.Thermal.HeatTransfer. The compiled model is translated into a system of heat conduction equations between adjacent soil cells. Explicit formulations are applied for heat fluxes at the collector pipe outer wall and the boundary conditions expressed at the ground surface and the maximum depth.

Heat Transport in the Soil Regime

Each soil cell consists of a serial connection of a thermal capacitance and thermal resistances. The general, transient 2-dimensional heat conduction equation without source term applied to the soil cells reads ([25], p. 125)

$${\displaystyle \frac{\partial }{\partial x}}\left({\lambda }_{eff}{\displaystyle \frac{\partial T}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\lambda }_{eff}{\displaystyle \frac{\partial T}{\partial z}}\right)={\displaystyle \frac{C_{eff,cell}}{V_{cell}}}{\displaystyle \frac{\partial T}{\partial \tau }}$$

(1)

with the effective heat conductivity ${\lambda }_{eff}$, the heat capacity $C_{eff,cell}$, and the volume $V_{cell}$ of the soil cell. For two neighbouring soil cells at depth $z_i$ and at x-locations $x_j$ and $x_{j+1}$, the heat flux (per unit depth in the y-direction) is expressed as

$${\displaystyle \frac{{\dot{Q}}_{x,j\to j+1}}{\Delta y}}={\displaystyle \frac{T_j-T_{j+1}}{\frac{\Delta x_j}{2{\lambda }_{eff,j}}+\frac{\Delta x_{j+1}}{2{\lambda }_{eff,j+1}}}}\Delta z_i.$$

(2)

Generally, a set of parameters is used to properly describe the soil characteristics in the presented model. These can be saved as records of fixed properties for a specific type of soil in the simulation library. The necessary fixed parameters are as follows:

• The dry bulk density of the soil, ${\rho }_d$ (as described for various soil types in [22]).
• The gravimetric composition of the mineral soil matrix with the fractions of sand, silt, and clay ($x_{sand}$, $x_{silt}$, and $x_{clay}$, respectively), which are provided by [26] or by the widely used USDA soil taxonomy triangle.
• The volumetric water content at saturation, ${\mathsf{\Theta }}_{sat}$, and the residual (minimal) volumetric water content, ${\mathsf{\Theta }}_{res}$ (provided for twelve reference soil types in [27]).
• The empirical van Genuchten parameters ${\alpha }_{VG}$, $n_{VG}$, and $m_{VG}$ (provided in [27]).
• Optionally, the volumetric content of organic material, $Y_{org}$.

Apart from these parameters, a time-dependent volumetric water content, $\mathsf{\Theta }$, as a function of depth can be specified for the soil zones.

Neglecting the influence of air captured within the soil matrix, the effective density of the soil, ${\rho }_{eff}$, is dependent on the dry bulk density of the soil type, ${\rho }_d$, the volumetric contents of unfrozen and frozen water, ${\mathsf{\Theta }}_u$ and ${\mathsf{\Theta }}_f$, and the volumetric content of organic material, $Y_{org}$, as well as its density, ${\rho }_{org}$ (see Equation (3)). ${\rho }_w$ and ${\rho }_I$ denote the densities of water and ice, respectively.

$${\rho }_{eff}={\rho }_d+{\mathsf{\Theta }}_u{\rho }_w+{\mathsf{\Theta }}_f{\rho }_I+Y_{org,vol}{\rho }_{org}$$

(3)

Volumetric changes as a consequence of the formation of ice are neglected. This simplification is also applied in the models described in [22,24]. The effects and ecological implications of frost heave are investigated separately in other studies; see [10,28].

The determination of the effective heat conductivity is based on the method introduced by Coté and Konrad [29]. From the gravimetric proportions of sand, silt, and clay in the mineral content of the soil, the thermal conductivity of the soil matrix, ${\lambda }_d$, can be derived by Equation (4).

$${\lambda }_d=(4,5·x_{sand}+3,5·x_{silt}+2,5·x_{clay})·1{\displaystyle \frac{\mathrm{W}}{\mathrm{m}\mathrm{K}}}$$

(4)

Subsequently, the thermal conductivity of the unfrozen soil, ${\lambda }_u\left(\mathsf{\Theta }\right)$, and of the completely frozen soil, ${\lambda }_f\left(\mathsf{\Theta }\right)$, can be determined—for further details on the formulae, refer to the cited work [29]. Finally, the effective thermal conductivity of the soil cell, ${\lambda }_{eff}$, is calculated by a linear interpolation between thermal conductivities for frozen and unfrozen soil following Equation (5).

$${\lambda }_{eff}\left(\mathsf{\Theta }\right)={\displaystyle \frac{{\mathsf{\Theta }}_u}{\mathsf{\Theta }}}{\lambda }_u\left(\mathsf{\Theta }\right)+(1-{\displaystyle \frac{{\mathsf{\Theta }}_u}{\mathsf{\Theta }}}){\lambda }_f\left(\mathsf{\Theta }\right)$$

(5)

Note that $\mathsf{\Theta }={\mathsf{\Theta }}_u+{\mathsf{\Theta }}_f$ and the water vapour phase are neglected in this approach. In the literature, alternative calculation approaches can be found for the determination of the thermal conductivity of soils. Approaches based on an equation by de Vries [30] directly associate the thermal conductivities of mineral constituents in the soil matrix and continuous media in the soil regime. An effective thermal conductivity is yielded by an average of the single thermal conductivities, weighted with estimated coefficients for liquid water, solid water, mineral content, etc. (see [31], p. 43, [32]). In contrast, Mottaghy et al. [33] calculate the effective thermal conductivity in unsaturated and partially frozen soil as the geometric mean of the thermal conductivity of each constituent and its volumetric fraction as exponent. However, as no procedure for an interpolation of effective thermal conductivities between unfrozen and unfrozen soil is found in the literature for the applied model by Coté and Konrad [29], the linear interpolation in Equation (5) is used. Significant deviations from this interpolation method within the small temperature range of changing ice content during freezing are not expected.

The volumetric unfrozen water content during ice-phase formation below $273.15$ K is calculated according to the van Genuchten–Wen model, which is described in [34] and compared against other methods in [35]. Soil-specific empirical parameters are used to describe the hydraulic traits. It directly links the unfrozen water content to an initial water content at the start of the freezing process and to the residual water content. The formulae can be seen in Equation (6). $h_{lat}$ denotes the latent heat of fusion of water, taken as 333.5 $\frac{\mathrm{kJ}}{\mathrm{kg}}$.

$${\mathsf{\Theta }}_u=\left\{\begin{array}{cc}{\mathsf{\Theta }}_{res}+(\mathsf{\Theta }-{\mathsf{\Theta }}_{res})·{(1+{(-{\alpha }_{VG}{h}_{lat}{\rho }_w{\displaystyle \frac{T-273.15K}{273.15K}})}^{n_{VG}})}^{-m_{VG}},& T\le 273.15K\\ \mathsf{\Theta },& T>273.15K\end{array}\right.$$

(6)

Hence, an apparent elevation of the volumetric specific heat capacity during the freezing process can be formulated to consider the latently stored heat on account of the water content in the soil [36], as expressed in Equation (7).

$$c_{vol,lat}={\displaystyle \frac{\partial {\mathsf{\Theta }}_u}{\partial T}}·{\rho }_w{h}_{lat}$$

(7)

Adding to the effective volumetric specific heat capacity $c_{vol,eff}$ of each soil cell are the contributions of the unfrozen and the frozen water, $c_{vol,u}$ and $c_{vol,I}$, respectively, and the contributions of organic and mineral materials, $c_{vol,org}$ and $c_{vol,solid}$; see Equation (8). The volumetric specific heat capacity of enclosed air is considered with the term $c_{vol,air}$. All these contributing volumetric specific heat capacities are determined as a product of their volumetric content, their density, and their gravimetric specific heat capacity [36].

$$c_{vol,eff}=c_{vol,lat}+c_{vol,w}+c_{vol,I}+c_{vol,org}+c_{vol,solid}+c_{vol,air}$$

(8)

The final equation for the heat capacity passed to the internal model of a thermal capacitance in each soil cell reads

$$C_{eff,cell}=c_{vol,eff}·V_{cell}.$$

(9)

Heat Transport into the Collector Pipe

A single collector pipe of a GHC usually consists of an uninsulated PE- or PE-Xa-pipe with common nominal diameters ranging from DA20 to DA32. Customary pipe data can be selected from pre-defined records in the simulation library. The length is commonly chosen to 70–100 m to limit pressure losses and multiple pipes are connected in parallel to a bigger distribution pipe to achieve uniform flow conditions. The heat transfer between the collector pipe and the surrounding soil cells is calculated in detail and depends on the actual flow conditions in the pipe. However, it is based on a quasi-stationary formulation of the heat transport, neglecting the thermal capacity of the pipe material itself. Heat conduction in the flow direction is neglected in the model. This assumption can also be found in other studies presenting simulation models for GHCs, such as [22].

A heat exchanger characteristic $\mathsf{\Phi }$ is used to determine the heat flux into the fluid circulating in the collector pipe and its outlet temperature. Following the commonly known denomination of transfer units N and the definition of the capacity flow ${\dot{C}}_{in}$ at the inlet of the heat exchanger, the basic relation to the heat transfer coefficient k and the heat exchanger surface A is given in ([25], p. 51) as

$$N={\displaystyle \frac{kA}{{\dot{C}}_{in}}}.$$

(10)

The heat transfer coefficient k is calculated as a combination of the convective heat transfer coefficient at the inner wall of the pipe and heat conduction in the pipe material; see Equation (11). $d_{out}$ and $d_i$ denote the outer and inner pipe diameters, respectively, whereas ${\lambda }_{pipe}$ describes the thermal conductivity of the pipe material ([25], p. 35).

$$k={\left({\displaystyle \frac{d_{out}}{k_{conv,i}{d}_i}}+{\displaystyle \frac{d_{out}ln\left(\frac{d_{out}}{d_i}\right)}{2{\lambda }_{pipe}}}\right)}^{-1}$$

(11)

The convective heat transfer at the inner pipe wall, $k_{conv,i}$, is calculated with Nusselt correlations dependent on the actual flow conditions in laminar or turbulent ranges; see, for example, [37] (p. 266). Based on the assumption of an external temperature profile at the outer pipe wall which runs in parallel to the temperature of the fluid, the following characteristic can be derived [38]:

$$\mathsf{\Phi }=min\left({\displaystyle \frac{2N}{N+2}};1\right)={\displaystyle \frac{T_{out}-T_{in}}{T_{wall}-T_{in}}},$$

(12)

where $T_{out}$ is calculated as the outlet temperature of the fluid. Its rise or decrease is limited to the temperature at the contact surface with neighbouring soil cells.

Heat Transfer at the Ground Surface

At the ground surface, the heat transfer is primarily described by convective heat exchange with the air and the irradiance, assuming a non-sealed surface, and by the influence of precipitation. The methodology implemented in the model is in analogy to the process described by Ramming [22]. Air temperatures ($T_{amb}$), and values for direct ($e_{dir,hor}$), diffuse ($e_{diff,hor}$), and long-wave ($e_{IR}$) irradiation and for precipitation are taken from the test–reference–year (TRY) datasets provided by the German meteorological service [39]. They can be adapted by the user to the specific location for which a simulation is performed.

Direct, indirect, and long-wave (infrared) irradiance sum up to a heat flux ${\dot{q}}_{irr}$ at the ground surface according to Equation (13). An optional factor representing a temporarily variable shielding of direct irradiance can be specified by the user. This enables taking temporary shielding below structures, trees, etc., into account. Its allowed values range from 0 (meaning that all direct irradiance is shielded) to 1 (direct irradiance is not shielded). By default, and for the analyses at hand, it is set to 1 so that incident irradiance is not attenuated. The absorption coefficient $a_{sw}$ for short-wave irradiation is set to 0.75 and the emission coefficient for long-wave radiation, ${ϵ}_{lw}$, is set to 0.98 by default ([22], p. 46):

$${\dot{q}}_{irr}=a_{sw}·(e_{dir,hor}·shielding+e_{diff,hor})+{ϵ}_{lw}·e_{IR}$$

(13)

Terrestrial radiant emittance of the ground surface is dependent on the surface temperature $T_{gs}$ [22]:

$${\dot{q}}_{terr}=-{ϵ}_{lw}·\sigma ·T_{gs}^4,$$

(14)

where $\sigma =5.67\mathrm{e}-08{\displaystyle \frac{\mathrm{W}}{{\mathrm{m}}^2{\mathrm{K}}^4}}$ denotes the Boltzmann constant.

As for the convective heat transfer, an approach from Duffie et al. [40] is adopted, where the heat transfer coefficient at the ground surface is related to the wind velocity, $v_{wind}$, so that the equation for the convection-driven heat flux reads

$${\dot{q}}_{conv}=(T_{amb}-T_{gs})·\left(2.3+3·v_{wind}·1{\displaystyle \frac{\mathrm{s}}{\mathrm{m}}}\right)·1{\displaystyle \frac{\mathrm{W}}{{\mathrm{m}}^{\mathrm{2}}\mathrm{K}}}.$$

(15)

As mentioned before, transport phenomena for moisture are not modelled explicitly in the soil regime. However, part of the influence of precipitation is considered at the ground surface by assessing the heat flux fed into or extracted from the top soil zones; see Equation (16). Precipitation is assumed to reach the ground surface with the wet bulb temperature $T_{wb}$ of the actual ambient air conditions ([22], p. 47).

$${\dot{q}}_{prec}={\dot{V}}_{prec}·{\rho }_w·c_w·(T_{gs}-T_{wb})$$

(16)

${\dot{V}}_{prec}$ denotes the volumetric flow of precipitation per square meter, taken from the TRY dataset, $c_w$ denotes the specific heat capacity of water. In total, the heat flux at the ground surface sums up according to Equation (17) (see [22], p. 48), which is used as the boundary condition in the GHC model.

$${\dot{q}}_{tot,gs}={\dot{q}}_{irr}+{\dot{q}}_{terr}-{\dot{q}}_{prec}+{\dot{q}}_{conv}$$

(17)

2.2.2. Hybrid Photovoltaic–Thermal (PVT) Collectors

In the presented simulation library, a comprehensive thermal and electrical model for uncovered collectors is implemented, which allows them to be used both as centralised and decentralised roof-mounted heat sources to feed into the network.

Description of Thermal Behaviour

According to the information in ISO 9806:2017 [41], the thermal characteristics of uncovered solar thermal collectors are commonly described by eight coefficients. Broetje [42] provides an example application of this standard, which is adopted for the modelling approach at hand. Similar to the model of the collector pipe in the GHC, the heat transfer to the circulating fluid in the PVT collector is described with a heat exchanger characteristic to make the calculation numerically stable for small and vanishing volumetric flows (see Equations (10) and (12)). For the PVT model $\mathsf{\Phi }$ is limited to a maximum value of 0.9, ensuring the outlet temperature of the fluid is always slightly below the actual temperature of the scalar capacitance representing the PVT collector in heating mode. This represents a conservative assumption as an over-estimation of the thermal power and the outlet temperature of the PVT modules should be avoided.

The specific thermal power per area provided by the PVT collector module is expressed in Equation (18), referred as to the aperture area $A_{apert,PVT}$ of the module [42].

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\dot{Q}}_{th,PVT}}{A_{apert,PVT}}}=& e_{dir,incl}·{\eta }_{dir}·IA{M}_{dir}\left({\varphi }_{incl}\right)\hfill \\ & +(e_{diff,incl}+e_{refl,incl})·{\eta }_{dir}·IA{M}_{diff}\hfill \\ & +(T_{hem}^4-T_{amb}^4)·{\displaystyle \frac{ϵ}{a}}·\sigma ·{\eta }_{hem}\hfill \\ & -{\eta }_{hem}·b_{wind}·v_{wind}·(e_{dir,incl}+e_{diff,incl})\hfill \\ & -a_1·({\overline{T}}_{mod}-T_{amb})\hfill \\ & -a_2·{({\overline{T}}_{mod}-T_{amb})}^2\hfill \\ & -a_3·v_{wind}·({\overline{T}}_{mod}-T_{amb})\hfill \\ & -a_5·{\displaystyle \frac{d{\overline{T}}_{mod}}{d\tau }}\hfill \end{array}$$

(18)

The notation is as follows:

• ${\varphi }_{incl}$: Incident angle of direct irradiation on module surface.
• ${\eta }_{dir}$: Conversion factor for direct irradiation at $({\overline{T}}_{mod}-T_{amb})=0$.
• $IA{M}_{dir}\left({\varphi }_{incl}\right)$: Incident angle modifier (IAM) for direct irradiation.
• $IA{M}_{diff}$: IAM for diffuse irradiation.
• ${\eta }_{hem}$: Conversion factor for hemispherical irradiation at $({\overline{T}}_{mod}-T_{amb})=0$
  ${\eta }_{hem}={\eta }_{dir}·IA{M}_{dir}\left({\varphi }_{incl}\right)·{\displaystyle \frac{e_{dir,incl}}{e_{tot,incl}}}+{\eta }_{dir}·IA{M}_{diff}·{\displaystyle \frac{e_{diff,incl}}{e_{tot,incl}}}$
• $T_{hem}$: Black body sky temperature; assumption $T_{hem}=T_{amb}-10\mathrm{K}$.
• $T_{amb}$: Ambient air temperature.
• $ϵ$: Hemispherical emissivity of PV glass cover.
• a: Absorption coefficient of PV glass cover.
• $\sigma$: Boltzmann constant $(5.67\mathrm{e}-08{\displaystyle \frac{\mathrm{W}}{{\mathrm{m}}^{\mathrm{2}}{\mathrm{K}}^{\mathrm{4}}}})$.
• $b_{wind}$: Wind-dependence coefficient of the collector. (design parameter)
• ${\overline{T}}_{mod}$: Average module temperature (temperature of the capacitance representing the collector module).
• $v_{wind}$: Wind velocity.
• $a_1...a_3$: Loss coefficients of thermal collector model according to ISO 9806:2017 [41] (design parameter).
• $a_5$: Overall specific heat capacity of the collector module (design parameter).

One further input to the model is the volumetric flow per module of the operational medium, ${\dot{V}}_{mod}$. Based on the control mode, this input can either be a fixed parameter taken from the manufacturer’s data or a varying value with an upper limit to apply a control strategy to achieve a desired output temperature. The inclination angles $\gamma$ and azimuth angles $\phi$ can be passed to the model by user inputs.

Description of Electrical Behaviour

The electrical model of the PVT collectors follows the approach given by Jonas et al. in [43] and adopted, for example, in [44]. It takes the temperature dependency of the electrical output power of the modules into account, with beneficial influence at lower module temperatures. The specific electrical power per area is expressed as

$${\displaystyle \frac{P_{el,PVT}}{A_{apert,PVT}}}={\eta }_{el,ref}·{\eta }_{el,loss}·e_{tot,mod}$$

(19)

in [43], with the total incident irradiation on the collector, $e_{tot,mod}$, as

$$e_{tot,mod}=e_{dir,incl}·IA{M}_{dir}\left({\varphi }_{incl}\right)+e_{diff,incl}·IA{M}_{diff}+e_{refl,incl}.$$

(20)

The incidence angle modifiers for electrical and thermal description are assumed to be identical if no explicit data are available. The electrical efficiency ${\eta }_{el,ref}$ at the reference temperature $T_{ref}=298.15K$ is altered by the factor ${\eta }_{el,loss}$, incorporating the influences described in the following. The module-specific parameters $a_{el,1}$, $a_{el,2}$, and $a_{el,3}$ summarise the dependency on the irradiation conditions. The parameter $a_{el,T}$ accounts for the influence of the temperature deviation from the reference conditions. Hence, ${\eta }_{el,loss}$ is expressed in [43] as

$$\begin{array}{cc}\hfill {\eta }_{el,loss}=& \left(a_{el,1}·e_{tot,mod}+a_{el,2}·ln(e_{tot,mod}+1)+a_{el,3}·[{\displaystyle \frac{{\left(ln(e_{tot,mod}+\mathrm{e})\right)}^2}{e_{tot,mod}+1}}-1]\right)\hfill \\ & ·\left(1-a_{el,T}·({\overline{T}}_{mod}-T_{ref})\right).\hfill \end{array}$$

(21)

Successive efficiency losses due to inverters are allocated individually. The net electrical power output of the PVT system is yielded after subtraction of the consumption in the circulation pump.

2.2.3. Solar Thermal (ST) Collectors

The modelling of ST collectors comprises their thermal and hydraulic description. In the presented simulation library, they are used as roof-mounted systems which can directly feed thermal energy into hot water storage in prosumer models and subsequently into the 5GDHC network itself.

Description of Thermal Behaviour

In the simulation library, a generic model of collector modules is implemented. The thermal efficiency is described by a quadratic dependency on the temperature difference between the temperature of the collector, summarised in a scalar capacitance like the PVT-collector model, and the ambient temperature.Both covered flat-plate collectors and evacuated tube collectors can be described with this approach. Relevant collector-specific coefficients are available in many manufacturers’ data sheets. Equation (22) gives the formula for the specific thermal output power of the solar thermal collector model per aperture area, taken from [45].

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\dot{Q}}_{th,ST}}{A_{apert,ST}}}& =(e_{dir,incl}·IA{M}_{dir}\left({\varphi }_{incl}\right)+[e_{diff,incl}+e_{refl,incl}]·IA{M}_{diff})·{\eta }_{dir}\hfill \\ & -a_1·({\overline{T}}_{mod}-T_{amb})\hfill \\ & -a_2·{({\overline{T}}_{mod}-T_{amb})}^2\hfill \\ & -a_5·{\displaystyle \frac{d{\overline{T}}_{mod}}{d\tau }}\hfill \end{array}$$

(22)

The denominations are analogous to those in Equation (18). Again, $a_1...a_5$ are coefficients specific to the chosen collector type.

2.3. Sub-Models for the Distribution Network

It is essential to the modelling approach for all components of the distribution network that a varying flow direction and bidirectional energy exchange with the surroundings is properly reproduced. The widely used Modelica standard library FluidHeatFlow [21], with its basic representations of pipe models supporting reverse flow, is therefore extended by user-friendly input parameters which can be ascribed to commonly known pipe parameters for insulated and uninsulated pipe types.

Dynamic Thermo-Hydraulic Model of Distribution Pipes

Pressure drops $\Delta p$ in straight pipes can generally be described with the following relation between the pressure loss coefficient $\xi$, the pipe length l, the inner pipe diameter $d_i$, the medium density $\rho$, and the average flow velocity $\overline{v}$ [46]:

$$\Delta p=\xi {\displaystyle \frac{l}{d_i}}·{\displaystyle \frac{\rho {\overline{v}}^2}{2}}.$$

(23)

For the mapping of pressure drops within pipes, three different and selectable calculation approaches are implemented. The first method, for hydraulically rough flow regimes, is based on the studies by Nikuradse; see [47]. A second method suited for turbulent flows in hydraulically rough and smooth flow regimes is the implicit Colebrook–Prandtl equation [48]. As this implicit formula enlarges the computational effort in network calculations, an explicit approximation (Swamee–Jain) given by [49] is implemented as a third approach; see Equation (24). Parameter s denotes the inner wall roughness. This is chosen as the default calculation method for pressure drops within the pipes. Additional pressure drops caused by armatures can be entered in the respective models by the user.

$$\xi ={\displaystyle \frac{0.25}{{\left(lo{g}_{10}(\frac{s}{3.71·d_i}+\frac{5.74}{R{e}^{0.9}})\right)}^2}}$$

(24)

All pipes in the distribution network are discretised into n segments of equal length along the pipe axis. Pressure drops are calculated for each segment, and are subsequently summed up for the entire length of the pipe.

As for the description of the thermal part, a dynamic model with consideration of the thermal capacity of the surrounding soil is implemented; an approach also adopted in [11]. The thermal energy exchange is calculated for each axial segment of the pipe and is summed up subsequently. The fluid temperature at the outlet of each axial segment and the temperature at the ground surface define the driving temperature difference for this calculation. In the described thermal model, the mutual influence of neighbouring pipes is not taken into account. A serial connection $R_{th,pipe}$ of thermal resistances between the inner pipe wall (diameter $d_i$, radius $r_i$) and its outer wall (diameter $d_{out}$, radius $r_{out}$) is calculated. The entered specifications of the pipe type, e.g., simple polyethylene (PE) pipe or insulated plastic jacket pipe, and its diameter define their value; see Equation (25), adapted from ([25], p. 35).

$$R_{th,pipe}={\displaystyle \frac{1}{2\pi r_i{k}_{conv,i}}}+{\displaystyle \frac{ln\left(\frac{d_{out}}{d_i}\right)}{2\pi {\lambda }_{eq}}}$$

(25)

Here, ${\lambda }_{eq}$ denotes the equivalent thermal conductivity of the pipe, incorporating all material layers. Equation (26) gives the relation for an example of three different layers, constituted by medium pipe, insulation, and casing, for example, with their respective inner and outer diameters and thermal conductivities (reformulation of the radial heat conduction through pipes in ([25], p. 36).

$${\lambda }_{eq}={\displaystyle \frac{ln\left(\frac{d_{out}}{d_i}\right){\lambda }_1{\lambda }_2{\lambda }_3}{{\lambda }_2{\lambda }_{3}ln\left(\frac{d_{out,1}}{d_{i,1}}\right)+{\lambda }_1{\lambda }_{3}ln\left(\frac{d_{out,2}}{d_{i,2}}\right)+{\lambda }_1{\lambda }_{2}ln\left(\frac{d_{out,3}}{d_{i,3}}\right)}}$$

(26)

Enclosing the pipe, ring-shaped elements representing the soil are modelled. These consist of serial connections of thermal resistances, $R_{th,j,a}$ and $R_{th,j,b}$, and thermal capacities, as shown in Figure 4. The total width of the modelled soil region, $s_g$, can be entered by the user, but is set to a total radial extension of $0.5$ m by default.

According to [50], the thermal resistance $R_{th,S_L}$ between a single buried pipe and the ground surface can be calculated with Equation (27). $S_L$ constitutes a shape factor for this geometric condition.

$$\begin{array}{cc}\hfill R_{th,S_L}& ={\displaystyle \frac{1}{{\lambda }_g{S}_L}}\hfill \\ \hfill S_L& ={\displaystyle \frac{2\pi }{arcosh\left(\frac{H}{d_{out,pipe}+s_g}\right)}}\hfill \end{array}$$

(27)

From the relations given in ([25], p. 36), the thermal resistances $R_{th,j,a}$ and $R_{th,j,b}$ within a single radial element j can be calculated:

$$\begin{array}{cc}\hfill R_{th,j,a}& ={\displaystyle \frac{ln\left(\frac{d_{out,j}}{d_{i,j}+\frac{d_{out,j}-d_{i,j}}{2}}\right)}{2\pi {\lambda }_g}}\hfill \\ \hfill R_{th,j,b}& ={\displaystyle \frac{ln\left(\frac{d_{i,j}+\frac{d_{out,j}-d_{i,j}}{2}}{d_{i,j}}\right)}{2\pi {\lambda }_g}}\hfill \end{array}$$

(28)

with ${\lambda }_g$ being the thermal conductivity of the soil.

2.4. Sub-Models for Prosumers

The implemented prosumer models incorporate functionalities for technical equipment such as roof-mounted systems and stratified hot water storage. Figure 5 shows an example of a prosumer model with both roof-mounted ST and photovoltaic (PV) systems and two storage models operated at different temperature levels. The most important components of the equipment are indicated in the icon itself, as shown in the upper left in Figure 5. The prosumer models feature connections to the external 5GDHC network as well as links to externally imported weather boundary conditions. An electricity connection enables the computation of balances of necessary imports from or possible exports to an external electricity grid. Explicitly modelled medium flows are differentiated visually from the routing of mere continuous or discrete signals between components. Their relative temperature levels are indicated from red to blue colouring for decreasing temperature.

In the zoomed-in part, the grid connections can be seen on the lower left side. The circulation pump, with the possibility to regulate reverse flow, can be seen upstream of the heat pump model, in the warm line indicated in red. Between these two components, the free-cooling module is located. In applications with prior heating operation for the heat pump, this configuration is used to slightly increase the evaporator temperature in the heat pump and therefore possibly contribute to raising its efficiency. Above the grid connections, sub-models for the calculation of heating and cooling demands from given time series or based on the weather conditions can be seen (see Section 2.4.1). On the upper left side, the roof-mounted systems are indicated. On the right, two storage units are displayed, the upper one for DHW and the lower one for heating water. They can be charged independently by the heat pump or by the ST system and operate at different temperature levels. The energy management system (EMS) is located in the centre. It contains all logical connections for the implemented control strategies and gathers the input parameters for the building equipment. The secondary side of the HP is not modelled explicitly and the supplied heat is transferred to the storages and the room heating system at the equation level.

In the depicted prosumer model, the ST system is used in a priority circuit: As long as its supply temperature is sufficiently high to feed into any of the storage units, the relevant valves are actuated and the heat pump operation is suspended. Storage temperatures are monitored and the feed-in of thermal energy by ST systems or the heat pump is limited by a two-step control hysteresis. If the return temperature of the ST system after releasing heat to the storage is still high enough to feed into the 5GDHC network circuit, the displayed heat exchanger is used. Electrical power consumption (negative sign) and production (positive sign) are linked at junctions; overshoot power and residual demand are exchanged with the external electricity grid.

Another prosumer model is shown in Figure 6. A heat exchanger connecting the roof-mounted PVT system upstream of the heat pump evaporator with the warm line is depicted. Only one storage unit for DHW is included in the prosumer model as the room heating system is realised as a floor heating system, constituting a storage capacity itself in practical applications. The model thus lends itself especially to single-family houses. The heat pump is implemented as an inverter-regulated model, directly providing the demanded heating power for the floor heating system and the DHW storage. The uncovered PVT system is only able to provide supply temperatures, which support the primary side of the heat pump. In order to elevate the evaporator entry temperature during heating operation, it is located upstream of the heat pump. If the PVT system can provide supply temperatures above the network temperature and the heat pump is not active, the circulation pump is switched on with a pre-defined volumetric flow. On the one hand, this operational mode is used for regeneration of other heat suppliers connected to the network, like the GHC. On the other hand, the shift of overshoot thermal energy into the grid makes it available to other prosumers with simultaneous heating demands in the network.

2.4.1. Heating and Cooling Demands

Room heating demand, DHW demand, and free cooling demand are modelled separately. Apart from reading direct time-series-based profiles from externally supplied files, a methodology for creating heating and cooling demands dependent on the weather conditions is implemented in the sub-models of prosumers.

Time series for the room heating demand are created with a methodology provided in [51], considering the actual weather conditions and the energetic efficiency of the building. The used approach translates an adjustable annual room heating demand, $Q_{dem,rh,ann}$ into an hourly distribution $Q_{dem,rh,h}$. Prosumer models can be assigned pre-defined sets of parameters for a sigmoid function, typical of its type of building usage such as single- or multi-family houses (SFH, MFH) or commercial usage. The sigmoid function gives an output value factor $F_{sigmoid}\left(T_{ref}\right)$, dependent on a multi-day weighted average of the ambient air temperature. This is further altered with coefficients dependent on the weekday for the building usage type, $F_{day}$, as well as daytime-specific coefficients, $F_{hour}$, distributing the room heating demand to the hours of the day. The final hourly room heating demand is calculated by Equation (29) and passed to the substation model of the prosumer. It encompasses the procedure provided in ([51], p. 44).

$$Q_{dem,rh,h}={\displaystyle \frac{Q_{dem,rh,ann}}{{\int }_0^{8760h}{F}_{sigmoid}\left(T_{ref}\right)d\tau }}·F_{sigmoid}\left(T_{ref}\right)·F_{day}·F_{hour}$$

(29)

Target temperatures for room heating, defining sink side temperatures for the heat pump, are calculated dependent on a gliding average value of the ambient temperature. DHW demand profiles can be provided to the model as text-based time series for tapping profiles. As a default, tapping profiles for households according to [52] are used. These time series are passed to the sub-models of DHW storage in the prosumer models.

Hourly time series for free cooling demand are generated prior to the simulations by scaling a defined annual cooling demand for the prosumer, $Q_{dem,fc,ann}$, with regard to the time series of ambient temperature. With a prescribed limit $T_{fc,lim}$ to the ambient temperature, above which a cooling demand may exist, the free cooling demand $Q_{dem,fc,i}$ at hour j is computed as follows:

$$Q_{dem,fc,j}=min\left(Q_{max,fc,h};Q_{dem,fc,ann}·{\displaystyle \frac{max(0;T_{amb,j}-T_{fc,lim})}{{\sum }_{j=0}^{8760}max(0;T_{amb,j}-T_{fc,lim})}}\right).$$

(30)

This methodology can similarly be found with the proposition of cooling degree days in [53,54]. The free cooling module is modelled as an ideal heat exchanger transferring heat to the network medium whenever a free cooling demand is present, and can be covered with respect to operational limitations. A maximum value for the cooling power can be set to limit the extractable power in free cooling operation. By default, a fixed value of 20 $\frac{\mathrm{W}}{{\mathrm{m}}^2}$, related to the total installed cooling surface, is chosen in the prosumer model. Another limitation to the use of free cooling is the definition of a maximum network temperature of 19 °C by default. This ensures a minimum temperature difference between the medium in cooling circuits in the building and the target room temperature during a free cooling operation.

If only a free cooling demand should be covered and the heat pump is not active, the circulation pump regulates a pre-defined volumetric flow to the free cooling module, which is set to half the nominal volumetric flow of the heat pump by default. In case the heat pump is active simultaneously, the circulation pump provides the nominal volumetric flow for heat pump operation.

2.4.2. Decentralised Substation Models

The decentralised substations primarily function as a brine/water-to-water heat pump. At the primary (evaporator) side, the heat pump model is connected to the network. A circulation pump provides user-specified nominal volumetric flow to the heat pump when it is active. Both the circulation pump model and the heat pump model can handle reverse flow on the evaporator side. If the residual possible pressure head is exceeded in an operational state, a warning is raised by the model. In this way, the heat pump model is compatible with both active and passive control strategies for the 5GDHC network operation, as it can be supported by an external network pump. Further input parameters are the nominal heating power, an internal pressure loss coefficient of the evaporator, and an efficiency factor ${\eta }_g$, which can either be fixed or variable with regard to the partial load ratio, following a definable functional relation. The sink-side supply temperature for the substation, $T_{HP,sink}$, is determined by the current set temperature of the DHW storage or the warm water storage for heating purposes depending on the operational mode. The coefficient of performance (COP) of the heat pump is computed using the temperatures on the evaporator side, $T_{HP,source}$, and on the supply side, $T_{HP,sink}$, combined with the efficiency factor, ${\eta }_{HP}$ ([55], p. 578). The electrical power $P_{HP,el}$ drawn by the compressor is calculated from the COP value and the demanded heating power, ${\dot{Q}}_{HP,heat}$.

$$CO{P}_{HP}={\eta }_{HP}{\displaystyle \frac{T_{HP,sink}}{T_{HP,sink}-T_{HP,source}}}={\displaystyle \frac{{\dot{Q}}_{HP,heat}}{P_{HP,el}}}$$

(31)

Note that the COP is prescribed a value of 1 if the minimum source-side temperature defined for the heat pump is undercut. In this way, the operation of an additional heating rod as a backup in the heat pump model can be indicated. Seasonal coefficients of performance (SCOP) for the heat pump operation are calculated by delivered heating energy and needed electrical energy in a given time span $\Delta \tau$; see Equation (32) from ([55], p. 589).

$$SCO{P}_{HP,\Delta \tau }={\displaystyle \frac{{\int }_{\Delta \tau }{\dot{Q}}_{HP,heat}d\tau }{{\int }_{\Delta \tau }{P}_{HP,el}d\tau }}$$

(32)

3. Results

This section provides results from the example simulations of a district. The input parameters to the sub-models are presented and results for two case studies are given. Case study 1 comprises results for variations in the GHC area and different environmental boundary conditions. Case study 2 gives a detailed insight into the energy balances and energetic system performance for a simulation run with fixed parameter settings.

3.1. Model Description, Demand and Supply Structure

The case studies comprise the simulation and analysis of a fictional district of new buildings and a ring network with warm and cold lines. It is designed as a passive 5GDHC network with uninsulated PE pipes, with each prosumer having an individual circulation pump. The network topology can be seen in Figure 7. For all further information on the design parameters and constant material properties, refer to Appendix A. Note that the material properties of the working media are kept constant as well.

In total, four different models for residential buildings are used in the case studies. Each of the models represents an aggregation of several single buildings to simplify the network topology for this example. The residential buildings are equipped with roof-mounted systems which provide heat either directly to the DHW and warm water storage and the network or solely to the network. A supermarket supplying excess heat from process cooling applications into the warm line is connected to the network as well. The total network length is 960 m. All pipes are modelled as dynamic thermo-hydraulic representations without coupling to neighbouring pipes, as shown in Figure 4. The pipe lengths are defined with respect to the connected aggregations of buildings, so that distances between single buildings are summed up to the length of a pipe section. On the left in Figure 7, the GHC can be seen. Above its icon, the sub-model for loading weather data for a location close to the city of Giessen, Germany, is shown. All components and prosumers are connected to the external electrical grid, depicted next to the weather data sub-model.

In

Table 1. Parameter settings and equipment of prosumer models for the annual simulations.

Prosumer	Aggregated Buildings	Room Heating (rh) Demand (MWh)	DHW Demand (MWh)	Cooling Demand (MWh)	Set Temperatures (DHW, rh)	Total Storage Volume (DHW, rh)	Roof-Mounted Systems ($\mathit{\gamma },\mathit{\phi }$)
MFH_001	3	70.20	21.94	7.02	58 °C, 28 °C to 40°C	$4.05$ m3, $9.24$ m3	$116.5$ m2 ST modules (a) ( 45 °,  15 °)
MFH_002	4	93.60	29.25	9.36	58 °C, 28 °C to 40 °C	$5.40$ m3, $11.83$ m3	$149.12$ m2 ST modules (a) ( 45 °,  15 °)
SFH_001	10	63.00	21.00	8.40	48 °C, 26 °C to 35 °C	$3.00$ m3, -	$65.4$ m2 PVT modules (a) ( 30 °,  10 °)
SFH_002	12	84.00	25.20	10.08	48 °C, 26 °C to 35 °C	$3.60$ m3, -	$78.48$ m2 PVT modules (a) ( 30 °,  10 °)
Supermarket_001	1	-	-	25.00	-	-	-

^(a) Aperture area.

, the parameters for the prosumer models, namely, the residential buildings and the supermarket, are summarised. The multi-family houses (MFHs) feature the same equipment, as explained in Figure 5, with two storage units, one for warm water, covering room heating demands, and one for DHW production. The roof-mounted ST systems feed directly into the storage or deliver excess heat into the network. A total of 16 ST modules are installed on each MFH. The PV systems are assumed to be turned off for all analyses. The single-family houses (SFHs) feature the equipment shown in Figure 6, with floor heating systems, one storage unit for DHW production, and roof-mounted PVT systems delivering heat to the evaporator side of the heat pump and the network. Here, three PVT modules are installed on each SFH. The storage unit sizes are determined dependent on an optional assistance by ST systems: Warm water tanks for room heating designed with a volume of $0.08$ m³ per m² ST collector area. This dimensioning is an established approximation found, e.g., in [56]. The DHW storage units are sized in such a way that they store 120% of the daily DHW demand without ST systems and 200% with assistance by ST systems. Example research dealing with an optimum design of solar-assisted storage can be found in [57,58].

The supermarket is assumed to deliver its excess heat from process cooling applications to the network based on a daily recurrent time series. For reasons of simplicity, the heating demand of the supermarket is not taken into account and is assumed to be covered by an individual system.

As the primary central heat source a one-dimensional model of a GHC is connected to the network. The installation depth is set to $1.5$ m. The used soil types are silty loam up to a depth of 3.3 m, followed by loamy sand up to 8 m, and basalt from there on up to the boundary at 25 m depth. The results are extracted from the last complete year out of a simulated three-year period. This allows for the complete model to reach steady-state conditions on an annual basis. It is necessary to extend the simulated time to this period since the modelled soil regime and its capacities in the GHC sub-model are equipped with fixed starting values for temperatures due to unknown temperature distributions at the starting time.

Specifics of Case Study 1

In case study 1, the focus is put on the behaviour of the GHC model within a district simulation. The following specifications and variations are applied in this case study:

• The GHC area is varied between 5500 m² and 1000 m². The minimum temperature in the warm network line at the connection points to the prosumers is logged as an indicator for shortfalls below common limits of heat pump source-side temperature values of $-5$ °C.
• All decentralised, roof-mounted systems in the prosumer models are turned off. This isolates the effect of the free cooling operation on the performance of the decentralised heat pumps. For this case, the supermarket is also detached from the network and only the residential buildings remain connected.
• The input for soil moisture is varied between a reference condition and a dry condition.
• Weather data are varied between a reference dataset and an alternative dataset with a particularly cold winter period, both provided by DWD [39]. This is used in order to assess the influence of temporarily rough climatic conditions on the performance of the GHC.

Figure 8 shows the set values for the water content to the soil zones, varying temporally and with depth. They are expressed in relation to the maximum water content at the saturated state, ${\displaystyle \frac{\mathsf{\Theta }}{{\mathsf{\Theta }}_{sat}}}$. The volume-weighted average of the effective field capacity (efc) up to a depth of 1 m is also given. The data for the efc at near-surface depths and at reference soil conditions are taken from [59]. The reference conditions are based on a constant ratio of ${\displaystyle \frac{\mathsf{\Theta }}{{\mathsf{\Theta }}_{sat}}}=60\%$ for depths greater than 5 m. The dry soil conditions, in contrast, are based upon a diminished ratio of 45% at these depths, and moisture values in more shallow regions are reduced by approximately 15% as well. A damping of temporal variations in the water content is applied between a depth of 1 m up to this constant value in 5 m, which considers the decreasing influence of precipitation and evaporation with growing distance from the ground surface.

Figure 9 juxtaposes the reference TRY dataset and the alternative TRY dataset with a distinct reduction in ambient air temperatures in winter time. The displayed temperature values of the datasets are equivalent to the multi-day average $T_{ref}$ relevant for the calculation of heating demands; see Section 2.4.1.

The colder temperatures of the winter time in the alternative dataset start in October of the second simulated year and last until the end of the third simulated year.

Specifics of Case Study 2

Case study 2 focuses on the performance of decentralised heat pumps and the overall system performances for a single simulation run. Multiple decentralised heat sources and a heterogeneous prosumer configuration are incorporated.The following specifications are applied in this case study:

• The GHC area is fixed at 3000 m².
• All prosumer models, including the supermarket, as well as roof-mounted PVT and ST systems Table 1 are in service. Free cooling operation is activated in every residential building.
• The input for soil moisture remains at the reference condition provided in Figure 8.
• Weather data from the reference TRY dataset are used as boundary conditions.

3.2. Results from Case Study 1

3.2.1. Variation in GHC Area at Different Environmental Boundary Conditions

In Figure 10, the simulation results from case study 1 are shown. For varying GHC areas, $A_{GHC}$, the district-wide annual performance factors of decentralised heat pumps, $SCO{P}_{HP}$, are given on the left axis. The data from simulations with reference soil conditions and activated free cooling are displayed in blue, those with reference soil conditions and deactivated free cooling are plotted in dashed orange lines. The results with dry soil conditions and activated free cooling are pictured in green dotted lines. The analysis with reference soil conditions and activated free cooling but two subsequent cold winter periods from the alternative TRY dataset is plotted in red dashed lines. Corresponding values for area-specific net extracted energy from the GHC are displayed on the right axis with the same line styles and in gray colour. These are calculated as the integral of discharging (positive) and charging (negative) power to the GHC over the span of a year in relation to the projected GHC area.

In general, starting at greater GHC areas, the courses of $SCO{P}_{HP}$ values feature a flat gradient before a distinct decline is observable for GHC areas below a certain limit. This qualitative course is also reported by Hüsing et al. [60]. For the three variants with the reference TRY dataset, this limit can be spotted between 2500 and 2000 m² GHC area. Between 2000 m² and 1500 m², the $SCO{P}_{HP}$ values drop from 4.1 to 3.65 at reference soil conditions and from 3.95 to 3.48 at dry soil conditions, respectively. At 2500 m² GHC area, the net extracted energy reaches values of 54 $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$ for activated free cooling and reference soil conditions, 58 $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$ for deactivated free cooling and reference soil conditions, and $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$ at activated free cooling and dry soil conditions.

When network temperatures at the connections to prosumers fall below a limit of $-5$ °C, a heating rod supports the decentralised heat pump operation, consequently reducing the overall performance of the heat pump. These points are indicated as big red triangles. A sufficient dimensioning of the GHC can only be found for GHC areas greater than 2500 m² without violation of this limitation for simulations with the reference TRY dataset. As for the influence of the free cooling operation, it is obvious that although the area-specific extraction decreases for activated free cooling at a fixed GHC area, the change in $SCO{P}_{HP}$ values is hardly notable. This result implies that a temporal resolution of actual cooling demands is crucial when their influence on possible downsizing of the GHC as the primary heat source is examined. If the cooling energy is not distributed evenly across the course of a year, as is the case for many room cooling applications, a simplified annual reduction in necessary source heat does not reflect the effective load on a GHC in the heating period correctly. Furthermore, it can be derived that energy from a pure free cooling operation in residential buildings on a practical scale is not an appropriate measure for reducing the capacity of a primary heat source in 5GDHC networks.

Examining the alternative TRY dataset with two consecutive cold winter periods as the boundary condition to the simulation, the limit for a minimum reasonable dimensioning of the GHC is shifted to greater GHC areas, above 3500 m² In the mentioned region, with comparatively small gradients in $SCO{P}_{HP}$ values and a sufficiently dimensioned GHC area, $SCO{P}_{HP}$ values are from 0.1 to 0.2 lower than for the reference TRY dataset. This can be explained by the diminished efficiency of the heat pump operation for higher sink-side temperatures. As the ambient temperatures in the heating period are lower, the target temperature for the heat pumps increases, thus reducing their efficiency. The course of the area-specific net extraction of this variant begins to fall below the values obtained for variants with the reference TRY dataset at GHC areas smaller than 2500 m² despite lower performance factors. This fact suggests an increasing proportion of heating demand is covered by heating rod operations in the decentralised heat pumps at smaller GHC areas.

Black crosses with reference to the left axis show data derived from VDI 4640-2 [61]. For an assumed $SCO{P}_{HP}$ of 4.25 in climate zone 7, the cited methodology reports a maximum value of 51.5 $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$ (values interpolated between the soil types loam and silt). The necessary source heat is calculated from the annual total heating demand of the district and the $SCO{P}_{HP}$ value. Reducing the source heat energy by the annual demand for free cooling and then further by the thermal network gains varied between 0% and 20%, yields the displayed results. The dimensioning rules following VDI 4640-2 do not take into account interactions of multiple heat sources or the temporal course of thermal gains in the distribution network and are based on coarse assumptions when applied by a user. The results of this case study acknowledge a further potential for downsizing GHC areas compared to VDI 4640-2. Choosing a GHC area of 3500 m² as a proper dimensioning, the reduction in GHC area can amount to 18%, for instance.

A prioritisation of individual, simulation-assisted design validation of the configuration of heat sources, consumers, and the distribution network over simplified dimensioning procedures is suggested when planning 5GDHC networks.

3.3. Results from Case Study 2

3.3.1. Simulation Performance

The simulated model of the district in OpenModelica consisted of 9602 equations in total in case study 2. A total of 7357 parameters were included in the mathematical model description. The total time for the simulation of the three-year period was 105 min. Simulations were performed on a Windows 11 system with specifications of i7-9700, 8 cores at 3 GHz, and 16 GB RAM. The solver DASSL in OpenModelica was chosen for the simulations. The biggest allowed computational time step was set to 20 min.

Abugabbara et al. [14] state a CPU time for annual simulation of their 5GDHC network of $4.56$ min. The used simulation framework is Dymola and the specifications of the used hardware are not mentioned. The number of building models is quoted as nine, all of which are connected within a ring network and have a balancing unit serving as the heat and cold supply. Sub-models of stratified storage tanks and building equipment other than heat pumps, like in the present study, are not included. Although the presented number of equations of 23,544 is higher than in the study at hand, this can be attributed to their use of the Modelica.Fluid library. This library generally features a higher amount of internal equations in sub-models for pipes, valves, and for the calculation of medium properties.

Heissler [11] quotes a simulation time of 7 to 9 days for the conducted co-simulation of 2 years of simulated time. A constant simulation time step of one minute is chosen in the presented approach. As the simulated district features a borehole heat exchanger field as the primary heat source, the comparability to the present study is further diminished.

The higher simulation time of the present study than in the example from [14] could originate from the detailed GHC model and other sub-models of higher discretisation order like the water storage tanks. Especially, the coupling of the GHC model to the network—and indirectly to the building equipment—leads to increased computational effort. The fine discretisation of soil cells near the collector pipe makes small internal time steps necessary for stable numerical computations. However, the achieved simulation time of 105 min for a simulated three-year period can still be seen as an efficient result. It enables us to conduct parameter studies and obtain access to detailed time-series results from individual buildings in a reasonable time span. Compared to other works on holistic simulation approaches for 5GDHC networks, there is a vast span of quoted simulation times noticeable, mostly dependent on the level of detail of the employed sub-models.

3.3.2. Network Temperatures and Thermal Network Gains

In Figure 11, the hourly course of temperatures in the network and the ambient air temperature are shown. In general, supply temperatures in the warm network line differ slightly between the prosumer models SFH_001 and MFH_001, which can be accounted for by the varying flow distance through the distribution network to these building representations. As these differences are very small, $T_{warmline,SFH001}$ lies behind the course of $T_{warmline,MFH001}$ in the displayed course of temperatures and cannot be spotted. It can be seen that the supply temperature of the GHC to the warm network line during the heating period falls below 0 °C, reaching a minimum of $-2.52$ °C in hour 980. The highest temperature at the connection of the GHC to the warm network line reaches $20.33$ °C in hour 5846. The maximum temperature difference between the supply and return temperature to the GHC is at $1.44$ K. This small temperature difference is due to an elevation of network temperatures close to the ground temperature along the routing of the distribution pipes.At the beginning of April (around hour 2160), there is a two-week period with relatively warm ambient temperatures of up to 22 °C. With a slight temporal shift, the elevation of the network temperatures can be stated to be in this range as well. Especially, the rise in supply temperature from the warm line for prosumer MFH_001 from 0 °C to 4 °C is notable. A parallel rise in the temperature in the cold network line above the temperature of the warm network line can be attributed to the increased thermal yield from decentralised ST and PVT systems feeding into the cold line. Equivalent increases in the cold network line temperature can be seen during summer time, where amplitudes of up to 5 °C above the warm network line temperature occur on an hourly basis. Generally, network temperatures tend to harmonize in the summer period as hardly any room heating demand is present. Note that there is no temperature set point or a fixed design temperature difference for any of the components in the network.

Figure 12 shows the monthly transfer of thermal energy between the GHC, and the decentralised ST and PVT systems, as well as the heat exchange between the distribution pipes and the surrounding soil. The net thermal yield of the GHC transferred to the network can be interpreted as the sum of the (positive) amount of discharged energy and the (negative) amount of charged energy. The charging contributes to the thermal regeneration of the soil regime in the GHC model while the discharging leads to continuous heat extraction and potential soil freezing. The thermal gains of the distribution pipes in the network hold positive values, indicating an inflow of thermal energy into the network. The same principle applies to the energy supplied by the decentralised ST and PVT systems, marked in yellow colours. Especially, in the spring time, starting in April, the GHC experiences a net regeneration in this monthly balance. In this period, the thermal gains of the distribution pipes in the network reach their highest monthly value, at more than 14 MWh. As stated in Section 3.3.3, the room heating demand subsides drastically in this time period. In parallel, ambient air temperatures rise after the heating period and the temperature in the upper soil regions is elevated. As the circulating network medium is still relatively cold, the driving temperature difference between the medium and the surroundings increases and the thermal network gains are raised. In the summer months, from June until late September, the network generates thermal losses as the medium temperature is often increased above the temperature of the surrounding soil due to energy fed in by the solar-powered systems.

The decentralised PVT systems feature a peak in thermal energy yield in spring time. In this period, their potential to transfer heat into the network is highest due to the elevated ambient air temperatures and rising solar irradiance values. Combined with the low network temperatures, the PVT systems can operate at low module temperatures, reducing their inherent thermal losses to the ambient surroundings and consequently increasing the overall thermal efficiency. The excess energy from cooling operations reaches a peak in June, when demand for free cooling in residential buildings is already present. At the same time, network temperatures are still low enough to fulfil all requirements to activate free cooling. The share of cooling energy fed in by the supermarket is constant throughout the year by definition. The annual sum values are reported in Section 3.3.3.

3.3.3. Energy Balances and System Performance

In Figure 13, the course of heating demands for room heating and DHW production is shown. The demand for free cooling, which is covered within the operational limits, is also displayed. Values for $SCO{P}_{HP}$ as a performance factor for the heat pump heating operation are given. These are complemented by a system performance factor $SCO{P}_{sys}$ for heating and free cooling purposes in the prosumer buildings. The combined $SCO{P}_{sys}$ at the building level takes into account all benefits constituted by useful heat for room heating and DHW production, $Q_{heat,rh,DHW}$, and for covered free cooling demand, $Q_{fc}$. These are related to the expenses in the form of electricity imported, $W_{el,import}$, at a building level according to Equation (33). The electrical energy for decentralised circulation pumps transferring excess heat to the network during free cooling operation, and the operation of ST systems, is considered in the expenses as well. It can be interpreted as an overall efficiency of the heating and cooling system. A similar definition and naming for a heating system with a GHC is found in [62].

Owing to the hydraulic setting where heat pump operation, free cooling, and thermal network regeneration by ST systems may occur simultaneously, a dedicated attribution of electrical efforts for the circulation pump to heating, cooling, or ST operation is not possible (see Section 2.4.1).

$$SCO{P}_{sys}={\displaystyle \frac{Q_{heat,rh,DHW}+Q_{fc}}{W_{el,import}}}$$

(33)

On the district scale, $SCO{P}_{sys}$ sums up all net electrical energy imports into the district and relates them to all benefits stemming from useful heating and free cooling energy combined, as can be read in the equation. In this way, all residual electrical demands for heat pumps and circulation pumps are properly considered.

An electrical coverage factor ${\eta }_{el,heat,fc}$ assesses the amount of energy consumed by the complete heating and cooling system, which is covered by the electrical yield from PVT systems. It is calculated from the imported amount of energy and the total electrical consumption at the building and district levels according to Equation (34) [63].

$${\eta }_{el,heat,fc}=1-{\displaystyle \frac{W_{el,import}}{W_{el,dem}}}$$

(34)

It can be seen that this district-wide system $SCO{P}_{sys}$ reaches values of up to 12.8 in June, exceeding the highest combined $SCO{P}_{sys}$ at the building level, which is slightly below 10. The annual value of the district-wide system $SCO{P}_{sys}$ is 5.98. This can be attributed to the production of electrical energy by the decentralised PVT systems, which is shared with all prosumers connected to the district electrical grid.

Figure 14 shows the partitioning of heating demand and its proportionate coverage by the different heat sources in the district investigated in case study 2. It can be seen that the total heating demand of the prosumers, shown on the left side, is either directly provided for by the decentralised ST systems or by the heat pumps. Thermal losses in the storage and occasional charging of the warm water storage in the MFH models during summer time (when hardly any heating demand is present) lead to the sum of heat supply by the ST systems and the heat pumps exceeding 100% of the total heating demand. The decentralised heat pumps feature an average annual $SCO{P}_{HP}$ of 4.43. The total amount of source-side heat demand taken from the network is $302.8$ MWh. This sum can further be divided into the contributions by the feed-in of decentralised ST systems (21.2%), the feed-in of decentralised PVT systems (16.7%), the net thermal gains collected in the distribution network (18.9%), and the feed-in by free cooling (3.8%). Whereas the supermarket provides 8.3% of the source heat demand (denoted as misc. prosumers), the central GHC, as the primary heat source, only needs to contribute a further $95.2$ MWh or 31.5%. The area-specific net extraction of the GHC throughout the year adds up to a value of 31.7 $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$. The electrical efforts of the decentralised circulation pumps add up to $8.98$ MWh, approximately 2.2% of the total heat demand in the district.

Table 2. Results for prosumer models and complete district from case study 2.

Component	Heating Demand, $Q_{\mathit{heat},\mathit{rh},\mathit{DHW}}$ (MWh)	Th. Yield ST System (MWh)	Th. Yield PVT System (MWh)	El. Import $W_{\mathit{el},\mathit{import}}$ (MWh)	El. Export $W_{\mathit{el},\mathit{export}}$ (MWh)	El. Coverage ${\eta }_{\mathit{el},\mathit{heat},\mathit{fc}}$ (%)	${\mathit{SCOP}}_{\mathit{HP}}$	${\mathit{SCOP}}_{\mathit{sys}}$	${\mathit{SCOP}}_{\mathit{HP}}$ Reference AWHP
MFH_001	92.14	41.34	-	22.43	0	0	4.02	4.22	4.01
MFH_002	122.85	51.79	-	29.84	0	0	4.03	4.19	4.01
SFH_001	84.00	-	22.33	16.64	12.01	17.4	4.87	5.23	4.75
SFH_002	109.20	-	28.34	20.41	14.03	18.3	4.96	5.53	4.89
Supermarket_001	-	-	-	0.03	0	0	-	-	-
District	408.10	93.13	50.67	86.29	22.98	3.4	4.43	5.98	4.38

summarises the further results and energy balances for the prosumers in the presented case study. The values of the electrical coverage factor show that at the district level, 3.4% of the electrical demand is covered by interchange between the prosumers.The $SCO{P}_{HP}$ for the operation of a reference AWHP is given in the data. Its operational state is identical to the conditions the regular heat pump connected to the network meets. As source-side temperature, however, the reference AWHP accesses the current ambient air temperature.

4. Discussion

The developed simulation model for 5GDHC networks and the contained sub-models for a GHC and prosumers described in more detail in this work are used in two simulation case studies of an example district. In the presented simulation runs, the capabilities of the integrated simulation model are laid out. Notable points for discussion from the presented results are addressed and recommendations for the scope of applicability of the developed simulation model are given.

The detailed model of the GHC used in the case studies enables the user to enter specifics of the soil regime at the location investigated with simulation assistance. Its default automatic discretisation control supports users in finding a justifiable trade-off between the granularity of the results and computational effort. This sub-model increases the simulation duration, especially when coupled with detailed prosumer models and meshed network topology. Nonetheless, this disadvantage is seen as a reasonable compromise as there is no need for externally created or approximated time series of supply temperatures of very shallow geothermal heat sources. The consideration of latent heat released during freezing periods comprises widely used methods based on van Genuchten parameters, available for a variety of common soils and extensible to user-specific parameter records. Beneficial future upgrades to the model might encompass a more detailed modelling of moisture transport within the soil regime and the consideration of thermal asymmetries, leading to further thermal gains at the collector boundaries. In addition, an ongoing monitoring of a real GHC installation close to the town of Giessen, Germany, is expected to support the researchers with further insights into the actual development of supply temperatures of a GHC to a 5GDHC network. With information from this monitoring, characteristics of the heat exchange described in Section 2.2.1 might be tuned for a more accurate representation of outlet temperatures from the necessarily simplified GHC model.

Case study 1 puts the focus on variations in the (environmental) boundary conditions and their effect on the performance of a GHC as the primary heat source in a 5GDHC network. Both annual, district-wide performance factors of the decentralised heat pumps and the net extracted amount of energy from the GHC are examined for this purpose. Settings for typical and dry soil conditions are juxtaposed and an alternative dataset for weather conditions with two consecutive cold winter periods is applied as a boundary condition.

In ([10], p. 78), a study on extractable energy amounts from a GHC as ground-coupled ice storage is described. The results from this examination show values for area-specific extractable energy ranging between 35 $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$ and 66.5 $\frac{\mathrm{kWh}}{{\mathrm{m}}^2}$ for a single-layer GHC with the soil type silty loam and within a typical range of moisture values. In [64], the performance of a shallow horizontal GHC under agricultural land as a single central heat source is investigated. A value for area-specific extractable energy from the GHC of 40 $\frac{\mathrm{k}\mathrm{W}\mathrm{h}}{{\mathrm{m}}^{\mathrm{2}}}$ to 70 $\frac{\mathrm{k}\mathrm{W}\mathrm{h}}{{\mathrm{m}}^{\mathrm{2}}}$ is reported there. The findings from the present study with reference soil conditions fit well in these quoted ranges, as stated in Section 3.2.1.

The general course of SCOP values of decentralised heat pumps shown in Section 3.2.1 features a range of relatively stable values for greater GHC areas; i.e., a range with a sufficiently dimensioned GHC area. For smaller GHC areas, the efficiency of decentralised heat pumps drops significantly from some point. This is attributed to the augmented operation of heating rods. In [60], not only is this behaviour also reported, but also a possible reduction in GHC area compared to VDI 4640-2 [61], by as much as 25%. The findings from the present study in Section 3.2.1, with possible reductions of approximately 18%, support these values. However, disadvantageous conditions in winter periods can increase the necessary GHC area and suggest individual simulation assistance.

As for the dry soil conditions, values for effective field capacity at depths of up to 1 m of 67% in winter and 19% in summer are assumed. The results for the performance of decentralised heat pumps show only slight deviations between reference and dry soil conditions for sufficiently dimensioned GHC areas. However, it is expected that the lower limit of this range will shift to greater GHC areas for even drier soil conditions, which are disadvantageous for geothermal usage of the ground. In drought periods, the effective field capacity can commonly drop to zero and the absolute water content of the soil is reduced even more. As only one soil type in the near-surface region is examined in the case study (silty loam), it has to be stated that a greater impact on extractable energy from the GHC by reduced water content can also be expected for other soil types with different thermo-hydraulic characteristics. This can be derived from findings in [10]; especially, soil types with higher proportions of coarse pores like sand tend to perform worse in dry conditions than in wet or saturated conditions. This is attributed to their low residual water content and simultaneously low adhesive forces preventing water from freezing. This makes the water more easily accessible for energetic usage of latent heat below the freezing point and raises the performance in wet conditions.

The results of case study 2 show that 5GDHC networks with connected decentralised heat pumps can be advantageous in terms of energy efficiency over air-to-water heat pump solutions. The seasonal performance factors of the heat pumps connected to the network are found to be slightly elevated, with 4.43, compared to those of reference air-to-water heat pumps with 4.38. Though this is only investigated for a single location of the example district in the simulation and may vary with the parametrisation of the heat pump models and different climatic boundary conditions, more advantages are worth noting: The possibility to connect a variety of low-temperature heat sources to the network makes it possible to share excess heat, for example, which otherwise remains unused and is released to the environment.

Network temperatures fluctuate freely in the presented case studies as they are determined primarily by the resulting soil temperatures in the GHC and along the distribution pipes. Other studies on simulation frameworks for 5GDHC networks employing set or minimum/maximum temperatures for the network lines therefore feature different performance factors for heat pump operation and system efficiency (see [11,14]). In [64], GHC supply temperatures between $2.1$ °C and $15.9$ °C are quoted. The minimum temperatures occur in February and March, which can also be found in the results of the present study; see Figure 11. Lower variations between the minimum and maximum supply temperatures of the GHC in [64] could be explained by the deeper installation depth of 2 m and the lack of additional decentralised heat sources active in summer. A minimum of the efficiency values for the operation of decentralised heat pumps is identified at the end of the heating period between March and April. Maxima are found in September and October. The locations of maximum and minimum values are in agreement with the findings in the study at hand, as can be seen in Figure 13, especially for SFH models. The quantitative differences of the efficiency values of the cited work can be attributed to the assisting decentralised heat sources employed in case study 2.

As stated in Section 3.3.3, the electrical energy drawn by the circulation pumps in the district constitutes 2.2% of the total heat demand. A similar value of 2% is reported by Gross ([65], p. 107) for a case study of a 5GDHC network with heat extraction from a mine water system. This proportion seems reasonable, as the temperature difference between the warm and cold lines is smaller than in conventional heating networks. For these, efforts for circulation usually add up to 0.5% to 1% of the transported heat ([66], p. 135).

The free cooling from residential buildings constitutes 3.8% of the source heat or 2.8% of the total heat demand in the district; see Figure 14. In ([65], p. 112), the author states that only 1% of the total heat demand can be covered by free cooling within the simulated district. There, it is mentioned that the effects of pure free cooling applications do not lead to a size reduction in the primary heat source in the network. This implication is also found in the results from case study 1 in Section 3.2. However, temperatures are predicted to rise as an impact of climate change and cooling demands in residential buildings are expected to increase as well [67,68]. This constitutes a strong encouragement for sustainable designs of energy supply systems on a district scale to always make use of this source of excess heat.

Excess heat from a supermarket is analysed in [65], too. Its available excess heat is 7.6% of the total heat demand of the district. Compared to 6.1% in the present study (8.3% of the source heat), this is a similar size of excess heat source. A reduction in the primary heat source by the feed-in from this steady excess heat source, as suggested in ([65], p. 116), may be supported by the results from case study 2. This would be part of further investigations on reductions in the GHC area.

The resulting feed-in of free cooling does not meet the initially defined annual demand values for the prosumer models in case study 2. Although these initial demands are set to approximately 10% of the total heat demand of the prosumers, the final contribution of $11.4$ MWh fed into the network constitutes only 2.8% of the total heat demand, which is a relatively low value for the used type of energy-efficient new buildings. Firstly, this can be attributed to the limiting factors regarding network temperatures for extracting the desired cooling power, described in Section 2.4.1. The location of the feed-in of thermal energy from roof-mounted PVT systems in the prosumer models is chosen to be upstream of the free cooling module in the dedicated warm line. Therefore, a compromise between maximising the utility of the PVT systems or the free cooling module has to be made. Secondly, in times with simultaneous operation of the PVT system and heating demands covered by the heat pump, the increase in the heat pump evaporator entry temperature is beneficial to the efficiency of the heat pump. In order to reduce the complexity in the prosumer models, the implementation of separated parallel hydraulic loops for free cooling and for the coupling with the PVT systems is avoided. The potential for energy feed-in from decentralised PVT systems, especially in spring time, suggests further investigations into the use of air-to-brine heat exchangers as a more cost-efficient alternative. This suggestion is supported by the findings in [9], where the most profitable time for operation of air-to-brine heat exchangers is seen to be in spring as well. Nonetheless, the following barriers for their implementation need to be considered in 5GDHC networks: Commonly used air heat exchangers are not suited for installation on pitched roofs and can cause sound emissions, which makes them unusable for use in residential areas or means that they require design adjustments [69]. Therefore, other dedicated installation areas possibly have to be found and prepared. What is more, the electrical yield of PVT modules is advantageous when covering simultaneous electrical power demands in the district. Overshoot power can be sold externally, raising economic returns for the operator. On top of that, covered PVT systems with the possibility to directly feed into warm water storage can generate additional benefits to the operation of the decentralised heat pumps and the total network efficiency, accordingly.

Roof-mounted ST systems contribute approximately 21% to the annual source to the network in the setting of case study 2. As these units usually provide higher supply temperatures than PVT systems, the exergetic losses of a feed-in of their heating energy to a network operating close to surrounding temperatures are higher. Common ST systems may not be fit for operating below limits around 25 °C, which makes it arguable whether the 5GDHC grid is a suitable heat sink. Nonetheless, in the chosen setup of the prosumer models, a cascaded usage of the ST systems for storage at different temperature levels at a building level is already implemented. A feed-in into the warm network line is used in this setting as a means for regenerating the geothermal heat source. Alternative ways of storing and releasing their overshoot heating power at the network level, as investigated in [70], could be considered as well.

The implemented thermo-hydraulic models for distribution pipes are developed based on reviews from modelling approaches described in several other works, as stated in Section 2.3. Sub-models for quasi-stationary and dynamic thermal behaviour are contained in the simulation library and can be chosen according to the desired degree of detail for simulations. In case study 2, the contribution of thermal gains from the distribution network accounts for 18.9% of the total source heat. Related to the total heat demand, it adds up to 14%. This contribution is plausible and in accordance with findings from other works [8,71]. The thermal gains may even add up to 20% to 45% of the total heating demand, as reported in [9] for various applications in Germany. The relative contributions of thermal network gains are especially high if a GHC is the only heat source in a network. Here, the freely fluctuating network temperatures provided by a GHC are pointed out as the main advantage for harvesting a considerable amount of heat from the distribution pipes in 5GDHC networks. These thermal gains are often seen as another advantage over heating networks of the 3rd or 4th generations. AGFW, a German national association of heating network operators, states in its annual report of 2022 that on average thermal losses in German heating networks amount to 11% of the supplied energy [72].

5. Conclusions

The simulation model presented in this work comprises the possibility for extensions of the network equipment due to the usage of interfaces from Modelica standard libraries. In this way, it can be used in the conceptual and advanced planning phase of 5GDHC networks and for comparisons with more traditional heating networks of former generations. Furthermore, the individual technical equipment of single buildings with roof-mounted thermal and electrical systems, heat pump models, storage capacities, or dedicated time-series for heating and cooling demands can be modelled in detail.

The results from case study 1 show that the influence of pure free cooling of residential buildings on the overall efficiency of the decentralised heat pumps is negligible. As the feed-in of heat from free cooling operation occurs primarily in the summer time, the impact on network temperatures seems to subside until the beginning of the following heating period. This is supported by the results from case study 2, where net thermal gains from the distribution network turn into thermal losses during summer time. What is more, a reduction in the minimum necessary GHC area by pure free cooling from residential buildings cannot be suggested from the results in case study 1. As the courses of district-wide SCOP values for decentralised heat pump operation at varied GHC areas are more or less identical for the settings with and without free cooling, it can only be promoted as a favourable functionality in 5GDHC networks, but not as a heat source of significant extent. Simulation-assisted calculations are advisable when seasonally shifted heating and cooling demands are present in order to adequately reflect the seasonal loads on environmental heat sources like ground heat collectors. The provided results also show that worst-case assessments in terms of (temporary) disadvantageous climatic boundary conditions should be taken into account for more detailed design phases, especially in networks with predominant heating applications.

Decentralised, roof-mounted PVT systems can contribute a share to the total source heat fed into the network on a comparable level to the entire distribution network and approximately half the share of a horizontal ground heat collector as a primary heat source, as shown in case study 2. This is attained with a number of PVT modules on a practical scale, where 22 SFHs are equipped with three modules each. Air-to-brine heat exchangers might be equally beneficial to the efficiency of decentralised heat pumps in spring and autumn due to elevated air temperatures and simultaneous low network temperatures after the heating period. However, the economical advantages of air heat exchangers over PVT modules need to be weighed up in more comprehensive analyses. The results from both case studies presented in this paper also suggest that a temporally resolved analysis of thermal gains from the distribution network and additional decentralised heat sources can lead to significant reductions in the GHC area compared to well-established, but simplified procedures based on tabular data. With the assistance of a highly customised simulation, both the demand-tailored dimensioning of heat sources like the GHC and the management of heat source operation can lead to reduced investment and operational costs in 5GDHC networks.