Assessment of the E ﬀ ective Variants Leading to Higher E ﬃ ciency for the Geothermal Doublet, Using Numerical Analysis-Case Study from Poland (Szczecin Trough)

: Numerical models of geothermal doublet allows us to reduce the high risk associated with the selection of the most e ﬀ ective location of a production well. Furthermore, modeling is a suitable tool to verify possible changes in operational geothermal parameters, which guarantees liveliness of the system. An appropriate selection of software as well as the methodology used to generate numerical models signiﬁcantly a ﬀ ects the quality of the obtained results. In this paper, the authors discuss the inﬂuence of such parameters as grid density and distance between wells on the e ﬃ ciency of geothermal heating plant. The last stage of the analysis was connected with estimation of geothermal power potential for a hypothetical geothermal doublet. Numerical simulations were carried out using the TOUGH2 code, which applies the ﬁnite-di ﬀ erence method. The research was conducted in the Szczecin Trough area (NW Poland), based on archival data from Choszczno IG-1 well. The results demonstrated that in the studied case of the Choszczno region, the changes in the distance of boreholes can have a visible inﬂuence on obtained results; however the grid density of the numerical model did not achieve a signiﬁcant impact on it. The results show the signiﬁcant importance of numerical modeling aimed at increasing the e ﬃ ciency of a potential geothermal heating plant.


Introduction
Geothermal waters and energy utilization have a beneficial impact that improves the living conditions in any region where existing geothermal resources enable its effective use [1,2]. It is the kind of renewable energy whose resources are available regardless of weather conditions throughout the year and around all day. However, the basic condition for successful investment is the existence of appropriate hydrogeothermal conditions at the place of its implementation [3][4][5]. In Poland, low-temperature geothermal resources occur in many parts of the country. Research carried out so far confirms that the main resources of geothermal waters in the Polish Lowlands are accumulated in the of the computation process depending on the studied case. Generally, higher density of calculation grid increases the accuracy of modeling, but on the other hand, it extends the time of computations. As long as this time increases from several to some tens of minutes, the grid densification does not cause problems. If, however, such an extension prolongs to some tens of hours (or even more), the computation time turns into a serious obstacle, which is the common case of more complicated models. For simple models, in which the succeeding rock layers are almost parallel, we can select quite a sparse calculation grid without declining the quality of modeling. The factor usually analyzed with the numerical modeling is the mutual influence of geothermal wells. Such modeling enables us to assess the distance between production and injection boreholes to extend the lifetime of the system but also to decrease the energy transfer loss and to reduce the costs of connection pipelines. Considering the above, TOUGH2 code was chosen for numerical analysis, which is described in numerous scientific publications [49][50][51][52][53].
The main aim of the study is to evaluate two parameters, namely, calculation grid density and distance between wells on the results of simulation carried out by using TOUGH2 code. Another objective is to determine the geothermal potential in the Szczecin Trough, based on the archival Choszczno IG-1 well. Numerical simulations are a common tool for geothermal investments especially in the initial stage of the design process, but also by improving the operating ones. The article proposes a methodology for the effective improvement of the simulation by selecting a computational grid density that will guarantee precise results while reducing simulation time. Attention was also paid to simulations related to the distance between production and injection wells that has an impact on the system's work but also on investment costs. In addition, based on the derived real hydrogeological model of the Szczecin Trough area, a numerical analysis of geothermal utilization was carried out in order to assess the thermal power of the source and the amount of energy obtained during operation. The temperature and pressure of extracted and injected fluids into the reservoir was estimated for 50 years of system activity. The obtained results are essential for the development of a sustainable geothermal heating plant and constitute a unique case study for exploitation of similar geothermal resources.
The paper is organized as follows. Summaries of a geological background that focuses on Szczecin Trough (NW part of Poland) are presented in Section 2. The numerical model is given in Section 3. The exposition of the results of numerical simulation is shown in Section 4. A brief estimation of the thermal power is presented in Section 5. Discussion considering three important issues, namely: (1) the changes of grid density; (2) the impact of distance between the geothermal wells; (3) estimation of geothermal power potential for a hypothetical geothermal doublet, is demonstrated in Section 6. Concluding remarks are formulated in Section 7.

