Finite Element Modeling of Geothermal Source of Heat Pump in Long-Term Operation †

: Model simulation allows to present the time-varying temperature distribution of the ground source for heat pumps. A system of 25 double U-shape borehole heat exchangers (BHEs) in long-term operation and three scenarios were created. In these scenarios, the di ﬀ erence between balanced and non-balanced energy load was considered as well as the inﬂuence of the hydrogeological factors on the temperature of the ground source. The aim of the study was to compare di ﬀ erent thermal regimes of BHEs operation and examine the inﬂuence of small-scale and short-time thermal energy storage on ground source thermal balance. To present the performance of the system according to geological and hydrogeological factors, a Feﬂow ® software (MIKE Powered by DHI Software) was used. The temperature for the scenarios was visualized after 10 and 30 years of the system’s operation. In this paper, a case is presented in which waste thermal energy from space cooling applications during summer months was used to upgrade thermal performance of the ground (geothermal) source of a heat pump. The study shows di ﬀ erences in the temperature in the ground around di ﬀ erent Borehole Heat Exchangers. The cold plume from the not-balanced energy scenario is the most developed and might inﬂuence the future installations in the vicinity. Moreover, seasonal storage can partially overcome the negative inﬂuence of the travel of a cold plume. The most exposed to freezing were BHEs located in the core of the cold plumes. Moreover, the inﬂuence of the groundwater ﬂow on the thermal recovery of the several BHEs is visible. The proper energy load of the geothermal source heat pump installation is crucial and it can beneﬁt from small-scale storage. After 30 years of operation, the minimum average temperature at 50 m depth in the system with waste heat from space cooling was 2.1 ◦ C higher than in the system without storage and 1.6 ◦ C higher than in the layered model in which storage was not applied.


Introduction
Increasing climate change awareness and the popularity of renewable energy sources are connected with the increasing attention to geothermal energy and devices such as heat pumps. According to [1], in 2050, electric heat pumps will become more common in most parts of the world. In 2017, there were 10,572,395 heat pumps units in operation in Europe [2]. This number is growing constantly. The main advantage of a heat pump solution is the possibility of providing not only heating during cold periods, but also space cooling during summer hot days. Ground source heat pumps (GSHP) are a popular solution, especially for bigger installations. Several regional studies of the shallow geothermal potential of heat pumps in different geological conditions have been conducted, for example, in Japan [3], Figure 1. Non-uniform grid model of the study area and 25 borehole heat exchangers (2D), crosssection lines (green). Spatial discretization of BHE nodes (in frame) according to [27].
To minimalize the influence on the result and calculation time, the size of the model domain was chosen after pre-calculations. Moreover, the second type of boundary conditions were applied to the model's northern and southern edges. A larger domain size does not affect the results. All boundary conditions applied to the models' domain are given in Figure 2. There were 3 types of boundary conditions applied in this model: the 1st boundary condition (Dirichlet) is temperature applied at the top surface; the 2nd boundary condition (Neumann) is heat-flux at the bottom surface and fluid-flux applied at the edges or in the aquifer elements (in the 3rd scenario); the 4th boundary conditions are borehole heat exchangers.  Non-uniform grid model of the study area and 25 borehole heat exchangers (2D), cross-section lines (green). Spatial discretization of BHE nodes (in frame) according to [27].
To minimalize the influence on the result and calculation time, the size of the model domain was chosen after pre-calculations. Moreover, the second type of boundary conditions were applied to the model's northern and southern edges. A larger domain size does not affect the results. All boundary conditions applied to the models' domain are given in Figure 2. There were 3 types of boundary conditions applied in this model: the 1st boundary condition (Dirichlet) is temperature applied at the top surface; the 2nd boundary condition (Neumann) is heat-flux at the bottom surface and fluid-flux applied at the edges or in the aquifer elements (in the 3rd scenario); the 4th boundary conditions are borehole heat exchangers.
Energies 2020, 13 A non-uniform grid of the area was used. In the vicinity of wells, grid points were refined to obtain a better accuracy (Figure 1). Optimal mesh refinement for nodes around BHEs were calculated according to the Nillert method [27]. The total model dimensions were set to 240 × 240 × 105 m. The mesh consists of 169,477 elements and 92,940 nodes.
To minimalize the influence on the result and calculation time, the size of the model domain was chosen after pre-calculations. Moreover, the second type of boundary conditions were applied to the model's northern and southern edges. A larger domain size does not affect the results. All boundary conditions applied to the models' domain are given in Figure 2. There were 3 types of boundary conditions applied in this model: the 1st boundary condition (Dirichlet) is temperature applied at the top surface; the 2nd boundary condition (Neumann) is heat-flux at the bottom surface and fluid-flux applied at the edges or in the aquifer elements (in the 3rd scenario); the 4th boundary conditions are borehole heat exchangers.   The multiple system considered consists of 25 borehole heat exchangers (BHEs). They are in 5 × 5 layout with 10-m spacing. The spacing was chosen to fulfill the heat pump's branch recommendations [30]. The BHEs are not connected in arrays, all (equal) parameters are set for single BHE. Each BHE is 100 m in length, which equates to 2500 m of heat exchangers' total length. The exchangers were of the double U-shaped type, which consist of two inlet pipes and two outlet pipes surrounded by grout (Figure 3). A single BHE setting is provided in Table 1. exchangers were of the double U-shaped type, which consist of two inlet pipes and two outlet pipes surrounded by grout ( Figure 3). A single BHE setting is provided in Table 1. In the simulation, an Eskilson and Claesson [31] (quasi-stationary) method, which assumes the local thermal equilibrium between all the elements of BHE during the whole time of the simulation, was applied. The difference in thermal energy per time between inlet and outlet of BHE is variable. It was set according to energy demand on a given day. The BHEs did not work longer than 6 hours a day. Their power was variable, from 0 to 3 kW. Each simulation was run for 30 years (262,800 h).

