A numerical analysis of the thermohydraulics of an EGS project in Turkey

A numerical study of the thermohydraulics of an enhanced geothermal system project in Turkey is presented. The solid structures are modelled as porous media, using the numerically determined hydraulic fracturing data of other authors. The influence of several numerical modelling aspects such as the domain size, grid resolution, temporal resolution as well as the discretization scheme are investigated and assessed to obtain highly accurate numerical solutions under the applied modelling assumptions. Using the suggested mathematical and numerical model, different production scenarios are investigated.


Introduction
Geothermal energy is becoming increasingly important in the energy research due to its large capacity and independence on weather conditions, unlike the solar and wind energy. Turkey is positioned as the seventh richest country in the world in its potential on geothermal energy. The estimated total geothermal energy potential of Turkey is about 60 GW [1]. Currently, there are approximately forty installed geothermal power plants in Turkey, and their total installed capacity is approximately 1100 MWe. The existing, conventional geothermal energy systems are based on natural cracks and reservoirs that are water-dominated. In cases, where the natural cracks and pores do not permit an economic operation, the permeability can be improved by hydraulic fracturing of the rock structures. Such systems are designated as Enhanced Geothermal Systems (EGS). Indeed, applying hydraulic fracturing, i.e. by EGS, the geothermal energy can be gained without a need for a natural water reserve. In Turkey, in a region near Dikili of the Izmir province, with considerable geothermal properties such as a geothermal gradient of approximately 7K/100m, an EGS project is planned. Hou et al. [2] published an investigation on several aspects of this project, including geology and geochemistry.
In the present analysis, a further numerical investigation of the thermohydraulics of this EGS project is presented, based on the boundary conditions and the numerically obtained hydraulic fracturing data published by Hou et al. [2]. The presently applied modelling is similar to that used by Hou et al. [2], which was performed by the software FLAC3D [3], which is a code specialized on geotechnical analysis. However, in [2] essentially no details on the mathematical and numerical modelling aspect were provided. In the present analysis, modelling details are provided and discussed. For guaranteeing high numerical accuracy, the influences of the domain size, grid resolution, temporal resolution as well as the discretization schemes on the numerical prediction are investigated. Thus, an accurate numerical model of the problem within the framework of the applied mathematical model is reached. The developed modeling strategy is, then, used to predict different production scenarios.

The considered EGS configuration
A triplet system with a single injection well, and two production wells are investigated in the current analysis. The configuration of the wells with the considered hexagonal region of the rock structures with its depth, thickness and width are is in Figure 1.
The considered geometrical arrangement has two symmetry planes. These symmetry planes are considered in the modelling. Thus, a ¼ th of the region shown in Figure 1 is defined to be the solution domain. The latter, which is based on the definitions provided in [2] is illustrated in Figure 2, where the different rock zones, dimensions (numbers indicate length in meters) as well as boundary conditions are also indicated.
The considered block of the rocks starts at a depth of 1710 m measured from the ground surface and reaches a depth of 3010 m. Each of the zones denoted as "upper" and "lower" has a depth of 150 m. The fractured zone, which is designated as "frac" is positioned between the "upper" and "lower" zones and has the dimensions of 1200x33x1000 (all dimensions in meters). The solid structure around the fractured zone is designated as "rock". The depths of the injection and extraction ports, measured from the upper border of the "frac" zone, are assumed to be the same and denoted by H. The separation of the injection and production wells is denoted by L. The extensions of the "rock" region in x and y directions are indicated by X and Y.

