Connecting Dynamic Heat Demands of Buildings with Borehole Heat Exchanger Simulations for Realistic Monitoring and Forecast

Space heating is a major contributor to the average energy consumption of private households, where the energy standard of a building is a controlling parameter for its heating energy demand. Vertical Ground Source Heat Pumps (vGSHP) present one possibility for a low-emission heating solution. In this paper, we present results of building performance simulations (BPS) coupled with vGSHP simulations for modelling the response of vGSHP-fields to varying heating power demands, i.e. different building types. Based on multi-year outdoor temperature data, our simulation results show that the cooling effect of the vGSHPs in the subsurface is about 2 K lower for retrofitted buildings. Further, a layout with one borehole heat exchanger per building can be efficiently operated over a time frame of 15 years, even if the vGSHP-field layout is parallel to regional groundwater flow in the reservoir body. Due to northward groundwater flow, thermal plumes of reduced temperatures develop at each vGSHP, showing that vGSHPs in the southern part of the model affect their northern neighbors. Considering groundwater flow in designing the layout of the vGSHP-field is conclusively important. Combining realistic estimates of the energy demand of buildings by BPS with subsurface reservoir simulations thus presents a tool for monitoring and managing the temperature field of the subsurface, affected by Borehole Heat Exchanger (BHE) installations.

1 Introduction

With around 435 TWh, space heating made up around 67 % of the average energy consumption of private households in Germany in 2018 (Bundesumweltamt, 2020). While this number is continuously decreasing, it still comprises the lion's share in private energy consumption. Gas or oil heating still constitute the primary heating systems (Bundesumweltamt, 2020), yet gradually decrease in usage and are partly substituted by direct-heat from shallow geothermal resources using borehole heat exchangers (BHE), also called vertical ground-source heat pumps (vGHSP). Heat from anthropogenic sources increasingly affects the shallow subsurface, especially in urban areas, creating so called heat islands (Ferguson and Woodbury, 2007; Taniguchi et al., 2007; Epting and Huggenberger, 2013; Mueller et al., 2018). Such elevated temperatures in these urban heat islands pose an increased potential for shallow geothermal exploitation (Epting and Huggenberger, 2013; Zhu et al., 2010), as higher groundwater temperatures increase the heat extraction rate (Tissen et al., 2019). However, the increasing spatial density of vGSHP-installations necessarily leads to the question of how strongly these installations affect each other and the overall subsurface, e.g. the groundwater. For instance, Herbert et al. (2013) conclude that regulation is needed when it comes to a sustainable management of the Chalk aquifer energy resource beneath the city of London. Similarly, García-Gil et al. (2014) show that thermal management of an exploited aquifer system is mandatory for a sustainable usage. Further, they show that factors like the presence and exploitation rate of vGSHP or river-level variations have to be considered to model the realistic state of a thermal aquifer resource. Mueller et al. (2018) as well, argue that modelling of groundwater flows calibrated by data from monitoring tools is important and valuable for managing groundwater resources and their temperatures. That is to say, the combination of groundwater modelling and modelling of realistic vGSHP operation is important for monitoring the state of a subsurface thermal resource and to forecast its long-term behaviour.

An essential input quantity in such coupled models of vGSHP operation and subsurface heat and fluid flow is the heat extraction via the GSHP system. Required heat extraction rates are determined by the energy demand of the respective buildings. In some cases, vGSHP field simulations use estimated values of the energy demand averaged over longer time frames as model input parameter (e.g. Mottaghy and Dijkshoorn, 2012). For instance, García-Gil et al. (2014) use annual extraction rates for the majority of systems implemented in their model. It should be noted that they study groundwater heat pumps, i.e. open systems using the groundwater directly. The ability to model the load of a user on the subsurface on finer timescales, be it days or even hours, is important for improving thermal management for a geothermal resource, since it increases the precision of the performance prediction and long-term system evaluation (Lamarche, 2009). This study assesses such a finely resolved coupling between user demand and response of the energy resource in a de-centralised system with consumers (houses) coupled to single geothermal vGSHPs. The energy demand itself is the result of a building performance simulation (BPS), which models the heating demand of a prescribed building in user-defined time frames. Results of the scalable BPS are thermal power curves, which can be used as input for modelling a realistic mix of energy input for heating. This approach allows for a realistic building-specific and temporally varying representation of heat extraction rates in the simulations when direct, location-specific heat demand data are not available. While using realistic BPS for vGSHP-simulations is not novel, the presented tool chain allows for a rapid set up of combined simulations of BPS and vGSHPs, focusing on all relevant properties and major drivers (e.g. building material, building shape, type of usage) without the need to define details of minor importance. It allows thus for evaluating various scenarios like retrofitting options and their impact on vGSHPs within relatively short time.

