Geofluids
Volume 2019 (2019), Article ID 3420142, 19 pages
https://doi.org/10.1155/2019/3420142

Research Article

Mineral Filling Pattern in Complex Fracture System of Carbonate Reservoirs: Implications from Geochemical Modeling of Water-Rock Interaction

Jianwei Feng,¹ Xu Zhang,² Peng Luo,³ Xizhe Li,⁴ and He Du¹

¹School of Geosciences, China University of Petroleum (East China), Qingdao, China
²School of Petroleum and Natural Gas Engineering, Chongqing University of Science & Technology, Chongqing, China
³Mining and Energy Division, Saskatchewan Research Council, Regina, Canada
⁴Langfang Branch, China Petroleum Exploration and Development Institute, Langfang, China

Correspondence should be addressed to Peng Luo; peng.luo@src.sk.ca

Received 19 February 2019; Revised 10 June 2019; Accepted 6 September 2019; Published 15 November 2019

Guest Editor: Victor Vilarrasa

Copyright © 2019 Jianwei Feng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Much research has been conducted on physical and numerical modeling that focus on stress state and structural controls on subsurface geofluid flow, yet very few attempts have been made to discover and quantify the mineral precipitation/dissolution kinetics in complex fracture system such as Tarim Basin of China. We conducted a geochemical simulation study using the outcrop fracture networks in Ordovician carbonate rocks in Tabei Outcrop Area of Tarim Basin. Structural analysis, filling analysis within the fracture networks and surrounding rocks were used to constrain the generation and geochemical evolution of the geofluids. Using an advanced reactive transport simulation platform TOUGHREACT, a pertinent thermodynamic system was applied to establish the geological model of the fracture-surrounding rock, where the corresponding calcium carbonate (CaCO₃) solution was configured to replace the deep saturated hydrothermal fluids. Different types of mineral parameters were considered with material balance and phase equilibrium calculation to perform numerical simulation of multi-field, e.g., pressure field, temperature field, seepage field and chemical field under formation conditions. The simulation results were consistent with field observations. The major findings of this simulation study include: (1) Along with fluid injection, local dissolution occurred within the fractures and matrix, but with the gradual saturation of calcium ions and the increasing pH value, considerable calcite precipitation occurred. (2) The dissolution/precipitation in different fractures was mainly affected by their structure and physical properties, resulting in changes in fluid flow rate, temperature, pressure and ion concentration over time. (3) In the same group, the degree of mineral filling of small-aperture fractures, low-angle fractures and shallow fractures was significantly higher than other types of fractures. (4) The better the connectivity between reticular fractures and the higher the linear density of fractures, the lower the mineral filling degree. (5) Dissolution phenomenon strengthened within large-aperture conjugated fractures gradually along the flow direction. The proposed methodologies in this study can be applied to model effective fracture filling of other deep reservoirs.

1. Introduction

In recent years, a new round of oil and gas resources evaluation in China shows that the middle-depth and deep carbonate rocks have become one of the important areas to increase reserves. Significant progress in exploration and development has been seen in Sichuan, Tarim and Bohai Bay Basins [1–4]. Due to the low porosity and permeability characteristics, widely developed fractures and dissolution pores are the most important reservoir space and seepage channel for carbonate rocks [5, 6]. However, fracture conductivity is largely affected by precipitation and dissolution of minerals transported by geothermal geofluids, which have suffered from multi-cycle karstification and tectonic activity [7–10]. Therefore, how to effectively evaluate the spatial distribution of filling and dissolution degree has become a key constraint for fracture prediction and modeling of carbonate rocks. Commonly practiced in fractured reservoir exploration, modeling or prediction of effective fractures (i.e., filling degree) is typically based on geometric models, such as qualitative description/analysis of field outcrops, drilled cores and thin sections [11–14], seismic techniques or logging methods [15], and various mathematical methods (e.g., Kriging interpolation, sequential Gaussian conditional simulation). These approaches, aiming to acquire the attributes of inter-well fracture networks, are useful as they utilize actual measured data and relate deformation to corresponding structural position. While most of them applie ideal geometric models or simplified assumptions, they cannot completely reflect the multi-phase precipitation of hydrothermal minerals and of fluids properties of underground rocks.

Physical model experiments and numerical simulation on water-rock geochemical interaction are powerful tools for investigating underground fluid flow in fault and fracture systems and exploring the mechanisms of hydrothermal fluid mineralization and dissolution, which can be induced by earthquakes, volcanism, deformation, and sediment remobilization [13, 16–19]. As one of the pioneer studies, Ramsay [20] discovered that the crystals making up the extension vein filling in many naturally deformed crustal rocks often show a fibrous structure and seem to be built up by a succession of “crack-seal” increments: the elastically deforming rock fails by fractures, and the walls of the open micro-crack are sealed together by crystalline material derived by pressure solution in the rock matrix. Applying advanced techniques (e.g., Scanning electron microscope-cathodoluminescence analysis, CT imaging technique, fluid-inclusion observations, and carbon and oxygen isotope records), some researchers [21–26] have analyzed cryptic textures, paleofluid compositions and migration path in fracture filling veins to gain insight into timing and physical conditions of mineral growth, especially in environments with multiple hydrothermal precipitation events. Further, many laboratories [27–31] have conducted batch water-rock interaction experiments to study the dissolution and filling mechanism of underground hydrothermal fluid, which further reveal the growth law of various minerals in faults and fractures and establish the reaction rate models for deep hydrothermal fluids.

Numerous studies have shown that the different flow velocity, saturation, pathways and temperatures of mineralizing fluid lead to different patterns of mineral dissolution/crystallization in different structural systems and environments [32–35]. Also, through the finite/discrete element modelling simulations of multi-phase fluid fields, researchers [12, 13, 36, 37] concluded that underground fluids tend to flow along fractures that are oriented parallel or nearly parallel to current maximum horizontal compressive stress under hydrothermal conditions of<3 km depth. On the other hand, at depth of more than 3 km where effective stress is compressive and fractures are expected to be closed, chemical alteration dictates the location of open conduits, either preserving or destroying fracture flow pathways regardless of their orientation. In a fluid flow model of crack-sealing, intergranular pressure solution and compaction take place around active faults in partially sealed fault zones. Dobson et al. [38] and Zhang et al. (2011) focused on structural controls (i.e., structural deformation, connectivity, complexity of ault and fracture network) on hydrothermal fluid flowand mineralization. Based on a thermodynamic model for H₂O-CO₂-NaCl-CaCO₃ system, Zhu et al. [35] simulated the interactions between CO₂-rich fluid and carbonate rocks. They found that the changes of CO₂ concentration, brine salinity, temperature and pressure affect actual characteristics of dissolution-filling during hydrothermal fluid migration from deep to shallow stratum of Tarim Basin.

As summarized above, the flow and mineralization of geothermal fluids in subsurface structures have been investigated in various studies, the majority of which agreed with strong connection between structural complexity (e.g., fault and fracture networks) and flow, precipitation, and dissolution of geothermal fluids. Nevertheless, there has been very limited research on spatial flow pattern for precipitation and dissolution of deep hydrothermal fluid in complex fractured reservoirs. In this study, we adopted a numerical modelling approach to analyze the controlling factors for precipitation and dissolution and to further track the 4-D evolution of mineral filling within complex fractured reservoirs in long term. A 3-D discrete conceptual model was built to include rock matrix and fractures at different scale. Using the coupled fluid flow and multi-mineral chemical reaction model, we studied the effects of several structural controls (i.e., fracture density, occurrence, aperture and connectivity) on spatial filling pattern in fracture networks, and further quantified this relationship for fractured carbonate reservoirs.

2. Geochemical Modeling

2.1. Geochemical Model

TOUGHREACT is a numerical simulation program for chemically reactive nonisothermal flows of multiphase fluids in porous and fractured media [39–42]. The program was written in Fortran 77 and developed by introducing reactive geochemistry into the multiphase fluid and heat flow simulator TOUGH2 (Pruess, 2004). The program can be applied to one-, two- or three-dimensional porous and fractured media with physical and chemical heterogeneity. After more than a decade of hard work, Lawrence Berkeley National Laboratory has developed many new capabilities, where Xu et al. [43] has incorporated these new capabilities into Version 2.0 of TOUGHREACT.