An outline of the modelling
The computational investigation is conducted within the framework of the general purpose Finite Volume Method [4] based Computational Fluid Dynamics software ANSYS Fluent 18.0 [5]. The rock regions are modelled as porous regions. Since the time scales of the process is quite long, thermal equilibrium between fluid and solid zones is assumed, which means that a single energy equation is solved [4,5] accounting for the heat transfer by conduction and convection, while the mechanical energies are neglected. Apparently, the radiative heat transfer [6] does not play a role in the present problem.
The introduced fluid is water at moderate temperatures. The water pressures that are encountered are expectedly quite high, also due to the hydrostatic accumulation of the pressure. A boiling of the injected water in the low, hot regions of the rock, followed by a condensation in the upper rock layers before, or during the production is principally possible, as it can happen in some EGS systems. However, in the present case, a study of the prevailing temperatures and pressures together with the properties of water has implied that it is unlikely that the liquid water would undergoes a phase change. Thus, the fluid considered in the system is modelled to be the liquid, incompressible water. For the fluid, i.e. liquid water, and solids, constant material properties are assumed, corresponding approximately to the prevailing thermodynamic conditions, on the average. The material properties of solids are assumed to be the same for all rock zones.
Turbulence modelling in porous media is a challenge [7]. In the current formulation, the turbulence effects are omitted as a first approximation. This neglect may at least partially be reasoned by the expectedly quite low superficial, and also physical, velocities in most parts of the solution domain. Thus, unlike the majority of engineering heat transfer applications [8][9][10][11][12][13][14][15][16][17], no turbulence model is applied. The influence of the porous solid region on the hydrodynamics of the fluid is considered by a sink of momentum in the Navier-Stokes equations that corresponds to a drop of static pressure according to law of Darcy [18]. The inertial resistance is omitted based on the above arguments and lack of data.
For the parameters that describe the rock zones, the hydraulic fracturing data provided by Hu et al. [2] is used as guidance in the current investigation. The currently applied rock zones' porosities and permeabilities (perm) [2] are provided in Table 1 (permeabilities in m 2 ). As one can see in the table, the porosities and permeabilities are, in general, quite small. Different from the remaining rock regions, the permeability of the fractured region is non-isotropic (Table 1). In the fractured region, the permeability is several orders of magnitude greater in the main plane of fracturing, i.e. in the x-z plane of Fig. 2), in comparison to the perpendicular direction (y, Fig.2). The four vertical boundaries of the solution domain (Fig. 2) are prescribed to be symmetry boundaries. The two horizontal boundaries are defined as isothermal walls as indicated in Figure 2 (TU: upper wall temp., TL: lower wall temp.).
According to the available information [2], the determined upper and lower wall temperatures have been: TU=137°C, TL=240°C. In the current, time dependent problem, initial conditions are also needed, together with the boundary conditions. As the initial condition, a stationary fluid with a linearly varying temperature distribution between the upper (TU) and lower (TL) isothermal planes (Fig. 2) is assumed.
Between the diameters of the boreholes and the size of the solution domain, there is a very large difference of scales, which is very difficult to resolve by a computational mesh, if a modelling of the pipes coupled with the rest of the solution domain is aimed at.
Therefore, currently, it is found convenient not to calculate the flow in the pipes and the heat transfer between the fluid in the pipe and the solid structures. This may be considered as a reasonable approach for the present purpose, since the diameters of the pipes are very small (approximately 0.2 m) in comparison to the size of the domain and the pipes are usually coated by cement with a low thermal conductivity. Consequently, the modelled outer walls of the pipes are assumed to be adiabatic walls, and the inlet and outlet boundaries of the solution domain are positioned at the location (Fig. 2) of the outlet surface of the injection pipe, and the inlet surface of the production pipe, respectively.
Moreover, a fine resolution of the processes in the near-fields of the inlet and exit planes of the pipes is not attempted, either. The prescribed inlet and outlet boundaries are planes with much greater dimensions (20 m x 20 m) in comparison to the cross-sectional areas of the real pipes. Since the low porosity and permeability (Table 1) tend to produce a fast dispersion of the flow, the currently adopted modelling may be considered to be convenient for the present purpose. At the inlet plane, a homogeneous inlet velocity normal to the plane, and a constant inlet temperature are prescribed as boundary conditions, while a uniform static pressure and the zerogradient condition for the velocities are applied at the outlet. The inlet temperature is prescribed to be 60°C.
In the numerical formulation, for the time discretization, a backward differencing scheme [4,5] that is second order accurate is mainly used. A first order scheme is also applied, for comparison. For the spatial discretization, the formally third order accurate QUICK [19] discretization scheme is mainly used for the convective terms. For the purpose of comparison, the second order upwind [20], power law [21], and first order upwind schemes [21] are also applied for the convective terms. The velocity-pressure coupling is treated by the PISO scheme [22].
Within the iterative procedure of a time step, the components of the velocity vector are under-relaxed by an under-relaxation coefficient of 0.5. No underrelaxation is applied to pressure and temperature. As the convergence criteria within each time step, the threshold values for the normalized residuals are fixed to 10 -12 for the energy equation and to 10 -6 for the continuity and momentum equations. The applied backward time differencing scheme is, essentially, unconditionally stable. Nevertheless, the time-step sizes (Δt) are prescribed as such that the cell Courant numbers (Co) [4] remain below unity, as it would be required for a convection problem by an explicit scheme.