The objective of this work is to use modeled heating demand data resulting from BPS of buildings in a residential settlement in order to assess the thermal impact of a vGSHP field operated for space heating. Using numerical simulations we aim at investigating the thermal effects on the subsurface given prescribed heat loads that are defined by the heating power demand of the houses. Hereby, the characterization of the regional hydrogeological setting and thermal state of the subsurface is essential for a proper assessment of the BHE field operation (e.g. Perego et al., 2021). A practical approach is presented combining two open-source codes: the Tool for Energy Analysis and Simulation for Efficient Retrofit (TEASER) (https://github.com/RWTH-EBC/TEASER, last access: 1 October 2021) (Remmen et al., 2018) for modeling the heating power demand of buildings and the code SHEMAT-Suite (https://git.rwth-aachen.de/SHEMAT-Suite/SHEMAT-Suite-open/-/wikis/home, last access: 1 October 2021) (Keller et al., 2020) for simulating groundwater flow and heat transport in porous media. The approach allows for studying the possibility of sustainable geothermal energy use in a residential neighborhood. Specifically, this study addresses the case of equipping an existing residential neighborhood with vGSHPs, where one vGSHP is placed next to every single-family house. This limits the arrangement of the BHEs compared to a newly planned neighborhood where BHEs could be placed at favourable or even optimal locations. In addition, it investigates the effects of energetic renovation of the houses (i.e. decreasing heat demand) on the thermal state of the subsurface. This case study is relevant in context of the decarbonization of the heating and cooling sector, as, for instance, the majority of heating systems is powered by hydrocarbons in Germany (Bundesumweltamt, 2020).

The houses in the studied settlement, the “Fliegerhorstsiedlung Teveren”, share a similar layout and building material parameters, making it a suited, well-controlled study environment for our approach. In addition, the study area reveals results of more general interest for vGSHP installations in sedimentary basins facing medium to high groundwater flow rates. In the following sections, we will first briefly describe the two codes used in this work. Then we will present an overview over the study area and the numerical model, before presenting results of transient simulations over different timescales.

2 Methods and workflow

For modelling the realistic transient vGSHP behaviour, we successively perform building performance simulation with the tool TEASER and subsurface heat and fluid flow simulation with the numerical code SHEMAT-Suite.

2.1 Simulation codes

Simulating the long-term operation of vGSHPs for direct heat use in residential buildings requires knowledge of the transient thermal energy demand of the buildings. TEASER (Tool for Energy Analysis and Simulation for Efficient Retrofit) (Remmen et al., 2018) is an open-source code designed for quickly generating Building Performance Simulation (BPS) models. It focuses on Reduced Order Models (ROMs), gradually known for reducing the computational effort and simplifying the parameterization of BPS models. This comes with the costs of reduced spatial and physical resolution. The required level of model complexity thus heavily depends on the use case of the simulation. While the investigation of architectural elements forces a high resolution, using thermal energy demand of a group of buildings as input for system simulations allows for reduced models focusing on the major driving forces (e.g. building material, building shape, type of usage). A common approach for such reduced order models are thermal networks, representing heat transfer and storage via chains of thermal resistors and capacitors. The building model of TEASER is based on the German standard VDI 6007-1 and is implemented in the modeling language Modelica as part of the open source library AixLib (Müller et al., 2016) within the IEA EBC Annex 60 (Wetter and van Treeck, 2017) and IBPSA Project 1 framework. The model has been used in various projects on district and urban scale (Lauster, 2018) and has been validated on the basis of the ASHRAE 140 test cases, proving competitive results compared to common BPS tools (Lauster et al., 2017). Besides the building model itself, the parameterization of such models is often a time-consuming and error-prone process. Keeping in mind the aim of this study, only a limited effort can be spend for creating the buildings' thermal demands and the information is sparse on the actual buildings. On the other hand, the thermal demands need to reflect the typical demand of the given building types, since the study aims at simulating the long-term behavior of vGSHPs. In addition, the reduced order model focuses on the aforementioned main driving forces and properties and defining all details is of no further use. Thus, TEASER comes with statistical enrichment processes to infer required properties out of acquired base data. By using those processes, the necessary user-input is reduced to a minimum, as given in Table 1. With the help of these parameters (Usage type, year of construction, number of floors, floor height and net lease area), TEASER is able to infer the building geometry, construction properties (insulation and thermal mass), zoning and usage conditions. For this, TEASER incorporates the data of several studies and standards aiming at characterising the German building stock (DIN EN 12524; Loga et al., 2005; Bundesministerium für Verkehr, Bau und Stadtentwicklung, 2007, DIN V 18599-10). The year of construction is for example used to determine the building standard applicable in this period (e.g. Bundesministerium für Wirtschaft und Technologie and Bundesministerium für Verkehr, Bau und Stadtentwicklung, 2014). A full description of TEASER and all related data sources can be found in Remmen et al. (2018) and Lauster (2018). Of course, models based on such basic parameters can only reflect the typical behavior of a given building and major deviations are to be expected if the building does not fit into statistics.

Figure 1 (a) Location of the study area in Germany. (b) Topography of the study area, its location relative to the city of Geilenkirchen and outline of regional flow model. (c) Layout of the settlement “Fliegerhorstsiedlung Teveren”, the grey boxes represent buildings. The digital elevation model seen in (b) and (c) was obtained from CGIAR (Jarvis et al., 2008).

In addition to building-specific parameters, TEASER requires transient input data, i.e. time series of outdoor temperature, for realistically simulating the heat demand of a specified building. After incorporating input data, time series, and generating the BPS model, dynamic building simulations can be run, and provide outputs, e.g. the net-heat demand used in this study. As the layout of the buildings in the Fliegerhorstsiedlung is standardized, the output of the dynamic building simulation is representative for all buildings incorporated in the simulation.

For modelling heat and fluid flow in a porous medium, we use coupled conservation equations for energy, mass, and momentum. These conservation equations are implemented in the open-source code SHEMAT-Suite (Keller et al., 2020). The transient flow equation used in our simulations is based on Darcy's law and the mass conservation equation. Similarly, the heat transport equation is calculated using an energy conservation equation building on Fourier's law of heat conduction. The two transport equations are coupled via the specific discharge (Darcy velocity), and solved numerically in a finite difference scheme (for more details see Clauser, 2003).

Figure 2 (a) Outdoor temperature curve of a weather station in the vicinity of the study area. (b) Heat demand curves for one house in the western model in pre-renovated (blue line) and post-renovated (green line) conditions.

2.2 Implementation of vGSHPs in SHEMAT-Suite and combination with TEASER

This study focuses on investigating the thermal impact of multiple vGSHPs in a residential neighborhood on the subsurface depending on a prescribed realistic heat load. Therefore, the fluid flow within the single vGSHPs is not computed explicitly. Instead, we apply the simplification for BHE fields developed by Mottaghy and Dijkshoorn (2012), which approximates each vGSHP within a subsurface heat and fluid flow model by transferring its heating load into an effective heat generation along the vGSHP. Thus, the inlet and outlet temperatures of the individual vGSHPs are not computed explicitly, but the temperature field around each vGSHP is modeled and its interaction with the natural heat and fluid flow regime can be studied. This approach is computationally feasible for whole vGSHP fields and long-term simulations without the need for simplifications in operation time or temporal resolution of the vGSHPs. Additionally, code developments towards High-Performance Computing (HPC) and the use of modern HPC architectures enable simulating larger model areas or finer discretization in order to resolve heterogeneous or complex model geometries in more complex geological study areas. The required input for the BHE field simulation with SHEMAT-Suite are the vGSHP locations and the transient heating load. The latter is defined by the heat demand of a building modeled with TEASER.

3 Model setup

The study area is located in the Lower Rhine Embayment, near the city of Geilenkirchen, North-Rhine Westfalia, Germany (Fig. 1a). It comprises the settlement “Fliegerhorstsiedlung Teveren” (blue polygon in Fig. 1b, c), located adjacent to the NATO-Airbase in Geilenkirchen-Teveren. Buildings in this settlement were constructed in the 1950s and are planned to be retrofitted to newest energy building standards. Their locations are depicted as grey rectangles in Fig. 1c. The following subsections briefly describe building parameters and the calculated heat demand of buildings, followed by a description of the subsurface vGSHP model.

3.1 Building parameters

We investigate a submodel around a neighborhood of semi-detached houses from the 1950s that are going to be retrofitted (western red rectangle in Fig. 1c). We assume each house to be equipped with one vGSHP to deliver its full heat demand throughout the year. Since the buildings share a very similar layout, each building is assigned the same heat demand curve, considering two scenarios: (i) buildings in original condition, (ii) retrofitted buildings.

Table 1 Building parameters used as simulation input for TEASER.

This is a relevant information for modelling the transient thermal energy demand of buildings, as the average demand correlates with the building type, as well as with the year of construction or year of renovation respectively (McKenna et al., 2013; Remmen et al., 2018). Table 1 summarizes all relevant building parameters for simulating the thermal energy demand of a building using the code TEASER (Remmen et al., 2018), described in Sect. 2.1. With these building parameters, and data enrichment functions used in TEASER, we calculate the heat demand curves for the specified buildings based on measured outdoor temperatures from a station near the study area (Fig. 2a). The bottom plot in Fig. 2 shows the modelled heat demand curves for buildings in original state (blue line) and for renovated buildings (green line) in the western part of the settlement.

Figure 3 (a) Positions of borehole heat exchangers (blue circles) in the local numerical model (see Fig. 1 for location). Each borehole heat exchanger is assigned to one housing unit. The vertical black line marks the location for the model cross-section shown in (b). (b) The S-N cross-section through the local model domain shows the different lithological model units (see Table 2 for their petrophysical properties). The horizontal model slices presented in the following figures are taken at a depth of 50 m, i.e. at the center of the BHE and within the Hauptkies Aquifer.

3.2 Geological and numerical model

The Lower Rhine Embayment is the northern part of the European Cenozoic rift system (Woolderink et al., 2021), characterised by NW-SE striking normal faults. This normal fault system yielded the development of a Horst-Graben-system, which created accommodation space for siliciclastic sediments from the Tertiary onwards (Schäfer et al., 2005). These sedimentary successions represent the encountered rock formations in our vGSHP-simulations. They consist of alternating sequences of relatively unconsolidated marine and continental deposits, with the fine-grained, shale-rich marine sediments being generally less permeable than the continental deposits. Thus, in terms of the hydrogeological system, several aquifers (sandy continental deposits) are divided by aquitards (shaly marine deposits and lignite seams). A 3D geological model, focusing on the detailed representation of Tertiary sediments in the study area, is the base for a regional heat and groundwater flow model (see Fig. 1b) that was calibrated against available groundwater level data (Erftverband, 2015) and temperature data from nearby deep wells, obtained via the Fachinformationssystem Geophysik (Kühne, 2006, https://fis-geophysik.de/, last access: 1 October 2021). We refer to Niederau et al. (2017) for a detailed description of the calibration process. The general gradient in hydraulic potential is from south to north, implying a northward groundwater flow. For evaluating the interaction of vGSHPs, knowledge of the natural groundwater flow direction and velocity is essential. Subsequently, this calibrated regional numerical heat and fluid flow simulation provides the initial and boundary conditions for the local vGSHP-field model of a representative neighborhood in the settlement Neuteveren (western red rectangle in Fig. 1c). This local model is discretised in an equidistant hexahedral grid with a cell size of 2 m × 2 m × 5 m (x × y × z), covering an area of 240 m × 160 m, with a vertical extent of 200 m. Table 2 lists the main petrophysical properties of geological units in the model which are inherited from the calibrated regional heat and fluid flow model and Fig. 3b displays a representative cross-section through the local model.

In the chosen vGSHP-layout, one ground source heat pump is assumed to be installed per housing unit, yielding 26 vGSHPs in total (Fig. 3a). All vGSHP units extend to a depth of 100 m below ground surface. This is a regular depth for geothermal borehole heat exchangers, as depths greater than 100 m require more extensive reviews based on the German mining law.

As the houses have the same design and approximate net-lease area, we assign the same heating load to each vGSHP, considering two scenarios: (i) heat demand of a house in original state (built 1953, Fig. 2 blue line) and (ii) heat demand of a house renovated to new energy standards (Fig. 2 green line).

4 Results

Operation of the vGSHP-field in a residential neighborhood was simulated over the timeframe of one year in order to assess the impact of seasonal heating demand and groundwater flow on the vGSHP-field. Additionally, we present a long-term simulation over 15 years for the scenario with renovated buildings for evaluating the sustainability of the vGSHP-field.

Figure 4 Map views of temperature slices at 50 m depth below surface. (a, c, e) Temperature slices for the vGSHP field assuming renovated buildings. (b, d, f) Temperature slices assuming original building structure.

Table 2 Geological units and ascribed petrophysical properties used for the vGSHP Models. The properties are inherited from the calibrated regional heat and fluid flow model.

For assessing the influence of retrofitting buildings on the subsurface temperature, we compare simulations considering heat demand of refitted buildings with simulations considering heat demand of buildings in original state. Figure 4 shows the resulting temperature solutions at a depth of 50 m below ground surface for refitted buildings (left) and for buildings in original state (right). We chose slices at 50 m depth, as this is at the center of the BHE, but also within the Hauptkies Aquifer, the main aquifer unit in the model next to quaternary deposits. Due to the significantly higher heat demand of the non-renovated houses, temperatures in and around the vGSHPs are around 2 ^∘C lower than for the renovated scenario. Depending on aquifer permeability, groundwater flow shows a significant impact on temperatures around the vGSHPs, as plumes of relatively lower temperatures are transported northward along the local groundwater flow direction. The evolution of such hot or cold plumes is a known effect when groundwater flow is present and is investigated numerically and analytically for example by Pannike et al. (2006) and Hähnlein et al. (2010). Since the heat demand for buildings in the original scenario is higher, the thermal plumes are more pronounced. Consequently, the interaction between single heat exchangers is stronger, as vGSHPs in the northern part of the model are affected by thermal plumes deriving from neighbouring installations. For determining the horizontal extend of the cooling plumes, we look at the temperature difference between simulated temperatures in the model cells and initial undisturbed ground temperatures. We define the maximal extend of a cooling plume at the point where this difference is zero. The cooling plumes reach up to 50 m north of the respective vGSHP, which does not favor a north-south alignment of vGSHPs (e.g. southern pairs of vGSHPs in Fig. 4).

Figure 5 (a) Time series for one year of two vGSHP in the southern part of the model. The influence of the southern vGSHP (blue) is expressed as a reduced temperature in the northern vGSHP (red). The difference between the two vGSHP is plotted and reaches its maximum around end of June. (b) Time series of temperature at 50 m depth over the same time period at a location (triangle) affected by the thermal plumes of the two vGSHPs marked in A. The time series are plotted for the original (red) and renovated (blue) scenario. A reference temperature at the same depth at a location (marked by the circle) unaffected by thermal plumes is plotted in gray. (c) North-south trending temperature profile at day 180 at a depth of 50 m below surface. Pronounced peaks are temperatures within the vGSHPs. Between the peaks, temperatures are lower than the undisturbed temperature of 11 ^∘C. The black line in the temperature field on the right shows the position of the temperature profile.

Figure 6 North South trending temperature profile for each year over a total simulation time of 15 years. Small subplot shows the location of the temperature profile.

For assessing the effect of a thermal plume on the neighbouring vGSHP, we study the temperature development over time at two vGSHPs in the southern part of the model (Fig. 5a). The simulation reveals that temperature differences between the two neighboring vGSHPs peak around day 180, i.e. at the end of June. Figure 5b presents the transient temperature development at one point (see the triangle in the cut-out map) north of these two vGSHPs for the original (red dashed line) and for the renovated scenario (blue dashed-dotted line), respectively. Other than the renovated scenario, the original scenario shows a plateau at around 10.1 ^∘C in March/April 2014 after steep temperature decrease. This is due to a constant temperature of around 2.3 ^∘C at the locations of the BHE installations. Figure 5c shows a north-south temperature profile at 50 m below ground level at day 180. Undisturbed simulated temperatures at that depth are around 10.9 ^∘C. The pronounced temperature peaks mark fluid temperatures directly in a vGSHP. The temperature signal of the created thermal plume extends from the most southern vGSHP until around 80 m along the profile. The temperature difference between the two vGSHPs, however, is relatively small with around 0.3 ^∘C. In contrast, the scenario with buildings in original conditions reveals significantly larger temperature differences. Here, the plume of the most southern vGSHP causes a temperature drop of about 1 ^∘C. This in turn affects the thermal energy rate that can be provided by a vGSHP. A slight temperature decrease at around 100 m in both scenarios reflects the thermal signal of a vGSHP adjacent to the profile (Fig. 5c). For evaluating the effective impact of the different building scenarios on a vGSHP over one year, we also calculated the root mean square difference (RMSD) of the temperature difference compared to an undisturbed subsurface at a vGSHP location for both scenarios. Averaged over the whole year, the influence of the renovated scenario on the temperature field in the subsurface of 2.31 ^∘C is about 40 % of the influence of the original scenario (5.78 ^∘C).

Seasonal variations in heat demand are reflected by northward migrating thermal plumes. When the heat demand decreases or is equal to zero (as in July/August in the presented simulations), temperatures in the vGSHPs almost re-equilibrate with the surrounding temperature for the scenario with renovated buildings. Contrary, the scenario with buildings in original state still preserves cooling plumes in August (Fig. 4d). Thermal plumes are only observed in more permeable units, namely the Hauptkies aquifer and the Quaternary aquifer. Maximum simulated groundwater velocities in these aquifers are in the range of 70 m a⁻¹. In deeper sedimentary units, such as the Inden and Ville Layers, no thermal plumes are observed and the heat transport towards vGSHP installations is dominated by conduction.

4.1 Model simulation over 15 years

Simulations over a timeframe of 1 year do not suffice for estimating the sustainability of a vGSHP-field. Mottaghy and Dijkshoorn (2012) simulate a vGSHP-field for 15 years of operation in order to analyse sustainable operation of the simulated vGSHP-field. Further, using geothermal heat supply just for space heating will cool the subsurface over a longer timeframe. To evaluate the sustainability of the vGSHP-field, we extended the simulation time to 15 years for the western model, assuming refitted buildings. For assessing the magnitude of cooling in more detail, we plot, similarly to Fig. 5, the temperature profile in June for each year.

Temperatures along the profile significantly decrease over the first three years, due to the development of thermal plumes, i.e. the effect of advective heat transport. However, temperatures stabilize at around 10 ^∘C for the rest of the simulation time. Accordingly, impact of thermal plumes also stabilizes. The possible influence of the prescibed heat load on this observation is discussed in Sect. 5.

5 Discussion

The presented results underline the importance of considering regional groundwater flow when modeling the transient behaviour of vGSHPs, as also concluded by other authors (e.g. García-Gil et al., 2014; Mueller et al., 2018). Knowledge of magnitude and direction of the regional groundwater flow facilitates an optimised layout of multiple vGSHP-units. For instance, by offsetting the southern vGSHPs in an EW direction, i.e. perpendicular to the groundwater flow direction, the performance of vGSHPs in the southern part of the western model improves slightly, as thermal plumes do not directly cross neighbouring installations. With the spatial distribution of vGSHPs likely becoming denser in the future, such information is important when planning new installations. We exemplarily presented the simulation of one vGSHP field layout with one 100 m deep vGSHP per housing unit and equally distributed heat loads. In subsequent studies, a systematic optimization approach could be applied for optimizing the BHE field operation in order to reduce the thermal impact on the subsurface and on neighboring BHEs. For instance Hecht-Méndez et al. (2013) present such a simulation-optimization procedure that aims at regulating the BHE operation in a field such that the thermal impact on the ground is minimized. They specifically investigate the influence of groundwater flow on the optimal operation strategy. Their results suggest that for high groundwater velocity as present in the Hauptkies aquifer it is beneficial to first operate the downgradient BHEs and then successively distribute the load on the upgradient BHEs. However, it might be more complex for our study area since the BHEs penetrate several layers with different groundwater flow velocities (i.e., aquifers and aquitards; see Fig. 3b). Furthermore, Beck et al. (2013) present a methodology for optimizing the arrangement of the vGSHP installations in a field. The combination of such optimization strategies with realistic hydrogeological models as well as building performance simulations could contribute to planning vGSHP-arrays, optimised for cost and energy output, while at the same time reducing the thermal impact on the ground.

It should be noted, that the temperature curves underlying the heating power demand for the vGSHP-field do not comprise extreme events, such as long cold spells during winter. Instead, we assume a more regular, sinusoidal seasonal temperature pattern. Accordingly for the 15 year simulations, temperatures stabilize at certain values, due to the periodic heating power demands, which enter the numerical model as effective heat generation along the vGSHPs. A stabilization of temperatures under the assumed conditions suggests, that the simulated vGSHP-field could be operated sustainably. However, with dry periods due to climate change becoming longer and more frequent, their long term effect on groundwater flow, or aquifer recharge could be considered. The lateral boundary conditions of the numerical model are constant throughout the simulation duration. That is, water with the same temperature enters the model at the same velocity, and may contribute to the stabilization of temperatures at the vGSHP positions. This may not be realistic, considering, for instance, the influence of climate change on the groundwater table. This in turn may yield changes in the hydraulic potential, and consequently, changes in direction and velocity of regional groundwater flow. Using historical data to model such complex, transient boundary conditions goes beyond the scope of this work. Future simulations, however, may include systematic sensitivity analysis of boundary conditions.

In our simulations, we only assumed heating power, not cooling during summer. Incorporating cooling would likely improve the vGSHP-field performance. A combination of heating and cooling, where surplus heat gets stored in the subsurface during summer to be available in winter, can optimize the operation of a vGSHP-field (e.g. Bayer et al., 2014). Again, an increasing influence of climate change, i.e. hotter summers, milder winters, can factor into such an optimized operation. However, the increasing energy efficiency of buildings due to better insulation, newly built or retrofitted, can drastically reduce the heating or cooling demand of buildings, and may make up for the increasing change of outdoor temperature.

The presented approach is restricted to evaluating the thermal influence of the BHE field operation on the subsurface, i.e. approximating the vGSHPs in the model and not calculating the inlet and outlet temperatures of the BHE installations explicitly. Therefore, no direct calculation of the performance or COP of the vGSHPs is possible. Investigating the influence of temperature changes on the operation and performance of single vGSHPs in the field would require the explicit simulation of single BHEs, which could be done in an extension of this study.

6 Conclusions

Based on results of a calibrated regional flow model around a settlement in the Lower Rhine Embayment (western Germany), we generated two representative models for simulating the operation of borehole heat exchangers using realistic heat demand curves from building performance simulations. Simulation results show the pronounced development of thermal anomalies at BHE-locations. Groundwater flow in aquifer layers advectively transports the thermal anomalies northward, resulting in cooling plumes that affect other BHEs in that direction. Simulations over a 15 years time period show that temperature change stabilises around 1 ^∘C below initial conditions at a depth of 50 m below surface. The incorporation of building performance simulations in our workflow allow for direct comparison of the temperature evolution across the BHE-field for old houses from the 1950s with renovated houses. Higher heat demand of old buildings directly impacts the temperature decrease in the subsurface during phases of high heat demand. Less heat demand of renovated houses favors the regression of thermal plumes during the summer months.

The simulation results support the conclusion that thermal management of the subsurface is necessary for its sustainable usage. Our work shows that coupling realistic heat demands of buildings on fine timescales with subsurface heat transport simulations provides a tool suited for such management. The presented workflow can help addressing economical issues as well as regulatory requirements considered during the planning phase of vGSHP installations, such as:

• Do the planned vGSHPs influence each other or neighbouring fields?

• What is the expected temperature change in the subsurface, and does it comply with existing regulations?

• Is a sustainable long-term operation possible, given specified building parameters?

Thus, the approach contributes to the optimal planning of vGSHP systems. Further, the considerable difference of heat demand between renovated older buildings and ones in their original state underlines the importance of up-to-date energy concepts for building, as this can significantly facilitate a sustainable management of the heat resource in the shallow subsurface.