Compared to previous versions, TOUGHREACT 2.0 adds some Fortran-90 extensions to the programming. The code can accommodate any number of chemical species present in liquid, gas, and solid phases. A variety of subsurface thermal, physical, chemical, and biological processes are considered under a wide range of conditions of pressure, temperature, water saturation, ionic strength, pH, and Eh [43]. TOUGHREACT 2.0, which is based on TOUGH2, incorporates a chemical reaction section that still uses a modular design. Due to the complex and varied geological conditions during diagenesis, EOS1, EOS7 and ECO2N modules are commonly used in numerous studies [43]. The EOS1 module is the most basic module, considering only the state of the water-steam two phases, and can simulate the basic water-rock reaction process. The EOS7 module considers the three components of water-salt-air, taking into account the gravity infiltration and diffusion, and is suitable for the simulation of water-rock reaction in marine sedimentation process. ECO2N is based on Spycher and Pruess’s research specifically for CO₂ geological storage, which can accurately describe the multiphase fluid migration process of water and CO₂ mixture in salt water [44]. Combined with the karst development in Tabei Outcrop Area, the ECO2N module was selected for this numerical simulation of water-rock interaction. During spatial discrete process, TOUGHREACT uses the IFD (Integral Finite Difference) method [45] which can handle irregular meshing and multiphase fluids in homogeneity/heterogeneous isotropic/anisotropic pores, migration in fractured media and water-rock reaction. The solution method in TOUGHREACT is implicit time weighting, and the coupling method between water flow, solute migration and chemical reaction module is sequential iterative method.

As observed from the carbonate formations of the northern outcrop area of Tarim Basin in northwestern China (i.e., Tabei Outcrop Area, Figure 1), a typical fractured vuggy reservoir contains a large number of different-scale fractures and lower-permeability rock matrix, as well as many scattered cavities or vugs. Tabei is located in the northwestern part of the Tarim Basin, including the two major tectonic units of Keping Fault Uplift and Bachu Uplift. The Ordovician carbonate strata are widely distributed in the area.. At the Keping Fault Uplift, the outcrop area extends in the south and southeast direction of the multi-row mountain system, and it is discontinuously distributed in the northwest part of the Bachu Uplift, mainly in the high parts [27]. As shown in outcrops of the reservoir formations in Figure 2, the fractures and vugs are irregular in pattern with very complicated in distribution, whereas the fractures are filled with various types of calcite, less commonly, dolomite and gypsum. Statistics of field observation showed that many of the small-aperture and low-angle fractures appear to be more easily filled with sparry carbonate cements, moreover, those fractures approximately parallel to the direction of the fluid are almost completely cemented. Several conceptual models for fracture fillings are shown in Figures 2(a)–2(c): (1) large-scale fracture network (model 1) are mainly half-filled or unfilled and directly connected to dissolved vugs; (2) small-scale fracture network (model 2) are mainly full-filled or half-filled and well connected; (3) low-angle fractures are usually filled significantly higher than the first two; (4) secondary associated fractures are the end of mineral-containing fluid channels and filled to a half degree; and (5) small-scale scattered fractures are isolated or separated from other fractures by pores of rock matrix, and commonly unfilled or slightly filled with sryptocrystalline calcite. Obviously, fractures cannot be modeled with uniform size distribution patterns in this paper.

Figure 1: Generalized map showing the location and the distribution of tectonic units in the Tabei Outcrop Area.

(a)

(b)

(c)

Figure 2: Schematic of conceptualizing fracture filling formation as a fracture-cavity-matrix system with well different-level mineral veins. (a₁-a₆) Profile outcrop pictures, (b₁-b₂) ground outcrop pictures and (c) conceptual model.

Different from the conventional dual-porosity and double-permeability concept models [46], the large-scale fracture or fractured space embedded in the rock matrix is conceptualized as main global flow channel. While low-porosity and low-permeability matrix continua, still mainly providing storage space as sinks, are indirectly or directly interacting with widely connecting fractures. The conceptual discrete model utilized in this study assumes that approximate thermodynamic equilibrium exists locally in each of the two independent units at all times [47]. On basis of the local equilibrium assumption, the thermodynamic variables, such as fluid compositions, temperature, pressures and saturation, for each continuum can be defined successfully.

As shown in Figure 3, to quantify the structural controls on subsurface geofluid flow and variations of mineral precipitation/dissolution within complex fracture networks, we designed four sets of discrete models (i.e., models with variable aperture, dip, depth, direction fractures) to simulate water-rock interactions. According to the outcrop observation, the degree of filling of the fractures changed greatly within a few meters, as a result, we set the conceptual model at scale of meter. Take single-fracture models with different aperture as an example, as shown in Figure 3, we defined a 3-D 1000 mm × 200 mm × 200 mm flow region containing three different-aperture fractures and one matrix segment. Each fracture was assumed to have an average aperture of either 5, 15, or 20 mm, parallel to the x-axis or fluid injection direction.

Figure 3: Schematic diagram of three-dimensional multiple geological model.

2.2. Fluid System

The extensively developed intra-fracture fillings (e.g., calcite and dolomite) and dissolution vugs/cavities indicate that the paleo-fluid system in the northern outcrop area of Tarim Basin is mainly characterized by thermal fluids caused by deep magma [27, 48, 49]. The ascending hydrothermal system generally consists of high-temperature water (H₂O). Carbon dioxide, hydrogen sulfide and possible methane released from hydrocarbon reservoirs (mainly considered as a gas reservoir), as well as of rock fragments, are dissolved into the water phase. To ensure the reliability of subsequent simulation results, in this study we applied a simplified fluid thermodynamic model for carbonate reservoirs, i.e., CO₂-H₂O hydrothermal system.

From the observed aqueous fluid inclusions in calcite veins of fractures of northern outcrop area of Tarim Basin, the estimated homogenization paleotemperature mainly range from 105–120°C and 125–160°C [27, 35], which indicates the characteristics of formation temperature changes since the Ordovician stratigraphic deposition. At the same time, based on the liquid inclusion test results in the study area, the estimated homogenization paleotemperature of the deep CO₂-rich hydrothermal fluid mainly concentrate in two intervals: 180–235°C and 250–330°C [27, 35, 48], which indicates two-phase hydrothermal mineral filling process since the Late Caledonian Period. Simultaneously, the entrapment paleopressures have been estimated from published PVT simulation results for aqueous inclusions [50, 51], which average 80 MPa and 150 MPa, respectively, recording two-stage two-phase tectonic hydrothermal episodes. According to previous studies, the paleotemperature gradient of the Paleozoic in Tarim Basin is about 3°C/100 m, and the paleosurface temperature is 20°C [27, 52]. Combined with the burial history curves of northern area of Tarim Basin [52], the two depth intervals of the inclusions are 2800–3300 m and 3500–4100 m, therefore the corresponding average formation paleopressure was approximately 75 MPa and 135 MPa, respectively.

As a result, initial reservoir temperature and pressure were set to 125°C and 105 MPa, respectively. An over-pressure of 10 MPa was applied to the injection (left) side, and the injected water and gas (CO₂) temperature was taken as 255°C (Figure 3). In TOUGHREACT simulator, the gas (CO₂) was set to be dissolved in water, i.e. water-soluble. It could be considered that the fluid type in this simulation belonged to single-phase flow. The boundary condition is given on the boundary of a seepage zone to express the physical condition at the boundary of the seepage zone, that is, the condition that the water head (or seepage flow) should meet on the boundary of the seepage zone. According to the current geological features, for the horizontal flow model (Figure 3), its upper and bottom surfaces (i.e. surfaces ADEF and BCHG) were both water-blocking layers, which belonged to second type of boundary conditions (i.e. Neumam condition). For the second type, the corresponding governing equation was: , . Here, K was permeability coefficient of the study area; H was water head on the top and bottom, and here referred to the underground heat flow; q (x,y,z,t) was flow into or out of the boundary per unit area and unit time, as a given quantity; Since this water-blocking boundary belonged to an isotropic medium rock, the expression could be simplified to: , indicated that the fluid flow was zero at this boundary. The south and north boundaries (surfaces CDEH and ABGF) were the first type of boundary condition, i.e., the constant pressure boundary (Dirichlet condition) (Figure 3). For this condition, the corresponding governing equation was: ,. Here, was a function of boundary flow as a function of time. Although the water point at each point on this boundary was constant at each moment, the boundary was in communication with the outer formation, and the fluid could still flow as long as the pressure was large enough. The east and west boundaries (surfaces ABCD and EFGH) were the flow boundary, that was, the second type of boundary (Figure 3). In the model, we further set the deep hydrothermal fluid to flow from the left (west) side to the right (east) side, and the CO₂ injection rate was 0.005 kg/s.

