The influence of variability of calculation grids on the results of numerical modeling of geothermal doublets - an example from the Choszczno area, north-western Poland

The numerical modeling enables us to reduce the risk related to the selection of best localization of wells. Moreover, at the stage of production, modeling is a suitable tool for optimization of well operational parameters, which guarantees the long life of doublets. The thorough selection of software together with relevant methodology applied to generation of numerical models significantly improve the quality of obtained results. In the following paper, we discuss the impact of density of calculation grid on the results of geothermal doublet simulation with the TOUGH2 code, which applies the finite-difference method. The study area is located between the Szczecin Trough and the Fore-sudetic Monocline, where the Choszczno IG-1 well has been completed. Our research was divided into the two stages. At the first stage, we examined the changes of density of polygon calculation grids used in computations of operational parameters of geothermal doublets. At the second stage, we analyzed the influence of distance between the production and the injection wells on variability in time of operational parameters. The results demonstrated that in both studied cases, the largest differences occurred in pressures measured in production and injection wells whereas the differences in temperatures were less pronounced.


Introduction
The numerical modeling is a common tool applied by various branches of science including the geothermics. Recently, modeling plays a crucial role even at the early stages of planning and designing of geothermal investments [1]. The numerical models enable us to recognize both the geological structure of rock formations and the processes which operate within them [2; 3, 4; 5; 6; 7]. This, in turn, facilitates the prediction of possible unwanted processes accompanying the exploitation of groundwater, as e.g., the clogging of wells [8] or corrosion of casings and surface installations. Such obstacles cause serious problems in many geothermal investments in Poland, particularly in those located in the Polish Lowlands due to (among others) high Total Dissolved Solids (TDS) of produced geothermal waters [9]. Moreover, numerical methods can be applied also to some related problems, e.g., mass exchange within the rock formations, which is crucial for prognosing the migration of contaminants [10]. In the case of mass and heat transfer, two methods are mostly applied in numerical modeling: finite element method and finite difference method [11]. Both methods are widely used by specialized simulators. Such tools facilitate the computations and enable us to determine the potential of explored area and/or to select optimal operational parameters of given geothermal installation in order to ensure its long life and failureless running [12;13;14]. Taking into account high costs of drillings, a reasonable idea is to run a series of simulations before any critical investment decision is made. This may eliminate many mistakes resulting from unwary solutions accepted at the initial stages of investments and may also facilitate the selection of optimal parameters of future installations, as e.g., pattern of wells (location and number), their spacing and wellbore diameters. Numerical computations always procure errors resulting, among others, from simplistic assumptions we apply in the modelings. However, the credible models must provide quantitative and qualitative results which are in agreement with the nature of processes they describe. The errors unrelated to modeling methodology originate from e.g., the quality (correctness) of geological, geophysical, hydrogeological and other data acquired in an exploration area. The modeling errors depend on selected numerical method, i.e. on the size of selected area of modeling (which cannot be either too small or too large -the latter unnecessarily extends the computation time) as well as on the selection of boundary conditions and mode of their input. Important is also the experience of staff members who run the simulations, precisely, how much they are capable of interpreting the results and evaluating their credibility. The crucial aspect of modeling is the proper selection of calculation grid, which enables us to optimize 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 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 optimize the distance between production and injection wells in order to extend the lifetime of the system but also to decrease the energy transfer loss and to reduce the costs of connection pipelines.

Conceptual geological model and hydrothermal parameters
The numerical modeling has been carried out for geological structure named the Szczecin Trough located in the northwestern Poland. More precisely, the study area was the Choszczno Anticline ( Fig.  1) in which the Choszczno IG-1 well was localized. For the purpose of modeling, it was considered as the production well of the geothermal waters. The Choszczno IG-1 well was located in the vicinity of Zam cin village, south from small Choszczno town. It was one of the wells drilled at the first stage of recognition of geological structure of the Polish Lowlands [15]. The main goal of this drilling was to identify possible stratigraphic gaps in both the Cretaceous and Jurassic successions as well as to provide data for interpretation of seismic sections completed in this area in 1957 [15].   In the study area, the most promising geothermal reservoirs occur in the Lower Jurassic succession [6; 7]. Lower Jurassic rocks occur in the whole Szczecin Trough, as alternating, lacustrine and alluvial sediments. The groundwater horizons are sandstone beds among which the best reservoir parameters were found in the The results of laboratory studies carried on at the State Geological Institute for samples from the Choszczno IG-1 well provided the mean values of effective porosity and density for particular stratigraphic units [18]. The best reservoir properties were found in both the Upper and Lower Sinemurian sediments, for which the effective porosity exceeds 20% [18]. The perspective Lower Jurassic geothermal reservoir comprises mostly white, grey and greyish-green, fragile, fine-grained sandstones accompanied by grey, compact mudstones and grey or greyish-green, sometimes fissile claystones [15]. During the drilling, the well-logging was carried on providing data for quantitative interpretation of Lower Jurassic interval.