Geological Background
Numerical modeling was performed for the Choszczno Anticline located within the geological structure named the Szczecin Trough (NW part of Poland) (Figure 1). For modeling purposes, the Choszczno IG-1 well had been selected as a geothermal production well. The Choszczno IG-1 well is located in the vicinity of Szczecin town.
The lithostratigraphic profile of the Szczecin Through is built of partially eroded and folded Carboniferous-Devonian deposits covered with Rotliegend, Zechstein, and Mesozoic deposits. The most promising groundwater horizons occurs in the Lower Jurassic sandstones, mostly Hettangian and Sinemurian. Pliensbachian and Upper Toarcian deposits have less favorable reservoir parameters [20,21].
In the Choszczno area, the Lower Jurassic succession occurs at a depth interval from 1164.5 to 1468.0 m below sea level (bsl) and reaches a thickness of 303.5 m (Figure 2). Beneath the Lower Jurassic reservoir, the Upper Triassic clayey-sandy sediments occur. These are underlain by Middle Triassic (Muschelkalk) carbonates. Total thickness of Upper-Middle Triassic strata reaches 444.5 m (as revealed by data from the Pławno-1 well, see [55]). The Lower Jurassic succession is covered by Middle Jurassic sediments, 114.0 m thick, followed by Upper Jurassic, alternating sandstones and mudstones of thickness 53 m. The Jurassic formation is overlain by Cretaceous sediments, 842.8 m thick, followed by Tertiary (6 m thick) and Quaternary (148.7 m thick) stratum [56,57]. In the Choszczno area, the Lower Jurassic succession occurs at a depth interval from 1164.5 to 1468.0 m below sea level (bsl) and reaches a thickness of 303.5 m (Figure 2). Beneath the Lower Jurassic reservoir, the Upper Triassic clayey-sandy sediments occur. These are underlain by Middle Triassic (Muschelkalk) carbonates. Total thickness of Upper-Middle Triassic strata reaches 444.5 m (as revealed by data from the Pławno-1 well, see [55]). The Lower Jurassic succession is covered by Middle Jurassic sediments, 114.0 m thick, followed by Upper Jurassic, alternating sandstones and mudstones of thickness 53 m. The Jurassic formation is overlain by Cretaceous sediments, 842.8 m thick, followed by Tertiary (6 m thick) and Quaternary (148.7 m thick) stratum [56,57].   In the Choszczno area, the Lower Jurassic succession occurs at a depth interval from 1164.5 to 1468.0 m below sea level (bsl) and reaches a thickness of 303.5 m ( Figure 2). Beneath the Lower Jurassic reservoir, the Upper Triassic clayey-sandy sediments occur. These are underlain by Middle Triassic (Muschelkalk) carbonates. Total thickness of Upper-Middle Triassic strata reaches 444.5 m (as revealed by data from the Pławno-1 well, see [55]). The Lower Jurassic succession is covered by Middle Jurassic sediments, 114.0 m thick, followed by Upper Jurassic, alternating sandstones and mudstones of thickness 53 m. The Jurassic formation is overlain by Cretaceous sediments, 842.8 m thick, followed by Tertiary (6 m thick) and Quaternary (148.7 m thick) stratum [56,57]. Geological cross-section through Szczecin Trough according to [36], modified. The black rectangle shows the range of the numerical model for which simulations of the geothermal doublet were carried out. Geological cross-section through Szczecin Trough according to [36], modified. The black rectangle shows the range of the numerical model for which simulations of the geothermal doublet were carried out.
In the Szczecin Trough, the density of heat flux to the Earth's surface is the highest in Poland and varies from almost 70 to over 100 mW/m 2 . The minimal values (about 70-75 mW/m 2 ) are recorded along the eastern margin of the trough, whereas the highest ones (over 90 mW/m 2 ) occur in its southern part. In the Choszczno area, heat flux density is about 80-85 mW/m 2 [58]. Such high flux corresponds to the pattern of subsurface temperatures, which vary from 50 • C in the top part to about 60 • C in  [20,21,59]. Such values correlate well with the temperature log recorded in the Choszczno IG-1 well from which both the geothermal step and geothermal gradient were calculated [60]. The geothermal step is 52.4 m/ • C, hence, temperatures of the groundwaters at about 1390 m bsl reach 45.5 • C. Water mineralization is around 100-125 g/dm 3 , while the potential discharge of wells is about 200 m 3 /h [20,21].
The results of laboratory studies carried out at the Polish Geological Institute for samples from the Choszczno IG-1 well provided the mean values of effective porosity and density for particular stratigraphic units. The best reservoir properties were found in both the Upper and Lower Sinemurian sediments, for which the effective porosity exceeds 20% [61].