For the vertical model (Figure 4(a)), four horizontal boundaries (surfaces ADHE, BCGF, ABCD and EFGH) were the first type of boundary, that was, when the fluid pressure was large enough, the fluid could communicate with surrounding rocks. However, the top and bottom boundaries were set to the second type or flow boundaries, that was, the fluid moved from deep to shallow at an initial injection rate of 0.005 kg/s. Table 1 shows initial water chemical elements. The geothermal fluid injection was assumed to last for a period of 100 years. In order to observe the chemical reaction during the fluid flow, the simulation of fluid flow and geochemical transportation was also set to be a period of 100 years.

(a)

(b)

(c)

(d)

Figure 4: Spatial distribution of mineral abundance, CO₂ concentration (aqueous phase) and formation pressure from simulations. (a) Three-dimensional conceptual model of deep hydrothermal fluid; (b) profile changes of redistributed calcite abundance at different injection stages; (c) profile changes of CO₂ concentration (aqueous phase) in aqueous phase over time; (d) profile changes of formation pressure over time.

Table 1: Initial water chemical characteristic.

Water type
Average of pH
Chemical characteristic
Concentration/(Mol·kg⁻¹)
Chemical characteristic
Concentration/(Mol·kg⁻¹) /(Mol·kg⁻¹)

CaCl₂-type
5.8
K⁺
1.577 × 10⁻⁴
Al³⁺
3 + 6.490 × 10⁻⁶

Na+
1.950 × 10⁻²
SO2⁴⁻
7.651 × 10⁻¹¹

Ca²⁺
9.270 × 10⁻¹
HCO⁻
9.149 × 10⁻³

Mg²⁺
2.170 × 10⁻¹⁴
Cl⁻
9.897 × 10⁻¹

Fe²⁺
3.015 × 10⁻³
SiO₂(aq)
1.378 × 10⁻³

2.3. Reservoir Properties

Core and thin section observations indicated that the main lithologies of the Ordovician strata in the study area included pelsparite, biogenic limestone, and some dolomitic limestone, mixed with other small amounts of clay and clastic rocks. In this study, for simplicity purpose, we considered the carbonate rock is composed of 10 vol% Calcite (CaCO₃) + 90 vol% CaMg(CO₃)₂. The PetraSim-2 platform of TOUGHREACT simulator requires that values for intrinsic seepage parameters (e.g., rock density ρ, porosity ϕ, permeability K) and wet heat conductivity (CWET) and specific heat (SPHT) be assigned or calculated at nodes located at the center of each computational cell [53]. A total of 26 porosity and 39 permeability measurements were conducted on carbonate core samples from the outcrops, as listed in Table 2. Porosity measurements of rock matrix ranged from 0.92 to 7.48%, with a mean value of 4.2%. As shown in Table 2, the horizontal permeability values of these samples (fracture samples and matrix samples) were all significantly higher than the vertical permeability, but still by an order of magnitude. Default values of the thermal conductivity of rocks set by the simulator were used, as the focus of this study is fluid movement and water-rock interactions. In addition, default values in the software were used for some hydrogeological and thermodynamic parameters (e.g., krl, Slr, Slg, and Pcap/Pa) that are difficult to obtain by conventional methods (Table 2).

Table 2: Hydrogeological and thermodynamical parameters used in the simulation (Note: Sl in the krl, krg, and Pcap columns refers to liquid saturation).

Parameter
Magnitude
Unit

Matrix density
2610
Kg/m³

Matrix porosity
4.2
%

Matrix vertical permeability
1.5 × 10⁻¹⁴
m²

Matrix horizontal permeability
5 × 10⁻¹³
m²

Matrix CWET
2.51
W/(m·K)

Matrix SPHT
920
J/(kg·K)

Fracture density
2200
Kg/m³

Fracture porosity
10
%

Fracture vertical permeability
1.5 × 10⁻¹²
m²

Fracture horizontal permeability
5 × 10⁻¹¹
m²

Fracture CWET
2.10
W/(m·K)

Fracture SPHT
800
J/(kg·K)

Relative permeability

Liquid [71]:

Irreducible water saturation Slr
Slr = 0.3 Slr = 0.3

Exponent
m = 0.457

Gas [72]:

krg = (1-S)²(1-S²)
S = (Sl-Slr)/(Sl-Slr-Sgr)

Irreducible gas saturation
Sgr = 0.05

Capillary pressure

Van Genuchten(1980)

Irreducible water saturation
Slr = 0.03

Exponent
m = 0.457

Strength coefficient
P₀ = 19.61 kPa

2.4. Water Chemistry

Minerals filter the groundwater and change the content of their ionic components during ion exchange with minerals [33]. Many researches have been carried out to investigate the geochemical composition and genesis of formation water of the representative Tahe oil field of northern outcrop area in Tarim Basin [54–56]. According to the classical scheme of Surin, the groundwater of Ordovician stratum is mainly the calcium chloride-type (CaCl₂) of oil field water, with average pH of 5.8. As a result, in this modelling, a simplified CaCl₂-type water model including calcium ions (Ca²⁺) and chloride ion (Cl⁻) was set as the initial injected hydrothermal fluid. The native reservoir water has a Ca²⁺ ionic concentration equal to 0.927 kg/mol, and a Cl⁻ ionic concentration of 0.9897 kg/mol, with a pH of 5.8 (Table 1). Finally, after generalize the aquifer hydraulic characteristic, vertical boundary and lateral boundary, the groundwater mathematical model has been established on the basis of the hydrogeology conceptual model.

2.5. Reaction Kinetics Model

Since hydrothermal geofluids generally have very high temperatures and carry a lot of acid gases such as CO₂ and CH₄, this acidic environment will undergo significant water-rock dissolution/precipitation reactions with the surrounding carbonate rocks (Zhu et al., 2012); [57, 58]. The kinetic rate expression used by the software in mineral dissolution and precipitation processes is from Lasaga [59, 60]:

where n is mineral code for dissolution/precipitation reaction; rn is dissolution/precipitation reaction rate, positive value indicates dissolution and negative value indicates precipitation; K is reaction rate constant; An is reactive specific surface area for minerals in 1 kg of water; Ω is saturation index; θ and ŋ are two parameters that are generally determined experimentally and usually set to unity by default. The calculation formula of the reaction rate constant used in the software is as follows [61]:

where K₂₅ is the reaction rate constant at 25°C (Table 3); Ea is reaction activation energy; R is gas constant; α is reaction activity; nu indicates neutral mechanism; H indicates acidic mechanism; OH indicates n alkaline mechanism. The chemical reaction kinetic parameters refer to research results of Xu et al. [61, 62].

Table 3: Parameters for calculating kinetic rate constants of mineral.

Minerals
A/(cm²/g)
Neutral mechanism
Acid mechanism
Base mechanism

k₂₅ (Mol/m²/s)
Ea (KJ/Mol)
k₂₅
Ea
n(H⁺)
k₂₅
Ea
n(H⁺)

Quartz
9.8
1.023 × 10⁻¹⁴
87.7

Low-albite
9.8
2.754 × 10⁻¹³
69.8
6.918 × 10⁻¹¹
65
0.457
2.512 × 10⁻¹⁶
71
-0.572

Anorthite
9.8
7.586 × 10⁻¹³
17.8
3.162 × 10⁻⁴
16.6
1.411