Thermal and Hydrogeological Settings
The model is considered as fully saturated. The initial temperature condition across the model was set to 9 °C (referenced ground temperature in the central Poland area). To show the thermal regime, heat flux from the ground during the year was set to the top surface according to the change of the earth's surface monthly average temperature (Table 2). Earth's temperature for the central Poland area is provided. These values were set as the first type of boundary conditions. Moreover, a geothermal heat flux of 0.06 W/m 2 assessed for this area was applied to the base surface of the model. Table 2. Earth's average temperature in the study area [32] and energy load during a year.  In the simulation, an Eskilson and Claesson [31] (quasi-stationary) method, which assumes the local thermal equilibrium between all the elements of BHE during the whole time of the simulation, was applied. The difference in thermal energy per time between inlet and outlet of BHE is variable. It was set according to energy demand on a given day. The BHEs did not work longer than 6 hours a day. Their power was variable, from 0 to 3 kW. Each simulation was run for 30 years (262,800 h).

Thermal and Hydrogeological Settings
The model is considered as fully saturated. The initial temperature condition across the model was set to 9 • C (referenced ground temperature in the central Poland area). To show the thermal regime, heat flux from the ground during the year was set to the top surface according to the change of the earth's surface monthly average temperature ( Table 2). Earth's temperature for the central Poland area is provided. These values were set as the first type of boundary conditions. Moreover, a geothermal heat flux of 0.06 W/m 2 assessed for this area was applied to the base surface of the model.
The BHE system was operating according to the assessed monthly load, as shown in Table 2.
In order to compare different operating regimes, three scenarios were taken into consideration: • 1st scenario: homogeneous parameters of the ground source, not balanced energy load (only heating in winter) • 2nd scenario: homogeneous parameters of the ground source, partially balanced energy load (heating in winter and cooling in summer) • 3rd scenario: parameters according to different layers of the model, not balanced energy (only heating in winter).
The hydrogeological and thermal parameters are shown in Table 3. The initial conditions of the models were obtain by applying the simulation in steady-state conditions first.  The area was poorly recognized in the field of shallow geothermal energy and this study reflects some assumptions according to thermal and hydrogeological parameters in the area. The model is considered as the confined one. Aquifers have good flow conditions. The fluid flow direction is from the south to the north, as the fluid flux speed is constant, being limited by groundwater conductivity.
Other borders of the model are defined as impermeable. The geological settings of the model are based on variability observed in boreholes located in the vicinity of Poddebice town, central Poland [33].
In the 1st and 2nd scenarios, the whole model had good porosity and permeability conditions and the groundwater flowed throughout the model. In the 3rd scenario, the permeability conditions worsened and the water flow was mostly limited to an aquifer of 10 m. The aquifer layer has better conditions of thermal and fluid conductivity ( Figure 4).