Numerical accuracy studies 4.1 Domain size
In the current time-dependent problem, the temperature front, which propagates in time, shall not interact with the spatial boundaries, i.e. be falsified by the prescribed boundary conditions, within the considered time domain. It has been aimed to analyze a time period of twenty years. The vertical boundary locations with the prescribed temperatures (Fig, 2), were acquired form [2]. Based on this definition [2] it can be assumed that the locations of the boundaries, each with a distance of 150 m from the fractured region allow sufficient room. This can similarly be assumed for the dimensions of the domain in the remaining directions (X, Y, Fig. 2).
Nevertheless, a course estimation was carried out, based on the analytical solution of a one-dimensional, transient heat conduction problem in a semi-infinite domain [23]. This implied a penetration depth of approximately 90 m, supporting the convenience of the assumed distance of 150 m. Initial numerical predictions were also carried out with varying X and Y (Fig. 2), which also supported the convenience of the abovementioned dimensions. Therefore, for the main calculations the domain size is prescribed by the dimensions X=150m and Y=200m (Fig. 2).

Grid independence study
The solution domain is discretized by essentially equidistant, orthogonal, block-structured, hexahedral grids. For ensuring grid independence three different grid sizes are investigated. The total number of grid nodes of these three grids are given as: The Coarse Grid: 54,000, The Medium Grid: 185,000, The Fine Grid: 410,000. The predictions are obtained for H = 590 m, for the same injection rate, using the QUICK discretization scheme, the second-order backward time differencing and the same time-step size for all three calculations.
The predicted variations of the temperature by the three grids along a horizontal line between the production and extraction ports, for a moment 90 days after the injection start are presented in Figure 3. As can be seen, the predictions of Medium and Fine grids are nearly identical. Therefore, for the Medium grid with a total number of grid nodes 185,000 is taken as the basis for further analysis.

Effect of spatial discretization schemes
The predictions are obtained applying the same grid (Medium), for H = 590 m, applying the QUICK, second order upwind, power law and first order upwind schemes, applying the second order backward differencing in time, as well as the same time-step size.
The temporal variations of the production temperature, i.e. the temperature at the exit of the solution domain, predicted for an injection volume flow rate of 250 l/s by the QUICK and power law schemes, are presented in Figure 4, for a period of two years.
The predictions of the second order upwind and the first order upwind schemes are not indicated separately in Figure 4, since they are practically identical with those of the QUICK and power law schemes.
As can be seen, the production temperature is constant for a certain period of time and starts to decay after a production time of approximately 1.5 years. One can also see that the different discretization schemes agree with each other quite well, except for the transition period around 1.5 -2 years, where the power law scheme (as well as the first order upwind scheme) predicts an earlier decline of the production temperature compared to the QUICK scheme (as well as the second order upwind scheme). In the present, main predictions, the principally more accurate QUICK scheme is used.

Effect of temporal discretization
Mainly, a second order backward time differencing is used, along with a time step (Δt) that ensures a cell Courant number (Co) not exceeding the unity. Because the current problem is dominated by the conduction heat transfer, especially in the non-fractured parts, the stability condition for sole conduction may also be checked, which requires the cell Fourier number (Fo) to remain below ½. A comparison of the occurring Co and Fo indicates that the ratio of Co to Fo is approximately 20 and, the satisfaction of the condition Co < 1 would directly satisfy Fo < ½ tenfold. For the Medium grid and a volume flow rate of 250 l/s, a time step size of Δt=6h, sufficiently fulfills the latter condition. For testing the accuracy of the temporal resolution, a computation is carried out using a time-step size, which is halved (Δt=3h). In the computations, the QUICK scheme is applied for the spatial discretization. All predictions are observed to be essentially identical, confirming the applied strategy in temporal discretization.