K-feldspar
9.8
3.890 × 10⁻¹³
38
8.710 × 10⁻¹¹
51.7
0.5
6.310 × 10⁻²²
94.1
-0.823

Calcite
9.8
Assumed at equilibrium

Chlorite
9.8
3.020 × 10⁻¹³
88
7.762 × 10⁻¹²
88
0.5

Dolomite
9.8
2.951 × 10⁻⁸
52.2
6.457 × 10⁻⁴
36.1
0.5

Siderite
9.8
1.260 × 10⁻⁹
62.76
6.457 × 10⁻⁴
36.1
0.5

Ankerite
9.8
1.260 × 10⁻⁹
62.76
6.457 × 10⁻⁴
36.1
0.5

Kaolinite
151.6
6.918 × 10⁻¹⁴
22.2
4.898 × 10⁻¹²
65.9
0.777
8.913 × 10⁻¹⁸
17.9
-0.472

Illite
151.6
1.660 × 10⁻¹³
35
1.047 × 10⁻¹¹
23.6
0.34
3.020 × 10⁻¹⁷
58.9
-0.4

Na-smectite
151.6
1.660 × 10⁻¹³
35
1.047 × 10⁻¹¹
23.6
0.34
3.020 × 10⁻¹⁷
58.9
-0.4

Ca-smectite
151.6
1.660 × 10⁻¹³
35
1.047 × 10⁻¹¹
23.6
0.34
3.020 × 10⁻¹⁷
58.9
-0.4

2.6. Porosity–Permeability Correlation Model

The aquifer porosity and permeability are constantly changing due to the dissolution of the initial minerals and the precipitation of secondary minerals [63]. The change in porosity is directly calculated from variation in volume fraction of minerals, and the change in permeability due to the variation in porosity can be calculated by the Kozeny-Caeman sphere particle model [41]. The specific expression is:

where K₀ is the initial permeability of rock; K is the permeability after chemical dissolution or filling; is the initial porosity of rock; is the porosity after chemical dissolution or filling. The model of porosity-permeability selected for the simulation process is relatively simple [63]. Though it does not consider the influence of various factors such as the size, shape, and connectivity of the pores, its complexity is sufficient for the purpose of this simulation study.

3. Results and Discussion

3.1. Filling Pattern in Fractures with Different Apertures

The flow rates and mineral filling/corrosion degree over time are presented in Figure 5. By increasing the fracture aperture at the flow direction from 5 mm to 30 mm, the rate of fluid within the model is increased compared to the initial injection rate. After 100 years of continuous injection, the formation pressure within the entire model gradually increased and reached a new equilibrium, which was slightly higher than the initial formation pressure of 105 MPa. However, when deep fluids migrated into fracture zones, the pressure dropped dramatically along direction of flow, approaching 108 MPa (Figure 5(a)). Further reviewing Figure 5(b), we could see that flow rates and temperature increased faster within the fractures compare to the reservoir matrix site, due to a larger permeability. Many previous studies ([64–67]; Duan et al., 2008) showed that CaCO₃ solubility decreases with fluid temperature increases, which means the increase of CaCO₃ precipitation around the fractures. After 10 years, the 255°C injection fluid was over-saturated with calcite minerals and under-saturated with respect to dolomite minerals because this fluid had not yet reached equilibrium with the reservoir rock. Therefore, at this moment, calcite slightly dissolved adjacent to the injection side because calcite solubility increased with relatively low temperature and gradually increasing flow rate (Figure 5(c)), which was different from dolomite. After 50 years, as the reservoir temperature rose and pressure declined, the chemical reaction in fractures changed from dissolution to precipitation gradually. Calcite has a higher filling tendency (volume fraction) at the fracture within aperture of 5 mm. As the fracture aperture increases, the degree of mineral filling was significantly reduced, but the filling range was relatively increased (Figure 5(d)). Interestingly, after 100 years, the difference of mineral precipitation in this fracture system became more apparent (Figure 5(e)). A calcite precipitation peak could be observed in the fracture with an aperture of 5 mm because the precipitation rate increases with temperature. However, in the early stage, due to the different reaction kinetics [33, 68], along the flow path, the fluid conditions changed continuously and the dolomite was slightly dissolved in the matrix, but a slight precipitation occurred in the large-aperture fracture, causing competing effects on calcite solubility (Figure 5(f)). In comparison, with the continuous injection of fluid, the dissolution of dolomite in the fracture was more obvious, and the larger the aperture, the stronger the degree of dissolution (Figure 5(g), 5(h)).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 5: Spatial distribution of geofluid flow and mineral abundance during simulations of different-aperture fracture networks (negative values indicate dissolution and positive precipitation). (a) Plans of reservoir pressure variation at different times; (b) plans of flow rate and temperature variations (length of the arrow representing flow rate, that is, the longer the arrow, the faster the flow rate); (c) changes of abundance of calcite in plans after 10 years; (d) changes of abundance of calcite in plans after 50 years; (e) changes of abundance of calcite in plans after 100 years; (f) changes of abundance of dolomite in plans after 10 years; (h) changes of abundance of dolomite in plans after 50 years; (h) changes of abundance of dolomite in plans after 100 years.

In order to further explore the filling mechanism in fractures, we extracted the solution ion curves in different time and different distance range (Figure 6). When the deep CO₂-rich dissolved fluid (i.e., acid solution) entered the CaCl₂-type formation water (i.e., weak acid fluid), the native chemical equilibrium would be broken, resulting in a linear decrease in Ca²⁺ concentration and a gradual increase in Mg²⁺ concentration, and continuously rising of pH in fractures (Figures 5(a) and 5(b)). Therefore, at the early stage, due to the high temperature and low pH of the injected fluid, the calcite and dolomite dissolution occurred constantly throughout the matrix and fractures, until Ca²⁺ and Mg²⁺ was saturated or over-saturated. At this time, the reactions of calcite and dolomite dissolution were expressed as.

(a)

(b)

(c)

(d)

Figure 6: Distribution of ion concentration, mineral abundance, and permeability obtained from six different simulations along section A-B. (a) Top view of 3-D model; (b) ion concentration and permeability change in points of fractures over time (C-D and E-F); (c) contrast of ion concentration and permeability changes between fractures after 100 years; (d) pH value of solution changes along section A-B over time.

As time went by, due to under-saturation of Mg²⁺, continuous dissolution of dolomite at the fractures promoted the continuous precipitation of calcite (Figure 6(b)). Zones of dolomite dissolution and calcite/dolomite precipitation moved gradually away from the injection point due to changes in temperature, flow rate and pH along the flow path. Moreover, key areas of calcite/dolomite precipitation gradually moved from rock matrix into fractures. Thus, the corresponding reaction of calcite/dolomite precipitation (new mineral) after heated was taken to be

In the geochemical modeling, the discrete pore-fluid flow was considered mainly taking place within the 3-D geological fracture zone. Undoubtedly, mineral dissolution and precipitation lead to changes in the porosity and permeability of the reservoir matrix and fractures (Xu et al., 2006); [69]. Permeability increases within fractures indicated that mineral dissolution was dominant at the early stage (Figure 6(b)), while permeability reduced when precipitation dominated in the later stage. The dissolution of dolomite and precipitation of dolomite played a certain role than other factors for calcite, giving rise to calcite redistribution in fracture zones. Nevertheless, the initial dolomite volume in carbonate reservoir did not exceed 10%, so its effect on the degree of fracture filling and permeability change was almost negligible.

Figure 6(c) shows Ca²⁺ and Mg²⁺ ion concentration and permeability changes perpendicularly across fractures after 100 years. The high-concentration zones of calcium ion (Ca²⁺) tended to accumulate in the reservoir matrix and the edges of fractures, indicating that the fracture zone was a favorable location for calcite precipitation. The reason for this phenomenon was that the surrounding rock of the carbonate rock was less impermeable than the fracture zone, and the vicinity of the contact surface became a typical flow transition boundary, such as flow rate, pressure and temperature. Moreover, the smaller the fracture aperture, the easier the calcite precipitation and the lower the total permeability of the fracture zone. In contrast, the temporal and spatial variation of magnesium ion (Mg²⁺) was exactly the opposite of the former, which of course had little effect on the change of fracture permeability. All of these indicated that due to the 3-D convective pore-fluid flow, the chemical reactions between the chemical species were finally promoted and accelerated within the 3-D geological fracture zones.

