# Numerical Modeling of Reactive Transport and Self-Sealing Processes in the Fault-Controlled Geothermal System of the Guide Basin, China

Authors: Zhaoyun Hou; Tianfu Xu; Guanhong Feng; Bo Feng; Yilong Yuan; Zhenjiao Jiang
Journal: Geofluids
 (2019)

---

## 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.


---

## Body

## 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 (Figure1(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].Figure 1
(a) Landform of the Guide Basin and (b) the geological map of the Zhacang geothermal field and sampling locations. A total of 6 water samples were collected with the numbers of ZR1, G1, G2, G3, G4, and G5. The solid red lines on (b) correspond to the faults, and the W-E trending fault of F2 is conductive and the N-S trending fault of F5 is impermeable. Hot springs discharge at the junction of F2 and F5.


### 2.2. Sampling and Analysis

Water samples were collected from the Zhacang geothermal field in 2015 to 2017 (Figure1(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.Table 1
Water chemistry and temperatures in recharge and discharge waters from the Zhacang geothermal field of Guide Basin (mg/L). The temperature of ZR1 refers to the sampling temperature (90°C); the bottom temperature is measured at 151.34°C.
Sample no.
ZR1
G1
G2
G3
G4
G5
Location
Borehole ZR1
Zhacang Valley
Zhacang Valley
Zhacang Valley
Hotspring Village
Gangyi Spring
Date
2017.09.09
2016.08.14
2016.08.14
2015.08.15
2015.07.27
2017.09.09
Description
Borehole
Hot spring
Hot well
Hot well
Hot well
Cold spring
Depth (m)
3050.7
Artesian
100.0
100.0
Unknown
Artesian
Elevation (m)
2514
2500
2500
2505
2475
2781
T (°C)
90.0
70.4
78.3
80.0
66.2
11.0
pH
8.33
8.20
7.99
8.40
8.36
8.21
K+
23.44
17.14
19.18
19.90
15.90
6.96
Na+
369.23
389.19
389.82
448.00
429.00
441.18
Ca2+
55.26
63.01
59.14
58.20
72.80
157.25
Mg2+
1.59
3.94
0.55
0.98
3.92
124.38
Al3+
0.096
0.038
0.068
—
—
0.005
TFe
0.0050
0.0026
0.0020
—
—
0.0120
Cl-
359.76
334.11
334.26
334.00
351.00
696.63
SO42-
383.43
532.54
533.09
562.00
538.00
474.65
HCO3-
39.92
55.79
53.80
39.70
73.30
244.09
CO32-
11.32
0.00
0.00
18.00
3.00
0.00
SiO2
125.61
93.40
87.79
56.69
76.35
11.76
TDS
1243.95
1401.05
1391.45
1480.78
1489.12
2339.38
Type
Cl·SO4-Na
SO4·Cl-Na
SO4·Cl-Na
SO4·Cl-Na
SO4·Cl-Na
Cl·SO4-Na·Mg
—: not detected.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 (Figure1). 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.Figure 2
Photos and SEM images of three rock samples (Y-3, Y-4, and Y-5) collected from the calcite veins in the discharge area of Zhacang geothermal field.
(a)
(b)
(c)
(d)
(e)
(f)Table 2
Mineralogical compositions (wt%) in calcite veins detected by X-ray diffraction.
Sample no.
Y-3
Y-4
Y-5
Description
Calcite vein
Calcite vein
Calcite vein
Depth (m)
Outcrop
Outcrop
Outcrop
Calcite
52
62
55
Quartz
4
12
0
K-feldspar
2
0
0
Albite
3
0
0
Illite/smectite
39
26
43
Gypsum
0
0
2Moreover, the hosting rocks of granite in the fault F2 were collected from borehole ZR1 (Table3). 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.%).Table 3
Summary of primary and secondary minerals considered in the numerical simulation and their initial volume fractions, kinetic constants (kn,25), activation energies (Ea), and specific reactive areas.
Mineral
Chemical formula
Volume fraction (vol.%)
kn,25 (mol m-2 s-1)
Ea (kJ/mol)
Specific reactive area (cm2/g)
Primary minerals
Quartz
SiO2
0.24
1.0e - 14
87.70
9.8
Albite
NaAlSi3O8
0.30
2.8e - 13
69.80
9.8
K-feldspar
KAlSi3O8
0.30
3.9e - 13
38.00
9.8
Calcite
CaCO3
0.04
Equilibrium
—
—
Muscovite
KAl2(AlSi3O10)(OH)2
0.12
2.8e - 14
22.00
9.8
Secondary minerals
Illite
K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2
0
1.7e - 13
35.00
151.6
Chlorite
Mg2.5Fe2.5Al2Si3O10(OH)8
0
3.0e - 13
88.00
9.8
Ca-smectite
Ca0.145Mg0.26Al1.77Si3.97O10(OH)2
0
1.7e - 13
35.00
151.6
Dolomite
CaMg(CO3)2
0
3.0e - 8
52.20
9.8
Kaolinite
Al2Si2O5(OH)4
0
6.9e - 14
22.20
151.6
Equilibrium: aqueous solution was forced to be at equilibrium with respect to the corresponding minerals as long as these are present or supersaturated (logQ/K>0) within a specific grid block.


## 2.1. Geologic Setting

The Zhacang geothermal field is located in the Guide Basin of the northeastern Tibetan Plateau, China (Figure1(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].Figure 1
(a) Landform of the Guide Basin and (b) the geological map of the Zhacang geothermal field and sampling locations. A total of 6 water samples were collected with the numbers of ZR1, G1, G2, G3, G4, and G5. The solid red lines on (b) correspond to the faults, and the W-E trending fault of F2 is conductive and the N-S trending fault of F5 is impermeable. Hot springs discharge at the junction of F2 and F5.


## 2.2. Sampling and Analysis

Water samples were collected from the Zhacang geothermal field in 2015 to 2017 (Figure1(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.Table 1
Water chemistry and temperatures in recharge and discharge waters from the Zhacang geothermal field of Guide Basin (mg/L). The temperature of ZR1 refers to the sampling temperature (90°C); the bottom temperature is measured at 151.34°C.
Sample no.
ZR1
G1
G2
G3
G4
G5
Location
Borehole ZR1
Zhacang Valley
Zhacang Valley
Zhacang Valley
Hotspring Village
Gangyi Spring
Date
2017.09.09
2016.08.14
2016.08.14
2015.08.15
2015.07.27
2017.09.09
Description
Borehole
Hot spring
Hot well
Hot well
Hot well
Cold spring
Depth (m)
3050.7
Artesian
100.0
100.0
Unknown
Artesian
Elevation (m)
2514
2500
2500
2505
2475
2781
T (°C)
90.0
70.4
78.3
80.0
66.2
11.0
pH
8.33
8.20
7.99
8.40
8.36
8.21
K+
23.44
17.14
19.18
19.90
15.90
6.96
Na+
369.23
389.19
389.82
448.00
429.00
441.18
Ca2+
55.26
63.01
59.14
58.20
72.80
157.25
Mg2+
1.59
3.94
0.55
0.98
3.92
124.38
Al3+
0.096
0.038
0.068
—
—
0.005
TFe
0.0050
0.0026
0.0020
—
—
0.0120
Cl-
359.76
334.11
334.26
334.00
351.00
696.63
SO42-
383.43
532.54
533.09
562.00
538.00
474.65
HCO3-
39.92
55.79
53.80
39.70
73.30
244.09
CO32-
11.32
0.00
0.00
18.00
3.00
0.00
SiO2
125.61
93.40
87.79
56.69
76.35
11.76
TDS
1243.95
1401.05
1391.45
1480.78
1489.12
2339.38
Type
Cl·SO4-Na
SO4·Cl-Na
SO4·Cl-Na
SO4·Cl-Na
SO4·Cl-Na
Cl·SO4-Na·Mg
—: not detected.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 (Figure1). 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.Figure 2
Photos and SEM images of three rock samples (Y-3, Y-4, and Y-5) collected from the calcite veins in the discharge area of Zhacang geothermal field.
(a)
(b)
(c)
(d)
(e)
(f)Table 2
Mineralogical compositions (wt%) in calcite veins detected by X-ray diffraction.
Sample no.
Y-3
Y-4
Y-5
Description
Calcite vein
Calcite vein
Calcite vein
Depth (m)
Outcrop
Outcrop
Outcrop
Calcite
52
62
55
Quartz
4
12
0
K-feldspar
2
0
0
Albite
3
0
0
Illite/smectite
39
26
43
Gypsum
0
0
2Moreover, the hosting rocks of granite in the fault F2 were collected from borehole ZR1 (Table3). 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.%).Table 3
Summary of primary and secondary minerals considered in the numerical simulation and their initial volume fractions, kinetic constants (kn,25), activation energies (Ea), and specific reactive areas.
Mineral
Chemical formula
Volume fraction (vol.%)
kn,25 (mol m-2 s-1)
Ea (kJ/mol)
Specific reactive area (cm2/g)
Primary minerals
Quartz
SiO2
0.24
1.0e - 14
87.70
9.8
Albite
NaAlSi3O8
0.30
2.8e - 13
69.80
9.8
K-feldspar
KAlSi3O8
0.30
3.9e - 13
38.00
9.8
Calcite
CaCO3
0.04
Equilibrium
—
—
Muscovite
KAl2(AlSi3O10)(OH)2
0.12
2.8e - 14
22.00
9.8
Secondary minerals
Illite
K0.6Mg0.25Al1.8(Al0.5Si3.5O10)(OH)2
0
1.7e - 13
35.00
151.6
Chlorite
Mg2.5Fe2.5Al2Si3O10(OH)8
0
3.0e - 13
88.00
9.8
Ca-smectite
Ca0.145Mg0.26Al1.77Si3.97O10(OH)2
0
1.7e - 13
35.00
151.6
Dolomite
CaMg(CO3)2
0
3.0e - 8
52.20
9.8
Kaolinite
Al2Si2O5(OH)4
0
6.9e - 14
22.20
151.6
Equilibrium: aqueous solution was forced to be at equilibrium with respect to the corresponding minerals as long as these are present or supersaturated (logQ/K>0) within a specific grid block.


## 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]:
(1)ddt∫VnMκdVn=∫ΓnFκ⋅ndΓn+∫VnqκdVn,where Vn is the volume of discretized elements, Γn represents the boundary of the model domain, M is the fluid mass when calculating the fluid transport and is the energy when modeling the heat transport, κ=1 and 2 represent the liquid and vapor phase of water, respectively, and 3 for heat transport, F denotes mass or heat flux, q denotes sinks and sources, n is a normal vector on surface element dΓn, and Fβ is for water transport following the Darcy’s law:
(2)Fβ=ρβuβ=−kkrβρβμβ∇Pβ−ρβg,where uβ is the Darcy velocity (volume flux) in phase β (liquid or vapor), k is absolute permeability, krβ is relative permeability to phase β, μβ is viscosity, ∇Pβ is the effective pressure gradient in phase β, ρβ is the density of phase β, and g is the vector of gravitational acceleration.The chemical transport is governed by the classical advection diffusion equation as follows:(3)∂ϕCi∂t=∇⋅ϕ⋅Di⋅τ⋅∇Ci−ϕ⋅u⋅∇Ci+R,where Ci refers to the concentration of aqueous component i, ϕ is the porosity of the porous medium, Di refers to the diffusion coefficient of species i, τ refers to the tortuosity, and u is the Darcy flux. R refers to the total reaction rate of any reaction affecting the component concentration i.Changes in porosity and permeability due to mineral dissolution and precipitation are described by the simplified cubic law as follows:(4)k=kiϕϕi3,where ki and ϕi 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 Table1, 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]:
(5)r=Ar⋅k25exp−EaR1T−1298.15⋅1−QKmn,where r refers to the reaction rate, Ar is the reactive surface area, k25 refers to the reaction rate constant at 25°C, Ea is the activation energy, T and R are the temperature (K) and ideal gas constant, Q and K refer to the ion activity product and equilibrium constant of a specific mineral reaction, and m and n 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 (Figure3(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.Figure 3
(a) Model domain of the conductive fault F2 and physical parameters, where land surface elevation extracted from digital elevation model, and the depth of F2 from the electromagnetic and gravity prospecting [19]; (b) irregular mesh generated by Gambit software (Fluent Inc.) with mesh refinement in the recharge and discharge zones; (c) boundary conditions for flow and heat in the fault F2 in the numerical model.
(a)
(b)
(c)As shown in Figure3, 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]:
(6)Tx,t−Ti=Tf−Ti+px+qx2exp−xd,where x is the distance from the boundary, Ti is initial temperature in surrounding rock, Tf is the time-varying temperature at the boundary, p and q are time-varying fit parameters, and d 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 (Figure4(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).Figure 4
(a) Steady-state temperature distribution and fluid flow field. (b) Calculated temperatures fit well with the observed temperatures in ZR1. (c) Gas saturation (Sg) distribution illustrates the existence of unsaturated zone at the top of the model. (d) The calculated chemical components in discharge area fit well with observations (ZR1, G1, G2, G3, and G4). The value of “HCO3-” is the sum of HCO3- and CO32-.
(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 (Figure4(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 Figure4(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.


## 3.1. Theory

The governing equations for multiphase fluid and heat transport in subsurface geothermal reservoirs can be expressed in a uniform formation [12]:
(1)ddt∫VnMκdVn=∫ΓnFκ⋅ndΓn+∫VnqκdVn,where Vn is the volume of discretized elements, Γn represents the boundary of the model domain, M is the fluid mass when calculating the fluid transport and is the energy when modeling the heat transport, κ=1 and 2 represent the liquid and vapor phase of water, respectively, and 3 for heat transport, F denotes mass or heat flux, q denotes sinks and sources, n is a normal vector on surface element dΓn, and Fβ is for water transport following the Darcy’s law:
(2)Fβ=ρβuβ=−kkrβρβμβ∇Pβ−ρβg,where uβ is the Darcy velocity (volume flux) in phase β (liquid or vapor), k is absolute permeability, krβ is relative permeability to phase β, μβ is viscosity, ∇Pβ is the effective pressure gradient in phase β, ρβ is the density of phase β, and g is the vector of gravitational acceleration.The chemical transport is governed by the classical advection diffusion equation as follows:(3)∂ϕCi∂t=∇⋅ϕ⋅Di⋅τ⋅∇Ci−ϕ⋅u⋅∇Ci+R,where Ci refers to the concentration of aqueous component i, ϕ is the porosity of the porous medium, Di refers to the diffusion coefficient of species i, τ refers to the tortuosity, and u is the Darcy flux. R refers to the total reaction rate of any reaction affecting the component concentration i.Changes in porosity and permeability due to mineral dissolution and precipitation are described by the simplified cubic law as follows:(4)k=kiϕϕi3,where ki and ϕi 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 Table1, 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]:
(5)r=Ar⋅k25exp−EaR1T−1298.15⋅1−QKmn,where r refers to the reaction rate, Ar is the reactive surface area, k25 refers to the reaction rate constant at 25°C, Ea is the activation energy, T and R are the temperature (K) and ideal gas constant, Q and K refer to the ion activity product and equilibrium constant of a specific mineral reaction, and m and n 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 (Figure3(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.Figure 3
(a) Model domain of the conductive fault F2 and physical parameters, where land surface elevation extracted from digital elevation model, and the depth of F2 from the electromagnetic and gravity prospecting [19]; (b) irregular mesh generated by Gambit software (Fluent Inc.) with mesh refinement in the recharge and discharge zones; (c) boundary conditions for flow and heat in the fault F2 in the numerical model.
(a)
(b)
(c)As shown in Figure3, 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]:
(6)Tx,t−Ti=Tf−Ti+px+qx2exp−xd,where x is the distance from the boundary, Ti is initial temperature in surrounding rock, Tf is the time-varying temperature at the boundary, p and q are time-varying fit parameters, and d 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 (Figure4(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).Figure 4
(a) Steady-state temperature distribution and fluid flow field. (b) Calculated temperatures fit well with the observed temperatures in ZR1. (c) Gas saturation (Sg) distribution illustrates the existence of unsaturated zone at the top of the model. (d) The calculated chemical components in discharge area fit well with observations (ZR1, G1, G2, G3, and G4). The value of “HCO3-” is the sum of HCO3- and CO32-.
(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 (Figure4(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 Figure4(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 Figure5 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:
(7)NaAlSi3O8albite+H++0.5H2O=0.5Al2Si2O5OH4kaolinite+Na++2SiO2aqKAlSi3O8K−feldspar+H++0.5H2O=0.5Al2Si2O5OH4kaolinite+K++2SiO2aqFigure 5
Changes of mineral abundance in volume fraction (dimensionless) after a simulated period of 8000 years, with positive values indicating precipitation and negative values indicating dissolution.
(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:
(8)1.77Al2Si2O5OH4kaolinite+4.4SiO2aq+0.29Ca2++0.52Mg2+=2Ca0.145Mg0.26Al1.77Si3.97O10OH2smectite−Ca+1.62H++0.73H2O,(9)2.3Al2Si2O5OH4kaolinite+2.4SiO2aq+1.2K++0.5Mg2+=2K0.6Mg0.25Al2.3Si3.5O10OH2illite+2.2H++1.5H2OFigure 6
Chemical components in fluid (mg/L) and pH value after a simulated time period of 8000 years.
(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 Figure6(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 H+:
(10)KAl3Si3O10OH2muscovite+10H+=K++3Al3++3SiO2aq+6H2Owhich 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
(11)Ca2++HCO3−+OH−=CaCO3calcite+H2OCa2++Mg2++4HCO3−=CaMgCO32dolomite+2H2O+2CO2calcite 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 Figure7.Figure 7
(a) Calculated temperatures along fast upflow paths at the position of ZR1, (b) the variation of mineral saturation indices, (c) Na+, K+, Ca2+, Mg2+, Al3+, HCO3-, Cl-, SO42-, and SiO2 concentrations, and (d) mineral abundance. The mineral abundance of kaolinite, chlorite, smectite-Ca, and dolomite maintains unchanged which is not shown in (d).
(a)
(b)
(c)
(d)The SI values of illite, K-feldspar, and quartz are larger than “0” (Figure7(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 depths>150m. However, calcite precipitates near the land surface (with depth<150m) 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 (Figures7(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.Figure7(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 (Figures8(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)).Figure 8
Relative variation of porosity and permeability within the fault F2 after 8000 years. Negative value means the decrease of porosity and permeability due to mineral precipitation.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.


## 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 Figure5 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:
(7)NaAlSi3O8albite+H++0.5H2O=0.5Al2Si2O5OH4kaolinite+Na++2SiO2aqKAlSi3O8K−feldspar+H++0.5H2O=0.5Al2Si2O5OH4kaolinite+K++2SiO2aqFigure 5
Changes of mineral abundance in volume fraction (dimensionless) after a simulated period of 8000 years, with positive values indicating precipitation and negative values indicating dissolution.
(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:
(8)1.77Al2Si2O5OH4kaolinite+4.4SiO2aq+0.29Ca2++0.52Mg2+=2Ca0.145Mg0.26Al1.77Si3.97O10OH2smectite−Ca+1.62H++0.73H2O,(9)2.3Al2Si2O5OH4kaolinite+2.4SiO2aq+1.2K++0.5Mg2+=2K0.6Mg0.25Al2.3Si3.5O10OH2illite+2.2H++1.5H2OFigure 6
Chemical components in fluid (mg/L) and pH value after a simulated time period of 8000 years.
(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 Figure6(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 H+:
(10)KAl3Si3O10OH2muscovite+10H+=K++3Al3++3SiO2aq+6H2Owhich 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
(11)Ca2++HCO3−+OH−=CaCO3calcite+H2OCa2++Mg2++4HCO3−=CaMgCO32dolomite+2H2O+2CO2calcite 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 Figure7.Figure 7
(a) Calculated temperatures along fast upflow paths at the position of ZR1, (b) the variation of mineral saturation indices, (c) Na+, K+, Ca2+, Mg2+, Al3+, HCO3-, Cl-, SO42-, and SiO2 concentrations, and (d) mineral abundance. The mineral abundance of kaolinite, chlorite, smectite-Ca, and dolomite maintains unchanged which is not shown in (d).
(a)
(b)
(c)
(d)The SI values of illite, K-feldspar, and quartz are larger than “0” (Figure7(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 depths>150m. However, calcite precipitates near the land surface (with depth<150m) 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 (Figures7(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.Figure7(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 (Figures8(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)).Figure 8
Relative variation of porosity and permeability within the fault F2 after 8000 years. Negative value means the decrease of porosity and permeability due to mineral precipitation.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 mThe 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.



---

*Source: 1853068-2019-04-23.xml*