Results
Using the described numerical model (X=150m, Y=200m, QUICK and Power Law schemes, second order backward time discretization, Δt ensuring Co < 1), the case with H = 590 m, L = 1000 m is analyzed for a period of 20 years, for a volume flow rate of 250 l/s. The contours of predicted temperature contours in plane y-z through the injection wells (Fig. 2) 20 years after production are presented in Figure 5, as predicted by the QUICK and Power Law schemes.
In this plane, with low permeability the conductive heat transfer is dominating. One can see that the temperature front has a rather low penetration in the y direction and remains quite confined to the fractured zone, as estimated before (Section 4.1). It is interesting to note that the results obtained by the QUICK (Fig. 5a) and Power Law (Fig. 5b) scheme for the temperature distribution in this plane, after 20 years are very similar.
The predicted temperature contours in plane x-z through the wells (Fig. 2) 20 years after production are presented in Figure 6. One can see that although the low temperature front has reached the production well in 20 years, the production temperature is still above the injection temperature, due to the mixing of the warmer water, which is heated up by the hot rocks that surround the fractured zone and convectively transported to the production well by a recirculating flow behind the production well (Fig. 6). Comparing the predictions provided by the QUICK (Fig. 6a) and Power Law (Fig.  6b) schemes, once can see that both discretization schemes deliver very similar predictions for the temperature distribution in this plane, after 20 years of production. However, some differences can still be observed, especially in the vicinity of the production well, where it can be observed that the QUICK scheme (Fig. 6a) predicts a slightly higher temperature level for the recirculating flow into the production port, compared to the Power Law scheme (Fig. 6b).
The influence of the injection rate is investigated by considering four values of the water volume flow rates, namely 250 l/s, 200 l/s, 150 l/s and 100 l/s. In all these calculations, the QUICK scheme is used. The temporal variations of the production temperature, i.e. the static temperature at the exit of the solution domain (i.e. inlet of the production well) for a period of 20 years, predicted for injection water volume flow rates of 250 l/s, 200 l/s, 150 l/s and 100 l/s are presented in Figure 7. At the start of the process, the production temperature is only a little above 200 °C and has the same value for all cases. As a common tendency for all cases, it can be seen that the production temperature stays constant and then declines in time. The period before the moment of decline increases with decreasing water volume flow rate. For a water volume flow rate of 250 l/s and 100 l/s this moment is reached after approximately 1.5 and approximately 4 years, respectively, with values in between for the remaining water volume flow rates.
One can observe that the production temperature decays more quickly with increasing water volume flow rate, as it would be expected. Because the gained thermal energy is aimed to be used for power generation, the production temperature should be high enough. The temperature after a production period of 20 years is near and a little below 110°C for the cases with water volume flow rates 200 l/s and 250 l/s, respectively. For the cases with water volume flow rates of 150 l/s and 100 l/s, the production temperatures predicted after 20 years are substantially higher, i.e. approximately 120°C, and 140 °C, respectively. The favorable injection rate should be defined by accounting for the geothermal capacity, the temperature draw-down and the economic aspects. The current information implies that the preferred injection rate should be within the range 100 l/s -150 l/s.

Conclusions
A numerical investigation of an intended EGS project in Turkey is carried out. Special attention is devoted to the issues pertaining to the numerical accuracy. For a given configuration, the velocity, temperature and pressure fields are predicted for 20 years, for different injection rates. The predictions show that the injection rates 150 l/s and 100 l/s are adequate, since they can ensure a sufficiently high production temperature level over a period of 20 years. The investigation will be extended in the future work. This includes the investigation of the possibilities of accounting for local inertial and turbulence effects in modelling the porous media. On the other hand, alternative concepts for the configuration of the injection and production wells, as well as for intensifying the heat exchange between the working fluid and the hot rocks will be investigated.