The curve of the pH value of the formation water was shown in Figure 6(e). It could be seen that in 5 × 10⁷ seconds, the pH of the formation water slowly increased to about 6.2. Before 1 × 10⁸ seconds, the pH value increased rapidly, reaching 8.0. However, after 1 × 10⁸ seconds, the pH value increased slowly and remained between 8 and 9. After 100 years, the pH of the formation water was basically stable. According to the A-B section, the pH value in the fracturewas significantly lower than that of the surrounding rock. For fractures with different apertures, when the aperture was larger, the lower pH value, the higher calcium ion concentration and the lower permeability could be observed.

3.2. Filling Pattern in Fractures with Different Dip Angles

According to the simulations of fractures with different dip angles, the difference of this mineral filling was more obvious (Figure 7). Calcite abundance had almost no difference in three fractures at the early injection stage. Native calcite also dissolved close to the injection side but precipitated later in farther region (Figure 7(a)). As temperatures increased away from the injection side, ions like HCO₃⁻ and Ca²⁺ were over-saturated andcalcite precipitation occurred, especially in the fracture zone. Zones of calcite precipitation move gradually away from the injection point due to changes in temperature along the flow path. A maximum of new 8% volume of new calcite had been precipitated in the low-angle fracture after 100 years, which was about 10% larger than the amount of vertical fracture deposition (Figure 7(b)). Notice that amounts of calcite precipitation and their redistribution depended on their precipitation kinetics. As a result, it could be seen that after a long geological history evolution, low-angle fractures were easily full-filled by minerals, and high-angle oblique fractures were easily half-filled with minerals, however, near-vertical fractures were prone to no filling.

(a)

(b)

(c)

Figure 7: Spatial distribution of geofluid flow and calcite abundance during simulations of fractures with different dip angles (negative values indicate dissolution and positive precipitation). (a) Changes of calcite abundance in cross section during fluid injection; (b) changes of calcite abundance of vertical fracture deposition for 200 years; (c) changes of calcite abundance of low-angle fracture deposition for 200 years.

3.3. Variation of Fracture Filling with Formation Depth

Generally, as the deep hydrothermal fluid migrates to the shallower depth, the corresponding temperature pressure gradually decreases as the depth decreases [70]. To investigate the effects of deep fluid transport on the calcite solubility, this study first assumed an ideal situation where fluid temperature and pressure were consistent with formation temperature and hydrostatic pressure in the formation. For simplicity, a three-dimensional 200 mm × 200 mm × 400 mm flow region containing continuum matrix and two different-scale fractures was defined (Figure 4).

As shown in Figure 4(b), interactive precipitation and dissolution occurred with depth due to the injection of CO₂ and H₂O. It was clearly observed that all variations in the abundance of calcite were focused within the fracture zones. Simultaneously, within the first 20 years, a chemical reaction boundary of precipitation/dissolution (approximately 80 mm from the top) appeared in the upper part of the model, above which strong dissolution dominated the process, while below which widely minor dissolution/precipitation dominated. Interestingly, since the calcite solubility is proportional to pressure (Duan and Li, 2008) and CO₂ (aqueous phase) concentration and inversely proportional to temperature, the distribution pattern of calcite was similar to that of these factors, which could be seen from Figures 4(c) and 4(d). It should be noted that by default TOUGHREACT does not consider pressure dependence of mineral solubilities. In this simulation study, we used the thermodynamic phase equilibrium model of Duan and Li (2008), from which the pressure versus solubility curve was derived by regression and added to the TOUGHREACT software. Mainly due to gravity segregation, the initial CO₂ concentration (aqueous phase) rapidly rose along the steep fracture zone to the top outlet, resulting in large area dissolution of calcite in the rock. However, it was not until 40 years that this dissolution/precipitation phenomenon changed significantly. As the CO₂ continued to accumulate in the aqueous phase, i.e., CO₂ concentration in the aqueous phase continued to increase, to the top, HCO₃⁻ and Ca²⁺ quickly reached saturated or over-saturated status, and the formation temperature gradually rose to a certain value. Therefore calcite began to precipitate. Nevertheless, the boundary at this time still existed and began to migrate down until it was 135 mm from the top, which was different from the boundary of the fracture zone (Figure 4(b)). This phenomenon became more distinct even after 100 years, in which the upper part of the formation was always the precipitation zone and the lower part was always the dissolution zone. As a final consequence, in the fracture region, maximum calcite precipitated at upper part, while maximum calcite dissolution took place within lower part but dissolved later (Figure 4(b)). Of course, from the depth-based simulation results, the size of the fracture at this time had slight effect on the dissolution and filling.

3.4. Filling Pattern in Fracture Network

Figure 8 shows the spatial redistribution of calcite volume fraction within the flow region for the mixed fracture networks for prescribed aperture of 5 mm and 15 m and uniform injection rate. Analysis of the dissolution/precipitation patterns further clarified the mechanism leading to differential filling patterns within fractures. In Patterns 1, 2, 3, and 4 (Figures 9(a), 9(b), 9(c), and 9(d)), the concentration distribution of precipitated calcite at the first 20 years were very similar to those obtained at 80 and 100 years, although the locations of maximum values of mineral abundance were slightly different. As shown in Figure 9(a), in complex fracture networks with constant apertures, calcite first precipitated in the facture parallel to the direction of flow path near the inlet site and then gradually migrated to the similar fracture in the farther zone. Moreover, as time went on, the high-value zones of mineral precipitation gradually migrated to or localized in the intersections of the fractures, which just showed that the degree of mineral filling was proportional to the fracture intensity (measurement for characterizing degree of development). In contrast, for different fracture combination, the high precipitation was focused in a few low permeability channels in variable-aperture fractures.

(a)

(b)

Figure 8: Spatial distribution of mineral composition varied with formation depth. (a) Changes of calcite abundance in profile during fluid injection (negative values indicate dissolution and positive precipitation); (b) changes of calcium concentration in cross section in profile during fluid injection.

Figure 9: Spatial distribution of mineral abundance for prescribed fracture patterns after simulations. (a) Changes of calcite abundance in cross section of Pattern 1 (net-shaped fractures) during fluid injection; (b) changes of calcite abundance in cross section of Pattern 2 (variable aperture); (c) changes of calcite abundance in cross section of Pattern 3 (variable angle); (d) changes of calcite abundance in cross section of Pattern 4 (variable connectivity).

In summary, the mineral filling degree could be calibrated to a full-filled type in fractures parallel to flow, the oblique small-aperture fractures and the intersections. Accordingly, the degree of filling at other locations might be defined as a half-filled or unfilled type. This indicated that small-aperture fractures that intersected at a low angle to the flow path not only transmitted sufficient ion concentration (i.e., HCO₃⁻ and Ca²⁺), but also maintained a higher temperature. Conversely, connecting large-aperture fractures transmitted water at sufficiently large rates, and it was difficult to provide a stable sedimentation environment, thereby mainly causing dissolution. In the same way, the secondary fracture, that was, the fracture in same direction as the main fracture, was more likely to precipitate minerals than the fracture diagonally inclined to the main fracture (Figure 9(b)). The reason was mainly that the secondary fracture along the direction of flow was easier to pass through more fluids to maintain higher temperatures and carried more precipitated ions than the reversed fracture.

For a fixed-angle combination, it might thus be concluded that those fractures with a small angle to the fluid transport direction (e.g., less than 30°) and near right angle were more likely to be filled, whereas those with medium angles had less filling. Generally, as time went on, more calcite was firstly precipitated in the region of up-temperature flow or major fractures (i.e., the yellow color region), then gradually transferred to secondary forward fractures and less calcite precipitated in the region of down-temperature flow or secondary reverse fractures or isolated fractures (i.e., the light green color region) (Figures 9(c), 9(d)).

3.5. Summary of Fracture Filling Patterns for Carbonate Reservoirs