Results and Discussion
In this project, three sets of simulations were run to verify different scenarios for the performance of heat pump geothermal heat exchangers. For each simulations, two check-points were selected: after 10 and 30 years of the system operation, at 87,600 and 262,800 h, respectively.

First Scenario: Homogeneous Parameters of the Ground Source, Not Balanced Energy Load
According to the simulation, the temperature of the ground decreased with time. At the beginning stage, the drop of temperature was most visible in the centre between BHEs. The difference in temperature between the points of the minimum temperature and the surrounding domain after the first year of operation is ca. 4.5 °C. Later, the influence of the water flow was more visible in the model and since the second year of the system operation, the points with minimum temperature shifted towards the north. At the same time, different BHEs started to work on different ground parameters. In Figure 5, the temperature model is presented for a depth of 50 m (the mid-section of the BHE). After 10 years of system operation (A), an increase in the cold plume can be observed. The range of the model zone where the temperature dropped to 5 °C was ca. 45 m from the centre of the BHE system and ca. 24 m from the BHE on the edge. During this time, the minimum ground temperature dropped to 0.1 °C.

Results and Discussion
In this project, three sets of simulations were run to verify different scenarios for the performance of heat pump geothermal heat exchangers. For each simulations, two check-points were selected: after 10 and 30 years of the system operation, at 87,600 and 262,800 h, respectively.

First Scenario: Homogeneous Parameters of the Ground Source, Not Balanced Energy Load
According to the simulation, the temperature of the ground decreased with time. At the beginning stage, the drop of temperature was most visible in the centre between BHEs. The difference in temperature between the points of the minimum temperature and the surrounding domain after the first year of operation is ca. 4.5 • C. Later, the influence of the water flow was more visible in the model and since the second year of the system operation, the points with minimum temperature shifted towards the north. At the same time, different BHEs started to work on different ground parameters. In Figure 5, the temperature model is presented for a depth of 50 m (the mid-section of the BHE). After 10 years of system operation (A), an increase in the cold plume can be observed. The range of the model zone where the temperature dropped to 5 • C was ca. 45 m from the centre of the BHE system and ca. 24 m from the BHE on the edge. During this time, the minimum ground temperature dropped to 0.1 • C. After 30 years of system operation (B), the decrease of temperature was even more visible and the cold plume travelled further north due to the fluid flow conditions. The minimum average temperature at the same surface of 50 m depth was -0.7 °C, which implies the risk of ground freezing. Particularly, the refrigeration fluid temperature in several BHEs can even drop to -3.5 °C. After 30 years of system operation (B), the range of the model zone in which the temperature dropped to 5 °C was ca. 100 m from the centre of the BHE system and ca. 75 m from the BHE on the north edge. The distribution of temperature in 3D view is shown in Figure 6. In the most affected BHEs, the average working fluid's temperature decreased significantly during the first 15 years of operation. Then, the drop of the temperature was rather small as time passed (Figure 7). The difference between the warmest and the coldest BHE's average temperature in this layout is 4.5 °C after 30 years of operation. After 30 years of system operation (B), the decrease of temperature was even more visible and the cold plume travelled further north due to the fluid flow conditions. The minimum average temperature at the same surface of 50 m depth was −0.7 • C, which implies the risk of ground freezing. Particularly, the refrigeration fluid temperature in several BHEs can even drop to −3.5 • C. After 30 years of system operation (B), the range of the model zone in which the temperature dropped to 5 • C was ca. 100 m from the centre of the BHE system and ca. 75 m from the BHE on the north edge. The distribution of temperature in 3D view is shown in Figure 6. After 30 years of system operation (B), the decrease of temperature was even more visible and the cold plume travelled further north due to the fluid flow conditions. The minimum average temperature at the same surface of 50 m depth was -0.7 °C, which implies the risk of ground freezing. Particularly, the refrigeration fluid temperature in several BHEs can even drop to -3.5 °C. After 30 years of system operation (B), the range of the model zone in which the temperature dropped to 5 °C was ca. 100 m from the centre of the BHE system and ca. 75 m from the BHE on the north edge. The distribution of temperature in 3D view is shown in Figure 6. In the most affected BHEs, the average working fluid's temperature decreased significantly during the first 15 years of operation. Then, the drop of the temperature was rather small as time passed (Figure 7). The difference between the warmest and the coldest BHE's average temperature in this layout is 4.5 °C after 30 years of operation. In the most affected BHEs, the average working fluid's temperature decreased significantly during the first 15 years of operation. Then, the drop of the temperature was rather small as time passed (Figure 7). The difference between the warmest and the coldest BHE's average temperature in this layout is 4.5 • C after 30 years of operation. The ground water flow provides the differences in temperature. The temperature in some BHEs (in the centre and north) dropped due the cold plume travelling. However, water of a higher temperature moved from the south to the BHE system. The southern edges of BHE system, which recover faster, have a higher temperature. In fact, the average temperature in the southern BHEs was almost constant after the initial first few years. For example, the minimal average temperature of the BHE on the southern edge stayed at level of ca. 4 °C after the fifth year.
The vertical distribution of temperature is given in Figure 8. In cross-sections A-A' and B-B', it is clearly visible that the northern side of the BHE system has a significantly lower temperature than the southern one. The BHEs' had an influence on each other and there was a decrease in temperature in the central BHEs. The ground water flow provides the differences in temperature. The temperature in some BHEs (in the centre and north) dropped due the cold plume travelling. However, water of a higher temperature moved from the south to the BHE system. The southern edges of BHE system, which recover faster, have a higher temperature. In fact, the average temperature in the southern BHEs was almost constant after the initial first few years. For example, the minimal average temperature of the BHE on the southern edge stayed at level of ca. 4 • C after the fifth year.
The vertical distribution of temperature is given in Figure 8. In cross-sections A-A' and B-B', it is clearly visible that the northern side of the BHE system has a significantly lower temperature than the southern one. The BHEs' had an influence on each other and there was a decrease in temperature in the central BHEs.