Numerical model
The TOUGH2 numerical simulator have been used. TOUGH2 solves mass and energy balance equations that describe fluid and heat flow in multiphase, multicomponent systems [19]. Fluid advection is described with a multiphase extension of Darcy's law; in addition there is diffusive mass transport in all phases. Heat flow occurs by conduction and convection, the latter including sensible as well as latent heat effects. The simulator base on the integral finite difference method. Our numerical model covered the area 42.6 x 29.3 km (Fig. 3A). The geothermal doublet was based upon the archival Choszczno IG-1 well, which was defined in the model as the production well. The hypothetical injection well was located about 2,000 m apart. After implementation of both wells into the model, the polygonal calculation grid was superimposed and calibration was run again. The production well was designed to 1,480 m depth in order to reach the groundwater horizons located in the Radowo Beds. The injection well was also designed to the same 1,480 m depth, which enabled us to generate the closed groundwater circulation system within the Lower Jurassic geothermal reservoir. The production   The regional-scale model (Fig. 3A) was verified by comparison of the results of computations with the recorded values (Fig. 4). The obtained model was recognized as satisfactory representation of the reality. However, in order to save the time of computations, the local-scale model was used in further modeling (Fig. 3B). In this model, we rejected the calculation elements located outside the potential influence zone of operating geothermal doublet, thus, the size of study area was reduced to 10 x 12 km. At the margins of that area the first type (Dirichlet) boundary conditions were applied to determine pressure and temperature values obtained for each node of calculation grid of the regionalscale model.

Polygonal calculation grid
The software widely used in calculations of mass and heat transfer includes many various types of calculation grids. Also the TOUGH2 code provides regular, polygonal and concentric (radially symmetric) grids [3]. Relevant grid is selected depending on the studied problem together with pertinent density of calculation grid. The proper selection of calculation grid influences both the correctness and the accuracy of generated models, and the time of computations. At the first stage of modeling we compared the results obtained for various polygonal grids, which are widely used for evaluation of well operation including the geothermal doublets. The polygonal grid comprises a set of polygons (grid blocks) of various dimensions. These blocks are densified around the wells, which provides more accurate results of modeling in the rock volume influenced by productioninjection cycles -an information deadly important from the point of view of the investors. At the first stage, we tested the impact of changes in density of grid nodes clustered around the wells. We considered four variants expressed by the areas of grid blocks (Fig. 5): 2,500 m 2 , 5,000 m 2 , 10,000 m 2 and 100,000 m 2 . We run the modelings 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 and 50-years-long lifetime of the plant. As a result, we observed changes in time of pressures and temperatures in production and injection wells (Fig. 6). Our results demonstrated (Fig. 6) that 5,000 m 2 density ensures sufficient precision because the differences of values of modeled parameters obtained for 5,000 and 2,500 m 2 grids appeared to be unimportant.  Figure 6. Changes of pressure and temperature in wells of geothermal doublet at various densities of calculation grid.

Changes of spacing of production and injection wells
The crucial factor in optimization of geothermal doublet operation is the proper spacing of production and injection wells. Too short distance results in too fast cooling of reservoir formation between the wells, which decreases the temperature and, thus, the thermal power capacity of geothermal system. On the contrary, too long distance between the wells generates excessive heat loss and requires more powerful pumps. High urbanization of the area also influences the localization of wells and limits the opportunities of new investments. At the second stage of modeling, we tested the influence of various spacing between production and injection wells (Fig. 7): 3,000 m, 2,000 m, 1,000 m and 500 m. In this modeling, we applied calculation grid density 5,000 m 2 and the same operational parameters of doublet as for modeling of grid densities. The results of modeling were represented by changes in time of pressures and temperatures in production and injection wells (Fig. 8).  Figure 8. Changes of pressure and temperature in wells of geothermal doublet at various spacing of production and injection wells.

Conclusions
Our paper aimed to turn attention of the Readers to the problem of credibility of numerical modeling. We also attempted to evaluate quantitatively the discrepancies between the results of modeling of the study area (case study). Generally, our studies demonstrated the importance of density of calculation grid and spacing between the wells of geothermal doublet for the results of numerical modeling.
The changes of density of calculation grid affected mostly the calculated pressures in both the production and injection wells referred to the initial pressure. Precisely, pressures obtained for the lowest grid density (100,000 m 2 ) differed significantly from those obtained for higher grid densities. On the contrary, differences in pressures calculated for two high-density grids (2,500 and 5,000 m 2 ) were rather low, which enabled us to conclude that 5,000 m 2 density of calculation grid is sufficient for credible results of modeling. The remarkable differences of temperatures appeared after about one year. For two high-density grids, modeled temperatures were similar although higher differences were found in the injection well. Finally, our results demonstrated that minimum density of calculation grid in the zone around the well is 5,000 m 2 .
In any analyzed variant of spacing between the wells, we did not observe abrupt temperature drop in production well below the initial value. Hence, the effect named "the breakthrough of the cold water front" did not occur, which enabled us to conclude that all the modeled spacing variants are safe from the point of view of energy transfer in the analyzed time span and in particular model.