当前位置：
X-MOL 学术
›
Lithosphere
›
论文详情
Our official English website, www.x-mol.net, welcomes your feedback! (Note: you will need to create a separate account there.)
Development Characteristics and Quantitative Prediction of Multiperiod Fractures in Superdeep Thrust-Fold Belt
Lithosphere ( IF 2.4 ) Pub Date : 2021-01-08 , DOI: 10.2113/2021/8895823 Jianwei Feng 1 , Junxiao Qu 1 , Pengfei Zhang 2 , Feng Qin 2
Lithosphere ( IF 2.4 ) Pub Date : 2021-01-08 , DOI: 10.2113/2021/8895823 Jianwei Feng 1 , Junxiao Qu 1 , Pengfei Zhang 2 , Feng Qin 2
Affiliation
For a superdeep imbricated thrust-fold belt in the Kuqa depression of Tarim Basin, NW China, the structural fractures have a great impact on the tight gas reservoir productivity. In this research, structural fractures were characterized from core data and imaging logging data. Numerically, the finite element (FE) method was applied to simulate the 3-D paleotectonic stress field of the key fracture-generating period as well as the present-day stress field in the Keshen gas field. Based on previously developed geomechanical models, we further derived the models of fracture aperture and porosity under the palaeostress field. A fracture permeability model considering fracture filling degree and stress was developed based on the Fracture Seepage Theory. Finally, we obtained a series of calculation models for the present-day fracture aperture and permeability. The predictions based on these models agreed well with actual measurement results, with most of the relative errors less than 20%. The developed 3-D FE geomechanical model and fracture parameter prediction method hold great promise for characterizing fractures in other deep low-permeability reservoirs.Fracture is the most widely developed fault structure in the upper crust, and its key difference from the fault is the indistinct displacement [1, 2]. All fractures/joints in a region are generally the product of long-term multiple tectonic activities and various stress states. Therefore, the study of fracture development, geometry, and distribution is important for exploring the history of deep tectonic deformation and restoring the evolution of the paleotectonic stress field in this area [3, 4]. In detail, the apertures, densities, and orientations of healed fractures/joints from borehole measurements are applied to investigate the polyphase tectonics and indicate paleostress intensity and orientations [5–7]. At the same time, the existence of fractures provides structural conditions for the ascending, dispersing, and infiltration of mineral fluids and is also a transportation channel and storage place for oil, natural gas, and groundwater. In particular, for deep/superdeep tight sandstones, structural fractures are widely developed due to their distinct characteristics from other conventional reservoirs, such as large burial depth, large pressure coefficient, high compaction, high reservoir heterogeneity, super-low permeability, and developed structural fractures [8–12]. This extensively developed structural fracture system is not only an important seepage channel and space of the tight sandstone reservoir but also seriously affects the productivity of a single well [13–17]. Therefore, understanding where the subsurface structural fractures develop and when they were formed, including their development degree, physical properties, and seepage characteristics, is very important for the effective exploration and development of fractured tight reservoirs. Accurate prediction on the spatial distribution and effectiveness of fracture networks has always been a universal challenge. Spatial quantitative characterization and modeling of complex fracture networks in deep or superdeep tight reservoirs within complex tectonic zones presents more difficulty [18–20].Many methods and techniques have been used to identify, characterize, and predict natural fractures in subsurface reservoirs, which can be generally divided into qualitative, semiquantitative, and quantitative categories. First, the approximate estimate of structural fracture parameters is usually based on geostatistical methods such as outcrop measurement, core observation, and microscopic thin section observation (Fossen et al., 1996; Golf-Racht,1996). Based on these basic means, fracture parameters such as density, length, and filling degree can be statistically analyzed and empirically determined. In particular, X-ray computed tomography (CT) scanning can quantitatively visualize the dynamic evolution process of microfractures from initiation, expansion, to destruction and is, therefore, a favorable technique for the characterization of multiscale fracture parameters ([21]; Zhu et al., 2011; [22] [23]). Second, the relationship among fracture development and distribution as well as principal curvature of fold has been mainly acquired from structural features. Also, geomechanical modeling has been used to reveal the relationship between fracture attributes and folding [24]. However, this curvature method has certain limitations, and the plastic deformation and creep of the brittle rock layer under long-term stress field are not fully considered.Third, based on a large number of geological and production data, discrete fracture network (DFN) modeling for low-permeability reservoirs has been established by using multivariate statistical method, deterministic method, or random interpolation method. This DFN model can describe the characteristics of the fracture system more realistically and finely from the collection form to complex seepage behavior. However, the accuracy of this model is generally influenced by the reliability of basic geological data and the analysis of the main controlling factors of fracture development. Fourth, the use of geophysical data, such as logging and seismic data, to characterize the density, opening, length, porosity, and permeability of fractures is still a common method. However, this method is still in the quantitative or semiquantitative stage due to the limited spatial resolution of original data, and the reliability and accuracy need to be further verified and improved. Finally, based on the theories of fracture mechanics and rock mechanics, the rock rupture criterion and dynamic numerical simulation method have been applied to establish a quantitative characterization model between strain energy, rupture rate, and volume density, which is an effective and most promising method ([25–27]; Liu et al., 2015; [22, 28–30]).In this paper, taking Keshen tight gas reservoir, a typical deep/superdeep imbricated thrust-fold belt as a case study, based on restoration of structural evolution history, paleostress tests, and rock mechanics experiments, finite element (FE) method was applied to simulate paleostress field in the key tectonic period. Based on comprehensive analysis of borehole breakouts, drilling-induced fractures, and logging data, the magnitude and direction of three principal stresses were restored and determined, and the distributions of in situ geostress field were simulated. The distributions of natural fracture parameters, such as density, aperture, porosity, and permeability, were predicted by the establishment of composite rock failure criterion and derivation of a superimposed fracture geomechanical model. At the same time, the actual fracture data measured from wells was used to verify the numerical simulation results. Finally, the combination of geomechanical modeling and finite element simulation effectively predicted the fracture distribution in tight reservoirs. This fracture prediction method is particularly useful for the in-depth exploration and fine development of deep gas reservoirs in imbricated thrust-fold belt.As an enrichment area for natural gas, the Kuqa depression is distributed along the northern edge of Tarim Basin, with the South Tianshan Mountains bounded in the north and the Tabei Uplift in the south. The total coverage of the Kuqa depression is 28,500 km2, spreading 550 km from east to west and 35–80 km from north to south [31, 32] (Figure 1). Structurally, the Kuqa depression is a typical reactivated foreland basin controlled by the South Tianshan Mountains. It has a large number of foreland basin sediments and structures, mainly formed under strong intraplate orogeny with large scale detachment since the Cenozoic [33–35]. A series of large-scale thrust faults and fault-related folds are widely developed in the Kuqa depression, mainly formed in the Cenozoic Himalayan movement. Laterally or from the Tianshan Mountains to the basin, the Kuqa depression is divided into five secondary tectonic units, including two sags and three structural belts, as shown in Figure 1(a) [36, 37]. As the second row of thrust-fold belts before the South Tianshan Mountains, the Kelasu belt is further divided into the Kela zone and the Keshen zone from north to south, based on the difference of fault structures. The Keshen zone/gas field is located in the footwall of the Kelasu belt, south of the Kelasu fault and north of Baicheng fault, and displays a series of imbricated thrust structures, with pop-up structure and fault-related fold in each imbricate fan. In recent years, many geologists believe that the main hydrocarbon-bearing formations in the folds are likely to be closely related to the location of the neutral surface of the folds, suggesting that folding can significantly influence/control the development and distribution of structural fractures [38–40]. Since degree of fold deformation can be reflected by the angle between two limbs of the fold (i.e., interlimb angle) [35], it can be found that there are only two types of folds in the Keshen area, which are gentle fold with interlimb angle within the range of 180°–120° and open fold with interlimb in range of 120°–70° (Figure 1(b)). The Keshen anticline is cut transversely by a series of gently dipping thrust faults striking nearly E-W, and dip of faults mainly ranges from 16° to 23°.The key gas-bearing interval in the Keshen gas field/anticline is the Cretaceous Bashkiqik Formation (K1bs), which is mainly a set of delta front sediments. The Bashkiqik Formation is mainly composed of medium-fine sandstones, siltstones, and mudstones interbedded with limited conglomerates. The buried depth of K1bs in the study area is generally in the range of 6,500 to 8,000 m, and the average total formation thickness can reach 290 m. According to the petrographic characteristics, from shallow to deep, the K1bs can be divided into lithological units, i.e., K1bs1, K1bs2, and K1bs3, which are interspersed with thin mudstone interlayers [11, 41]. The porosity of K1bs reservoir ranges from 2% to 7% through core tests and permeability lies in the range of (0.05–0.5) md, while total fracture permeability reaches (1.0–10.0) md. As a result, all the above evidence proves that K1bs reservoir in the Keshen gas field belongs to a superdeep tight reservoir with low porosity and low permeability.In this study, rock mechanics experiments were carried out with drilled cores from wells Ks1, Ks201, Ks202, Ks205, and Ks7 in the tight sandstone formations. All the samples were processed to a diameter of 25 mm and a height of 50 mm, of which 16 samples were performed to analyze the uniaxial and triaxial mechanical parameters of rocks. Respectively, the surface roughness was generally greater than class 5 (i.e., Ra<2.5 mm, Gb1031-83), and the errors of nonparallel and nonperpendicular were <5 mm and <0.11°.All the mechanical experiments were carried out in the Provincial Key Laboratory of Geotechnical Engineering in Sichuan University. Four sets of triaxial experiments were taken to acquire the compressive strength data, such as rocks’ Poisson’s ratio and Young’s modulus (Table 1). The experimental rules and procedures are mainly based on the “Hydropower Engineering Rock Test Code” (DLJ204-81) and “International Society of Rock Mechanics Recommended Test Methods” (ISRM-2014). The triaxial experiment was carried out on a TAW-100 microelectromechanical servo-controlled rock triaxial stress testing machine with a measurement accuracy of about 1.0%. Internal friction angles and cohesion values of rocks were determined by a TAW-100 triaxial testing machine with a deviation of approximately 1.0%. The confining pressure during experiments was set to 40 and 60 MPa, respectively (Table 1). Twelve uniaxial compression experiments were also performed on the above testing machine, which could automatically calculate the rocks’ uniaxial compression strength, Young’s modulus, and Poisson’s ratio under zero confining pressure (Table 2).Drilled cores are the most direct and intuitive source for evaluating subsurface fractures, and the filling composition of fractures (e.g., calcite and quartz) and the change in aperture [41] can be finely measured. However, due to the high coring cost and tedious processing, drilled cores are not always available. At present, logging methods widely used to identify wellbore fractures mainly include acoustic borehole televiewer logging (ABT), resistivity logging (RD, RS), density logging (DEN), compensated neutron logging (CAL), gamma ray logging (GR), and formation microimage logging (FMI) ([42]; Tao et al., 2010; [43–45]). Past experience has shown that some conventional logging curves are not sensitive to fracture characteristics of the Keshen gas field, but through continuous trial and analysis, some log curve combinations with relatively sensitive responses can still be selected. As a result, deep tight sandstone rocks in K1bs of Keshen gas field could be identified by RS, RD, GR, and FMI logging methods, and the relationship between gas production and fractures development was studied through gas experiments and log interpretation. From log curves, information on the change of GR values from low to high in the intervals of fractures could be obtained, that is, the RD and the RS did not coincide, and both had lower values. For GR logging, the enrichment of the uranium element indicates the existence of fractures, so the difference between uranium content and thorium uranium content can be used to identify fractures. Similarly, due to the different detection depths, the existence of natural fractures will lead to differences in the resistivity values of deep and shallow lateral logging. Generally, the difference in resistivity of steeply dipping fractures (>75°) is “positive anomaly,” and the higher the angle, the larger the difference, the larger the fracture aperture, and the larger the extension length, the better its effectiveness. In comparison, in general, the difference in resistivity of gently dipping fractures (<45°) is “negative anomaly.” The larger the difference, the larger the fracture aperture, the larger the extension length, and the better its effectiveness. However, the resistivity difference of oblique fracture is not obvious. Therefore, the high “positive anomaly” between RD and RS logs and the relatively low GR value of sandstone interval indicate the existence of large-scale, steeply dipping, low-filled fractures, most of which belong to the Set II and some conjugate fractures in Set IV. At the same time, the high “positive anomaly” and the relatively high GR value of the sandstone section indicate the existence of Set I and the Set V fractures with large scale, high angle, and high-filled degree. However, the high “negative anomaly” of between RD and RS logs and the relatively low GR value of sandstone interval indicate the existence of large aperture, gently dipping, and low-filled fractures. These fractures mainly belong to Set III and decollement type (as described later).Especially, formation microimage (FMI) logging compensates for this deficiency and enables continuous imaging of wellbore formations and fractures. And manual interpretation can be used to effectively measure fracture parameters, such as linear density, aperture, and porosity under in situ conditions without paleomagnetic orientation correction [46, 47]. In the K1bs formation of the Keshen gas field, a large number of FMI logs are carried out, which provided a critical basis for the statistical analysis and development characteristics of fracture parameters. The drilling core of the local well section is a good verification of the results of FMI interpretation, and the effective filling degree of fractures is obtained by manual identification. The main basis for FMI logging to identify fractures is that during drilling, the open effective fractures are invaded by drilling fluid, and the resistivity of the fractures is different from that of the surrounding rock. Fractures identifiable by FMI logs include drilling-induced fractures, high-conductivity fractures, and few high-resistance fractures, the latter two being natural fractures [46, 47]. In FMI logs, the induced fractures are generally arranged neatly, in the form of en-echelon combination, while the natural fractures are irregularly arranged (Figure 2). Morphologically, the natural fracture surface is irregular, and the opening degree (i.e., aperture) changes greatly, but the induced fracture surface is relatively regular, and the opening degree changes little [48–50].Taking well As201 as an example, fractures identified in the FMI logging have high dip angles from 60° to 90°, high linear density from 1 m-1 to 9.6 m-1, high fracture length from 0 to 8 m, and low fracture porosity from 0.015% to 0.15% (Figure 2). The gas content in tight sandstones was analyzed by logging data, and the profile of well As201 was compared. The high gas content of the interval with high fracture density indicated that fractures provided a large permeability or reservoir space for the low-permeability tight sandstone reservoir. In contrast, the lower gas content of the interval with lower fracture density indicated that the reservoir space was relatively poor. At the same time, statistical results showed that the characteristics of fracture parameters were closely related to increase of total gas content, that is, when the linear density of fractures exceeded 3 m-1, the length exceeded 2.7 m, the dip angle exceeded 65°, and the porosity was 0.055% or more.In Keshen area, core observation and FMI analysis show that the fracture dip angle is mainly in the range of 75° to 90° (i.e., vertical/subvertical fractures), followed by ranges 45°–75°and 0–45° (Figure 3(a)). In plane, the proportion of vertical/subvertical fractures is increasing from north to south, while the proportion of steeply dipping fractures is just the opposite. The proportion of gently dipping fractures is small, and they are only distributed in the east section of Block 1 and Block 2, indicating that the Keshen area is dominated by vertical and steeply dipping fractures. The statistical results also show that most fractures in K1bs were half-filled, with dolomite and calcite as the main filling minerals. On the whole, the fractures in northern blocks are filled to a higher degree than that in the middle and south, and the partially filled fractures account for 65% of the total [51]. The proportion of filled fractures is 70% in the east section of Block 1, 68% in the west section of Block 1, and 63% in Block 2, but only 27% and 26% in Block 3 and Block 4, respectively (Figure 3(b)). Since fractures are exposed to the surface, the aperture and porosity parameters will change slightly from the in situ state. From FMI interpretation and core observation, statistical results show that the fracture aperture is mainly distributed in two intervals, dominantly 0.2–0.6 mm and sparsely 0.6–1 mm (Figure 3(c)).In general, linear density (fractures/m or m-1) is a simple, effective, and intuitive parameter for characterizing the intensity of fracture development [52, 53]. The results from cores and FMI show that the average linear density of fractures within three lithologic members ranges mainly from 0.7 m-1 to 7.3 m-1 and has obvious spatial heterogeneity. In general, the fracture density has the following distribution trend, with the northern density higher than the south and the eastern density higher than the west. Structurally, the highest values of fracture density are mainly concentrated in the cores of anticlines and the fault zones, especially near the major faults. In contrast, the fracture density appears to be relatively low in the limbs and saddles of anticlines (Figure 3(d)). Based on core observations, the structural fractures in the Keshen gas field can be divided into two main types, namely, shear fractures and tensile fractures, respectively, accounting for 72% and 28% of the total fracture number (Figure 3(e)). Generally, shear fractures, or small fault with slip amount from tens of millimeters to tens of centimeters, can extend up to several tens of meters on the cores and image log (Table 3). Due to smooth straight surfaces and consistent orientations (e.g., cores of wells Ks2, Ks206, and Ks802), these shear fractures, in some positions, can form conjugate patterns under strong compression condition. Typically, a single tension fracture, or opening-mode fracture, extending only a few centimeters on the cores with irregular and rough surfaces, is usually developed in wells at the top of anticline. The tensile fracture has a dendritic structure, a relatively shorter distance, and frequently bypasses rock grains, wherever they are observable by drill core and FMI (e.g., cores of wells Ks201, Ks8004, and Ks9).In general, structural fractures in K1bs can be divided into five distinct fracture sets that are oriented NW-SE and NNE-SSW (set I), NW-SE and nearly S-N (set II), NE-SW and NWW-SEE (set III), and nearly E-W (set IV) and nearly S-N (set V), in which the first two are the main developmental types (Figure 4). Each fracture set is characterized by different specific dip angles, mechanical properties, and fracture density. To summarize, (1) core observations reveal that some first and second sets of fractures mechanically form planar and sectional conjugate fracture systems, while the third set of fractures also form planar conjugated shear fracture system, which are developed in the later tectonic movement period along with anticline formation. (2) The vertical and steeply dipping fractures in the first two sets are formed under a regional compression environment, mainly including large-scale shear fractures, with few or medium mineral filling, large aperture (0.4 to 1 mm), high linear density, and favorable spatial connectivity. Commonly, these fractures are mainly developed in the early horizontal strata and the limbs of late anticlines. (3) The third set of conjugate tensile fractures is formed under a local extensional environment at the top of the anticline, mainly including subvertical fractures with a wide aperture range of 0.4–2 mm, low linear density, different filling degree, and small scale. Because the tensile stress is perpendicular to the hinge at the top of anticline, this type of conjugated fractures has the characteristics of both tensile and shear fractures and occasionally appears as en-echelon arrangement. (4) The nearly E-W striking tensile fractures (set IV) closely related to set III are a kind of longitudinal fractures, which are mainly developed in the top after strong fold uplift. These fractures are perpendicular to the stratum and are especially developed in brittle tight sandstones. They have the characteristics of large aperture (most >0.8 mm), short extension, vary-widely degree of mineral filling, and high linear density, but medium connectivity. (5) In the nearly E-W striking fractures, those gently dipping oblique shear fractures constitute a set of sectional conjugate fracture system, which are mainly formed in the local stress field after the folds are uplifted. Commonly, these fractures have a medium mineral filling degree and aperture (0.4 to 0.8 mm), low linear density, large scale and good connectivity, symmetrically distributed on the top, and hinge zone of anticlines. (6) Other nearly S-N striking expansion fractures are mainly formed in a strong compression environment and are passive tensile fractures generated under triaxial compressive stress state. These fractures are generally scattered, with low linear density, low filling degree, but large aperture (some >1.2 mm). Additionally, some of the variably oriented net-shaped fractures are scattered around the main faults, with neither a fixed direction nor a large aperture, and most of which have been filled to cause poor effectiveness (Table 3, Figures 3 and 4).Since the late Cretaceous, with southward thrust of the Tianshan Mountains and northward movement of the Indian subcontinent, the Kuqa depression experienced a complex evolutionary process [54–56]. The rock acoustic emission and structural trace analysis are commonly the effective technologies for the determination of paleotectonic stress values and formation periods of fractures [57–59] According to previous research results, four regional tectonic movements in Kuqa depression are recorded, i.e., the late Yanshanian movement and the early, middle, and late (i.e., Kuqa Period) Himalayan movements, which control the formation of a series of faults, folds, and fractures [60–62]. Among them, the late Yanshanian movement in Kuqa depression is relatively weakest with stress magnitude of 28–35 MPa during the late Cretaceous (test data comes from the lower Jurassic drilling cores, with a depth range of 3618–3624.2 m, and the Triassic formation in the outcrop area), followed by the early Himalayan movement with lower magnitude of 48–60 MPa at the end of Paleogene (test data comes from the Lower Cretaceous and Paleogene in outcrop area). According to Geomechanics Theory [1] and Andersonian faulting modes [63], some nearly E-W striking full-filled tensile fractures were inferred to be formed under late Yanshanian palaeostress field, while the NE-SW, NWW-SEE, and nearly E-W striking tensile fractures are mainly the product of local extension stress environment at the top of anticline due to a rapid uplift during middle and late Himalayan periods. In other words, the development of fractures in the Keshen area was mainly controlled by the regional tectonic stress field and boundary faults in the early stage and was mainly controlled by the fold uplift and local stress field in the later stage. As the compression force of the Himalayan movement in the middle period gradually increased, the Keshen anticline started to form, the middle principal stress at the top of the anticline was transferred to the nearly E-W direction, and the minimum principal stress was transferred to the vertical direction. In this way, a sectional conjugated shear system striking E-W with a large interval was formed. In contrast, influenced by the following three-stage Himalayan stress fields (~NNW350° compression), a large number of NW-SE, NNE-SSW, and nearly S-N striking shear fractures were generated, with widely developed cutting or limiting phenomena at different timing. Due to the difference of tectonic stress intensity and strength of fluid cementation (Figure 5), the earlier the fracture formation, the lower the density, the smaller the scale, and the higher the filling degree, mainly characterized by either half-filled or full-filled with calcite minerals [60–62]. Therefore, the widely developed net-shaped fractures under the late Himalayan stress field (77.9 MPa), characterized by steep dips and lower filling degree/high permeability, are extremely important to future gas reservoir exploration and development in the Keshen anticline (Figure 5).Stress field simulation is usually influenced and restricted by many factors such as buried depth, structure (e.g., folds and faults), strata thickness, rock mechanical properties, and boundary conditions [15, 56, 64–66]. In this simulation, based on log interpretation and rock mechanics experiments, the mechanical parameters of different elements in K1bs formation were calculated by dynamic-static corrections. According to the well-to-seismic integration, the thickness of three lithologic members was obtained by 3-D seismic event tracking and time-depth conversion. After establishing geomechanical models through C++ programming, the 3-D paleotectonic stress field and the present-day geostress field simulation were carried out in six blocks of the Keshen area. Then, the relationship models between stress-strain and fracture parameters were implanted into the FE platform to predict the spatial distribution of fractures.The ANSYS software based on the FE simulation platform was used to simulate paleo- and present-day stress fields. The construction of the initial 3-D geologic model was based on the restoration of paleostructure in the late Himalayan/Kuqa period and incorporated the mechanical stratigraphy of K1bs1, K1bs2, and K1bs3 and six isolated blocks (Figures 6(a) and 6(b)). The average thickness of members K1bs1, K1bs2, and K1bs3 within the geologic model was set to 50 m, 160 m, and 70 m, respectively. In order to eliminate the boundary effects in the simulation process, it is necessary to nest a larger rectangular parallelepiped outside the Keshen geologic model with a planar area of about 1.5 times (Figures 6(c)–6(e)). Considering that the Kuqa period after Pliocene deposition was a critical period for fracture formation in K1bs reservoir, the thickness of upper cover layers was set to about 5000 m.The setting of mechanical parameters of fault zone directly affects simulation results, as this is the foundation of modeling [22]. Before formation of faults or folds, the overall Young’s modulus was about 115–135% of sediment values in the model, but the faults became soft/weak zones after formation. According to the evolution history of structures in the Keshen gas field (Figure 4), the number of faults in the late Himalayan period is basically the same as that in present-day structure. In FE numerical simulation, the static Young’s modulus, Poisson’s ratio, rock cohesion, and internal friction angle of each geologic unit were summarized (Table 4). The geological model was segmented using a triangular grid so that a large number of elements and nodes could be subdivided. It is worth noting that the meshing of target layers and fault zones needs to be much thinner than the mesh of surrounding rocks. The geologic model of the late Himalayan period was finally divided into 48,256 nodes and 125,677 elements, and the elements and nodes of the present-day model were 28,456 and 125,677, respectively (Figures 6(c) and 6(d)).Stress direction, magnitude, and displacement constraints are the primary mechanical boundary conditions for FE numerical simulation. Based on the orientations of the active faults and fractures in the different tectonic movements, stress directions of the fracture formation period were determined (Figure 4). Based on previous research results, the parallel or nearly parallel faults striking NEE80° and conjugate shear fractures indicate a compression of approximately NNW350° at the late Himalayan period [37, 55]. However, the paleostress values (81–85 MPa) of the late Himalayan period obtained by acoustic emission are too small to cause significant macroscopic fracture of the rock and could not be directly used for FE numerical simulation ([61], Zeng et al., 2005, [37]). Therefore, based on the assumption that structural fracture density is basically constant under the ancient and modern tectonic stress fields [28], we adopt the concept of “equivalent paleotectonic stress” in this study. The so-called equivalent method is based on the evolutionary history of tectonics and stratigraphic deformation in the Keshen area, and the paleostress is roughly estimated based on the theory of elastic deformation of rocks. Definitely, the effect of formation pore pressure and hydrostatic pressure should also be taken into account in the equivalent stress calculation, but such large-scale calculation is difficult to be realized in FE geological model, and the range and saturation of paleofluid are difficult to be determined. Sun et al. [67] believe that the maximum principal stress can reach a high value during the formation of reverse fault in terms of Andersonian faulting and Mohr-Coulomb theory (e.g., [68]), and it has a linear proportional relationship with the vertical principal stress. And, at the depth of 3000 m in situ, the maximum principal stress is about 7 ~ 8 times of the vertical principal stress. According to the burial history of the Keshen area, it is concluded that the average paleoburial depth of K1bs in the late Himalayan movement/Kuqa period is about 4500–5500 m (Figure 8). According to the conclusion of Sun et al. [67], the maximum principal stress of the Keshen gas field during this period is calculated as 548–650 MPa; after removing the effective formation pressure (calculated according to the Biot coefficient of 0.38, about 45 MPa), the maximum effective principal stress is between 503–605 MPa.After many attempts, the reasonable boundary condition scheme was finally confirmed. In detail, the magnitudes of 608 MPa and 120 MPa were applied to the four boundaries in the NNW-SSE direction and NEE-SWW direction as maximum principal stress (σ1) and minimum principal stress (σ3) (Figure 6(c)). The southern edge was set to a fixed boundary so that the blocking effect of the southern Indian plate could be reflected appropriately. Simultaneously, the bottom boundary of the model was also fixed/pinned, and gravity load could be automatically applied by software so that the system was in equilibrium.The in situ stress state is generally characterized by stress tensor, which consists of three orthogonal stresses, maximum horizontal principal stress (σH), vertical principal stress (σv), and minimum horizontal principal stress (σh), each with a magnitude and orientation [69, 70]. The horizontal stress orientation in the study area was determined from drilling-induced fractures (DIF) and borehole breakouts (BO). As shown in Table 5 and Figure 3, the results suggested that the orientation of σH in wells of Keshen gas field was mainly concentrated in two ranges, namely, NW330° and NE 35°. It is obvious that the orientation of σH was nearly S-N in the central region, which was significantly different from the eastern and western regions. Moreover, in each block, the stress direction of σH is mainly in the nearly S-N direction in the middle of anticline and is mainly the NNW or NW direction in the western part of anticline and the NNE or NE direction in the eastern part. As can be seen from Figure 5, the magnitudes of σH at different depths of K1bs formation were significantly larger than that of σv, indicating a widely distributed Type-III stress environment (i.e., σH>σv>σh, compressional Andersonian stress state). In the Keshen gas field, the depth range was in the range of 5500~9900 m, and the drop reached 4400 m. In this case, the determination of horizontal principal stress must be considered in the vertical direction to provide a basis for the subsequent fine stress field simulation. Combined with multiwell logging interpretation, vertically, the average gradient of horizontal principal stress was 0.02548 MPa/m, indicating that the σH gradually decreased from north to south and gradually increased from top to bottom. It can be said that the vertical gradient of the horizontal principal stress is also sensitive to the depth in both the paleostress and present-day stress fields.Finally, a reasonable boundary-setting scheme was determined. A compressive stress of 215 MPa was applied to the east and west sides of the boundary model, as initial horizontal σh (Figure 6(d)). Simultaneously, a compressive stress of 270 MPa was applied to the south and north sides, as initial horizontal σH. And a stress gradient of -0.02548 MPa/m was applied from bottom to top of the model, that is, the horizontal principal stress decreased by 0.02548 MPa for every 1 m of depth reduction. The vertical stress was automatically generated by software, with a gravitational acceleration of 9.8 m/s2.According to the theory of rock mechanics, fracture mechanics and the law of conservation of energy, starting from the results of CT scanning of rock mechanics, the mechanical relationship model between stress-strain or strain energy and the fracture volume density is established, and the fracture linear density is derived and established as a link. The mechanical model of density and aperture is in good agreement with the core experiment results to over 95%.Otherwise, if σ1+3σ3≤0, and θ=0, there was Dlf=Dvf.In the first and third lithologic members (K1bs1 and K1bs3), the maximum principal stress was mainly ranged from 531 MPa to 575 MPa, indicating compression environment (Figures 9(c) and 9(f)). The intermediate or vertical principal stress was mainly in the range of 58–114 MPa (Figures 9(a) and 9(d)). The minimum principal stress was mainly ranged from 87.2 MPa to 12 MPa, indicating both compression and extension tectonic environment (Figures 9(a) and 9(d)). As shown in Figures 6(a) and 6(d), tensile stress (i.e., negative values) was only distributed in the upper part of northeastern blocks or along the long axis of anticlines. The distribution characteristics of the three principal stresses were very similar, and their distribution trends indicated that in the Keshen area, the fold was the primary factor controlling stress distribution, and the fault was the second main controlling factor. For example, tectonic movement caused a portion of accumulated energy in strata to be released rapidly, the faults appeared as low-stress zones, and a stress transition zone was formed around the main faults. In general, the paleostress values in the Keshen area gradually decreased from west to east and from south to north. Vertically, the paleostress values gradually increased from shallow to deep, with the highest value appearing in the superdeep strata in the southwest. The results also showed that the directions of minimum and intermediate principal stresses had obvious zoning, with boundary slanting across the entire Keshen area (Figure 9(g)). Specifically, the minimum principal stress was in a vertical direction in the northeast and in a horizontal direction in the southwest, but the characteristics of intermediate principal stress were just the opposite (Figure 9(h)). It could be shown that when the paleoburial depth increased to a certain value (about 6050–6200 m), the stress state would undergo a significant transition, which was basically consistent with logging interpretation and test results. Differently, the direction of maximum principal stress was nearly S-N, which was basically consistent with that of regional tectonic stress (Figure 9(i)). The application of the Coulomb failure criterion and Griffith propagation criterion to geologic materials as well as stress in the Earth to predict newly created fault/fracture orientations is understood by Anderson’s Theory of faulting [63]. It could be seen from Figure 9 that for the paleostress state, the maximum principal stress (σ1) in the study area is always nearly S-N. In the lower depth region, the distributions of principal stresses indicated an Andersonian II-type stress state, which tended to form thrust faults or gently dipping conjugate shear fractures. However, in the deeper areas, the distributions of principal stresses indicated an Andersonian II-type strike-slip stress state, that is, tending to form strike-slip faults or en-echelon shear fractures. The reason was that as the depth increases, the vertical stress increased beyond the intermediate principal stress.Coincidentally, the distributions of maximum, minimum, and intermediate principal stresses were also similar within the Keshen area. However, the orientation of minimum principal stress under present-day conditions was horizontal throughout the Keshen area, with values varying between 92 MPa and 132 MPa (Figures 10(a) and 10(b)). The values of maximum principal stress were mainly in the range of 160~196 MPa, with direction horizontal in most areas and vertical only in parts of the southwest (Figures 10(c) and 10(d)). In contrast, the distribution of intermediate principal stress was very similar to maximum principal stress, except that dominant orientation gradually changed to S-N in the southwestern region (Figures 10(e) and 10(f)). In general, the highest values of the three principal stresses were concentrated in deeper layers of the southwestern part of the study area, while the lowest values were located in the north, northeast, and fault zones and at the top of the anticline. As in the above, under current conditions, in the southwestern part of the study area, when depth exceeded a certain depth (about 7950–8100 m), the state of geostress rapidly changed from type III (i.e., σH>σv>σh, σH is the horizontal maximum principal stress, σh is the horizontal minimum principal stress, and σv is the vertical principal stress) to type I (i.e., σv>σH>σh) [83]. Furthermore, according to fine simulations of Keshen 2 anticline/Block 2, from the west to the east, the direction of σH gradually changed from NW to near S-N and then to NE, which was in good agreement with the test results of wells (Figure 10(g)). According to Anderson’s Theory of faulting, the distributions of present-day principal stresses indicated an Andersonian II-type strike-slip stress state, that is, tending to form strike-slip faults or en-echelon shear fractures.According to calculation results, the distribution of fracture parameters was similar to that of the stress field. There were also two major fracture development zones in the northeast and northwest, with a small transition zone in the middle. In the whole area, the highest value of fracture linear density reached 7.3 fractures/m, and the lowest value was less than 0.05 fractures/m. Moreover, the degree of fracture development gradually decreased from east to west, from north to south, and from deep to shallow (Figure 11(a)). The high-value areas of fracture density were mainly concentrated at the top of the anticline and hanging wall of fault but were moderate in limbs and fault zones. From the Keshen 2 anticline or the Block 2, the fracture density was the highest in the wells along anticline hinge, such as wells Ks101, Ks202, Ks2, and Ks201, exceeding 6.2 fractures/m (Figure 11(b)). In general, the distribution characteristics of paleoaperture of fractures were very similar to that of linear density, and the high values were mainly located at the anticline core, hanging wall of fault, and northeastern area (Figure 11(c)). As primary controlling factor for the quality and productivity of tight gas reservoirs, the high-value areas of fracture porosity were mainly located in (1) the north and northeast, (2) around hinge of anticlines, (3) around faults, and (4) the upper sandstones (Figures 11(d) and 11(e)). Moreover, the high-value areas of fracture porosity (0.075–0.145%) basically coincided with the areas with larger aperture (0.24–0.5 mm), such as on top of the fold and fault zones, mainly due to low degree of fracture filling (10–20%) during late folding period. It was also found that faults at structural high point controlled the distribution of fracture porosity more obviously than faults at low point.After superposition between the present-day and palaeostress fields in the Keshen area, the overall values of fracture aperture under present-day conditions decreased significantly, mainly in the range of 0.04–0.14 mm (Figures 12(a) and 12(c)). At the same time, the present-day fracture porosity under the current conditions was also significantly lower than the original, mainly in the range of 0.004–0.04%, but the overall distribution trend had not changed much (Figures 12(b) and 12(d)). As described above, the distribution of fracture aperture under current conditions was very similar to that of fracture porosity. The high-value areas of fracture aperture and porosity were mainly located near the cores of anticlines and fault zones. However, due to the difference in shape and amplitude of folds, the control effect on the development of fractures was not greatly different. Folding of open anticlines in the northern part of the study area (Blocks 1, 2, and 3) had significant influence and controlling on the effectiveness of fractures, but were not apparent in the gentle anticlines in Blocks 4, 5, and 6 (Figures 13(a) and 13(b)). In detail, the positive effect of folding on fracture apertures was about 0–370 m to the axis. Similarly, the maximum influence distance of faulting on fracture apertures in the Keshen area reached about 400 m, and the maximum distance on fracture porosity was about 380 m (Figures 13(c) and 13(d)).The permeability anisotropy has a noticeable influence on the deliverability and the physical parameters of whole reservoir and gas reservoir capacity [84, 85]. Under current conditions, obviously, as a whole, the fracture permeability had obvious regional characteristics. Generally, the permeability in the north was higher than that in the south, and the east was higher than the west. For example, for the first lithology unit (K1bs1), in the northern part of the study area (Blocks 1, 2, and 3), the average fracture permeability in the east-west direction was greater than 180 md, the average permeability in the north-south direction was greater than 90 md, and the average in the vertical permeability was only greater than 30 md (Figures 12(e)–12(g)). In contrast, in the southern blocks (Blocks 4, 5, and 6), the average east-west fracture permeability was lower than 8 md, the north-south average was lower than 10 md, and the vertical average was lower than 16 md. It could be seen that in deep tight sandstone reservoirs, fracture permeability decreased with increasing depth. On one hand, for the same unit of strata (e.g., K1bs1), the northeast blocks had a shallow burial depth, and the average fracture permeability was generally high, while the southwest blocks had a large depth, and the average permeability was generally low. On the other hand, for the three lithological units, the average permeability of K1bs1 was generally higher than that of K1bs2 (Figures 12(h)–12(j)), and the average value of K1bs2 was generally higher than the K1bs3 (Figures 12(k)–12(m)). Obviously, from northeast to southwest, the direction of maximum fracture permeability changed, that is, from east−west E−W>north−south S−N>vertical to vertical>north−south S−N>east−west E−W. The reason for this overall change in fracture permeability could be explained simply by the change in present-day stress state, as shown in Figure 14. In detail, in the shallow part of the Keshen area, the σ1 was in the S-N direction, the σ3 was in a vertical direction, and the σ2 was in the E-W direction. Accordingly, permeability in the E-W direction was larger than that in the S-N direction, and permeability in the S-N direction was larger than the vertical permeability. In deep strata, when buried depth exceeded 7950 m, the stress state changed, that is, the σv became σ2, and the σ1 and σ3 were transferred to the horizontal direction. More interestingly, if the shallow stress state was rotated by 90° along σ1 axis, it became a deep stress state, and fracture permeability values in each direction changed accordingly [68]. As a result, the highest fracture permeability then shifted to the vertical, the lowest permeability shifted to the E-W direction, and the S-N permeability was the intermediate (Figure 14).In detail, the high values of fracture permeability in the Keshen area were mainly located in high parts and near fault zones, that is, structure also played an important role in controlling the distributions of fracture permeability [86]. As shown in Figures 4 and 5, during the middle and late stages of the Himalayan movement, as the compression stress continued, a series of anticlines rose rapidly and was finalized. At this time, due to the occurrence of tensile stress in cores of the anticlines, a large number of nearly E-W striking tensile fractures (e.g., Set I and Set II) developed. Compared with a series of early-stage conjugate shear fractures (e.g., Set III and Set IV) formed under regional stress field and late-stage conjugate fractures in the anticline limbs, these tensile fractures are either in the form of en-echelon or conjugated, and because of the late formation time, and high linear density and large aperture (Figures 11 and 12), they still had high permeability values.According to Table 6, the overall numerical simulations in wells were in good agreement with core observations, FMI, and CT results. Moreover, according to data analysis, the average relative error of simulation results for fracture apertures was mostly less than 20%, and the individual exceeded 35%. Coincidentally, the simulation results of fracture porosity of most wells were 80% coincident with the actual measured values, and the relative error was less than 20%. These were sufficient to show that the quantitative prediction method based on superimposed geomechanical models was credible for the study of structural fractures in deep tight gas reservoirs. The reason for the large error was probably due to local stress concentrations caused by frequent lithological changes and fault turns, which are currently difficult to reflect in large-scale FE model. Also, the description/prediction for present-day spatial distribution of the filling degree in fractures was still a problem, mainly determined by interpolation method based on well points, which brought a lot of uncertainty to the prediction of fracture porosity and permeability.For a further verification of the simulation results, we provided analogous results from another well Ks2-1-8, which was located on the southern limb of the Keshen 2 anticline. Well Ks2-1-8 was a well-documented cored well with a variety of geological data such as continuous cores, FMI, and CT scanning, approximately 85 m long. As shown in Figure 15, the overall modeling results were consistent with in situ core observations, FMI, and full-sized CT scanning. The average fracture aperture and porosity in K1bs2 was as high as 0.085 mm and 0.028%, respectively, according to the core and FMI measurements, while modeling results yielded 0.072 mm and 0.022%. Both the actual measurement and simulation results showed that the fracture development of well Ks2-1-8 decreased gradually from shallow to deep, and the fracture dip angle also became smaller as depth increased. Located in the southern steep limb of the anticline and close to the boundary fault, a set of nearly S-N striking shear fractures were developed in the k1bs2 of well Ks2-1-8, and the fracture density gradually decreased from top to bottom but showed a higher Gas production capacity. According to field test data, this well was high-yield, especially in upper thick sandstones, and the total daily gas production exceeded 50×104 m3. The lower production interval (6820–6835 m) originally had a daily gas production of less than 10×104 m3. After reservoir acidification reconstruction, the oil pressure reached 88.2 MPa, and the daily gas production exceeded 83.2×104 m3. The numerical simulation and actual measurement results showed that not only the fracture linear density was high but also the fracture length was large and the connectivity was high. In contrast, gently dipping fractures (<60°) developed in the lower-capacity interval, and their length was small, especially the vertical connectivity was poor. It was certain that most of the conjugated net-shaped fractures would develop in thicker sandstone layers, mainly because the tectonic stress field in thick strata was more stable and closer to the regional stress field. Commonly, structural fractures are difficult to form in ductile mudstone sections with low Young’s modulus and high Poisson’s ratio, and mudstone intercalations will become barrier for gas migration [11, 71, 87, 88]. However, when the sandstone layer was thicker and the mudstone layer was thinner, the large-scale fractures generated in the sandstone would fully penetrate the mudstone, which indicated that lithology variations could change the local stress field and fracture development.As the most important gas-producing blocks in the Keshen gas field, Blocks 2 and 3 were typical long-axis anticlines, with low-permeability characteristics. The average porosity of the reservoir matrix in Blocks 2 and 3 was less than 0.1%, which had a limited contribution to the reservoir spaces. However, structural fractures could effectively increase the reservoir permeability by 2–3 orders of magnitude, as evidenced by well-testing results, such as testing permeability of 1–80 md (Figure 16). Therefore, fractures were not only the main seepage channel of the reservoir but also determined the high and stable yield of the gas well in Blocks 2 and 3. Through single-well numerical simulation, it was found that the average absolute open flow of gas well (AOF) in Blocks 2 and 3 was only 4–5×104 m3/d, depending on the seepage capacity of the reservoir matrix. However, the contribution rate of fractures to the gas well productivity reached 60–99%, and the more developed the fractures, the better the effectiveness, and the higher the contribution rate (Figure 16). Due to the low porosity and low permeability of superdeep tight sandstone reservoir-Block 2, fractures are the most important seepage channels for gas reservoirs, which control the productivity, high, and stable production of a gas well. According to the characteristics of gas reservoir test production, reservoir physical properties, and fracture distribution, Block 2 is divided into three basic dynamic types: high-stable production area, high-declining production area, and low-stable production area (Figure 17). As shown in Figure 17(d), the high and stable production areas of gas reservoirs are distributed along the hinge of anticline, but they are not continuous in space. However, these high production areas of gas reservoir are basically consistent with the high-value areas of fracture density, aperture, and permeability, especially in the eastern part of the block (Figures 17(a)–17(c)). For example, well As3 in the middle of block is located in the high-stable production area. Before reservoir fracturing improvement, the average absolute open flow (AOF) reached 236×104 m3/d. In comparison, the AOF of well As3-1 located in the transitional position was only 13×104 m3/d before reservoir fracturing, and it reached 175×104 m3/d after fracturing improvement. The reason can be inferred that the simulated fracture development and permeability of well As3 are higher than those of well As3-1, and the fractures of the former are mainly oriented NW-SE and nearly S-N (i.e., set II), while the latter are mainly oriented nearly E-W (i.e., set IV). According to the fracture stress sensitivity experiments [89, 90], the smaller the angle between the fracture surface and the direction of the maximum principal stress, the better the fracture permeability. Under normal circumstances, we all think that the fault zone is the area with high permeability, but the fracture zone around faults is not necessarily the area with high permeability. The wells As2-1-3 and As2-1-5 located in the east, although they are close to the major fault, their daily gas production (DGP) is lower than 10×104 m3/d, and the AOF before and after reservoir fracturing is lower than 10×104 m3/d, which is contrary to conventional understanding (Figure 17(d)). Based on fracture prediction and single-well measurement, it was realized that the fracture density at this location has reached 5.5 m-1 or more (Figure 17(a)), but the fracture aperture and permeability are only between 0.54−0.56×e−4 and 70-90 m, respectively, which is very low compared to the core of the anticline (Figures 17(b) and 17(c)). As mentioned above, some conjugate shear fractures developed in the limbs of anticline in the early Himalayan movement, their density is low, and the filling degree is also high. Although a series of derivative fractures were formed by the fault activity in the later period, they are all messy and small in scale, that is, their connectivity is not good. As mentioned above, all of this evidence was sufficient to demonstrate that the superimposed geomechanical models were reliable for predicting the 3-D fracture seepage parameters in the imbricated thrust-fold belt.A series of adapted superimposed geomechanical models coupling paleo- and present-day stress fields were developed to quantitatively predict fracture distributions in deep tight sandstone reservoirs of Tarim Basin. The models presented good applicability for fracture research in complex tectonic areas. The modeling results were verified by the actual measurements of many wells in the studied area. More in-depth scientific research was recommended here, including the establishment of composite failure criteria for interbedded strata and the exploration of new methods to characterize complex mineral filling laws in subsurface fracture systems. (1)By comparing tectonic movement and fracture distributions, the macroscopic developmental features of fractures were summarized based on core observation and FMI interpretation. The Keshen gas field was dominated by vertical and steeply dipping fractures, while gently dipping and horizontal fractures were less developed. In the coring section of K1bs, the average linear density of fractures was distributed mainly between 0.7 and 7.3 m-1. The formation of structural fractures in the Keshen area was mainly in the late Himalayan period (dominated by large-aperture, low-filled, and subvertical conjugated shear fractures). Most shear fractures striking in NNW-SSE and NE-SW were formed mainly under NNW350° compression during Himalayan periods, and most tensile fractures striking in NEE-SWW and NWW-SEE were formed under local extensional stress field induced by rapid fold uplift during the late Himalayan period(2)In the piedmont imbricated thrust-fold belt, for open folds in the north (e.g., Blocks 1, 2, and 3), curvature or uplifting of folds was the primary factor controlling development and distribution of seepage parameters, followed by fault activity. However, for gentle folds in the south (e.g., Blocks 4, 5, and 6), regional stress field and faulting played important roles in controlling the development and distribution of fractures, while fold uplifting had little effect on fractures. Vertically, the fracture density, aperture, and porosity decreased with increasing strata depth, which primarily coincided with increasing degree of compaction in this superdeep tight sandstone reservoir. In addition, the thickness of single-bed sandstone and frequency of mudstone interlayer also locally influenced the distribution of fractures(3)In the Keshen gas field, the areas of high fracture permeability were mainly located in the high structures and fault zones, and the permeability values decreased with depth. In most areas, the permeability was highest in the E-W direction, medium in the S-N direction, and lowest in the vertical direction. However, in the southwest due to transformation of stress state, the permeability was highest in vertical direction, medium in S-N direction, and lowest in E-W directionBy comparing tectonic movement and fracture distributions, the macroscopic developmental features of fractures were summarized based on core observation and FMI interpretation. The Keshen gas field was dominated by vertical and steeply dipping fractures, while gently dipping and horizontal fractures were less developed. In the coring section of K1bs, the average linear density of fractures was distributed mainly between 0.7 and 7.3 m-1. The formation of structural fractures in the Keshen area was mainly in the late Himalayan period (dominated by large-aperture, low-filled, and subvertical conjugated shear fractures). Most shear fractures striking in NNW-SSE and NE-SW were formed mainly under NNW350° compression during Himalayan periods, and most tensile fractures striking in NEE-SWW and NWW-SEE were formed under local extensional stress field induced by rapid fold uplift during the late Himalayan periodIn the piedmont imbricated thrust-fold belt, for open folds in the north (e.g., Blocks 1, 2, and 3), curvature or uplifting of folds was the primary factor controlling development and distribution of seepage parameters, followed by fault activity. However, for gentle folds in the south (e.g., Blocks 4, 5, and 6), regional stress field and faulting played important roles in controlling the development and distribution of fractures, while fold uplifting had little effect on fractures. Vertically, the fracture density, aperture, and porosity decreased with increasing strata depth, which primarily coincided with increasing degree of compaction in this superdeep tight sandstone reservoir. In addition, the thickness of single-bed sandstone and frequency of mudstone interlayer also locally influenced the distribution of fracturesIn the Keshen gas field, the areas of high fracture permeability were mainly located in the high structures and fault zones, and the permeability values decreased with depth. In most areas, the permeability was highest in the E-W direction, medium in the S-N direction, and lowest in the vertical direction. However, in the southwest due to transformation of stress state, the permeability was highest in vertical direction, medium in S-N direction, and lowest in E-W directionWhether in the field of structural geology or petroleum exploration and development, the identification and prediction of deep fractures are a very difficult task, which has attracted the interest and attention of many geologists. In this study, the established geomechanical models have achieved certain results in practical applications, and the accuracy is relatively reliable. However, these mechanical models and geological models have also undergone some assumptions and simplifications, so there may be errors in local locations. In particular, how to establish a complex mechanical model for multiperiod fractures and establish a geological model and mechanical parameter model closer to the real underground conditions is worthy of further discussion.For data stored in a repository, such as software and input files used to produce the results in this work are available at doi: 10.17632/76ssybpycr.1.The authors declare no conflict of interest regarding this publication.This research was first financially supported by the National Natural Science Foundation of China (42072234), then the National Oil and Gas Major Project (2016ZX05047-003, 2016ZX05014002-006, and 2017ZX05013006-003), and the Fundamental Research Funds for the Central Universities (17CX05010) and Shandong Province Key Research and Development Program (2018GGX109003).
更新日期：2021-01-08
京公网安备 11010802027423号