As the numerical simulation results demonstrated above, fluid precipitation/dissolution process in a fracture system was a complex and subtle reactive transport process. We standardized the data and plotted them into a uniform cross-section map to generalize mineral filling pattern for different fracture configurations. As shown in Figure 10(a), under the preferable fluid and temperature conditions for carbonate reservoirs, the effect of fracture aperture on the calcite dissolution/precipitation was obvious, that was, the smaller the aperture, the higher the calcite filling degree (volume fraction), and vice versa. Certainly, this inherent relationship between precipitation and structural complexity was not constant. The curve showed that as the fracture aperture decreased toless than 10 mm (which could be considered as a dimensionless relative value), the mineral filling no longer changed, that was, tended to be stable. For the second type of curve (Figure 10(b)), the trend was very similar to the first type. In detail, another effect that had an influence on filling degree in the conceptual model was the variable dip angles of fractures under complex stress state. Obviously, in the process of fluid injection, as the inclination angle of the fracture increased, the amount of calcite precipitation or filling was continuously reduced until the inclination angle reached 90°, at which the filling was close to zero. However, unlike the former two, the variation of the latter two types of calcite precipitation/dissolution curves generally had a tendency of decreasing first and then increasing rapidly (Figures 10(c) and 10(d)).

(a)

(b)

(c)

(d)

Figure 10: Typical model curves that characterize the behaviour of precipitation/dissolution when prescribed geofluid flow within various fracture networks. (a) Change curve of calcite abundance along with variable aperture of fractures (coordinates can be dimensionless); (b) change curve along with variable dip angle; (d) change curve along with variable distance; (c) change curve along with variable angle (coordinates can be dimensionless); (d) change curve along with variable angle.

For the third curve model, with different fluid injection directions (Figure 10(c)) from the other three (Figures 10(a), 10(b), 10(d)), there was a clear chemical reaction interface, roughly in the upper half of the reservoir. Furthermore, vertically, in the matrix and fractures above the interface, calcite dissolution or a small amount of precipitation was dominant, while a large amount of precipitation was dominant under the interface (Figure 10(c)). The fourth curve model presenting as a “U” shape showed that the angles between fracture and flow direction also played an important role during chemical precipitation/dissolution processes in fracture networks. In detail, within range of 45–75° from the direction of flow, the fluid in fractures was less likely to precipitate, and thus the degree of filling was relatively low (Figure 10(d)). In contrast, fractures within ranges of 0–45° and 75–90° from the flow direction were relatively more prone to be filled with calcite minerals, and thus the degree of filling was relatively high. In addition, it could be seen that as the linear density of connected fracture network increased, the precipitation volume fraction or the degree of mineral filling decreased. All of this was mostly due to the structural complexity of fracture networks, which controlled on flow rate, local temperature and pressure changes, and variations of ion concentration, but those were not the only factors that contributed to this effect. As a result, the precipitation/corrosion laws obtained by numerical simulation were basically consistent with the phenomena observed in the outcrop of Tabei Outcrop Area.

4. Conclusion

Based on the outcrop observation in the northern area of Tarim Basin, we summarized the principal components to be ideally considered when modelling hydrothermal fluid flow in deep low-permeability carbonate reservoirs: the transported pores and their properties, the geometry of the fissure system, and the host rocks. Four sets of reactive transport models were developed to simulate the influence of structural complexity on subsurface preferential flow and precipitation/dissolution processes within different fracture combinations. Several findings were drawn from the simulation results:

(1)For the planar model, the pH change had obvious “three-stage characteristic” throughout the process of groundwater rock reaction (i.e., slowly increasing, increasing rapidly, and slowly increasing to a stable stage). In the early stage (pH< 6.2), local dissolution occurred within the fractures and rock matrix. But with the gradual increasing of calcium ions and pH value (pH value <9.1), there was obvious calcite precipitation in the fractures, and this process lasted until the end(2)During the whole process, dolomite was always in a state of dissolution until it reached 0.00195 mol/L, especially in fractures, which accelerated the precipitation of new calcite minerals. The redistribution of minerals (e.g., calcite, dolomite) led to changes in both porosity and permeability of the fractures, which generally manifested as a rapid decrease of porosity and permeability first (about 0.1 of original value) and then slowly increased until it stabilized (about 0.8 of original value).(3)For isolated fracture group, although the amount of mineral precipitation in the small-aperture fracture was relatively small, the filling degree was significantly higher than that of the large-aperture fracture. Specifically, for fractures having aperture of less than 10 mm, the degree of filling did not change much, and the volume fraction was about 0.008. When the fracture aperture exceeded 10 mm, the filling degree rapidly decreased and gradually became stable, and its volume fraction was about 0.002(4)Due to the vertical gravity differentiation and temperature change, the fractures in shallow stratum were prone to be severely filled, while the deep fractures were either dissolution-dominated or filled with a small amount. There was a clear precipitation/dissolution reaction boundary between them. In early stages, the boundary was roughly one-fifth of the total height from the top surface, and in the later stage, the boundary gradually migrated to one-third(5)For fracture network, the angle between fracture surface and direction of fluid flow and the degree of connectivity had the most significant effects on the degree of filling. Within range of 45–75° from direction of flow, the fluid in fractures was less likely to precipitate, and thus the degree of filling was relatively low. Connectivity, to some extent, had affected the mineral precipitation/dissolution boundary during the water–rock interaction. And dissolution phenomenon strengthened within large-aperture conjugated fractures gradually along the flow direction. In addition, in complex system, fractures with higher linear density had a relatively lower filling degree, which could provide more space for fluid to migrate without easily precipitating the carried ions

In other words, the precipitation/dissolution of minerals in the fracture system was largely affected by structural complexity and its properties, which directly led to the coupling changes of fluid flow rate, temperature, pressure and ion concentration in different forms of fractures. Because the developmental characteristics of fault and fracture systems were similar in various petroliferous basins of the world, geological fluid migration, mineral dissolution/precipitation mechanisms and laws under such formation conditions were also applicable to structural geology and fluid mechanics studies in other basins.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was financially supported by the National Oil and Gas Major Project (2016ZX05047-003, 2016ZX05014002-006), the National Natural Science Foundation of China (41572124, 41628201), the Fundamental Research Funds for the Central Universities (17CX05010), and Pilot project of Chinese academy of sciences (XDA140102).

References

A. Q. Chen, M. G. Hou, H. D. Chen, S. L. Xu, Y. J. Zhong, and W. J. Liu, “Marine carbonate exploration history and the basic driving action from sedimentology in China,” Acta Sedimentologica Sinica, vol. 35, no. 5, pp. 1055–1062, 2017.

D. S. Li, “China's multicycle superimposed petroliferous basins: theory and explorative practices,” Xinjiang Petroleum Geology, vol. 34, no. 5, pp. 497–503, 2013.

N. Li, C. W. Xiao, L. H. Wu et al., “The innovation and development of log evaluation for complex carbonate reservoir in China,” Sino-Russian Logging International Academic Exchange Conference, vol. 38, no. 1, pp. 1–10, 2014.

X. Shi, J. X. Dai, and W. Z. Zhao, “Analysis of deep oil and gas reservoirs exploration prospect,” China Petroleum Exploration, vol. 1, pp. 1–11, 2005.

W. L. Ding, T. L. Fan, B. S. Yu, X. B. Huang, and C. Liu, “Ordovician carbonate reservoir fracture characteristics and fracture distribution forecasting in the Tazhong area of Tarim Basin, Northwest China,” Journal of Petroleum Science and Engineering, vol. 86-87, pp. 62–70, 2012.

Y. B. Li, W. F. Pu, B. Wei, Y. F. Chen, and B. J. Bai, “The feasibility of CO₂ and N₂ injection for the Tahe fracture-cavity carbonate extra-heavy oil reservoir: an experimental study,” Fuel, vol. 226, pp. 598–606, 2018.

J. S. Caine, J. P. Evans, and C. B. Forster, “Fault zone architecture and permeability structure,” Geology, vol. 24, pp. 1025–1028, 1996.

Q. Jin, J. J. Mao, Y. S. Du, and X. X. Huang, “Fracture filling mechanisms in the carbonate buried-hill of Futai oilfield in Bohai Bay basin, East China,” Petroleum Exploration and Development, vol. 42, no. 4, pp. 454–461, 2015.

