Modeling neighborhood-scale shallow geothermal energy utilization: a case study in Berlin

Nowadays, utilizing shallow geothermal energy for heating and cooling buildings has received increased interest in the building sector. Among different technologies, large borehole heat exchanger arrays are widely employed to supply heat to various types of buildings. Recently, a 16-borehole array was constructed to extract shallow geothermal energy to provide heat to a newly-developed public building in Berlin. To guarantee the quality of the numerical model and reveal its sensitivity to different subsurface conditions, model simulations were conducted for 25 years with two finite element simulators, namely the open-source code OpenGeoSys and the widely applied commercial software FEFLOW. Given proper numerical settings, the simulation results from OpenGeoSys and FEFLOW are in good agreement. However, further analysis reveals differences with respect to borehole inflow temperature calculation implemented in the two software. It is found that FEFLOW intrinsically uses the outflow temperature from the previous time step to determine the current inflow temperature, which makes it capable of much faster simulation by avoiding iterations within a single time step. In comparison, OpenGeoSys always updates the inflow and outflow temperature based on their current time step values. Because the updates are performed after each iteration, it delivers more accurate results with the expense of longer simulation time. Based on this case study, OpenGeoSys is a valid alternative to FEFLOW for modeling ground source heat pump systems. For modellers working in this field, it is thus recommended to adopt small enough time step size, so that potential numerical error can be avoided.