Second Scenario: Homogeneous Parameters of the Ground Source, Balanced Energy Load (Heating in Winter and Cooling in Summer)
A balanced energy load affects the drop in the system's temperature less. The average temperature of the system did not drop below 1 • C during the 30 years of operation.
At the beginning stage, the drop in temperature is the most visible in the centre, between BHEs. The difference in temperature between the points of the minimum temperature and the surrounding domain after the first year of operation is ca. 4.4 • C.
Starting from the second year, different BHEs started to work on different ground parameters due to the interaction between the BHEs and water flow. In Figure 9, the temperature model is presented for a depth of 50 m (the mid-section of the BHE). After 10 years of system operation (A), the zone with a temperature of 5 • C was concentrated near BHEs and was not present outside of the BHEs layout.

Second Scenario: Homogeneous Parameters of the Ground Source, Balanced Energy Load (Heating in Winter and Cooling in Summer)
A balanced energy load affects the drop in the system's temperature less. The average temperature of the system did not drop below 1 °C during the 30 years of operation. The distribution and layout of the temperature in the 3D model is provided in Figure 9. After 10 years of operation, the minimum ground temperature dropped to 1.9 • C at a 50-m depth.
After 30 years of system operation (Figure 10), despite the partly balanced energy load, there is a cold plume observed to be travelling north, according to the fluid flow direction. However, its shape is more slender than those observed in the first scenario. The range of the model zone in which the temperature dropped to 5 • C is ca. 50 m from the centre of the BHE system and ca. 30 m from the BHE on the north edge. The minimum average temperature at the surface of a 50 m depth was 1.4 • C, which dismiss the risk of ground freezing. Starting from the second year, different BHEs started to work on different ground parameters due to the interaction between the BHEs and water flow. In Figure 9, the temperature model is presented for a depth of 50 m (the mid-section of the BHE). After 10 years of system operation (A), the zone with a temperature of 5 °C was concentrated near BHEs and was not present outside of the BHEs layout.
The distribution and layout of the temperature in the 3D model is provided in Figure 9. After 10 years of operation, the minimum ground temperature dropped to 1.9 °C at a 50-m depth. After 30 years of system operation ( Figure 10), despite the partly balanced energy load, there is a cold plume observed to be travelling north, according to the fluid flow direction. However, its shape is more slender than those observed in the first scenario. The range of the model zone in which the Starting from the second year, different BHEs started to work on different ground parameters due to the interaction between the BHEs and water flow. In Figure 9, the temperature model is presented for a depth of 50 m (the mid-section of the BHE). After 10 years of system operation (A), the zone with a temperature of 5 °C was concentrated near BHEs and was not present outside of the BHEs layout.
The distribution and layout of the temperature in the 3D model is provided in Figure 9. After 10 years of operation, the minimum ground temperature dropped to 1.9 °C at a 50-m depth. After 30 years of system operation ( Figure 10), despite the partly balanced energy load, there is a cold plume observed to be travelling north, according to the fluid flow direction. However, its shape is more slender than those observed in the first scenario. The range of the model zone in which the Figure 10. Results of the 3D simulation of ground temperature after 10 years (a) and 30 years (b) of heat pump operation within the 2nd scenario (homogeneous model); the view from the north-west.
In the most affected BHEs, the temperature decreased significantly during the first 15 years of operation. Then, the drop in temperature was rather small over time ( Figure 11). The difference between the warmest and the coldest BHE's average temperature was ca. 3.1 • C after 30 years of operation. The vertical distribution of the temperature is given in Figure 12. In cross-sections A-A' and B-B', smaller difference between the northern and southern site can be observed compared to the first scenario. The influence of the system into the neighboring ground is lower than in the first scenario (C-C', D-D'). As the groundwater flow was set as homogenous, the BHEs in the centre were the most affected by the heat withdrawal.