S. J. Jolley, D. Barr, J. Walsh, and R. Knipe, “Structurally complex reservoirs: an introduction,” in Structurally Complex Reservoirs, S. J. Jolley, D. Barr, J. Walsh, and R. Knipe, Eds., pp. 1–24, Geological Society, London, 2007, Special Publications.

S. X. Wang, L. P. Guan, and H. L. Zhu, “Prediction of fracture-cavity system in carbonate reservoir: a case study in the Tahe oilfield,” Applied Geophysics, vol. 1, no. 1, pp. 56–62, 2004.

S. Becker, C. Hilgers, P. A. Kukla, and J. L. Urai, “Crack-seal microstructure evolution in bi-mineralic quartz-chlorite veins in shales and siltstones from the RWTH-1 well, Aachen, Germany,” Journal of Structural Geology, vol. 33, pp. 676–689, 2011.

S. E. Laubach, J. E. Olson, and J. F. W. Galea, “Are open fractures necessarily aligned with maximum horizontal stress?” Earth and Planetary Science Letters, vol. 222, pp. 191–195, 2004.

S. E. Laubach, R. M. Reed, J. E. Olson, R. H. Lander, and L. M. Bonnell, “Coevolution of crack-seal texture and fracture porosity in sedimentary rocks: cathodoluminescence observations of regional fractures,” Journal of Structural Geology, vol. 26, pp. 967–982, 2004.

M. Sirat, I. S. Al-Aasm, S. Morad et al., “Saddle dolomite and calcite cements as records of fluid flow during basin evolution: Paleogene carbonates, United Arab Emirates,” Marine and Petroleum Geology, vol. 74, pp. 71–91, 2016.

K. Zhao, Z. Gao, W. He, M. Wang, and Y. Hu, “Research on detection and identification of cracks in cemented backfill with ground penetrating radar,” Metal Mine, vol. 454, pp. 17–20, 2014.

R. B. Hanna and H. Rajaram, “Influence of aperture variability on dissolutional growth of fissures in karst formations,” Water Resources Research, vol. 34, no. 11, pp. 2843–2853, 1998.

M. Holland and J. L. Urai, “Evolution of anastomosing crack–seal vein networks in limestones: insight from an exhumed high-pressure cell, Jabal shams, Oman Mountains,” Journal of Structural Geology, vol. 32, no. 9, pp. 1279–1290, 2010.

K. Koegelenberg, A. F. M. Kisters, and C. Harris, “Structural controls of fluid flow and gold mineralization in the easternmost parts of the Karagwe-Ankole Belt of North-Western Tanzania,” Ore Geology Reviews, vol. 77, pp. 332–349, 2016.

C. M. Sayers, “Stress-induced fluid flow anisotropy in fractured rock,” Transport in Porous Media, vol. 5, pp. 287–297, 1990.

J. G. Ramsay, “The crack–seal mechanism of rock deformation,” Nature, vol. 284, pp. 135–139, 1980.

C. Bai, R. J. Spontaka, C. C. Kocha, C. K. Saw, and C. M. Balik, “Structural changes in poly (ethylene terephthalate) induced by mechanical milling,” Polymer, vol. 41, pp. 7147–7157, 2000.

T. Hirono, M. Takahashi, and S. Nakashim, “In situ visualization of fluid flow image within deformed rock by X-ray CT,” Engineering Geology, vol. 70, no. 1, pp. 37–46, 2003.

H. S. Mii, G. R. Shi, C. J. Cheng, and Y. Y. Chen, “Permian Gondwanaland paleoenvironment inferred from carbon and oxygen isotope records of brachiopod fossils from Sydney Basin, Southeast Australia,” Chemical Geology, vol. 291, pp. 87–103, 2012.

W. T. Parry and N. J. F. Blamey, “Fault fluid composition from fluid inclusion measurements, Laramide age Uinta thrust fault, Utah,” Chemical Geology, vol. 278, no. 1-2, pp. 105–119, 2010.

C. A. S. Teixeira, A. O. Sawakuchi, R. M. S. Bello, S. F. Nomura, D. J. Bertassoli, and M. A. C. Chamani, “Fluid inclusions in calcite filled opening fractures of the Serra Alta formation reveal paleotemperatures and composition of diagenetic fluids percolating Permian shales of the Paraná Basin,” Journal of South American Earth Sciences, vol. 84, pp. 242–254, 2018.

Z. M. Zhang, K. Shen, W. D. Sun et al., “Fluids in deeply subducted continental crust: petrology, mineral chemistry and fluid inclusion of UHP metamorphic veins from the Sulu orogen, eastern China,” Geochimica et Cosmochimica Acta, vol. 72, pp. 3200–3228, 2008.

P. Li, J. Y. Ren, H. Z. Yang, D. S. Hu, and Y. P. Zhang, “Hydrothermal activities in the fault system and their effects on carbonate rocks in the Bachu area,” Geotectonica et Metallogenia, vol. 35, no. 3, pp. 378–385, 2011.

Q. Liu, Y. R. Lu, F. Zhang, and K. N. Xiong, “Study of simulation experiment for carbonate rocks dissolution under hydrodynamic pressure,” Rock and Soil Mechanics, vol. 31, no. 1, pp. 91–102, 2010.

Q. Q. Meng, D. Y. Zhu, W. X. Hu, and Z. J. Jin, “Dissolution-filling mechanism of atmospheric precipitation controlled by both thermodynamics and kinetics,” Science China Earth Sciences, vol. 56, no. 12, pp. 2150–2159, 2013.

D. M. Shao, “Influence of temperature on dissolution rate in Ordovician carbonate rock in different water flow rate,” Coal Geology & Exploration, vol. 40, no. 3, pp. 62–66, 2012.

H. Yasuhara, A. Polak, A. Mitani, A. S. Grader, P. M. Halleck, and D. Elsworth, “Evolution of fracture permeability through fluid–rock reaction under hydrothermal conditions,” Earth and Planetary Science Letters, vol. 244, no. 1-2, pp. 186–200, 2006.

R. W. Lewis and W. K. S. Pao, “Numerical simulation of three-phase flow in deforming fractured reservoirs,” Oil & Gas Science and Technology, vol. 57, no. 5, pp. 499–514, 2002.

J. Y. Li, Y. S. Wang, C. H. Liu, D. W. Dong, and Z. Q. Gao, “Hydrothermal fluid activity and the quantitative evaluation of its impact on carbonate reservoirs: A case study of the Lower Paleozoic in the west of Dongying sag, Bohai Bay Basin,” Petroleum Exploration and Development, vol. 43, no. 3, pp. 359–366, 2016.

W. Shao, T. A. Bogaard, M. Bakker, and R. Greco, “Quantification of the influence of preferential flow on slope stability using a numerical modelling approach,” Hydrology and Earth System Sciences, vol. 19, pp. 2197–2212, 2015.

D. Y. Zhu, Q. Q. Meng, W. X. Hu, and Z. J. Jin, “Differences between fluid activities in the central and North Tarim Basin,” Geochimica, vol. 42, no. 1, pp. 82–94, 2013.

K. J. Heffer and J. Lean, “Earth stress orientation-a control on, and a guide to, flooding directionality in a majority of reservoirs,” in Reservoir Characterization III, W. Linville, Ed., pp. 799–822, PennWell Books, Tulsa, 1993.

J. E. Gaiser, “Stress direction hints and flow,” AAPG Exploration, vol. 24, pp. 36–39, 2003.

P. F. Dobson, S. Salah, S. Nicolas, and E. L. Sonnenthal, “Simulation of water-rock interaction in the Yellowstone geothermal system using TOUGHREACT,” Geothermics, vol. 33, no. 4, pp. 493–502, 2004.

E. Sonnenthal, A. Ito, N. Spycher et al., “Approaches to modeling coupled thermal, hydrological, and chemical processes in the drift scale heater test at Yucca Mountain,” International Journal of Rock Mechanics and Mining Sciences, vol. 42, no. 5-6, pp. 6719–6987, 2005.