Introduction
The current energy transition in the world demands increasing amount of renewable and sustainable energy sources. By 2020, almost half of the electricity in Germany is already generated with renewable sources. In comparison renewable heat source only accounts for 15.2% in the building heating and cooling sector (UBA 2021). When coupled with Borehole Heat Exchangers (BHE), Ground Source Heat Pump (GSHP) systems extract geothermal energy from the shallow subsurface by circulating fluid in the BHEs, and transfer the heat into the building through a heat pump. It is thus one of the low-carbon and emission-free technologies that satisfies the energy transition target of the building sector. There are different types of BHEs which vary in different pipe configurations. Besides coaxial pipes, commonly 2U pipes are used in urban building projects, where two U-pipes are installed in a single borehole. Compared to single U-tube BHEs they offer advantages such as smaller thermal resistances and the borehole can still function in case one pipe is leaking fluid (Zeng et al. 2003;Casasso and Sethi 2014). During summer, when no thermal energy is extracted, the subsurface can be recharged by different sources, which include heat conduction, groundwater flow and, heat flux from the surface due to elevated air temperature, as well as vertical heat flux from deep earth. Additionally, BHEs can store the surplus thermal energy from building cooling into the subsurface, which can increase the system performance . Ahmadi et al. (2017) revealed in their study that using GSHP systems for cooling purposes can be more efficient than traditional air conditioning systems.
There are varying standards regulating the subsurface temperature change. For example in Berlin it is not allowed to exceed 3 K after 25 years. In the State of Hesse, however, the induced temperature drop is not supposed to be more than 0.1 K at the property boundaries after 50 years. Therefore it is often legally required, that the designer of the GSHP system has to conduct proper numerical simulations to determine the temperature evolution in the subsurface, caused by the extraction and/or injection of thermal energy over the long term. For this purpose typically the commercial software FEFLOW is utilized. Besides FEFLOW, OpenGeoSys (OGS) is a finite element numerical and open-source software for various thermo-hydro-mechanical-chemical (THMC) processes in porous and fractured media (cf. OGS 2021b). One of the applications of OGS is to simulate heat transport process caused by the operation of BHEs. As a result the associated evolution of temperature distribution in the subsurface and inside the boreholes can be predicted and analyzed. Both FEFLOW and OGS are based on the dual continuum approach, posed by Al-Khoury et al. (2010) and extended by Diersch et al. (2011aDiersch et al. ( , 2011b and Diersch (2013). They provide a fully-transient computation method for longterm simulations of the temperature profile in the subsurface. The methods are based on validated capacitor-resistor models which provide a significantly shorter computation time compared to fully discretized models (cf. . The simplified thermal resistance and capacity models are also used by Morchio and Fossa (2020) to validate and predict the performance of U-pipe BHEs. FEFLOW can also apply a quasistationary method proposed by Eskilson and Claesson (1988) alternatively. Besides there are several analytical computation methods available for an approximate estimation of the GSHP system performance. They usually consider the BHE as a line heat source as used by Lamarche and Beauchamp (2007) and Zeng et al. (2003) for example. Furthermore, Beier (2014) provides an analytical solution for simple 1U-BHE temperature developments by approximating the U-pipe as a single one separated in two halfs. It uses the thermal resistance R b which is usually obtained from practical thermal response tests (TRT). Distributed TRTs can be conducted to extract detailed and depth-dependent information from the subsurface (cf. Acuña et al. 2009). Claesson and Bennet (1987) introduced the multipole method to calculate this value without field results which was extended by Hellström (1991). This method is offered by FEFLOW to calculate the thermal resistance inside the BHE. In this case the equations provided by Diersch (2013) would be simplified and calculated with the thermal resistance R b . Chen et al. (2020) found a correlation between the temperature difference of the maximum or the mean temperature respectively, referring to the initial state and the total heat injected into the subsurface. Based on this a linear relationship could be found between the working fluid temperature change and the amount of the extracted heat from the subsurface.
Considering the total thermal balance of the subsurface, groundwater has a measurable effect on the extractable heat and the performance of BHEs. There have been several investigations about the heat advection caused by groundwater on the temperature distribution of the subsurface and around BHEs (cf. Taniguchi et al. 2003;Molina Giraldo et al. 2011;Rivera et al. 2015). Choi et al. (2013) found that even small velocities and the direction of groundwater flow can have noticeable influences on the system's performance, wherefore this boundary condition must considered properly in advance. Meng et al. (2019) stated that with an increasing groundwater velocity, the controlling heat transport process in the subsurface changes from thermal conduction to convection. With high Darcy velocity, low temperature accumulation in the direct surrounding of BHEs becomes less strong. Groundwater also affects the results of thermal response tests by increasing the measured value of the thermal conductivity of the subsurface (cf. Witte 2001; VDI 2020b). Zhu et al. (2010) found that a temperature difference of 2K of a 20 m thick aquifer yields enough geothermal energy to satisfy the annual space heating demand of Cologne 2.5 times. Therefore the influence of groundwater can not be ignored when performing long-term simulations to predict the subsurface temperature.
In modern GSHP projects, BHE array is often utilized to increase the amount of energy exploited from the subsurface. Arola and Korkka-Niemi (2014) found that in urban areas 50-60% more heat can be extracted due to the effect of heat islands. You et al. (2018) introduced an analytical solution for BHE arrays which also considers the groundwater flow. Due to elevated heat extraction rate by the BHE arrays in comparison to thermal recharge, low temperature zones may occur in the subsurface. The German VDI 4640 guideline suggests a minimum distance of 6 m between BHEs in an array, so that the system's performance can be ensured (VDI 2020a). Li et al. (2012) came up with an arrangement strategy considering the influence of groundwater and stated that the distance between BHEs should be adjusted depending on the flow velocity. Moreover, they concluded that the optimal setup would be a single line of BHEs perpendicular to the groundwater flow or several shifted lines.
Considering the sustainable character of GSHP systems, a 16-BHE array has been planned to provide heating and cooling in a public building project in Berlin. As typically required by the regulations, a numerical model has been set up and long-term simulations are conducted with the software FEFLOW. With generous access of design data provided by the engineering office, the FEFLOW model has been transferred to OGS, preserving nearly the same geometric and numerical settings. This allows a detailed investigation on the comparability of the two software, and whether OGS can be an alternative for modeling and simulating heat transport processes related to shallow geothermal energy utilization. To our knowledge, it is the first time that such comparison is performed with a background of using GSHP system to extract geothermal energy.
For this purpose, the general methodology of the two software will first be reviewed in the subsequent chapter. Hereafter, the Berlin model will be configured to enable the software comparison, with a focus on diverging boundary conditions. The generated results will be summarized and investigated. This will lead to a discussion not only on the quality of the simulation results, but also the boundary conditions of the numerical simulation versus the required simulation. Based on this, advice will be given on how to use numerical simulations for a long-term prediction of the temperature development in the subsurface. The discussion also covers legal issues which may occur in different states of Germany.

Method
In this work, the open-source software OGS (cf. Kolditz et al. 2012;Shao et al. 2016), as well as the commercial software FEFLOW, are utilized to generate simulation results for a 16-BHE GSHP system project located in Berlin. Basically, both modeling software operate similarly and are based on the same assumptions and numerical approaches.

Dual-continuum finite element approach
Al- Khoury et al. (2010) developed the dual-continuum finite element approach that can be applied to simulate the subsurface heat transfer process induced by BHE operation. It was later-on further developed by Diersch et al. (2011aDiersch et al. ( , 2011b, Diersch (2013) and implemented in the FEFLOW and OGS software. Considering both the thermal convection and conduction in the soil, the governing equation for the heat transport in the subsurface reads as with as the tensor of thermal hydrodynamic dispersion. Moreover, the porosity ε , density ρ and heat capacity c are presented with the indices s for soil and f for fluid which is the groundwater in this case.
The commonly used 2U type of heat exchangers are divided into four areas (cf. Fig. 1). Within each area, it contains either an inlet (i1, i2) or an outlet pipe (o1, o2), surrounded by the corresponding grout zone (g1, g2, g3, g4). This leads to eight governing equations regulating the heat transport process on either the pipe or the grout. These equations can be summarized as follows for the domain and its boundary Ŵ: with the respective heat fluxes between pipe and grout where T k and T g denote the temperatures in the pipe and the surrounding grout zone. Interested readers are referred to Diersch et al. (2011a) for more details of the numerical methods applied to solve them.

Imposed thermal load
In common subsurface heat transport simulations, a monthly varying thermal load Q heating BHE (t) is imposed on the ground loop and acts as the driving force in the numerical model. The simulation then computes the pipe and soil temperature distribution that is needed to generate the respective energy to supply this load. Therefore, the amount of extracted heat at the entry of each BHE at a given time t must satisfy where the flow rate V r is equal to the sum of flow rates in two separate loops of the 2U-type BHE. Since the temperature difference T in Eq. (4) is the only parameter that has a degree of freedom, re-arranging the equation leads to how the T value can be obtained by Assuming an operation during January (cf. Table 1) with 31 days and using the fluid properties listed in Table 2 the resulting temperature difference amounts to T = 0.950 K exemplary. Since the thermal load is externally imposed on the numerical model, one can use the inlet and outlet temperatures from the BHEs to calculate the T value and check how accurate the model simulation is. Besides, there are different types of conditions for the flow and temperature control and both OGS and FEFLOW are capable of operating with various impositions (cf. OGS 2021a; FEFLOW 2021). As written before a month-dependent power demand is used in the model.

Iterative and non-iterative approaches
As part of the spacial and temporal discretization, the coupled governing equations [Eqs.
(1) to (3)] are linearized and transformed to the following linear algebra form The primary unknowns in Eq. (6) are the pipe T π and soil temperatures T s . In order to make Eq. (6) solvable, at least one temperature value has to be imposed as a Dirichlettype boundary condition. In both FEFLOW and OGS software, the inlet temperature T in is chosen to be the boundary value, which can be obtained by using Eq. (5) with the given thermal load Q heating BHE and other parameters. However, for which T out to be used for the boundary value calculation, FEFLOW and OGS have chosen different strategies. In FEFLOW, the outflow temperature from the previous time step T n−1 out is used for the calculation of boundary value T in at the current time step. The advantage of this choice is that, once the last time step T out value is known, the entries in the left and right hand side matrices, together with the boundary value T in , will not change within a single time step. This means that Eq. (6) can be solved directly without any further iteration. The requirement of taking this choice is that the time step size must be controlled to a small value, so that the change in T out is not too large to cause an inaccurate simulation result (cf. FEFLOW 2021).
In the contrary, the OGS software always choose the T out value from the current time step. The consequence is that after the solution of Eq. (6), the current time step value T out will be updated. A corresponding boundary value T in must be updated accordingly. With this choice, several Picard iterations are required and the solution of Eq. (6) must be repeated. So the accuracy of the model can be guaranteed with the expense of increasing computational resources. It will be shown in this paper, how the above choices will influence the modeling results, with regard to both the achievable precision as well as the simulation time.

Numerical models
In order to analyze the numerical simulation as precise as possible, a comprehensive model has been set up. Several previous investigations were performed to achieve characteristic parameters for the BHEs and the geological subsurface. All conditions and parameters are applied equally in both software, so that computational differences between FEFLOW and OGS can be analyzed. The GSHP system is located in Berlin and contains a total of 16 BHEs. Each of them is imposed with the same thermal load, so that they will extract the same amount of thermal energy. For the simulations it is assumed that there is no energy loss in the pipe network. Furthermore, the heat pump COP is already taken into account as shown in Table 1. Thanks to the detailed geological information based on the previous investigations, numerous parameters could be attached explicitly and the numerical model could be set up comparatively realistically.

Model validation
For model validation, OGS has already been compared against both analytical solution and also experimental data. Hein et al. (2016) proved that OGS obtains comparable results as the analytical solution for the Beier Sandbox benchmark (cf. Beier 2014). For deep coaxial BHEs, Chen et al. (2019) provided a verification of the numerical results of OGS against analytical calculations. In Chen et al. (2020) an array containing 56 BHEs which have been numerically simulated with OGS has been validated against real monitoring data. As a widely used commercial software, FEFLOW has been accepted in the industry as the standard tool for long-term simulation of GSHP system operation. It can therefore be assured that both numerical modeling software are already well validated.

Preliminary geological investigation
Several investigations were performed in advance to evaluate model-specific parameters. A thermal response test has been performed to obtain site-specific parameters including the thermal conductivity of the subsurface, which is a main factor influencing the performance of the system. This parameter can be calculated in the post process of the TRT by utilizing a sequential forward evaluation (cf. VDI 2020b). According to the analysis a mean thermal conductivity of 2.5 Wm −1 K −1 which is in good agreement with the depth-dependent values from Table 3. In a further analysis of the TRT the undisturbed underground temperature have been measured. As a result a mean value of 11.6 • C could be determined, which has been applied as the initial condition in the numerical model. Besides, a borehole resistance of 0.09 mKW −1 has been ascertained.
In FEFLOW this values could be used to determine the heat transfer coefficients alternatively. In OGS, however, this value plays an subordinate role and is not used for the numerical computation. Considering the results from the TRT as well as the preliminary BHE setup, the array has been designed using the Earth Energy Designer (EED), which is an analytical software tool. With an explicitly defined thermal load over the year, the EED was able to calculate a required total of 16 BHEs with a length of 99 m each. 1 Based on degree-day measurements, the month-dependent thermal load distribution could be imposed as shown in Table 1.

Model domain and mesh
The numerical model in FEFLOW and OGS has a dimension of approximately 1100 × 800 × 175 m 3 , with its two sides aligned with the groundwater flow direction. The finite element mesh consists of 1,605,660 nodes and 3,139,837 elements. Within these elements, 624 1D line elements are used to represent the 16 BHEs, and the remaining 3,139,213 prism elements are created to discretize the subsurface soil compartment. The entire mesh is originally created in FEFLOW and then imported into OGS, so that differences from the spacial discretization can be eliminated. Both the Péclet number as well as the Courant number has been considered when creating the mesh, so that a stable numerical simulation is guaranteed. The Randow et al. Geothermal Energy (2022)   numerical mesh is constructed with a bottom-up approach, meaning that the lowest bottom layer is firstly placed parallel to the x-y-plane, but subsequently non-horizontal layers has been added on top repeatedly, based on borehole drilling logs obtained from the site. In the FEFLOW simulation, the global absolute Universal Transverse Mercator (UTM) coordinate system is used in the model. In OGS, however, the coordinates are attached to a local system. Therefore, when importing the mesh into OGS, a coordinate transformation is necessary to obtain the coordinate system as shown in Fig. 2b.
The BHE locations are restricted to the utilisable space at the building site. Therefore, considering the limited space and design recommendations from literature (cf. Li et al. 2012), the design of the BHE array is roughly line-shaped and they are placed perpendicular to the groundwater flow direction. Whenever possible, a predetermined distance have been held constantly between the boreholes. However, due to the limited size of the building site, five BHEs had to be placed besides this regular line -three in the upper part and two at the lower end.
The numerical simulation is configured to start at the beginning of a heating season (1st September in both models). The original FEFLOW model used a time stepping with maximum time steps of 5 days which provides acceptably precise results. In further investigations, however, the time stepping has been modified in different simulations in order to analyze its influence on the performance of the utilized software.

Model parameters
Analogously to the mesh, the same model parameters have been applied in both FEFLOW and OGS software. For example, the total energy demand for the building Q Building amounts to 126 MWh annually. As the heat pump is working to shift geothermal energy into the building (cf. Hein et al. 2016), the amount of heat extracted from the subsurface can be quantified with the combination of building thermal load and the heat pump COP.
Assuming COP heating = 4.6 , Q heating BHE sums up to 98.6 MWh which is equally distributed to every BHE with 6.1625 MWh each. Moreover, the thermal load fluctuates during a year, largely influenced by the atmospheric temperatures in the Berlin region. In winter, more heat is typically needed. In summer, the BHE array is switched off, allowing the subsurface to thermally recharge. Based on the pre-investigation and calculation by EED software, this monthly-dependent distribution of thermal load can be found in Table 1. The values are imposed each year equally, so that in a specific month the corresponding amount of thermal load (in kW) is applied on each BHE. Table 2 shows a list of BHE related parameters applied in the numerical model. It summarizes known geological values and results from different preinvestigations. Each of the 16 BHEs are set up equally and the parameters pertain to each of them. The comprehensive model consists of 58 non-horizontal layers in the subsurface with six different hydraulic property groups. 2 Thanks to a high-resolution geological report in the Randow et al. Geothermal Energy (2022) 10:1 Berlin area more detailed and depth-specific values could be achieved and applied in the numerical model. Table 3 shows the layer-dependent parameters of BHE #2 for example, which is located in the north-west part of the array (cf. Fig. 2b). As a result of the non-horizontal layers, the BHEs are not situated at the same vertical height, but differ corresponding to their respective location in the model. In the case of BHE #2, its top z-coordinate is z = 55.28 m . Therefore and due to the length of 99 m the BHE is only partially located in the upper part of the fifth material group and does not reach the sixth layer set at all, which is the same setup for all BHEs in the model (cf. Fig. 3).

Initial and boundary conditions
At the beginning of the simulation, the subsurface temperature in the whole domain as well as in all BHEs is set to 11.6 °C which is the measured undisturbed underground temperature. This same value is also imposed on the upstream side of the model and serves as a constant Dirichlet-type boundary condition for the inflow groundwater temperature. Considering that the atmospheric temperature is naturally fluctuating, three different configurations have been simulated to analyze their impact on the results. It is either (i) a no-flux boundary condition at the surface of the domain, or (ii) a constant ground surface temperature value set to 11.6 • C, or (iii) a fluctuating temperature curve over the year. In case (iii) the values are defined according to the equation following Bechtel (2015) In the above equation, the annual mean temperature T mean = 11.6 • C and an amplitude A gs = 10 • C were applied. Additionally, the time shift ϕ adjusts the function so that the peak is during summer and lowest temperature occurs in winter. In the original simulation no ground surface boundary condition has been applied with the intention to decrease the computation complexity of the model. Its influence will be discussed in the results part. At the bottom surface of the model, a Neumann-type boundary condition is imposed to represent the geothermal heat flux from the earth's core.

Fig. 3 Hydraulic group layers and groundwater flow in the subsurface
It is calculated by the product of subsurface thermal conductivity and the local geothermal gradient. According to the thermal conductivity of layer group six (cf. Table 3) a bottom heat flux value of q geo = 0.042 Wm −2 has been imposed as a Neumann-type boundary condition. It is caused by the combined effect from natural radiogenic heat of the crust and heat flow from the mantle. Additionally, BHEs are implemented as a fourth order specific inner boundary condition in FEFLOW.

Groundwater implementation
Both OGS and FEFLOW are able to include groundwater flow in their simulations. Based on the hydrogeological survey of the Berlin site, it is known that a natural  Figure 4 states that only in the upper layers an appropriate hydraulic head difference occurs. From z ≈ −10 m and deeper, h and thus the groundwater flow becomes minor. The major groundwater flow effect, however, depends on the permeability k f values of the different underground rock materials. In general, there are two common methods to determine the k f values; either in lab experiments for a specific material or in tracer tests particular for the examined area. Only in hydraulic group two there is a high k f value combined with a hydraulic head difference of h = 4 m . Thus the heat transport due to convection in this layer is clearly larger than it in the other material groups. Although the k f value differs a lot from each other in the lower layers, the hydraulic gradients here are comparatively small. Therefore, the flow velocities are almost negligible in each of them. The resulting distribution of the Darcy velocity u is supplied to the heat transport model as a set of input parameters. As several investigation found, the hydrodynamic dispersion in the aquifer layers has a great influence on the temperature profile. For example, Gillbricht and Radmann (2017) mentioned that the dispersion is a unity of different hydro-geological parameters to describe the distribution of particles in the subsurface besides the influence of the groundwater flow which can not be determined by typical physical laws. There are several attempts to define this parameter (cf. Pearson and Hanshaw 1970;Beims 1983;Gelhar et al. 1985), but there is no unequivocal model for a valid definition of the dispersion in groundwater filled media. Zech et al. (2015) reevaluated data values of the dispersion gathered by several calculation attempts and stated that there is no unique scaling method. In general, higher dispersivity values in longitudinal and transversal direction increase the distribution of substances and also heat in the geological subsurface. FEFLOW uses α L = 5 m and α T = 0.5 m for the two dispersivity values by default. For the OGS simulation two scenarios have been set up. The first one have the same dispersivity values as in FEFLOW and the second one set α L = α T = 0 m. Fig. 5a shows the long-term development of the mean fluid temperature from all 16 BHEs simulated by OGS and FEFLOW. Here, a five-day time stepping OGS simulation is contrasted to a one-day time stepping FEFLOW setup. Further analyses showed that these two scenarios yield the best results for the respective software (cf .  Table 4 and Figs. 7, 8). In the first 5 years, a relatively rapid decrease of the fluid temperatures can be observed, which attributes to the imbalance between heat extraction and thermal recharge. Due to the heat exploitation, a quasi-steady state with a lower temperature range emerges in the later years which shows an equilibrium between the heat extraction during winter and recovery in the following summer. In Fig. 5b a more detailed profile is illustrated, indicating the close resemblance of both simulation results. It can be found, that the FEFLOW curves are deviating from the OGS results, since the two software take different T out values for the calculation. Nevertheless, the difference between the simulations is negligible minor and it can be stated that both software deliver equivalent results based on the respective numerical setups. Further analyses will investigate the origin of these differences.

Fluid temperature development
When looking at the big picture, a characteristic low temperature plume is developed in the subsurface after 25 years of operation. As shown in Fig. 6, the greatest temperature difference between the initial and the current state can be observed next to the BHEs. Further away from the array, the temperature difference decreases until the subsurface is no longer affected by the GSHP system operation. Additionally, the shape of the thermal plume is noticeably influenced by the groundwater flow in different layers. In the upper layers, where a high Darcy velocity is present due to corresponding hydraulic head differences and high k f values, the extend of the thermal plume reaches 324 m with a isotherm temperature of 11.2 • C . That also means that a convective heat transport occurs in the respective layers wherefore the subsurface surrounding the BHEs are continuously recharged, thus increasing the system's performance. The isotherm of 11.0 • C is comparatively smaller and does not exceed the property boundary in the downstream direction at all. In some areas in Germany this contour is a crucial factor for the installation of GSHP systems since the temperature difference in the subsurface at the property borders is not allowed to exceed a defined value after a specific operation time. In this case it could be proved by a numerical simulation, that the maximum difference is smaller than 0.3 K at the property boundaries (cf. Fig. 6). A more detailed discussion will be lead in the legal issue section. The minimum subsurface temperature after 25 years of operation amounts to 9.17 • C and can be found in the lower part right next to a BHE where the advection does not have a large influence and the heat transport mechanism is mainly driven by conduction. While above this depth the groundwater flow contributes to the heat dispersion, closer to the lower end of the BHE it can also extract from the subsurface below which benefits the temperature development positively.

Vertical temperature distribution
In addition to original FEFLOW and OGS scenario, an alternative OGS scenario has been set up. It recreates the computation method from FEFLOW using the previous outflow temperature to find the according fixed inflow temperature T in value in Eq. 6. Nevertheless, it still conducts multiple Picard iterations in each time step of simulation. Figure 7 illustrates vertical profiles of different simulation setups in the first January of operation. As a consequence the temperature difference between inlet and outlet of the BHE can be obtained from the graphs from the two temperature values at the top. As predetermined by Eq. (5) this gap should theoretically equal 0.950 K in January, based on the thermal load value from Table 1 accordingly. This target value thus indicates how precise the generated results of a simulation are. Due to the slightly different computational approaches implemented in the two software, this temperature difference also serves as a measure on how far the simulation results have deviated from the reality. Several simulations have thus been performed to analyze the behavior with different time stepping configurations.
In Fig. 7a all of the three scenarios employ a maximum time step size of five days. Since an automatic time stepping was applied in FEFLOW which also have been adopted in OGS consequently, the respective maximum time steps have not been used permanently, but for the majority of the simulation period. Especially at the beginning of months the interval between two discretized calculation points is decreased so that   additional advective heat transport. Comparing the temperature difference of the fluid in each simulation case, FEFLOW delivers a value of T FEFLOW = 0.883 K which is a relative deviation of 7.05 % to the target value 0.950 K . On the other hand, with an equivalent time stepping and utilizing the outflow temperature from the previous time step, OGS model 2 achieves a gap of T OGS,2 = 0.944 K which amounts to a relative difference of 0.63 % . The original OGS scenario (model 1) results in T OGS,1 = 0.929 K which deviates 2.21 % to the target temperature difference in January. As mentioned above, the simulations made with OGS are very near to the target value compared to FEFLOW. One option to increase the output precision is to shorten the utilized time stepping, especially the accuracy of FEFLOW results may benefit greatly from this adjustment. In Fig. 7b the time stepping was decreased to a maximum of one day, and simulations are repeated for both software. This leads to lower fluid temperatures which indicates that a greater time stepping may overestimate the values. All the simulated temperature profiles now converge closer with each other with a maximum difference of approximately 0.04 K , which is negligibly for long-term simulations. With smaller time step size, all scenarios could increase their precision compared to the target temperature difference of 0.950 K with a maximum relative deviation of 1.05 % achieved by FEFLOW. Additionally, a 15-days time stepping scenario has been simulated. All setup results are summarized in Table 4. In general for these scenarios OGS can deliver more accurate results compared to FEFLOW. The above comparison indicates that FEFLOW is more sensitive to increasing time step size, as the relative temperature difference deviations enhance more rapidly than OGS. As usual, a comparatively shorter time stepping can achieve more precise results which, however, requires extended simulation time.

Time stepping comparison
In Fig. 8 the sensitivity to the utilized time stepping is more clearly visible. In this case, the outflow temperatures are compared to represent the precision of the results. As stated in Table 4 the smaller the time stepping size, the corresponding simulation results will achieve higher accuracy. The solver for the linear equation system (cf. Eq. (6)) generally tries to determine a pair of in-and outflow temperature values which satisfies the thermal load given for the respective month. As found earlier, FEFLOW and the alternative OGS scenario using the outflow temperature of the previous time step are more sensitive. With larger time step size, the outflow temperature will deviate more from the previous step value, thus increasing the error. Additionally, since FEFLOW choose to solves Eq. (6) without iterations, the created numerical error will not be suppressed and tends to accumulate over time.
In Fig. 8 the outflow temperatures of both OGS model 1 (using T n out ) and OGS model 2 (using T n−1 out ) as well as FEFLOW (using T n−1 out ) are compared. According to Table 4, OGS model 1 shows a good resilience against a varyring time stepping scheme, which is also illustrated in Fig. 8c. Even with large time steps this OGS scenario has the ability to produce precise results. The larger time step size only leads to more Picard iterations, which increases the total simulation time. The maximum difference between 1-day and 15-days time step size is only 0.28 K during the first year of operation.
In Fig. 8b results from the adapted OGS model 2 are illustrated. Comparing the outputs to the original OGS simulation shows some minor differences that can not be neglected. The greatest difference between the three time stepping setups amounts to 0.87 K . In comparison to the original method some more deviations can be observed. It can be noticed that using the outflow temperature of the last time step will distort the simulation results. As the Picard iteration is present, the distortion is less in comparison to FEFLOW where no iterations are applied. In contrast to the results in Fig. 8c, the one-day time stepping scenario has only relatively small deviations. However, there are recognizable differences occurring for both the 5-days and 15-days time stepping compared to model 1. From this comparison it can be concluded that the choice of using the outflow temperature from the previous time step results in significant differences in the final simulation results.
This effect is actually enlarged by the bigger time step size. As present in Fig. 8a, where the FEFLOW results with different time stepping scenarios are shown. In this case, no iterations are conducted in FEFLOW during the simulation. Compared to other scenarios, the temperature profiles in Fig. 8a deviate the most from each other. The maximum occurring temperature difference between 1-day and 15-days time step settings reaches 1.21 K . The original time step size in the FEFLOW model was chosen to be 5 days. When reducing the time step size to one day, the relative deviation from the target values could be improved to 1.05 % (cf. Table 4). On the other side, however, a choice of 15-days time step size will lead to a relative difference only increased by 0.32 % , but the simulation time can be saved significantly. Users of FEFLOW may choose the corresponding time step size to achieve satisfying results. This phenomenon only occurs in the FEFLOW simulations, but not in OGS.

Ground surface boundary condition
As mentioned earlier in the model configurations, there are three different boundary conditions settings applied on the surface of the model domain. Figure 9 shows that the surface temperature has very minor influence on the resulting in-and outflow temperature profiles of the BHEs. There is even no recognizable difference between the cases whether it is no-flux boundary or a constant temperature is imposed on the model surface. Caused by the relatively large size of the model, especially Fig. 9b shows that a constant temperature boundary condition at the top layer has no effect on the fluid temperatures in the BHEs. Even after 25 years of simulation, most of the top layer still has the initial temperature of 11.6 • C . This suggests that a constant boundary condition would be redundant or even sometimes unflavoured, as the imposed surface temperature values may interfere with the BHE temperature values at the same location. Although the variable boundary condition with a fluctuating surface temperature better represents the reality, it only has negligible influence on the simulated temperature profile. The maximum temperature only deviates 1.5 % in comparison to the constant surface temperature setting. It can also be observed that the simulated mean BHE temperatures value equal 99.805 % , suggesting that the fluctuating surface boundary condition has one negligible effect. On the required computation resources, however, applying no boundary condition for the surface temperature the simulation time could be halved. It can be stated that decreasing the complexity of the model by disregarding any influence of the surface delivers sufficient results while shortening the computation time rapidly.

Influence of dispersion
In Fig. 10 the effect of hydrodynamic dispersivity on the simulation results is revealed. In this case the averaged fluid temperature of all BHEs are simulated in two configurations, one with and the other without the hydrodynamic dispersion effect. Already in the first year of simulation, there is a temperature difference up to 0.96 K between the two scenarios. Over time, this behaviour accumulates to a maximum divergence of 1.38 K and the influence of the dispersion becomes more obvious. The relative deviation is illustrated in the lower part (scaled to the right axis) of Fig. 10. Similar to the graphs of the averaged fluid temperature, it fluctuates throughout the years. While in summer there is a relatively small difference to the original simulation, a greater deviation is present in each winter. In the later years of the simulation it reaches a quasi-steady state analogously to the absolute temperature values. The relative deviation, however, accounts to more than 20 % based on whether dispersion is considered in the model or not. Therefore it can be concluded that hydrodynamic dispersion has a strong impact on the simulated temperature profile. It is thus necessary to predict the longitudinal and transverse dispersion length as precisely as possible, so that accurate temperature distribution can be predicted by the numerical model.

Software comparison
In this paper two different numerical software are employed for the simulation of subsurface heat transfer induced by GSHP system operation. Based on the detailed analysis of modeling results it is proven that the open-source software OGS can produce at least equal results as FEFLOW. When predicting long-term temperature development over 25 years, results generated with a five-day time step size in OGS have been compared to results from FEFLOW with a one-day size. FEFLOW simulations with a coarser temporal discretization could not generate comparative outputs.
Digging deeper into the detailed solution methods there are some differences in computational scheme implemented in FEFLOW and OGS. The most influential is the non-iterative method utilized by FEFLOW in its solution process. By using the outflow Fig. 11 Temperature difference comparison in FEFLOW temperature of the last time step as a fixed value, there is an intrinsic error existing in the solution of linear equation system. The size of the error, however, decreases along with the time step size but can not be fully eliminated. Since numerical modellers always need to find a balance between the accuracy of the result and the required computational time, a certain inaccuracy might be acceptable in respective cases. Figure 11 further illustrate the precision issue of FEFLOW in a quantative manner. The figure is split into two parts, where the upper bar chart illustrates the temperature difference between in-and outlet of the BHE. The values are compared to the target temperature difference calculated from the imposed thermal load in each specific month (cf. Table 1). Since there are only small divergences, the difference to the target value is also plotted in the lower part of Fig. 11. The green bar chart and graph, respectively, label the results made with the in-and outflow temperature values from the current time step. To obtain this pair of values, the previous outflow temperature T n−1 out is used in Eq. (5) to calculate the inflow temperature of the current time step, which is then imposed as a fixed value in Eq. (6) to compute the temperature values in the model including the outflow temperature for the current time step. Equation (5) can be calculated directly since it displays a linear dependency, while Eq. (6), however, underlies typical numerical inaccuracies. It would therefore be actually more accurate to take the previous outflow temperature T n−1 out into account when calculating the temperature difference from a specific time step, which is illustrated by the orange graph. It can be seen that there is no difference between this method and the blue target values, because the temperature difference is directly linked by Eq. (5). As mentioned above, an alternative OGS scenario has been set up to reproduce this computation manner and the corresponding simulation results show an analogous behavior, although the reduced precision is mitigated due to the iterative calculation in OGS. As a general suggestion, modelers need to be always aware of this computational setting embedded in FEFLOW and need to adjust their time stepping scheme accordingly.
Besides the precision issue discussed above, FEFLOW currently has an advantage regarding the required simulation time. Utilizing a five-days time stepping with recording 3 took 114 h while OGS required 244 h for the whole 25 years of simulation. Currently, OGS is still being further developed to reduce the computation time for BHE model simulations. The approach to accelerate the simulation and how to utilize a parallel solver will be topic of our next publication.

Legal issues
In the Federal Mining Act of Germany, geothermal energy is deemed as a freely mineable resource (BBergG 1980) which means that it can be exploited on a personal property. In general, the legal basis for installing GSHP systems is the water law. For the purpose of operating GSHP systems with a power of more than 30 kW as this project in Berlin, however, there might be special regulations during the permission progress. To identify which legal regulations apply, the temperature difference induced by the GSHP system at the property border after a specific period of operation time is crucial. If it inflicts a temperature loss of a neighboured land in a certain amount, the jurisdiction changes from the usually considered water law to mining law and an application of concession has to be requested to exploit the extended subsurface. The temperature distribution in the subsurface is therefore simulated with numerical models, wherefore it is necessary, that the simulation can predict the development as precise as possible. Due to the federal legal structure in Germany, however, there is no consistent legal situation throughout all German states and the regulations sometimes differ from each other considerably.
As this paper references to a real building project in Berlin, the state of Berlin has a less restrictive permission process compared to others. Here it is common law, that after 25 years of operation a maximum temperature change of 3 K is permitted at the property borders to be excluded from the mining law. As seen in Fig. 6 this value for the specific GSHP system in Berlin accounts to a minimum of 11.2 • C in the downstream direction, despite of a relatively large BHE array. The predicted temperature drop remains under the permitted temperature change by far. Considering the legend of Fig. 6 is continuous, it can be concluded that there is no point in the numerical model where a critical temperature change will occur.
Based on the simulation results presented here, a reader may conclude that there is only negligible impact on the groundwater temperature by the operation of GSHP system. However, a neighborhood-scaled GSHP project as presented in this paper would have to follow the more complex mining law in the State of Hesse. There must be no temperature change at the property borders within a period of 50 years. Such strict target value leads to the problem, that practically every GSHP installation falls into the regulation of the mining law, because numerical models will always predict some temperature changes in the subsurface, even the value may just be half a degree and falls well in the tolerance of the monitoring sensors.
Although the regulations in each German state bear on the same Federal Mining Act, their interpretation differs among the regions in Germany. This situation causes different level of difficulty when installing GSHP systems for building heating and cooling. We would ambitiously propose, to better facilitate the heat and energy transition in Germany, the legal standard and administrative requirements should be carefully reviewed and unified in different states. Such improvement will greatly simplify the installation of GSHP system and extend its application.

Conclusion
In this study, a neighborhood-scale GSHP system model has been setup based on a building project in Berlin. It involves 16 BHEs to extract geothermal energy from the shallow subsurface. To guarantee a proper long-term performance of the system in advance, two comprehensive numerical simulations, one with the commercial software FEFLOW and the other with open-source code OGS, have been built with parameters generated from several preinvestigations and geological analyses at the site. It is the first time that OGS has been compared to an established numerical simulation software as FEFLOW. The most important finding of this work is that FEFLOW always adopts the outflow temperature from the previous time step to calculate the temperature distribution for the current time step. This means modelers need to apply a smaller time step size in order to to keep a similar accuracy. Additional main conclusions are: • The analysis of the simulation results shows that both software predicts similar fluid temperature development over the long term. • The groundwater flow and the thermal dispersivity are found to be the strongest factors on the temperature development in the subsurface. • Accurate results were obtained in FEFLOW using a one-day time step size, while comparative precision can be achieved by OGS with 5-days time stepping. • OGS is more resistant against a coarser temporal discretization than FEFLOW and can even achieve comparative results with 15-days time steps. • Adopting a smaller time stepping and calculating iteratively will lead to more accurate results to the expense of increased computational resources.
Summarizing all achieved results in this paper it can be stated, although there are some computational differences in the software FEFLOW and OGS, both can generate sufficient and comparable predictions and that OGS is a valid alternative for longterm heat transport process simulations in the shallow subsurface. Currently, OGS is being further developed to accelerate the numerical simulations which is a topic for following publications.