Numerical Methodology and Mathematical Model
The numerical model covered the area of 11.5 × 10.0 × 2.3 km (Figure 3). The geothermal doublet was based upon the archival information from the Choszczno IG-1 well, which was defined in the model as the production well. In this section a mathematical model of the whole system is presented also. Governing equations for fluids are strictly connected with the elementary balance of mass, momentum, and energy. However, the governing equation in solids is focused on the balance of energy. For porous rocks in the governing equations, porosity should be considered. In addition, the porosity of rocks, as follows: ε = V f /V, is considered in governing equations. Typically, for each geological formation the values of rock porosity are the factors that take into account the amounts of fluid volume V f in respect to the whole volume of rocks formation V. This paper uses the traditional system of equations for numerical modelling, namely [62,63]: in which: ∂ ∂τ is local derivative, div means divergence, ρ; ρv; ρe represents the relevant conserved variable vector. These variables (ρ; ρv; ρe) are generated from the well-known equations: balance of mass (ρ), the balance of momentum (ρv-three equations), the balance of energy (ρe). However, ρ = ρ(x, τ) represents the density that can change in time τ and space x. Indices f and s refer to fluid and solid, respectively. v = v i e i is velocity vector including e i (versor in one direction defines by subscript i) and v i (scalars magnitude of vector v); ρv is the momentum density vector, p represents pressure; I = δ ij e i e j is unit tensor, δ ij means Kronecker's delta. Additionally, diffusive momentum flux t takes into account laminar viscous stress flux. ρS v is momentum sources which includes Darcy's force acting as an extra source in the momentum balance equation. The total specific energy e = u + 1 2 v 2 consists specific kinetic energy 1 2 v 2 and specific internal energy u = ct, where c is specific heat capacity and t is temperature. q is the diffusive heat flux. Last but not the least, S f e ; S s e represent the energy sources in fluid and solid body, respectively. Measuring fluid mass flow rates in porous materials may help to predict the sum of effects of both the thermal transpiration and surface velocity slip. An additional bulk influence of viscous flow (Poiseuille's type) is visible for the pores greater than 2-5 times the free length of molecules. In this specific case, the mass flow rate is the common effect of one bulk and two surface phenomena. Therefore, the resultant flow velocity of filtration is governed by the Poiseuille-Knudsen-Reynolds equation [64].
Based on the integrated finite difference method (IFDM), the simulator discretizes both space and time. In the case of space, it is discretized to form basic differential equations. The TOUGH code has been continuously developed and successively validated for more than 35 years, as evidenced by the work initiated by Pruess [65][66][67][68] and subsequent numerous tests confirming its reliability based on benchmark experiments and multi-variant validation [69][70][71][72][73][74][75][76][77][78]. Therefore, the authors of this paper have based their work on previously conducted validations and focused only on model calibration based on available geological parameters. Also, the issues concerning the generation of the discretization grid and its usefulness for modelling various physical phenomena have been developed in recent years, which is evidenced by a number of works [79][80][81][82]. Therefore, it was based on previous works on mesh for TOUGH2 code and prepared the grid as follows. The number of elements with defined volume is determined, as well as the method of connection between elements, enabling the creation of regular or irregular one-, two-, or three-dimensional interpolation meshes. Time discretization occurs by determining the number of simulation steps and can also apply to its total duration. Intervals for individual simulation steps were automatically selected by the simulator. Numerical simulations were conducted in PetraSim, which uses TOUGH code for the calculating process. Density of computational grid was changed by the maximum area near the wells located in the Choszczno Region, but without any changes in number of grid elements that equals 1616 elements.
At the first stage of the simulation, the hypothetical injection well was located about 1000 m from the production borehole. After implementation of both wells into the model, the polygonal calculation grid was superimposed, and calibration process was carried out again. The production well was designed at a depth of 1360 m. The injection borehole reached the depth of 1420 m. Both wells were designed at the level of Upper Synemurian Beds groundwater horizons, which guaranteed a closed groundwater circulation system within the Lower Jurassic geothermal reservoir. The yield of the geothermal doublet was planned as 120 m 3 /h and the lifetime of the installation was assumed as 50 years, at the temperature of injected water at the level of 25 • C.
Energies 2020, 13, x FOR PEER REVIEW 6 of 20 on benchmark experiments and multi-variant validation [69][70][71][72][73][74][75][76][77][78]. Therefore, the authors of this paper have based their work on previously conducted validations and focused only on model calibration based on available geological parameters. Also, the issues concerning the generation of the discretization grid and its usefulness for modelling various physical phenomena have been developed in recent years, which is evidenced by a number of works [79][80][81][82]. Therefore, it was based on previous works on mesh for TOUGH2 code and prepared the grid as follows. The number of elements with defined volume is determined, as well as the method of connection between elements, enabling the creation of regular or irregular one-, two-, or three-dimensional interpolation meshes. Time discretization occurs by determining the number of simulation steps and can also apply to its total duration. Intervals for individual simulation steps were automatically selected by the simulator. Numerical simulations were conducted in PetraSim, which uses TOUGH code for the calculating process. Density of computational grid was changed by the maximum area near the wells located in the Choszczno Region, but without any changes in number of grid elements that equals 1616 elements.
At the first stage of the simulation, the hypothetical injection well was located about 1000 m from the production borehole. After implementation of both wells into the model, the polygonal calculation grid was superimposed, and calibration process was carried out again. The production well was designed at a depth of 1360 m. The injection borehole reached the depth of 1420 m. Both wells were designed at the level of Upper Synemurian Beds groundwater horizons, which guaranteed a closed groundwater circulation system within the Lower Jurassic geothermal reservoir. The yield of the geothermal doublet was planned as 120 m 3 /h and the lifetime of the installation was assumed as 50 years, at the temperature of injected water at the level of 25 o C. The numerical model of the Choszczno region included 8 stratigraphic formations. For each geological formation the values of rock density, porosity, permeability, thermal conductivity, and specific heat value were defined individually as shown in Table 1. The values were introduced from well data [55], most of them were average values for subsequent stratigraphic formations. Values of thermal conductivity coefficient and specific heat were introduced with reference to lithological data [55]. Table 1. Parameters defined for the model of Choszczno region [55]. The numerical model of the Choszczno region included 8 stratigraphic formations. For each geological formation the values of rock density, porosity, permeability, thermal conductivity, and specific heat value were defined individually as shown in Table 1. The values were introduced from well data [55], most of them were average values for subsequent stratigraphic formations. Values of thermal conductivity coefficient and specific heat were introduced with reference to lithological data [55].
The simulations were started from a conceptual model for which a regular calculation grid was implemented. Constant values of pressure and temperature were attributed to the top surface (in contact with the atmosphere) and the bottom surface, using the first type (Dirichlet) boundary condition. Then, the model was calibrated with the measured values. The results of calibration are shown in Figure 4. Due to the availability of the temperature curve from the Choszczno IG-1 well, the curve was used for the calibration process. During calibration, the permeability and thermal conductivity values were mainly changed. The results obtained during the calibration process were compared with the Energies 2020, 13, 2174 7 of 20 temperature curve from the Choszczno IG-1 well. Figure 4 shows the final calibration result which was considered sufficient for further analysis. The simulations were started from a conceptual model for which a regular calculation grid was implemented. Constant values of pressure and temperature were attributed to the top surface (in contact with the atmosphere) and the bottom surface, using the first type (Dirichlet) boundary condition. Then, the model was calibrated with the measured values. The results of calibration are shown in Figure 4. Due to the availability of the temperature curve from the Choszczno IG-1 well, the curve was used for the calibration process. During calibration, the permeability and thermal conductivity values were mainly changed. The results obtained during the calibration process were compared with the temperature curve from the Choszczno IG-1 well. Figure 4 shows the final calibration result which was considered sufficient for further analysis. The obtained model was evaluated as a good representation of reality. Then, the regular grid was changed to a polygonal type, which allowed increasing density of the mesh near the boreholes. The calibration process was refined. Before starting numerical simulations, the model was recalibrated to eliminate any errors related to the distribution of initial conditions in the newly defined model grid. The estimated results were compared once again with temperature from the Choszczno IG-1 well. The process was completed after receiving the best match. It should be mentioned that in The obtained model was evaluated as a good representation of reality. Then, the regular grid was changed to a polygonal type, which allowed increasing density of the mesh near the boreholes. The calibration process was refined. Before starting numerical simulations, the model was re-calibrated to eliminate any errors related to the distribution of initial conditions in the newly defined model grid. The estimated results were compared once again with temperature from the Choszczno IG-1 well. The process was completed after receiving the best match. It should be mentioned that in the study [78], similar accuracies were obtained in the compatibility between the temperature distribution obtained experimentally and numerically. Additionally, the study [78] shows that the model is calibrated and achieves convergence between the measured and calculated pressure distribution. Therefore, it was assumed that a similar convergence would be obtained for pressure if only a pressure distribution for the Choszczno IG-1 borehole would be available.
Energies 2020, 13, 2174 8 of 20 In this article, the direct "optimization" method, consisting in the cyclical change of the input calculation parameters (according to the established scenario) in order to approximate more refined and higher resolution of the proposed cases, is presented. The evaluation process includes two goal functions. The first concerned issues related to the numerical modelling and included the choice of resolution-computational grid density. Dense grids improve the accuracy of the solution, but unfortunately their use extends the calculation time, which significantly affects the selection process. This is particularly important for direct "optimization" methods, which were used in the article to look for the extreme of the second objective function. The second objective function was to determine the smallest safe distance between the wells of the geothermal doublet. The measure of safety was analyzed by the time when breakthrough of the cold front between the injection and production well, was observed. The smallest distance for which no significant decrease in the temperature of the geothermal fluid was estimated, was considered as the most appropriate but addressable only for analyzed cases. Therefore, mathematically speaking, what has been proposed here is not a "global" optimization scheme, but it is an assessment of the effective variants leading to higher performance for the studied cases, using numerical analyses.
However, limitation of the distance between production and injection wells affects the costs of geothermal doublet by reducing the length of pipeline between wells [83] and therefore various types of optimization have been developed. Apart from the proposed direct "optimization" method, we also distinguish other techniques of searching for an optimal solution, among them: surrogate models [84], Bayesian functions examining the probability of occurrence of various conditions [85], optimizations taking into account the process economy [86], or discrete network modeling [87]. However, it is still important to take proper account of the quantity and quality of the energy produced, which has been described, among others, in the works: [88][89][90] in which the quality of the energy produced has been optimized. Another aspect is the optimization of the location of boreholes [91,92] and the deposit's lifetime [93,94]. Recently, articles have been published which consider many different criteria of optimization and put different weightings on every standard due to different types of conditions [95][96][97].

Results of Numerical Simulation
At the first stage of the simulations, the results obtained for various polygonal grids were compared. The polygonal grid comprises a set of polygons (grid blocks) of various dimensions ( Figure 5).
These blocks were densified around the wells, which provided more accurate results of modeling in the rock volume influenced by production-injection cycles. The impact of changes in the density of grid nodes clustered around the wells for values ( Figure 5): 1 m 2 , 10 m 2 , 100 m 2 , 1000 m 2 , 10,000 m 2 , was analyzed.
Simulations were carried out for each grid density variant taking into account the following parameters: designed discharge of production well 120 m 3 /h, temperature of injected water 25 • C, 1000 m distance between wells and 50-year-long lifetime of the heating plant. As a result, changes in time of pressure and temperatures values in production and injection wells ( Figure 6) were observed.
Based on results in Figure 6, it can be observed that for density from 1 to 100 m 2 estimated values are similar. The most visible differences are visible for the density between 1000-10,000 m 2 . Furthermore, taking into consideration the simulation time, which increases with increasing density, for 1 m 2 time is approximately 4 times longer than for 100 m 2 density. For further calculations, the 100 m 2 density was assumed as the most effective for the Choszczno region. Presented and preceding published [98] tests allowed to choose the numerical grid to ensure that further refinement insignificantly influence the temperatures and pressures variation in time. Since then, several simulations have been carried out, which affected the final density of the grid. Therefore, computational results for 100 m 2 density can be treated as mesh independent.
At the second stage of the simulation, the distance between production and injection wells in ranges of: 3000, 2500, 2000, 1500, 1000, 500 m (Figure 7), was analyzed. Pressure in the production Energies 2020, 13, 2174 9 of 20 borehole zone changes from 13.749 MPa (for 3 000 m) to 13.762 MPa (for 500 m), while in the injection well pressure increased from 14.001 MPa (for 500 m) to 15.018 MPa (for 3000 m). The temperature of the injected water for a period of 50 years in the distance range of 1000 to 3000 m between wells is stable, only for the distance range of 500 m the temperature in production well zone drops up to 1.6 • C (Figure 7). The obtained results were compared with the results of other analyses carried out for the research area in recent years [20,21,36]. The results were also referred to the parameters of geothermal heating plants located in Stargard and Pyrzyce, in a short distance from the research area. In both cases it was confirmed that they are correct. These blocks were densified around the wells, which provided more accurate results of modeling in the rock volume influenced by production-injection cycles. The impact of changes in the density of grid nodes clustered around the wells for values ( Figure 5): 1 m 2 , 10 m 2 , 100 m 2 , 1000 m 2 , 10,000 m 2 , was analyzed.
Simulations were carried out for each grid density variant taking into account the following parameters: designed discharge of production well 120 m 3 /h, temperature of injected water 25 o C, 1000 m distance between wells and 50-year-long lifetime of the heating plant. As a result, changes in time of pressure and temperatures values in production and injection wells ( Figure 6) were observed. These blocks were densified around the wells, which provided more accurate results of modeling in the rock volume influenced by production-injection cycles. The impact of changes in the density of grid nodes clustered around the wells for values ( Figure 5): 1 m 2 , 10 m 2 , 100 m 2 , 1000 m 2 , 10,000 m 2 , was analyzed.
Simulations were carried out for each grid density variant taking into account the following parameters: designed discharge of production well 120 m 3 /h, temperature of injected water 25 o C, 1000 m distance between wells and 50-year-long lifetime of the heating plant. As a result, changes in time of pressure and temperatures values in production and injection wells ( Figure 6) were observed. Based on results in Figure 6, it can be observed that for density from 1 to 100 m 2 estimated values are similar. The most visible differences are visible for the density between 1000-10,000 m 2 . Furthermore, taking into consideration the simulation time, which increases with increasing density, for 1 m 2 time is approximately 4 times longer than for 100 m 2 density. For further calculations, the 100 m 2 density was assumed as the most effective for the Choszczno region. Presented and preceding published [98] tests allowed to choose the numerical grid to ensure that further refinement insignificantly influence the temperatures and pressures variation in time. Since then, several simulations have been carried out, which affected the final density of the grid. Therefore, computational results for 100 m 2 density can be treated as mesh independent.
At the second stage of the simulation, the distance between production and injection wells in ranges of: 3000, 2500, 2000, 1500, 1000, 500 m (Figure 7), was analyzed. Pressure in the production borehole zone changes from 13.749 MPa (for 3 000 m) to 13.762 MPa (for 500 m), while in the injection well pressure increased from 14.001 MPa (for 500 m) to 15.018 MPa (for 3000 m). The temperature of the injected water for a period of 50 years in the distance range of 1000 to 3000 m between wells is stable, only for the distance range of 500 m the temperature in production well zone drops up to 1.6 °C (Figure 7). The obtained results were compared with the results of other analyses carried out for the research area in recent years [20,21,36]. The results were also referred to the parameters of geothermal heating plants located in Stargard and Pyrzyce, in a short distance from the research area. In both cases it was confirmed that they are correct.

(c) (d)
Grid density: Figure 6. Changes of the pressure (a) in production and (b) in injection well zones, and temperature (c) in production and (d) in injection well zones for various densities of grid calculation.
Based on results in Figure 6, it can be observed that for density from 1 to 100 m 2 estimated values are similar. The most visible differences are visible for the density between 1000-10,000 m 2 . Furthermore, taking into consideration the simulation time, which increases with increasing density, for 1 m 2 time is approximately 4 times longer than for 100 m 2 density. For further calculations, the 100 m 2 density was assumed as the most effective for the Choszczno region. Presented and preceding published [98] tests allowed to choose the numerical grid to ensure that further refinement insignificantly influence the temperatures and pressures variation in time. Since then, several simulations have been carried out, which affected the final density of the grid. Therefore, computational results for 100 m 2 density can be treated as mesh independent.
At the second stage of the simulation, the distance between production and injection wells in ranges of: 3000, 2500, 2000, 1500, 1000, 500 m (Figure 7), was analyzed. Pressure in the production borehole zone changes from 13.749 MPa (for 3 000 m) to 13.762 MPa (for 500 m), while in the injection well pressure increased from 14.001 MPa (for 500 m) to 15.018 MPa (for 3000 m). The temperature of the injected water for a period of 50 years in the distance range of 1000 to 3000 m between wells is stable, only for the distance range of 500 m the temperature in production well zone drops up to 1.6 °C (Figure 7). The obtained results were compared with the results of other analyses carried out for the research area in recent years [20,21,36]. The results were also referred to the parameters of geothermal heating plants located in Stargard and Pyrzyce, in a short distance from the research area. In both cases it was confirmed that they are correct.

(c) (d)
The well's distance: Figure 7. Variants of spacing between wells (a) in production and (b) in injection well zones, and temperature in (c) production and (d) injection well zones for var i ous well distance.

Thermal Power Estimations
Based on the results of numerical simulations, the thermal power for the geothermal doublet versus grid density ( Figure 8) and versus distance between production and injection wells (Figure 9), were estimated. However, the wellhead temperature doesn't differ very much from the influence on thermal power as to be noticeable. The thermal power of the geothermal doublet was calculated according to the formula:

Thermal Power Estimations
Based on the results of numerical simulations, the thermal power for the geothermal doublet versus grid density ( Figure 8) and versus distance between production and injection wells (Figure 9), were estimated. However, the wellhead temperature doesn't differ very much from the influence on thermal power as to be noticeable. The thermal power of the geothermal doublet was calculated according to the formula: where: P-thermal power of geothermal doublet (W), The brine outflow . V b and reference temperature t 0 were assumed at 120 m 3 /h and 25 • C. Brine temperature t b was estimated by numerical modeling. Thermal properties of the brine c b and ρ b were estimated based on the [99] assumed salinity as 120 g/L and mean temperature between t b and t 0 at the level of approximal 40 • C. The rest of the data was determined on the basis of a commonly accepted methodology [13,[100][101][102]. The brine density was estimated at 1090 kg/m 3 and specific heat at 3.65 (kJ/(kg K)). For specific variability of thermal power, it was possible to determine the amount of energy produced in the analyzed time period (i.e., 50 years). The amount of thermal energy generated by the installation and coming from geothermal energy was determined by the relationship: where: Q-thermal energy (J), P-thermal power of geothermal doublet (W), τ-time (s).
Based on energy, it is possible to determine the average annual heat output of the installation over a 50-year period: where: P a -annual average thermal power of geothermal doublet (W), Q-energy received during time period (J), ∆τ-analysed time period (s).
In addition, based on energy received for 50 years (Q), it can be estimated the average amount of energy produced annually, as follows: where:   The amount of average energy produced annually and average annual thermal power for the analyzed variants are listed in the Table 2. Figure 8 shows the results of numerical simulations of the thermal power for the geothermal doublet versus grid density in order to check mesh independence of resolution. However, the wellhead temperature does not differ very much in Figure 6., in which the influence on thermal power has been considered. In addition, Table 2 presents the average amount of energy produced annually (TJ/yr) and averaged power of the geothermal doublet to highlight that slight fluctuation of those parameters is observed.  The amount of average energy produced annually and average annual thermal power for the analyzed variants are listed in the Table 2. Figure 8 shows the results of numerical simulations of the thermal power for the geothermal doublet versus grid density in order to check mesh independence of resolution. However, the wellhead temperature does not differ very much in Figure 6., in which the influence on thermal power has been considered. In addition, Table 2 presents the average amount of energy produced annually (TJ/yr) and averaged power of the geothermal doublet to highlight that slight fluctuation of those parameters is observed.

Discussion
Conducted simulations confirmed the possibilities of geothermal energy utilization in the Choszczno Region, which was indicated as a prospective area in terms of geothermal energy utilization [20,21,36]. Taking into account the results of grid density simulations ( Figure 6), it is visible that 100 m 2 density ensures sufficient precision. The differences of values for estimated parameters between 1 and 100 m 2 density were not significant. The estimated pressure values for the production well zone vary in the range of 13.75-13.78 MPa, which gives changes up to 0.22%. In the case of the injection well, the observed changes are higher, in the range of 14. 15-14.20 MPa, which gives a difference up to 0.35%. Pressure values in both cases stabilize between first and second year of exploitation. The values of the estimated temperatures for the production well zone ranged between 52.55-53.20 • C (change at the level of 1.2%). Larger differences in the estimated values of temperature were observed for the injection borehole during the first year of doublet operation. After approximately 2 years of system operation, the estimated temperature in extreme cases ranged between 32-45 • C, which generates a large calculation error (approximately 30%). However, this error relates to the largest of the computational grid cells, which were considered in the well zones (1000 m 2 and 10,000 m 2 ).
For changes in the distance between the boreholes (Figure 7), in the case of pressure changes in the production borehole zone, the increase from 13.749 MPa (for 3000 m) to 13.762 MPa (for 500 m), was observed (change at the level of 0.1%). The opposite situation was observed in the injection borehole, where the pressure increased from 14.001 MPa (for 500 m) to 15.018 MPa (for 3000 m) (change at the level of 7%). It should be emphasized that the depth of the injection borehole also changes with the distance. Pressure stabilization is visible between first and second years of system operation. Taking into account the obtained results for the assumed yield (120 m 3 /h) and temperature of the injected waters (25 • C) for a period of 50 years in the distance range of 1000 to 3000 m between borehole, it is possible to exploit geothermal waters without the risk of cooling groundwater horizon. In the case of 500 m distance the temperature in production well zone drop up to 1.6 • C; however, it is still not large enough to cause significant cooling of the water horizons in the assumed time of exploitation.
It should be mentioned that similar studies were carried out in the Netherlands with a brine outflow 200 m 3 /h and 400 m 3 /h and with a distance between the boreholes of 800 m and 1600 m. However, there, due to better geological conditions, the temperature (for 200 m 3 /h) was kept at 60 • C for a longer time. In addition, the period of numerical simulation corresponded to 200 years and included from one double well to several working in the system [79]. A similar scope of the double well study was also conducted in Italy [103], but for use in heat pumps. In the article [103] the authors focused on a smaller brine outflow (about 85 m 3 /h) and a shallower double well, namely 125 m. The boreholes were also close to each other, only 25 m. In the work [103], the author referred to his own experiment, from which similar data as in this work were obtained, i.e., soil parameters, brine outflow and temperature at the inlet and outlet of the boreholes. The numerical simulations covered Energies 2020, 13, 2174 14 of 20 60 years and only during the first 10 years slight temperature fluctuations were observed, which is consistent with other work [104] on low-temperature geothermal energy sources.
In the case of power generation based on the obtained results only for the distance of 500 m, thermal power is reduced after the 10th year of exploitation ( Figure 9). For other variants (when the distance equals at least 1 km), it might be noticed that the distance between wells does not influence the thermal power of the doublet and the received amount of energy ( Figure 9 and Table 2).
The averaged power of the geothermal doublet is estimated at 3.7 MW. This is the typical heat capacity of installations that are built in Germany, for example. The thermal energy also corresponds to the capacity of existing units of installations that are built in Germany, for example. The calculated thermal energy produced annually also corresponds to the energy amount received from existing units [105]. Also, the expected service life of the installation and the pressure drops resulting from its operation are at a comparable level [106,107].
The influence of grid density is not noticeable (Figure 8). Higher computational grid density does not result in higher values of thermal power and energy production ( Figure 8 and Table 2). The differences are less than 1%.

Conclusions
The results of the research work proves that modeling of the geothermal doublet can be effectively improved by selecting the maximum cell surface of the grid density, which can give precise results (errors at the level of several percent without significant impact on the final results) and will evaluate the computational time necessary to perform the simulation. Also, the influence on the distance between production and injection well was described. Both results were referred to the potential amount of geothermal energy of the analyzed geothermal doublet, which is crucial for a potential investor.
The changes of grid density calculation affected mostly the estimated value of pressure in both production and injection well zones. The estimated values of pressure and temperature obtained for the lowest grid density (between 1-100 m 2 ) were rather low, which enabled us to conclude that 100 m 2 grid density calculation is the most effective for credible results of modeling for the studied case. On the contrary, differences in pressure and temperature values calculated for two high-density grids (1000 and 10,000 m 2 ), were higher and can generate a significant measurement error. In any analyzed variant of spacing between the wells in range 1000-3000 m, the temperature drop in production well below the initial value was not observed. Hence, the effect named "the breakthrough of the cold water front" did not occur, which enabled us to conclude that in the analyzed range of spacing between the boreholes, the required level of safety is maintained from the point of view of energy transfer-in the analyzed time for a particular model. The only drop of temperature after 50 years of exploitation was observed for the distance of 500 m, but reached only 1.6 • C, which does not yet indicate that the aquifer is cooling down. However, a visible decrease indicates that further shortening of the distance may lead to a situation in which the decrease in temperature in the production well zone will be observed. For power generation, the changes in the distance between production and injection wells showed effect only for the shortest of the analyzed distances (500 m). For the rest of the distances between boreholes (equal at least to one kilometer), the results obtained during the simulation do not influence thermal power. In the case of density changes, there was no visible effect (difference at the level of one percent) on thermal power and energy production.
The presented procedure involving a multi-variant simulation of the model's grid density can contribute to shrinking the simulation time. In the case of simulations of different distances between wells, it is possible to reduce investment costs associated with the construction of surface infrastructure, which can influence on reducing the cost of the entire investment.
It should be emphasized that all carried out analyzes refer to a research area. Due to the high variability of geological conditions, it is difficult to determine whether in other areas the analysis would give similar results. It is likely that the results for the computational grid density may be refined in other cases. However, in the case of distances between wells, they always should be analyzed individually for each considered region due to the high costs associated with the increase in distance between wells. The obtained results show the significant importance of numerical modeling aimed at increasing the efficiency of a potential geothermal heating plant. The possibility of choosing the best location of boreholes is important both before and during the exploration and can significantly increase the efficiency of the entire system.