T. Xu and K. Pruess, “Modeling multiphase non-isothermal fluid flow and reactive geochemical transport in variably saturated fractured rocks: 1. Methodology,” American Journal of Science, vol. 301, pp. 16–33, 2001.

T. Xu, E. Sonnenthal, N. Spycher, and K. Pruess, TOUGHREACT user’s Guide: A Simulation Program for Non-isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media, V1.2.1, Ernest Orlando Lawrence Berkeley National Laboratory, Berkeley, CA, USA, 2008.

L. Zheng, J. A. Apps, Y. Zhang, T. Xu, and J. Birkholzer, “On mobilization of lead and arsenic in ground water in response to CO₂ leakage from deep geological storage,” Chemical Geology, vol. 268, pp. 281–297, 2009.

T. Xu, N. Spycher, E. L. Sonnenthal, G. Zhang, L. Zheng, and K. Pruess, “TOUGHREACT version 2.0: a simulator for subsurface reactive transport under non-isothermal multiphase flow conditions,” Computers and Geosciences, vol. 37, no. 6, pp. 763–774, 2011.

P. Karsten and S. Nicolas, “ECO2N-A fluid property module for the TOUGH2 code for studies of CO₂ storage in saline aquifers,” Energy Conversion and Management, vol. 48, pp. 1761–1767, 2007.

T. N. Narasimhan and P. A. Witherspoon, “An integrated finite difference method for analyzing fluid flow in porous media,” Water Resource Research, vol. 12, no. 1, pp. 57–64, 1976.

J. E. Warren and P. J. Root, “The behavior of naturally fractured reservoirs,” Society of Petroleum Engineers Journal, vol. 3, pp. 245–255, 1963.

Y. S. Wu, Y. Di, Z. J. Kang, and P. Fakcharoenphol, “A multiple-continuum model for simulating single-phase and multiphase flow in naturally fractured vuggy reservoirs,” Journal of Petroleum Science and Engineering, vol. 78, pp. 13–22, 2011.

J. W. Cao, R. Y. Xia, Q. Y. Zhang, B. Liang, and Y. Dan, “Comparative study about the characteristics of inclusions in the paleokarst fillings in the northern outcrop and underground areas of Tarim Basin,” Carsologica Sinica, vol. 35, no. 1, pp. 114–124, 2016.

X. B. Lu, Y. Wang, F. Tian et al., “New insights into the carbonate karstic fault system and reservoir formation in the southern Tahe area of the Tarim Basin,” Marine and Petroleum Geology, vol. 86, pp. 587–605, 2017.

D. H. Liu, X. M. Xiao, J. K. Mi et al., “Determination of trapping pressure and temperature of petroleum inclusions using PVT simulation software—a case study of lower Ordovician carbonates from the Lunnan low uplift, Tarim Basin,” Marine and Petroleum Geology, vol. 20, no. 1, pp. 29–43, 2003.

J. W. Zhang, Y. H. Zou, P. P. Li, X. Y. Fu, and W. Wang, “A new PVT simulation method for hydrocarbon-containing inclusions and its application to reconstructing paleo-pressure of gas reservoirs,” Petroleum Geology and Experiment, vol. 37, no. 1, pp. 102–107, 2015.

H. L. Li, N. S. Qiu, Z. J. Jin, and Y. L. He, “Geothermal history of Tarim basin,” Oil&Gas Geology, vol. 26, no. 5, pp. 613–617, 2005.

E. M. Edward and R. W. Healy, “Numerical investigation of steady liquid water flow in a variably saturated fracture network,” Water Resources Research, vol. 29, no. 12, pp. 4091–4102, 1993.

H. Cui, P. Guan, and X. Jian, “Magmatic hydrothermal fluids-formation water compound system and diagenetic response of carbonate reservoir rocks in northern Tarim Basin,” Acta Scientiarum Naturalium Universitatis Pekinensis, vol. 48, no. 3, pp. 433–443, 2012.

C. S. Jia, X. J. Feng, D. Rao, and R. X. Gao, “Isotopic characteristics of oil field waters from Ordovician oil accumulations in Tahe oilfield and its geological significances,” Petroleum Geology & Experiment, vol. 29, no. 3, pp. 292–297, 2007.

Q. X. Yu, Z. Shi, Q. Li, and H. Guo, “Analysis of the genesis and accumulation conditions of the Ordovician in the Yubei 1 well area of the Tarim Basin,” Xinjiang Oil & Gas, vol. 1, pp. 1–5, 2016.

H. Kweon and M. Deo, “The impact of reactive surface area on brine-rock-carbon dioxide reactions in CO₂ sequestration,” Fuel, vol. 188, pp. 39–49, 2017.

G. P. D. Silva, P. G. Ranjith, and M. S. A. Perera, “Geochemical aspects of CO₂ sequestration in deep saline aquifers: a review,” Fuel, vol. 155, pp. 128–143, 2015.

A. Lasaga, “Chemical kinetics of water-rock interactions,” Journal of Geophysical Research: Solid Earth, vol. 89, no. B6, pp. 4009–4025, 1984.

A. C. Lasaga, J. M. Soler, J. Ganor, T. E. Burch, and K. L. Nagy, “Chemical weathering rate laws and global geochemical cycles,” Geochimica et Cosmochimica Acta, vol. 58, pp. 2361–2386, 1994.

T. F. Xu, P. Rose, S. Fayer, and K. Pruess, “On modeling of chemical stimulation of an enhanced geothermal system using a high pH solution with chelating agent,” Geofluids, vol. 9, no. 2, 177 pages, 2009.

T. F. Xu, E. Sonnenthal, N. Spycher, and K. Pruess, “TOUGHREACT: a new code of the TOUGH family of non-isothermal multiphase reactive geochemical transport in variably saturated geologic media,” in Proceedings of the TOUGH Symposium, pp. 12–14, Lawrence Berkeley National Laboratory, Berkeley, CA, USA, 2003.

J. Na, B. Feng, C. Y. Lan, T. F. Xu, and X. Bao, “Effectiveness of using supercritical CO2 as stimulation agent for enhanced geothermal systems,” Journal of Central South University (Science and Technology), vol. 45, no. 7, pp. 2447–2458, 2014.

D. Buhmann and W. Dreybrodt, “Calcite dissolution kinetics in the system H₂O-CO₂-CaCO₃ with participation of foreign ions,” Chemical Geology, vol. 64, no. 1-2, pp. 89–102, 1987.

P. M. Dove, “The dissolution kinetics of quartz in sodium chloride solutions at 25°C to 300°C,” American Journal of Science, vol. 294, pp. 66–71, 1994.

H. Guo, Y. B. Wang, X. M. Ni, C. Q. Tao, and W. L. Han, “Chemical kinetics of ca/mg elements from calcite / dolomite in carbonic acid solution,” Journal of China Coal Society, vol. 41, no. 7, pp. 1806–1812, 2016.

J. W. Morse and R. S. Arvidson, “The dissolution kinetics of major sedimentary carbonate minerals,” Earth Science Reviews, vol. 58, no. 1-2, pp. 51–84, 2002.

J. R. O'Neil and S. Epstein, “Oxygen isotope fractionation in the system dolomite-calcite-carbon dioxide,” Science, vol. 152, no. 3719, pp. 198–201, 1966.

C. B. Zhao, B. E. Hobbs, A. Ord, M. Kühn, H. B. Mühlhaus, and S. L. Peng, “Numerical simulation of double-diffusion driven convective flow and rock alteration in three-dimensional fluid-saturated geological fault zones,” Computer Methods Applied Mechanics Engineering, vol. 195, no. 19–22, pp. 2816–2840, 2006.

A. Schepers and H. Milsch, “Dissolution–precipitation reactions in hydrothermal experiments with quartz–feldspar aggregates,” Contributions to Mineralogy and Petrology, vol. 165, no. 1, pp. 83–101, 2013.

M. T. Van Genuchten, “A closed-form equation for predicting the hydraulic conductivity of unsaturated soils,” Journal of American Soil Science Society, vol. 44, no. 5, pp. 892–898, 1980.

A. T. Corey, “The interrelation between gas and oil relative permeabilities,” Producers Monthly, vol. 19, no. 1, pp. 38–41, 1954.