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 (CaCO3) 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 [14]. 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 [710]. 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 [1114], 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, 1619]. 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 [2126] 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 [2731] 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 [3235]. 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 H2O-CO2-NaCl-CaCO3 system, Zhu et al. [35] simulated the interactions between CO2-rich fluid and carbonate rocks. They found that the changes of CO2 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 [3942]. 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 CO2 geological storage, which can accurately describe the multiphase fluid migration process of water and CO2 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.

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.

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 (H2O). 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., CO2-H2O 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 CO2-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 (CO2) temperature was taken as 255°C (Figure 3). In TOUGHREACT simulator, the gas (CO2) 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 CO2 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.

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 (CaCO3) + 90 vol% CaMg(CO3)2. 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).

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 [5456]. According to the classical scheme of Surin, the groundwater of Ordovician stratum is mainly the calcium chloride-type (CaCl2) of oil field water, with average pH of 5.8. As a result, in this modelling, a simplified CaCl2-type water model including calcium ions (Ca2+) and chloride ion (Cl-) was set as the initial injected hydrothermal fluid. The native reservoir water has a Ca2+ 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 CO2 and CH4, 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 K25 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].

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 K0 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 ([6467]; Duan et al., 2008) showed that CaCO3 solubility decreases with fluid temperature increases, which means the increase of CaCO3 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)).

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 CO2-rich dissolved fluid (i.e., acid solution) entered the CaCl2-type formation water (i.e., weak acid fluid), the native chemical equilibrium would be broken, resulting in a linear decrease in Ca2+ concentration and a gradual increase in Mg2+ 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 Ca2+ and Mg2+ was saturated or over-saturated. At this time, the reactions of calcite and dolomite dissolution were expressed as.

As time went by, due to under-saturation of Mg2+, 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 Ca2+ and Mg2+ ion concentration and permeability changes perpendicularly across fractures after 100 years. The high-concentration zones of calcium ion (Ca2+) 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 (Mg2+) 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 × 107 seconds, the pH of the formation water slowly increased to about 6.2. Before 1 × 108 seconds, the pH value increased rapidly, reaching 8.0. However, after 1 × 108 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 HCO3- and Ca2+ 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.

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 CO2 and H2O. 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 CO2 (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 CO2 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 CO2 continued to accumulate in the aqueous phase, i.e., CO2 concentration in the aqueous phase continued to increase, to the top, HCO3- and Ca2+ 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.

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., HCO3- and Ca2+), 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)).

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