Abstract
Strong chemical reactions in the geothermal systems may cause sealing of fractures, which reduces the permeability in the reservoir and subsequently affects the heat production. However, it is difficult to reveal the sealing range in a deeply buried reservoir based on a limited number of downhole logs. This study recreated the sealing processes of the fault-controlled geothermal system in the Guide Basin, China, by reactive transport modeling. The modeling domain was discretized based on multiple interacting continua (MINC) approach, to address the nonequilibrium heat transport processes between the matrix and conduit in the fractured fault damage zone. Once the model was validated by observations of major ions in spring water and downhole temperature logs in the discharge area, it was used to determine the coupled processes of fluid, heat, and chemical transport in the reservoir and the resultant sealing ranges. It was found that the dissolution of albite and K-feldspar leads to the precipitation of smectite-Ca and illite in the middle and bottom of the fault under the condition of high concentration of Ca2+ and Mg2+ in the recharge water. Calcite veins were formed in discharge zone, because the horizontal fast flow in shallow subsurface zone supplied abundant Ca2+ and HCO3-. As a consequence, the permeability in the discharge zone reduced by 15% when compared to the original permeability of 100 mD. Moreover, another three self-sealing areas were formed near the recharge zone, the deep upgradient zone, and the downgradient area where the fast upward fluid flow occurred. Self-sealing subsequently prevented the deep circulation of the flow and heat absorption, which tends to make the fault-controlled geothermal system inactive.
1. Introduction
Water circulation in the geothermal systems accompanies with mineral deposition along the flow paths due to chemical reaction, which fills the fractures and resists water circulation. Such process is also referred to as self-sealing [1]. Delineating self-sealing process is critical to understand the history of geothermal processes and to locate the heat production wells.
Self-sealing processes were often investigated based on the mineral analysis of rock samples collected at outcrop zone and drillholes. Dubois et al. [2] used transmitted light microscopy and scanning electron microscope (SEM) to reveal the secondary mineralization processes in rock samples collected in the region of Soultz, France, and recovered the fluid and heat transport history in the Rhine continental rift basin. Ngwenya et al. [3] conducted laboratory experiments and mass balance calculations and found that diagenetic processes can induce up to two orders of magnitude of the total sealing observed in natural faults. Tian et al. [4] used the one-dimensional models to study the effects of mineral alteration on the self-sealing for CO2 geological storage. Based on the field outcrop observations and laboratory experiments [5], numerical modelling of reactive transport is necessary to upscale the local present-day observations to the long-term and field-scale self-sealing processes.
The Zhacang geothermal field in the Guide Basin of northeastern Tibetan Plateau is a typical fault-controlled geothermal system. Magma in the depth of 20 to 40 km induced a terrestrial heat flow of 76.5 mW/m2 in this region [6, 7]. Active faults allow the surface water to flow deeply and, after heated, discharge as geothermal springs with temperature exceeds 90°C. The calcite scaling in the discharge area and the calcite veins in the outcrop of reservoir rocks indicate the occurrence of self-sealing process in this geothermal field. However, the self-sealing range is still underdetermined, which potentially affects the efficiency of heat production.
In this study, a realistic reactive transport model is established by TOUGHREACT simulator [8–10], to assess the geochemical evolution of geothermal fluids in the Zhacang geothermal field. The model is validated by downhole temperature logs and chemical compositions in water collected from hot springs and borehole. Subsequently, the validated model is used to describe the chemical evolution and self-sealing processes and importantly to determine the zone where porosity and permeability decrease significantly.
2. Study Area and Sampling
2.1. Geologic Setting
The Zhacang geothermal field is located in the Guide Basin of the northeastern Tibetan Plateau, China (Figure 1(a)). It is surrounded by Laji Mountains to the north, Zhama Mountains to the East, Waligong Mountains to the west, and Luobuleng Mountains to the southwest. Water circulation in the Zhacang field is controlled by the conductive fault F2, which allows the water from Gangyi Stream seeped into the subsurface and, after heated, discharge as geothermal springs at the eastern end of the fault F2 (Figure 1(b)). The damage zone of the fault F2, with the depth of 3000 m, width of 40 m, and length of 2500 m, controls the fluid and heat transport in the Zhacang geothermal field [11].
2.2. Sampling and Analysis
Water samples were collected from the Zhacang geothermal field in 2015 to 2017 (Figure 1(b)). A total of six water samples were collected, including one borehole water sample (ZR1), one hot spring water sample (G1), three hot well water samples (G2, G3, and G4), and one cold spring water sample (G5). Temperature and pH were measured in situ with portable instruments with uncertainties of ±0.1°C and ±0.05, respectively, while the common ions were analyzed by ion chromatography, inductively coupled plasma atomic emission spectroscopy (ICP-AES), inductively coupled plasma mass spectrometry (ICP-MS), and acid titration (AT) in Jilin University, China. The geochemical analysis data at these sampling positions are summarized in Table 1. The hydrochemical characteristic of borehole ZR1 represents the deep fluid chemistry, presenting a high SiO2 concentration of 125 mg/L. All the hot waters (G1, G2, G3, and G4), discharged from the end of the fault F2 (Figure 1), have the similar chemical components with the type of SO4·Cl-Na, while the recharge water chemistry of the fault F2 (G5) is the type of Cl·SO4-Na·Mg. Ca2+ and Mg2+ concentrations in the recharge water are higher than those in the discharge water, which provide a source of cations for deep water-rock reactions. This indicates that the precipitation occurs along the flow path from recharge area to discharge area, causing self-sealing.
Calcite veins relating to the self-sealing processes were observed in the discharge area of the Zhacang geothermal field. Three rock samples (Y-3, Y-4, and Y-5) containing secondary mineral deposits were collected in the veins (Figure 1). The SEM images (Figure 2) and X-ray diffraction analysis (Table 2) show that these secondary minerals are mainly composed of calcite with minor hydrothermal alteration minerals of illite and smectite.
(a)
(b)
(c)
(d)
(e)
(f)
Moreover, the hosting rocks of granite in the fault F2 were collected from borehole ZR1 (Table 3). X-ray diffraction analysis shows that it mainly consists of quartz (24 vol.%), albite (30 vol.%), K-feldspar (30 vol.%), calcite (4 vol.%), and muscovite (12 vol.%).
3. Numerical Modeling
3.1. Theory
The governing equations for multiphase fluid and heat transport in subsurface geothermal reservoirs can be expressed in a uniform formation [12]: where is the volume of discretized elements, represents the boundary of the model domain, is the fluid mass when calculating the fluid transport and is the energy when modeling the heat transport, and 2 represent the liquid and vapor phase of water, respectively, and 3 for heat transport, denotes mass or heat flux, denotes sinks and sources, is a normal vector on surface element , and is for water transport following the Darcy’s law: where is the Darcy velocity (volume flux) in phase (liquid or vapor), is absolute permeability, is relative permeability to phase , is viscosity, is the effective pressure gradient in phase , is the density of phase , and is the vector of gravitational acceleration.
The chemical transport is governed by the classical advection diffusion equation as follows: where refers to the concentration of aqueous component , is the porosity of the porous medium, refers to the diffusion coefficient of species , refers to the tortuosity, and is the Darcy flux. refers to the total reaction rate of any reaction affecting the component concentration .
Changes in porosity and permeability due to mineral dissolution and precipitation are described by the simplified cubic law as follows: where and are the initial permeability and porosity, respectively.
Following the geochemical data in geothermal fluid and rock samples, the ions in fluid considering in the water-rock interaction in this work are listed in Table 1, with the initial concentration in the model domain following those detected in ZR1. The minerals accounted for include quartz, albite, K-feldspar, calcite, muscovite, illite, Ca-smectite, chlorite, dolomite, and kaolinite (Table 3).
The governing equations above are solved by nonisothermal reactive geochemical transport code TOUGHREACT [9, 10], where fluid and heat transport modeling are fully coupled, which are sequentially coupled with chemical transport processes. The thermodynamic data are drawn from the EQ3/6 database of Wolery [13]. The kinetic reaction processes between fluid and minerals are controlled as follows [14]: where refers to the reaction rate, is the reactive surface area, refers to the reaction rate constant at 25°C, is the activation energy, and are the temperature (K) and ideal gas constant, and refer to the ion activity product and equilibrium constant of a specific mineral reaction, and and are constants resulting from experimental data; they are usually but not always taken equal to 1. The coefficients in equation (5) are listed in Table 3 [15].
3.2. Model Geometry, Fluid Flow, and Heat Transfer Simulations
The modeling domain in the damage zone of F2 was discretized into 3099 elements, and the grids are refined at both recharge and discharge zones (Figure 3(b)). Moreover, to consider the nonequilibrium heat exchange between matrix and fractures, “multiple interacting continua (MINC)” is applied to further describe the grid blocks as dual continuum domain, with the outer layer representing the fracture with a volume fraction of 5%, and the inner layer representing the matrix with a volume fraction 95%, and the fracture spacing is set to 10 m determined according to field investigation. The porosity of the fractures and matrix is 0.5 and 0.025, respectively, and the permeability is 10-13 m2 and 10-17 m2, respectively [16]. Other physical parameters are summarized in Figure 3(a) [17], and the model domain is assumed to be isotropic and homogeneous.
(a)
(b)
(c)
As shown in Figure 3, water from the Gangyi Stream seeps into the F2 within 500 m in the west, with the inflow rate of 210 m3/d and recharge temperature of 7.2°C (following local average atmospheric temperature). Jiang et al. [11] simulated the transport of the helium ratio in this model domain and suggested that there is no external fluid input from the deep subsurface in the Zhacang field of the Guide Basin. When the water flow in the fault encounters the impermeable fault F5, the hot water is naturally discharged through the hot springs near the borehole ZR1, which is defined as a constant pressure boundary with the pressure of 1 bar and no heat flux boundary. Other part on the land surface is defined as no flow boundary with constant temperature of 7.2°C. The bottom of the model is defined as a no flow boundary with a temperature of 150°C according to downhole temperature logs in ZR1, while heat transfer from both lateral sides of model domain is calculated by semianalytical solution as follows [18]: where is the distance from the boundary, is initial temperature in surrounding rock, is the time-varying temperature at the boundary, and are time-varying fit parameters, and is the penetration depth for heat conduction.
3.3. Model Validation
When the temperature stabilizes in the model domain after a simulation period of 8000 years, the calculated temperatures in ZR1 compare well with the downhole temperature logs (Figure 4(b)). It is illustrated in Figure 4(a) that the steady-state temperature in this fault-controlled geothermal system is mainly controlled by the convection. The average flow velocity in the fault F2 is approximately between 0.01 and 0.02 m/d, and a fast upflow with a flow velocity of 0.04 m/d occurs when the flow is blocked by the impermeable fault F5 (Figure 4(a)). The water-rock interaction in this fast upflow path is of great significance to the hot springs (see further).
(a)
(b)
(c)
(d)
In addition to the deep fluid flow, part of recharge water can flow rapidly in the shallow subsurface zone without heating, driven by the hydraulic gradient between recharge area and discharge area (Figure 4(a)). These unheated horizontal fluids eventually meet the heated upflow fluids in the discharge zone and affect the water chemical composition of the hot springs. In addition, these unheated horizontal fluids are also related to the formation of hydrothermal alteration minerals, i.e., calcite veins outcropped on the surface (details in Section 4). The gas saturation in the model domain shows that an unsaturated zone exists within the depth of 50 m, which however does not affect the chemical reaction process in the modeling domain (Figure 4(c)).
Moreover, as shown in Figure 4(d), the calculated water chemistry compares well with observations in the discharge area. This verifies the rationality of the initial geochemical and mineral composition in this reactive transport modeling.
4. Results and Discussion
4.1. Water-Rock Interactions in Fault-Controlled Geothermal System
After a simulated time period of 8000 years, mineral dissolution and precipitation occurred with fluid and heat transport. The variation of mineral components is expressed by changes of mineral abundance in volume fraction (dimensionless). Mineral abundance changes are equal to the change in porosity. It is shown in Figure 5 that a chemical active zone is located at upgradient area of fluid flow within depth between 1300 and 3000 m, due to the reaction of low-temperature recharge water and high-temperature minerals in reservoir rocks. The volume fraction of quartz increases at depth between 2500 m and 3000 m, with the maximum increases reach 0.02 (Figure 5(a)). Under the condition of low SiO2 concentration (<12 mg/L) in recharge water (G5 in Table 1), the formation of quartz is related to the dissolution of albite, K-feldspar, and muscovite in the zones near the quartz formation (Figures 5(b), 5(c), and 5(f)) in the following:
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Such chemical reaction not only supplies SiO2 for the formation of the quartz but also releases K+ and Na+ (Figures 6(d) and 6(f)) and generates intermediate mineral of the kaolinite. However, no precipitation of kaolinite is observed in the simulation, because it reacts quickly with recharge water with high concentration of Ca2+ and Mg2+ (Table 1), forming smectite-Ca at depth of 1300 to 2300 m and illite at 2300 to 3000 m, as follows:
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Consequently, concentrations of Ca2+ and Mg2+ in water decrease along the flow path due to the consumption of smectite-Ca and illite in middle-bottom of the model.
It is also observed in Figure 6(i) that the pH values increase from 7.2 to 9.1 from recharge area to discharge area. This is because the dissolution of muscovite consumes a large amount of : which release Al3+ (~0.005 mg/L) in discharge water. This can explain the minor components of Al3+ detected in the discharge water under a condition of undetected Al3+ in recharge water. Subsequently, due to the chemical reaction of calcite is precipitated in the depth between 1300 and 2300 m and also in both recharge and discharge areas (Figure 5(g)), while dolomite is formed at the bottom of the fault zone (Figure 5(h)). However, the large number of calcite veins outcropping in the geothermal discharge area is not due to the upflow fluid along fractured path, as merely limited amount of Ca2+ and HCO3- (<3 mg/L and 17 mg/L, respectively) in this upflow fluid. It is rather attributed to the horizontal flow fluid in shallow subsurface zone that carries abundant Ca2+ and HCO3- (>100 mg/L and 150 mg/L, respectively) (Figures 6(c) and 6(e)). In addition, the concentration of Fe2+ in the recharge water and primary minerals is merely 0.012 mg/L (Figure 6(g)), there is almost no precipitation of chlorite (not listed in the figure).
Except for the deep reactive zone in the upgradient of fluid flow, another chemical active zone exists at the cross of the fault F2 and impermeable fault F5 (the fast upflow rising along the fractured path), where the fluid is forced to flow upward rapidly with an average flow velocity of 0.04 m/d. To understand the influence of upward fast flow on the chemical reaction, the mineral saturation index (SI), mineral abundance, and chemical component concentration along ZR1 are illustrated in Figure 7.
(a)
(b)
(c)
(d)
The SI values of illite, K-feldspar, and quartz are larger than “0” (Figure 7(b)), indicating the precipitation processes in the fast upflow path, especially for illite (Figure 7(d)). The SI value of calcite is equal to zero along the upflow path (Figure 7(b)), indicating no calcite precipitation along ZR1 within . However, calcite precipitates near the land surface (with ) in the discharge area (Figure 7(d)), because the upflow meets the horizontal flow in shallow subsurface zone carrying a large amount of Ca2+ and HCO3- (Figures 6(c) and 6(e)).
In contrast, the dissolution of muscovite and albite occurred (Figures 7(b) and 7(d)). Although the SI value for albite is only slightly smaller than 0 (~-0.0001), a large amount of dissolution has occurred. Smectite-Ca, kaolinite, and dolomite are secondary minerals, which are not dissolved.
Figure 7(c) shows the evolution of the concentration of each chemical component during deep geothermal fluid flowing upward. The concentration of most chemical components maintains unchanged, but the concentrations of Al3+ and K+ decrease significantly. The decrease of the Al3+ concentration is due to the decrease in the dissolution of muscovite under lowering temperature (Figure 7(d)), while the decrease of K+ concentration is simultaneously affected by three minerals: muscovite, K-feldspar, and illite. A small amount of muscovite dissolution can produce K+, and a small amount of K-feldspar precipitates can consume K+. However, a large amount of illite precipitates consumes a large amount of K+, which controls reduction of K+.
4.2. Self-Sealing Effect
As a result, four domains of mineral precipitation (self-sealing) were formed, including the recharge area, deep upgradient zone, fast upflow zone, and discharge area. The mineral precipitation further leads to the decrease of porosity and permeability (Figures 8(a) and 8(b)). The permeability of the recharge area drops by 20%, relating to the precipitation of calcite (with the mineral abundance increasing by 0.06, Figure 5(g)). The permeability of the discharge area drops by 15%; the precipitated minerals include calcite (0.04 of volume fraction), albite (0.01 of volume fraction), and K-feldspar (0.01 of volume fraction) (Figures 5(b), 5(c), and 5(g)). Moreover, the permeability along the fast upflow path and in the deep upgradient zone drops by 5%, due to the precipitation of illite (with the mineral abundance increasing by 0.04, Figure 5(e)).
The large amount of Ca2+ and Mg2+ contained in the recharge water of Gangyi Spring is a major reason for the self-sealing effect. The secondary precipitation in the deep reservoir can reduce the circulation depth of the fluid. Consequently, insufficient heat can be extracted by the fluid flow from recharge area to discharge area. For present-day exploitation, we expect to locate the heat extraction well at the zone with high temperature and permeability. As shown in Figure 8, at the position of drillhole ZR1, it is preferred to locate the well screen at the depth between 1800 and 2800 m.
5. Conclusions
This study delineated the zone of self-sealing in the Zhacang geothermal field in the Guide Basin by a 2D reactive transport model. The following conclusions can be drawn: (1)The Zhacang geothermal field in the Guide Basin is a fault-controlled geothermal system, and the temperature distribution is dominated by fluid convection. There are two main flow paths for recharge water flowing toward discharge area: (i) a heated deep circulation flow, including a fast upflow with velocity 0.04 m/d on the downgradient of fluid flow and (ii) an unheated horizontal flow in the shallow subsurface zone which eventually mixed the heated upflow in the discharge area, causing self-sealing(2)The recharge water carried a large amount of Ca2+ and Mg2+, which leads to the precipitation of calcite, smectite-Ca, and illite and dissolution of albite and K-feldspar. Muscovite was only slightly dissolved in high-temperature domains, and the dissolution of muscovite controls the distribution of Al3+(3)Self-sealing occurred in four major domains: recharge area, deep upgradient area, fast upflow zone at downgradient area, and discharge area. In the recharge area, the permeability dropped by 20%, while in the discharge area dropped by 15%. Along the fast upflow path, the permeability dropped by 5%, and in the deep upgradient zone, the permeability dropped by 5%. For present-day geothermal exploitation of the borehole ZR1 position, the best well screen depth is between 1800 and 2800 m
The range of problems concerning the processes in the geothermal system is very broad. The present modeling results are specific to the conditions and parameters considered. The “numerical experiments” do give a detailed understanding of the dynamic evolution and provide useful insight into fluid flow, geochemical, and mineralization processes in the 2D modeling domain.
Data Availability
All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (Grant Nos. 41572215 and 41502222), the China Geological Survey (Grant Nos. DD20179621 and DD20189113), and the Special Project for New Energy, Jilin Province (Grant No. SXGJSF2017-5). We are grateful to Yude Lei of Environmental Geological Prospecting Bureau of Qinghai Province for assisting our field sampling works.