- Email: [email protected]
A case study of the modeling of a hydrothermal reservoir: Khankala deposit of geothermal waters
Geothermics 59 (2016) 56–66
Contents lists available at ScienceDirect
Geothermics journal homepage: www.elsevier.com/locate/geothermics
A case study of the modeling of a hydrothermal reservoir: Khankala deposit of geothermal waters A. Farkhutdinov a,b,∗ , P. Goblet a , C. de Fouquet a , S. Cherkasov c a b c
Mines ParisTech, 35 rue Saint Honoré, 77305 Fontainebleau, France Bashkir State University, 32 Zaki Validi Str., 450074 Ufa, Russia Vernadsky State Geological Museum, 11 Mokchovaya Str., 125009 Moscow, Russia
a r t i c l e
i n f o
Article history: Received 24 October 2014 Received in revised form 17 August 2015 Accepted 12 October 2015 Keywords: Geothermal waters Khankala deposit Kriging Numerical simulation
a b s t r a c t In 2013 the Grozny State Oil Technical University together with members of the “Geothermal resources” consortium started a pilot project to build a geothermal plant on the basis of the most promising XIII layer of the Khankala thermal waters deposit. Planned capacity of the projected facility is 22.8 GJ/h, which will heat a greenhouse complex. Geostatistical analysis and numerical simulations are required for achieving sustainability in exploitation of the resource and in order to define new drilling sites. In this paper intrinsic random functions of k order kriging (IRF-k) is used for the XIII layer map creation and estimation of the temperature distribution within it. The results obtained by geostatistical approach are used for numerical simulation of the hydro-thermal process relied to reinjection of cold thermal water back into the productive layer. The work carried out allows highlighting the most promising zones for exploitation and shows that the geothermal water temperature in the production well starts to decrease after 6 years. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Geostatistical analysis and simulation is now widely applied in geosciences in order to study complex phenomena. In geothermy, initial data, as a rule, are obtained from irregularly located wells, which complicate interpolation and extrapolation of different parameters in order to model the behavior of a geothermal reservoir. Thus, a geostatistical approach is useful for mapping the reservoir itself, as well as describing permeability, porosity, salinity, transmissivity, temperature, etc., A major task for geostatistics here is to improve the reliability of the mapping of the reservoir’s geometry and the spatial distribution of its characteristics. These descriptions are extremely important for the design of the wells and their further exploitation. In the current paper, geostatistical analysis is applied to the existing data on the Khankala deposit of geothermal waters (Russia) in order to define new drilling sites. The Khankala deposit is located 10 km to the south–east from the city of Grozny (Fig.1). The deposit has a layered structure and contains multiple layers of Chocrack and Karagan horizons of mid-
∗ Corresponding author at: Mines ParisTech, 35 rue Saint Honoré, 77305 Fontainebleau, France. E-mail address: [email protected] (A. Farkhutdinov). http://dx.doi.org/10.1016/j.geothermics.2015.10.005 0375-6505/© 2015 Elsevier Ltd. All rights reserved.
dle Miocene sandstones, interlayered with clays. Tectonically, the deposit occupies the southeastern plunge of the Octyabrsk anticline bordered by two co-axial major faults. The field was actively developed from 1974 to 1994, but its exploitation was ended due to the beginning of the war in Chechnya. In 2013 the Grozny State Oil Technical University together with members of the “Geothermal resources” consortium started a pilot project to build a geothermal plant on the basis of the most promising XIII layer of the Khankala deposit. The water temperature is about 100 ◦ C. Planned capacity of the projected facility is 22.8 GJ/h, which will heat a greenhouse complex. The drilling of a “doublet” is planned (Fig. 2), i.e., one production and one injection well with reinjection of all the cooled thermal water. A numerical model is needed to predict the evolution of the resource in order to study how the production temperature would decrease as a result of the recycling of the reinjected cold water.
2. Theory and methods The aim of the work is to estimate the temperature distribution within the XIII layer of the Khankala deposit of thermal waters, develop a new structural map and make, by numerical simulation, a prognosis for deposit exploitation.
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
57
Fig. 1. General map of the location of the Khankala deposit.
2.1. Geostatistics A geostatistical analysis is used to estimate the geometry of the geothermal site. This analysis is carried out using the ISATIS program (Géovariances, 2009). Intrinsic random functions of order k kriging (IRF-k) was used because of the nonstationarity in data. The idea behind IRF-k, developed by Matheron (1973), is that it is possible to analyze the data considering some variations of the differences between data points to achieve a stationary covariance. Increments of order k are also called authorized linear combinations of order k (k = 0 indicates a difference Z (x) − Z (x + h), k = 1 is an increment of an increment, etc., Usually k is not higher than 2). The increments filter out the mean (k = 0) or “drift”, supposed to be expressed by a polynomial of order k of the coordinates. For example, in a two-dimensional space for filtering linear drift (for k ≥ 1) a0 + a1 x + a2 y, the conditions of authorization are (Chilès and Delfiner, 1999):
i
Fig. 2. “Doublet” used in geothermal water extraction-injection.
i = 0,
i
i xi = 0,
i yi = 0
i
A generalized covariance is used for the estimation. For an intrinsic random function Z(x) of order k there is a symmetrical
58
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
where |x˛ − xˇ | is the distance between 2 points in the horizontal plane. 2.2. Mathematical modeling The estimated temperature is used for the geothermal flux calculation and a structural map is used as a support for the creation of the mesh (Fig. 3). Numerical simulation was made using the Metis program, developed at the Ecole des mines de Paris (Goblet, 1980). Metis simulates liquid flow, mass and heat transport in porous or fractured aquifers. The equations are solved by the finite element method. The aquifer is assumed to respond almost instantaneously to changes in regime and the flow is supposed to be steady (Goblet, 2005a). The equation describing the steady flow is: div (k grad h) = 0
(2)
where K denotes the hydrogeological permeability and h the hydraulic head pressure. Heat transport is solved under transient conditions. It takes into account:
Fig. 3. Scheme of work with data.
- convection in the aquifer (flow of heat by water movement); - conduction in the aquifer and the surrounding formations (heat flux resulting from the temperature gradient); - kinematic thermal dispersion in the aquifer (heat flux resulting from the heterogeneity of the local velocity field). Convection is the dominant mechanism. The exchange with the surrounding formations by conduction delays progression of the cold front. The conduction and the dispersion in the aquifer have the effect of spreading the transition zone between cold water and hot water (Goblet, 2005b,c). The equation describing the heat transport in the aquifer is:
Fig. 4. Structural map of the XIII layer after georeferencing and digitization.
function K(h), called a generalized covariance, which for any pair of increments, Z() and Z() of order k is equal to:
E [Z () Z ()] = ˛ ˇ K |x˛ − xˇ | ˛ ˇ
(1)
div grad-w UD = a
∂ ∂t
(3)
where is the thermal dispersion tensor; the temperature;␥w the volumetric heat capacity of water; a the volumetric heat capacity of the aquifer (rock + water);UD the Darcy velocity. The thermal dispersion includes the pure conduction and heat dissipation kinematics according to the relationship: = 0 + w ˛|UD |
Fig. 5. Omnidirectional sample variogram of the XIII layers top elevation (left) and initial data location (right).
(4)
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
59
where 0 is the coefficient of the thermal conduction of the aquifer (rock + water);˛ the dispersivity tensor. The transport of heat in the low permeability formations is restricted to a mechanism of pure conduction:
div s grad = s
∂ ∂t
(5)
where s is the conduction coefficient; s the calorific capacity of the surrounding formations.
3. Initial data 3.1. Structural map of the XIII layer
Fig. 6. Cross-validation of model for altitude of the XIII layers top estimation.
A first step is to collect data for structural mapping of the XIII layer. A map from 1967 year is used for this purpose. It represents more than 100 wells crossing the XIII layer as the deposit is situated in the Octyabrsk anticline, which contains an oil deposit of the same name. Georeferencing was made in the Quantum GIS program (Qgis Development Team, 2014) using wells found within the Khankala deposit. Elevation at each well which was written on the map was copied and reinterpreted using geostatistical methods. There are 106 points of the top of the XIII layer after map digitization and georeferencing (Fig. 4).
Fig. 7. Estimation of the elevation (top) of the XIII layers and associated standard deviation of the estimation error (bottom).
60
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 8. Comparison of old and new structural maps of the top of the XIII layer.
Fig. 9. Reference plane choice for z coordinates.
3.2. Temperature measurements Geological observations and temperature measurements in the wells were carried out in 1968 and 1988 during the exploitation of the Khankala deposit of thermal waters (Krilov, 1988; Shpak, 1968). Measurements were carried out after temporary stoppage of pumping so as to establish a thermal equilibrium between the well and the surrounding rocks, in total about 100 measurements were made in 14 production wells, each 100–200 m in depth from the surface. As a result of the tragic events that occurred in the territory of the Chechen Republic, data on the inclination of the wells have been lost. The only information available, allowing an estimation of the in-depth position of boreholes in the area, are corrections to the inclination of the wells, which is equal to a maximum of 11 m for depths greater than 1 km, which means a relatively small
inclination. As a consequence, all wells are supposed to be vertical during the geostatistical analysis. The area of study was limited by two faults located in the south–west and north–east of the study area (Farkhutdinov et al., 2014). 3.3. Numerical simulation Conditions for the Metis numerical simulation are chosen according to the data obtained from hydrogeological tests and laboratory examinations (Ermolaev et al., 1954; Krilov, 1983,1988). Aquifer parameters are as follows: – permeability: 1.4e–12 m2 ; – longitudinal and transverse thermal dispersivity: 10 × 2 m; – volumetric heat capacity of water 4.18 MJ/m3 /K;
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
61
Fig. 10. Scatter-diagrams of the temperature versus depth in cases of the different reference plane.
– volumetric heat capacity of the aquifer 2.75 MJ/m3 /K; – volumetric heat capacity of the surrounding formations 2.09 MJ/m3 /K; – conductivity of the aquifer 1.45 W/(m K); – conductivity of the surrounding formations 1.3 W/(m K). – Two boundary conditions are entered: – injection flow of 200 m3 /h; – a constant heat flow of 10,051 MW at the injection well (corresponding to an injection temperature of 45 ◦ C). No boundary conditions were applied to the vertical sides of the model, assuming no liquid or heat flow through these limits. Geothermal flux which is imposed at the base of the model and unequally distributed (see Section 4.1.2) is calculated according to Fourier’s law, taking into account the temperature of the XIII layer of the reservoir, its depth, the temperature at the Earth’s surface in this area and the average coefficient of thermal conductivity of rocks lying between the producing formation and the surface: →
q = −k∇ T
(6)
Application of this formula gives a fairly high average heat flow (155 MW/m2 ), which can be explained by “overthrust-nappe theory” as caused by tectonic movements, considering the North Caucasus as a mobile tectonic zone (Kamaletdinov et al., 1991).
Table 1 Selection of the drift. Drift
Mean error 2
2
1 x y x xy y 1 x y x2 xy y2 x3 x2 y xy2 y3 1xy Without drift
7.116 11.54 42.41 96.11
The sample variogram confirms the known non-stationarity: the mean varies smoothly over the study area, as the fold sinks in a southeasterly direction. An IRF-k modeling is needed, indirect fitting is divided into 3 steps: a Drift finding. b Indirect fitting of the generalized covariance. c Variable estimation and kriging error variance.
IRF − k kriging system :
∗
Estimation : Z =
Z˛ 0
k˛ˇ f ˛
×
fˇ 0
t
×
˛
Estimation varience : Var( ) = K00 −
The first step is to create a structure map of the XIII layer. The omnidirectional sample variogram of the elevation of the top of the XIII layer is calculated (Fig. 5): 1 N(h) 2 [Z (xi ) − Z (xi + h)] 2N (h) i=1
(7)
where N(h) is the number of point pairs separated by a vector h.
=
K˛ˇ
f0
l
4.1.1. Structural map of the productive layer a) Variogram model
(h) =
l
4. Modeling and results 4.1. Geostatistical approach
˛
˛ l
t
×
K˛0
f0
Several possible drifts are first supposed (Table 1) and coefficients are fitted by the least squares method. A quadratic drift (local linear combination of the following functions: 1, x, y, x2 , xy and y2 ) is chosen as the most suitable. To find the optimal model of the generalized covariance, several possible models are considered (Table 2). They are listed by score, which is equal to the theoretical variance divided by experimental variance (it should be close to 1). The chosen model is spherical with a range of 1000 m and a sill of 236.2 m2 . Initial data one-by-one were “hidden” and estimated,
62
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 11. Coordinate system rotation.
• Estimation
Table 2 Selection of the generalized covariance model. Score
Nugget effect (Sill, m2 )
Spherical (Sill, m2 )
Linear (Sill, m2 )
Exponential (Sill, m2 )
1.031 1.055 1.166 1.206
0 0 0 30.8
236.2 0 0 0
0 453.2 0 0
0 0 236.5 216.6
then the difference between initial and estimated value was calculated. Unique kriging neighborhood was selected as the most appropriate, giving the best estimation by cross-validation tests with a coefficient of correlation equal to 0.99 (Fig. 6).
A structural map of the XIII layer is built using this model (Fig. 7). Faults are considered as screens for the estimation. Because there is only one well at the other side of the faults, estimations are not possible to the north and south of the anticline. A standard method of linear interpolation was used during the creation of the map in 1967, based on the measurement points. All data in the vicinity of a point has been taken into account in the new geostatistical interpretation (Fig. 8). Due to the high density of data in the most part of the map, it is not surprising that results are similar. The difference in the south–east is due to a lower density of wells, and to the introduction of a drift.
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
63
Fig. 12. Variograms calculation.
4.1.2. Initial temperature distribution The temperature distribution within the XIII layer and at the surface is supposed to be known, in order to define geothermal flux as one of the boundary conditions for the simulation. The three -dimensional structure of the temperature distribution is examined on directional sample variograms. There were two options for choice of the initial coordinate system for working with temperature measurements: a normal reference plane and the top of the XIII layer as a reference plane (Fig. 9), and both of them were tested. As we can see from scatter-diagrams (Fig. 10) variance of data is less in the case of the XIII layers top reference plane choice, especially in the vicinity of the layer itself. Before working with variograms, the coordinate system was rotated (Fig. 11) so that the anticline folds directly to the east with the depth increasing from west to east and from the axis to the limbs of the structure (i.e., in x and y directions). This facilitates the drift finding as the first step of the IRF-k modeling. 4 horizontal variograms for different directions were calculated to evaluate the benefits of choosing one or another variant for coordinate system (vertical variogram remains unchanged in both cases) (Fig. 12). Although the number of pairs is insufficient for reliable conclusions, especially with increasing distance, it is observed that the variation of temperature with the top of the XIII layer chosen as reference plane is reduced. Thus this variant is used as a selection for the coordinate system. It can be explained by the fact that the XIII layer represents anticline containing geothermal waters and the temperature follows its structure, with undergoing changes only with large differences in depth of formation. Vertical variogram is non-stationary, which corresponds to the temperature progressive increase with depth. It was the main reason for selecting the IRF-k modeling. Samples which are situated less than 300 m from the XIII layer were used and grid for block of 40 × 40 × 40 m size estimation was created. Next were steps of IRF-k modeling (see Section 4.1.1.) (Table 3).
Table 3 The XIII layers temperature estimation parameters. Drift
1yz
Variogram Neighborhood
Spherical (range = 1500 m, sill = 12.7 ◦ C2 ) Unique
Fig. 13. Cross-validation of the XIII layers temperature estimation.
Unique neighborhood was chosen as the most appropriate regarding to the small number of samples. Cross-validation was performed to evaluate chosen parameters (Fig. 13). Cross-validation based on 37 samples showed error variance equal to 6.8 ◦ C and IRF-k kriging was performed (Fig. 14). All components are now available for geothermal flux determination taking surface temperature as average for this region and equal to 24 ◦ C (Fig. 15).
64
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 14. Temperature of the XIII layer estimation (top) and error standard deviation (bottom).
Fig. 15. Geothermal flux estimation.
This geothermal flux was applied as one of the boundary condition for numerical simulation. For zones on the other sides of faults the average value was taken. 4.2. Simulation results The mesh, created for numerical simulation, geometrically represents the 3D structure of the deposit, taking into account the geometry of the XIII productive layer obtained by IRF-k estimation and consists of 141,550 prismatic (6 nodes) elements and 75,348 nodes. It is refined near production and injection wells, which are represented by one-dimensional 2-node elements (Fig. 16). The average size triangle covers an area of 2.5 km2 . Productive layer thickness is equal to 40 m, as it was mentioned as almost
constant in relevant geological reports. We do not have information on the thickness of XIII layer in all wells due to the consequences of two wars (in order to make more precise estimation). Elevation value for zones on the other sides of the faults was obtained by IRF-k estimation (see Section 4.1.2.) without taking these faults into account. It was done because of lack of certainty in shift range and real influence of faults on water flow. Simulation time is equivalent to 30 years of exploitation with a constant average regime (see paragraph initial data). Hydrogeological tests of 1981 and 1988 (Krilov, 1988) showed that the permeable/impermeable status of the north–east and south–west faults is uncertain, so both hypotheses were tested. The simulation
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
65
Fig. 16. The mesh refined near wells.
does not predict any change of temperature in production well for the first 6 years (Fig. 17). Only after this time period the cold plume expands into the production well and the temperature started to decrease (Figs. 17 and 18).
Fig. 18 shows that there would be a difference in well water temperature during production assuming faults are impermeable or not. Thus it was the correct decision to place injection and production wells almost parallel to the faults.
Fig. 17. Simulation of geothermal water reinjection in case of permeable (top) and impermeable (bottom) faults.
66
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 18. Temperature of geothermal water in production well.
There are a lot of uncertainties in temperature estimation especially for zones where the average value of flux was taken. But regarding the small number of samples it is one of the best possibilities for today which allows making initial forecasts. This work will be continued after drilling the first doublet and obtaining new, more precise data. 5. Summary and conclusions The map of the temperature distribution within the XIII layer of the Khankala deposit of geothermal waters allows, depending on the chosen site for the production well, the prediction of the water temperature in the production well, which in conjunction with a flow rate provides important information on the potential future of the geothermal power station. Numerical simulation shows that some corrections in the exploitation regime will be necessary after a certain period of time in order to make geothermal reservoir management sustainable for more than 30 years (Ungemach et al., 2011). One of the possible solutions could be Aquifer Thermal Energy Storage (Réveillère et al., 2013), which means reinjecting hot water during summer when there is not much need of heating. Acknowledgements The research has been executed with support of the Ministry of Education and Science of Russian Federation, contract no. 02.G25.31.0056. References Chilès, J.-P., Delfiner, P., 1999. Geostatistics: modeling spatial uncertainty. In: Wiley Series in Probability and Statistics. John Wiley & Sons Inc., New York, pp. 695.
Ermolaev, A., Lyalin, L., Napolsky, M., Feldman, S., 1954. Estimation of the Octyabsk deposit oil and gas reserves as of January 1, 1954 (in Russian), 338 p. Farkhutdinov, A., Farkhutdinov, I., Ismagilov, R., 2014. History of the discovery and exploitation of Khankala geothermal waters deposit (in Russian). Bull. Bashkir State Univ. 19, 93–96. Géovariances, 2009. Isatis Technical Ref., ver. 9.06. Geovariances & Ecole Des Mines De Paris: Avon Cedex, France. Goblet, P., Programme METIS –Simulation d’écoulement et de transport miscible en milieu poreux et fracturé Notice de conception mise à jour le 6/09/10. Rapport Mines ParisTech—Centre de Géosciences, R100907PGOB. Goblet, P., 2005. Programme METIS—Simulation d’écoulement et de transport miscible en milieu poreux et fracturé. Notice d’emploi. Rapport Ecole des Mines de Paris—Centre d’Informatique Géologique LHM/RD/05/24. Goblet, P., 2005. METIS —Tutoriel et notice d’utilisation. Rapport Ecole des Mines de Paris—Centre d’Informatique Géologique LHM/RD/05/ 2005; 25. Goblet, P., 1980. Modélisation des transferts de masse et d’énergie en aquifère [Modelling mass and energy transfers in aquifers]. In: PhD Thesis, Ecole Nationale Supérieure des Mines de Paris. University Pierre and Marie Curie, Paris, France, 200 p. Kamaletdinov, M., Kazantseva, T., Kazantsev, U., Postnikov, D., 1991. Overthrust-nappe tectonics of the lithosphere. Science, Moscow, 254 p. (in Russian). Krilov, V., 1988. Complex investigations of wells and evaluation of exploitation reserves of geothermal waters in terms of geothermal circulation system on Khankala deposit of CHIASSR (in Russian), 128 p. Krilov, V., 1983. Report on the subject Recommendations for the further exploitation of karagan-chokrak deposits at Khankala and Goity withdrawals of thermal waters of CHIASSR (in Russian);104. Matheron, G., 1973. The intrinsic random functions and their applications: Advances in Applied Probability, v. 5, 439–468. Qgis Development Team, 2014. QGIS 1.8.0 Geographic Information System User Guide. Open Source Geospatial Foundation Project, webpage consulted on September 10, http://docs.qgis.org/1.8/pdf/QGIS-1.8-UserGuide-en.pdf. Réveillère, A., Hamm, V., Lesueur, H., Cordier, E., Goblet, P., 2013. Geothermal contribution to the energy mix of a heating network when using aquifer thermal energy storage: modeling and application to the Paris basin. Geothermics 47, 69–79. Shpak, A., 1968. Report on the calculation of exploitation reserves of thermal water at Khankala deposit of CHIASSR (for heating and hot water supply) as of January 1,(in Russian);352. Ungemach, P., Antics, M., Lalos, P., 2011. Sustainable geothermal reservoir management—a modelling suite. In: Proc. Australian Geothermal Energy Conference, 16–18th November 2011, Melbourne, Geoscience Australia Record 2011/43, pp. 267–275.
Contents lists available at ScienceDirect
Geothermics journal homepage: www.elsevier.com/locate/geothermics
A case study of the modeling of a hydrothermal reservoir: Khankala deposit of geothermal waters A. Farkhutdinov a,b,∗ , P. Goblet a , C. de Fouquet a , S. Cherkasov c a b c
Mines ParisTech, 35 rue Saint Honoré, 77305 Fontainebleau, France Bashkir State University, 32 Zaki Validi Str., 450074 Ufa, Russia Vernadsky State Geological Museum, 11 Mokchovaya Str., 125009 Moscow, Russia
a r t i c l e
i n f o
Article history: Received 24 October 2014 Received in revised form 17 August 2015 Accepted 12 October 2015 Keywords: Geothermal waters Khankala deposit Kriging Numerical simulation
a b s t r a c t In 2013 the Grozny State Oil Technical University together with members of the “Geothermal resources” consortium started a pilot project to build a geothermal plant on the basis of the most promising XIII layer of the Khankala thermal waters deposit. Planned capacity of the projected facility is 22.8 GJ/h, which will heat a greenhouse complex. Geostatistical analysis and numerical simulations are required for achieving sustainability in exploitation of the resource and in order to define new drilling sites. In this paper intrinsic random functions of k order kriging (IRF-k) is used for the XIII layer map creation and estimation of the temperature distribution within it. The results obtained by geostatistical approach are used for numerical simulation of the hydro-thermal process relied to reinjection of cold thermal water back into the productive layer. The work carried out allows highlighting the most promising zones for exploitation and shows that the geothermal water temperature in the production well starts to decrease after 6 years. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Geostatistical analysis and simulation is now widely applied in geosciences in order to study complex phenomena. In geothermy, initial data, as a rule, are obtained from irregularly located wells, which complicate interpolation and extrapolation of different parameters in order to model the behavior of a geothermal reservoir. Thus, a geostatistical approach is useful for mapping the reservoir itself, as well as describing permeability, porosity, salinity, transmissivity, temperature, etc., A major task for geostatistics here is to improve the reliability of the mapping of the reservoir’s geometry and the spatial distribution of its characteristics. These descriptions are extremely important for the design of the wells and their further exploitation. In the current paper, geostatistical analysis is applied to the existing data on the Khankala deposit of geothermal waters (Russia) in order to define new drilling sites. The Khankala deposit is located 10 km to the south–east from the city of Grozny (Fig.1). The deposit has a layered structure and contains multiple layers of Chocrack and Karagan horizons of mid-
∗ Corresponding author at: Mines ParisTech, 35 rue Saint Honoré, 77305 Fontainebleau, France. E-mail address: [email protected] (A. Farkhutdinov). http://dx.doi.org/10.1016/j.geothermics.2015.10.005 0375-6505/© 2015 Elsevier Ltd. All rights reserved.
dle Miocene sandstones, interlayered with clays. Tectonically, the deposit occupies the southeastern plunge of the Octyabrsk anticline bordered by two co-axial major faults. The field was actively developed from 1974 to 1994, but its exploitation was ended due to the beginning of the war in Chechnya. In 2013 the Grozny State Oil Technical University together with members of the “Geothermal resources” consortium started a pilot project to build a geothermal plant on the basis of the most promising XIII layer of the Khankala deposit. The water temperature is about 100 ◦ C. Planned capacity of the projected facility is 22.8 GJ/h, which will heat a greenhouse complex. The drilling of a “doublet” is planned (Fig. 2), i.e., one production and one injection well with reinjection of all the cooled thermal water. A numerical model is needed to predict the evolution of the resource in order to study how the production temperature would decrease as a result of the recycling of the reinjected cold water.
2. Theory and methods The aim of the work is to estimate the temperature distribution within the XIII layer of the Khankala deposit of thermal waters, develop a new structural map and make, by numerical simulation, a prognosis for deposit exploitation.
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
57
Fig. 1. General map of the location of the Khankala deposit.
2.1. Geostatistics A geostatistical analysis is used to estimate the geometry of the geothermal site. This analysis is carried out using the ISATIS program (Géovariances, 2009). Intrinsic random functions of order k kriging (IRF-k) was used because of the nonstationarity in data. The idea behind IRF-k, developed by Matheron (1973), is that it is possible to analyze the data considering some variations of the differences between data points to achieve a stationary covariance. Increments of order k are also called authorized linear combinations of order k (k = 0 indicates a difference Z (x) − Z (x + h), k = 1 is an increment of an increment, etc., Usually k is not higher than 2). The increments filter out the mean (k = 0) or “drift”, supposed to be expressed by a polynomial of order k of the coordinates. For example, in a two-dimensional space for filtering linear drift (for k ≥ 1) a0 + a1 x + a2 y, the conditions of authorization are (Chilès and Delfiner, 1999):
i
Fig. 2. “Doublet” used in geothermal water extraction-injection.
i = 0,
i
i xi = 0,
i yi = 0
i
A generalized covariance is used for the estimation. For an intrinsic random function Z(x) of order k there is a symmetrical
58
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
where |x˛ − xˇ | is the distance between 2 points in the horizontal plane. 2.2. Mathematical modeling The estimated temperature is used for the geothermal flux calculation and a structural map is used as a support for the creation of the mesh (Fig. 3). Numerical simulation was made using the Metis program, developed at the Ecole des mines de Paris (Goblet, 1980). Metis simulates liquid flow, mass and heat transport in porous or fractured aquifers. The equations are solved by the finite element method. The aquifer is assumed to respond almost instantaneously to changes in regime and the flow is supposed to be steady (Goblet, 2005a). The equation describing the steady flow is: div (k grad h) = 0
(2)
where K denotes the hydrogeological permeability and h the hydraulic head pressure. Heat transport is solved under transient conditions. It takes into account:
Fig. 3. Scheme of work with data.
- convection in the aquifer (flow of heat by water movement); - conduction in the aquifer and the surrounding formations (heat flux resulting from the temperature gradient); - kinematic thermal dispersion in the aquifer (heat flux resulting from the heterogeneity of the local velocity field). Convection is the dominant mechanism. The exchange with the surrounding formations by conduction delays progression of the cold front. The conduction and the dispersion in the aquifer have the effect of spreading the transition zone between cold water and hot water (Goblet, 2005b,c). The equation describing the heat transport in the aquifer is:
Fig. 4. Structural map of the XIII layer after georeferencing and digitization.
function K(h), called a generalized covariance, which for any pair of increments, Z() and Z() of order k is equal to:
E [Z () Z ()] = ˛ ˇ K |x˛ − xˇ | ˛ ˇ
(1)
div grad-w UD = a
∂ ∂t
(3)
where is the thermal dispersion tensor; the temperature;␥w the volumetric heat capacity of water; a the volumetric heat capacity of the aquifer (rock + water);UD the Darcy velocity. The thermal dispersion includes the pure conduction and heat dissipation kinematics according to the relationship: = 0 + w ˛|UD |
Fig. 5. Omnidirectional sample variogram of the XIII layers top elevation (left) and initial data location (right).
(4)
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
59
where 0 is the coefficient of the thermal conduction of the aquifer (rock + water);˛ the dispersivity tensor. The transport of heat in the low permeability formations is restricted to a mechanism of pure conduction:
div s grad = s
∂ ∂t
(5)
where s is the conduction coefficient; s the calorific capacity of the surrounding formations.
3. Initial data 3.1. Structural map of the XIII layer
Fig. 6. Cross-validation of model for altitude of the XIII layers top estimation.
A first step is to collect data for structural mapping of the XIII layer. A map from 1967 year is used for this purpose. It represents more than 100 wells crossing the XIII layer as the deposit is situated in the Octyabrsk anticline, which contains an oil deposit of the same name. Georeferencing was made in the Quantum GIS program (Qgis Development Team, 2014) using wells found within the Khankala deposit. Elevation at each well which was written on the map was copied and reinterpreted using geostatistical methods. There are 106 points of the top of the XIII layer after map digitization and georeferencing (Fig. 4).
Fig. 7. Estimation of the elevation (top) of the XIII layers and associated standard deviation of the estimation error (bottom).
60
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 8. Comparison of old and new structural maps of the top of the XIII layer.
Fig. 9. Reference plane choice for z coordinates.
3.2. Temperature measurements Geological observations and temperature measurements in the wells were carried out in 1968 and 1988 during the exploitation of the Khankala deposit of thermal waters (Krilov, 1988; Shpak, 1968). Measurements were carried out after temporary stoppage of pumping so as to establish a thermal equilibrium between the well and the surrounding rocks, in total about 100 measurements were made in 14 production wells, each 100–200 m in depth from the surface. As a result of the tragic events that occurred in the territory of the Chechen Republic, data on the inclination of the wells have been lost. The only information available, allowing an estimation of the in-depth position of boreholes in the area, are corrections to the inclination of the wells, which is equal to a maximum of 11 m for depths greater than 1 km, which means a relatively small
inclination. As a consequence, all wells are supposed to be vertical during the geostatistical analysis. The area of study was limited by two faults located in the south–west and north–east of the study area (Farkhutdinov et al., 2014). 3.3. Numerical simulation Conditions for the Metis numerical simulation are chosen according to the data obtained from hydrogeological tests and laboratory examinations (Ermolaev et al., 1954; Krilov, 1983,1988). Aquifer parameters are as follows: – permeability: 1.4e–12 m2 ; – longitudinal and transverse thermal dispersivity: 10 × 2 m; – volumetric heat capacity of water 4.18 MJ/m3 /K;
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
61
Fig. 10. Scatter-diagrams of the temperature versus depth in cases of the different reference plane.
– volumetric heat capacity of the aquifer 2.75 MJ/m3 /K; – volumetric heat capacity of the surrounding formations 2.09 MJ/m3 /K; – conductivity of the aquifer 1.45 W/(m K); – conductivity of the surrounding formations 1.3 W/(m K). – Two boundary conditions are entered: – injection flow of 200 m3 /h; – a constant heat flow of 10,051 MW at the injection well (corresponding to an injection temperature of 45 ◦ C). No boundary conditions were applied to the vertical sides of the model, assuming no liquid or heat flow through these limits. Geothermal flux which is imposed at the base of the model and unequally distributed (see Section 4.1.2) is calculated according to Fourier’s law, taking into account the temperature of the XIII layer of the reservoir, its depth, the temperature at the Earth’s surface in this area and the average coefficient of thermal conductivity of rocks lying between the producing formation and the surface: →
q = −k∇ T
(6)
Application of this formula gives a fairly high average heat flow (155 MW/m2 ), which can be explained by “overthrust-nappe theory” as caused by tectonic movements, considering the North Caucasus as a mobile tectonic zone (Kamaletdinov et al., 1991).
Table 1 Selection of the drift. Drift
Mean error 2
2
1 x y x xy y 1 x y x2 xy y2 x3 x2 y xy2 y3 1xy Without drift
7.116 11.54 42.41 96.11
The sample variogram confirms the known non-stationarity: the mean varies smoothly over the study area, as the fold sinks in a southeasterly direction. An IRF-k modeling is needed, indirect fitting is divided into 3 steps: a Drift finding. b Indirect fitting of the generalized covariance. c Variable estimation and kriging error variance.
IRF − k kriging system :
∗
Estimation : Z =
Z˛ 0
k˛ˇ f ˛
×
fˇ 0
t
×
˛
Estimation varience : Var( ) = K00 −
The first step is to create a structure map of the XIII layer. The omnidirectional sample variogram of the elevation of the top of the XIII layer is calculated (Fig. 5): 1 N(h) 2 [Z (xi ) − Z (xi + h)] 2N (h) i=1
(7)
where N(h) is the number of point pairs separated by a vector h.
=
K˛ˇ
f0
l
4.1.1. Structural map of the productive layer a) Variogram model
(h) =
l
4. Modeling and results 4.1. Geostatistical approach
˛
˛ l
t
×
K˛0
f0
Several possible drifts are first supposed (Table 1) and coefficients are fitted by the least squares method. A quadratic drift (local linear combination of the following functions: 1, x, y, x2 , xy and y2 ) is chosen as the most suitable. To find the optimal model of the generalized covariance, several possible models are considered (Table 2). They are listed by score, which is equal to the theoretical variance divided by experimental variance (it should be close to 1). The chosen model is spherical with a range of 1000 m and a sill of 236.2 m2 . Initial data one-by-one were “hidden” and estimated,
62
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 11. Coordinate system rotation.
• Estimation
Table 2 Selection of the generalized covariance model. Score
Nugget effect (Sill, m2 )
Spherical (Sill, m2 )
Linear (Sill, m2 )
Exponential (Sill, m2 )
1.031 1.055 1.166 1.206
0 0 0 30.8
236.2 0 0 0
0 453.2 0 0
0 0 236.5 216.6
then the difference between initial and estimated value was calculated. Unique kriging neighborhood was selected as the most appropriate, giving the best estimation by cross-validation tests with a coefficient of correlation equal to 0.99 (Fig. 6).
A structural map of the XIII layer is built using this model (Fig. 7). Faults are considered as screens for the estimation. Because there is only one well at the other side of the faults, estimations are not possible to the north and south of the anticline. A standard method of linear interpolation was used during the creation of the map in 1967, based on the measurement points. All data in the vicinity of a point has been taken into account in the new geostatistical interpretation (Fig. 8). Due to the high density of data in the most part of the map, it is not surprising that results are similar. The difference in the south–east is due to a lower density of wells, and to the introduction of a drift.
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
63
Fig. 12. Variograms calculation.
4.1.2. Initial temperature distribution The temperature distribution within the XIII layer and at the surface is supposed to be known, in order to define geothermal flux as one of the boundary conditions for the simulation. The three -dimensional structure of the temperature distribution is examined on directional sample variograms. There were two options for choice of the initial coordinate system for working with temperature measurements: a normal reference plane and the top of the XIII layer as a reference plane (Fig. 9), and both of them were tested. As we can see from scatter-diagrams (Fig. 10) variance of data is less in the case of the XIII layers top reference plane choice, especially in the vicinity of the layer itself. Before working with variograms, the coordinate system was rotated (Fig. 11) so that the anticline folds directly to the east with the depth increasing from west to east and from the axis to the limbs of the structure (i.e., in x and y directions). This facilitates the drift finding as the first step of the IRF-k modeling. 4 horizontal variograms for different directions were calculated to evaluate the benefits of choosing one or another variant for coordinate system (vertical variogram remains unchanged in both cases) (Fig. 12). Although the number of pairs is insufficient for reliable conclusions, especially with increasing distance, it is observed that the variation of temperature with the top of the XIII layer chosen as reference plane is reduced. Thus this variant is used as a selection for the coordinate system. It can be explained by the fact that the XIII layer represents anticline containing geothermal waters and the temperature follows its structure, with undergoing changes only with large differences in depth of formation. Vertical variogram is non-stationary, which corresponds to the temperature progressive increase with depth. It was the main reason for selecting the IRF-k modeling. Samples which are situated less than 300 m from the XIII layer were used and grid for block of 40 × 40 × 40 m size estimation was created. Next were steps of IRF-k modeling (see Section 4.1.1.) (Table 3).
Table 3 The XIII layers temperature estimation parameters. Drift
1yz
Variogram Neighborhood
Spherical (range = 1500 m, sill = 12.7 ◦ C2 ) Unique
Fig. 13. Cross-validation of the XIII layers temperature estimation.
Unique neighborhood was chosen as the most appropriate regarding to the small number of samples. Cross-validation was performed to evaluate chosen parameters (Fig. 13). Cross-validation based on 37 samples showed error variance equal to 6.8 ◦ C and IRF-k kriging was performed (Fig. 14). All components are now available for geothermal flux determination taking surface temperature as average for this region and equal to 24 ◦ C (Fig. 15).
64
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 14. Temperature of the XIII layer estimation (top) and error standard deviation (bottom).
Fig. 15. Geothermal flux estimation.
This geothermal flux was applied as one of the boundary condition for numerical simulation. For zones on the other sides of faults the average value was taken. 4.2. Simulation results The mesh, created for numerical simulation, geometrically represents the 3D structure of the deposit, taking into account the geometry of the XIII productive layer obtained by IRF-k estimation and consists of 141,550 prismatic (6 nodes) elements and 75,348 nodes. It is refined near production and injection wells, which are represented by one-dimensional 2-node elements (Fig. 16). The average size triangle covers an area of 2.5 km2 . Productive layer thickness is equal to 40 m, as it was mentioned as almost
constant in relevant geological reports. We do not have information on the thickness of XIII layer in all wells due to the consequences of two wars (in order to make more precise estimation). Elevation value for zones on the other sides of the faults was obtained by IRF-k estimation (see Section 4.1.2.) without taking these faults into account. It was done because of lack of certainty in shift range and real influence of faults on water flow. Simulation time is equivalent to 30 years of exploitation with a constant average regime (see paragraph initial data). Hydrogeological tests of 1981 and 1988 (Krilov, 1988) showed that the permeable/impermeable status of the north–east and south–west faults is uncertain, so both hypotheses were tested. The simulation
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
65
Fig. 16. The mesh refined near wells.
does not predict any change of temperature in production well for the first 6 years (Fig. 17). Only after this time period the cold plume expands into the production well and the temperature started to decrease (Figs. 17 and 18).
Fig. 18 shows that there would be a difference in well water temperature during production assuming faults are impermeable or not. Thus it was the correct decision to place injection and production wells almost parallel to the faults.
Fig. 17. Simulation of geothermal water reinjection in case of permeable (top) and impermeable (bottom) faults.
66
A. Farkhutdinov et al. / Geothermics 59 (2016) 56–66
Fig. 18. Temperature of geothermal water in production well.
There are a lot of uncertainties in temperature estimation especially for zones where the average value of flux was taken. But regarding the small number of samples it is one of the best possibilities for today which allows making initial forecasts. This work will be continued after drilling the first doublet and obtaining new, more precise data. 5. Summary and conclusions The map of the temperature distribution within the XIII layer of the Khankala deposit of geothermal waters allows, depending on the chosen site for the production well, the prediction of the water temperature in the production well, which in conjunction with a flow rate provides important information on the potential future of the geothermal power station. Numerical simulation shows that some corrections in the exploitation regime will be necessary after a certain period of time in order to make geothermal reservoir management sustainable for more than 30 years (Ungemach et al., 2011). One of the possible solutions could be Aquifer Thermal Energy Storage (Réveillère et al., 2013), which means reinjecting hot water during summer when there is not much need of heating. Acknowledgements The research has been executed with support of the Ministry of Education and Science of Russian Federation, contract no. 02.G25.31.0056. References Chilès, J.-P., Delfiner, P., 1999. Geostatistics: modeling spatial uncertainty. In: Wiley Series in Probability and Statistics. John Wiley & Sons Inc., New York, pp. 695.
Ermolaev, A., Lyalin, L., Napolsky, M., Feldman, S., 1954. Estimation of the Octyabsk deposit oil and gas reserves as of January 1, 1954 (in Russian), 338 p. Farkhutdinov, A., Farkhutdinov, I., Ismagilov, R., 2014. History of the discovery and exploitation of Khankala geothermal waters deposit (in Russian). Bull. Bashkir State Univ. 19, 93–96. Géovariances, 2009. Isatis Technical Ref., ver. 9.06. Geovariances & Ecole Des Mines De Paris: Avon Cedex, France. Goblet, P., Programme METIS –Simulation d’écoulement et de transport miscible en milieu poreux et fracturé Notice de conception mise à jour le 6/09/10. Rapport Mines ParisTech—Centre de Géosciences, R100907PGOB. Goblet, P., 2005. Programme METIS—Simulation d’écoulement et de transport miscible en milieu poreux et fracturé. Notice d’emploi. Rapport Ecole des Mines de Paris—Centre d’Informatique Géologique LHM/RD/05/24. Goblet, P., 2005. METIS —Tutoriel et notice d’utilisation. Rapport Ecole des Mines de Paris—Centre d’Informatique Géologique LHM/RD/05/ 2005; 25. Goblet, P., 1980. Modélisation des transferts de masse et d’énergie en aquifère [Modelling mass and energy transfers in aquifers]. In: PhD Thesis, Ecole Nationale Supérieure des Mines de Paris. University Pierre and Marie Curie, Paris, France, 200 p. Kamaletdinov, M., Kazantseva, T., Kazantsev, U., Postnikov, D., 1991. Overthrust-nappe tectonics of the lithosphere. Science, Moscow, 254 p. (in Russian). Krilov, V., 1988. Complex investigations of wells and evaluation of exploitation reserves of geothermal waters in terms of geothermal circulation system on Khankala deposit of CHIASSR (in Russian), 128 p. Krilov, V., 1983. Report on the subject Recommendations for the further exploitation of karagan-chokrak deposits at Khankala and Goity withdrawals of thermal waters of CHIASSR (in Russian);104. Matheron, G., 1973. The intrinsic random functions and their applications: Advances in Applied Probability, v. 5, 439–468. Qgis Development Team, 2014. QGIS 1.8.0 Geographic Information System User Guide. Open Source Geospatial Foundation Project, webpage consulted on September 10, http://docs.qgis.org/1.8/pdf/QGIS-1.8-UserGuide-en.pdf. Réveillère, A., Hamm, V., Lesueur, H., Cordier, E., Goblet, P., 2013. Geothermal contribution to the energy mix of a heating network when using aquifer thermal energy storage: modeling and application to the Paris basin. Geothermics 47, 69–79. Shpak, A., 1968. Report on the calculation of exploitation reserves of thermal water at Khankala deposit of CHIASSR (for heating and hot water supply) as of January 1,(in Russian);352. Ungemach, P., Antics, M., Lalos, P., 2011. Sustainable geothermal reservoir management—a modelling suite. In: Proc. Australian Geothermal Energy Conference, 16–18th November 2011, Melbourne, Geoscience Australia Record 2011/43, pp. 267–275.