Third Scenario (Layered Model): Parameters According to Different Layers of the Model, Not Balanced Energy (Only Heating in Winter)
In Figure 13, the temperature of the ground source is shown after 10 and 30 years of the simulation at a 50-m depth. After 10 years of system operation, the temperature at a depth of 50 m was 2.7 • C. The cold plume is rather round-shaped and did not shift significantly towards the north (it was not affected by the ground water flow).
Energies 2020, 13, x FOR PEER REVIEW 13 of 19 In Figure 13, the temperature of the ground source is shown after 10 and 30 years of the simulation at a 50-m depth. After 10 years of system operation, the temperature at a depth of 50 m was 2.7 °C. The cold plume is rather round-shaped and did not shift significantly towards the north (it was not affected by the ground water flow). After 30 years of operation, the cold plume travelled slightly north and was slightly more visible in the permeable layers. At a depth of 50 m, the range of the model zone in which the temperature dropped to 5 °C was ca. 45 m from the centre of the BHE system and ca. 25 m from the BHE on the northern edge. The most visible decrease of the temperature is observed inside the BHEs' system layout ( Figure 14). After 30 years of operation, the cold plume travelled slightly north and was slightly more visible in the permeable layers. At a depth of 50 m, the range of the model zone in which the temperature dropped to 5 • C was ca. 45 m from the centre of the BHE system and ca. 25 m from the BHE on the northern edge. The most visible decrease of the temperature is observed inside the BHEs' system layout ( Figure 14). In Figure 13, the temperature of the ground source is shown after 10 and 30 years of the simulation at a 50-m depth. After 10 years of system operation, the temperature at a depth of 50 m was 2.7 °C. The cold plume is rather round-shaped and did not shift significantly towards the north (it was not affected by the ground water flow). In the northern BHEs, there is a risk of ground freezing, while the ground temperature can drop to -0.2 °C (layers with the worst thermal conditions). At the same time, the average temperature of the refrigerant can be as low as -2 °C at the coldest time of the season ( Figure 15). The figure shows the average temperature of BHE's working fluid temperature within the 3rd scenario (system without cooling) within a single year. This was for the last and 30th year considered in the simulation. In the northern BHEs, there is a risk of ground freezing, while the ground temperature can drop to −0.2 • C (layers with the worst thermal conditions). At the same time, the average temperature of the refrigerant can be as low as −2 • C at the coldest time of the season ( Figure 15). The figure shows the average temperature of BHE's working fluid temperature within the 3rd scenario (system without cooling) within a single year. This was for the last and 30th year considered in the simulation. In the northern BHEs, there is a risk of ground freezing, while the ground temperature can drop to -0.2 °C (layers with the worst thermal conditions). At the same time, the average temperature of the refrigerant can be as low as -2 °C at the coldest time of the season ( Figure 15). The figure shows the average temperature of BHE's working fluid temperature within the 3rd scenario (system without cooling) within a single year. This was for the last and 30th year considered in the simulation.  Figure 16 depicts differences in the vertical distribution of temperature around BHEs. Some sections or layers interfered more strongly than the others due to the different hydrogeological and thermal parameters of the layers. The influence of the groundwater flow in the aquifer is less significant to the system.  Figure 15. Average temperature of BHEs within the 3rd scenario (system without cooling) during the 30th year of operation. Figure 16 depicts differences in the vertical distribution of temperature around BHEs. Some sections or layers interfered more strongly than the others due to the different hydrogeological and thermal parameters of the layers. The influence of the groundwater flow in the aquifer is less significant to the system.  Figure 17 shows the comparison between the cold plumes that appeared after 30 years of operation for the different scenarios. A temperature value of 5 • C was arbitrary chosen as the cold plume boundary. Normally, heat pumps can work on temperatures below 5 • C; however, lower temperatures are no clearly visible in these 3D models.
Energies 2020, 13, x FOR PEER REVIEW 16 of 19 Figure 16. The temperature distribution as a cross-section view for the 3rd scenario. The cross-section lines are the same as those used in Figure 1. Figure 17 shows the comparison between the cold plumes that appeared after 30 years of operation for the different scenarios. A temperature value of 5 °C was arbitrary chosen as the cold plume boundary. Normally, heat pumps can work on temperatures below 5 °C ; however, lower temperatures are no clearly visible in these 3D models. Figure 17. Cold plumes in 3D models for 3 different scenarios (1st scenario, 2nd scenario and layered model (3rd scenario)) after 30 years of system operation. The cold plume from the 1st scenario is the most developed in length and width. It travelled further to the north according to the ground flow direction. This BHE system might threats future installations in the vicinity.
In the 2nd scenario, due to the balanced energy load and energy storage, the range of the cold plume is definitely smaller, even also directed to the north (groundwater flow). Moreover, thermal energy storage (understood here as cooling of space in summer) can partially overcome the negative influence of traveling of a cold plume. The most exposed to freezing are BHEs from the north edge of the system located at the core of the plumes. In the 1st and 2nd scenario, the influence of the groundwater flow on the thermal recovery of the southern BHEs is clearly visible. They worked with a higher temperature regime the whole time.
In the 3rd scenario, where the best groundwater conditions are limited to the aquifer, the heat conductions have a greater influence on the system (more rounded shape of the plume). There was no significant thermal recovery in the south. The cold plume travelled much more slowly to the north (similarly to the 2nd scenario with a balanced energy load) and the cold plume shift was slight.
The model allows to simulate the time-varying temperature of the ground source for heat pumps. The temperature regime is important for proper heat pump operation. Even the ground source heat pumps can work with a temperature below 0 • C; this temperature is associated with the freezing of the ground and water in the BHE vicinity, which can have a negative impact on the BHE itself (cracks, air presence, etc.) [30]. It was observed that after 15 years, the system's temperature stabilized and the system's elements worked with similar parameters. However, after that long period of exploitation, there is a risk of ground freezing in several locations. It is important to determine how the cold plume location changes, how it propagates and how BHEs interfere with each other and with the cold plume.
In the layered model, the influence of the groundwater flow was limited to the aquifer and provided lower temperatures in the centre BHEs. The thermal parameters' diversity in each section of the BHE have a lower influence on the thermal layout of the whole system than the balanced energy load. This balance can be obtained by small-scale and short-time thermal energy storage, which requires the use of waste heat from space cooling during the summer months. In the system with small-scale energy storage, the zone of a temperature of 5 • C is two times smaller and the cold plume does not travel as far as in the system without storage.
The proper energy load of the geothermal source heat pump installation is crucial, and it can benefit from small-scale storage. After 30 years of operation, the minimum average temperature at a 50-m of depth in the system with waste heat from space cooling was 2.1 • C higher than in the system without storage and 1.6 • C higher than in the layered model where storage was not applied.
The groundwater flow influences the improvement of parameters of BHEs located in the south (in the model). At the same time, it shifts the cold plume to the north, which negatively influences the northern BHEs. BHEs that are the most affected by the cold plume or by interference with other BHEs can be found thanks to the simulation.
Balanced energy load (including e.g., summer cooling) can increase system efficiency by providing higher temperatures. The study shows differences in the temperature in the ground around BHEs. Therefore a proper layout of the BHE system can be adopted prior to the system construction.

Conclusions
Geological and hydrogeological conditions have a great impact on BHEs' system performance. Thus, good recognition of the parameters and suitable long-term simulation are necessary prior to installing a sustainable ground source for heat pumps. Installations with seasonal storage (if possible) are preferred for long-time operation due to their better ground thermal recovery which causes a higher effectiveness of heat pumps' installation and prevents ground freezing. FEM modelling of such installations would answer the question of whether a given heat pump and BHE's installation can work efficiently in specific conditions and locations of installation.