3D geomechanical modeling and numerical simulation of in-situ stress fields in shale reservoirs: A case study of the lower Cambrian Niutitang formation in the Cen'gong block, South China

Accepted Manuscript 3D geomechanical modeling and numerical simulation of in-situ stress fields in shale reservoirs: A case study of the lower Cambrian Niutitang formation in the Cen'gong block, South China

Jingshou Liu, Wenlong Ding, Haimeng Yang, Ruyue Wang, Shuai Yin, Ang Li, Fuquan Fu PII: DOI: Reference:

S0040-1951(17)30276-7 doi: 10.1016/j.tecto.2017.06.030 TECTO 127542

To appear in:

Tectonophysics

Received date: Revised date: Accepted date:

19 January 2017 20 June 2017 26 June 2017

Please cite this article as: Jingshou Liu, Wenlong Ding, Haimeng Yang, Ruyue Wang, Shuai Yin, Ang Li, Fuquan Fu , 3D geomechanical modeling and numerical simulation of in-situ stress fields in shale reservoirs: A case study of the lower Cambrian Niutitang formation in the Cen'gong block, South China, Tectonophysics (2017), doi: 10.1016/ j.tecto.2017.06.030

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT 3D geomechanical modeling and numerical simulation of in-situ stress fields in shale reservoirs: A case study of the lower Cambrian Niutitang Formation in the Cen’gong block, South China Jingshou Liua,b,c,d, Wenlong Dinga,b,c,d*, Haimeng Yange, Ruyue Wangf,g, Shuai Yinh, Ang Lia, Fuquan Fua a. School of Energy Resources, China University of Geosciences, Beijing 100083, China

PT

b. Key Laboratory for Marine Reservoir Evolution and Hydrocarbon Abundance Mechanism, Ministry of Education, China University of Geosciences, Beijing 100083, China

RI

c. Beijing Key Laboratory of Unconventional Natural Gas Geology Evaluation and Development Engineering,

SC

China University of Geosciences, Beijing 100083, China

d. Key Laboratory for Shale Gas Exploitation and Assessment, Ministry of Land and Resources, China University

NU

of Geosciences, Beijing 100083, China

e. Oil Recovery Plant No.3, Zhongyuan Oilfield Co. Ltd., SINOPEC, Puyang 066004, Henan, China

MA

f. State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development, Beijing 100083, China

D

g. Petroleum Exploration and Production Research Institute, SINOPEC, Beijing 100083, China

PT E

h. School of Earth Science and Engineering, Xi’an Shiyou University, Xi’an 710065, China

Abstract: An analysis of the in-situ state of stress in a shale reservoir was performed based on

CE

comprehensive information about the subsurface properties from wellbores established during the development of an oil and gas field. Industrial-level shale gas production has occurred in the

AC

Niutitang formation of the lower Cambrian Cen'gong block, South China. In this study, data obtained from hydraulic fracturing, drilling-induced fractures, borehole breakout, global positioning system (GPS), and well deviation statistics have been used to determine the orientation of the maximum horizontal principal stress. Additionally, hydraulic fracturing and multi-pole array acoustic logging (XMAC) were used to determine the vertical variations in the in-situ stress magnitude. Based on logging interpretation and mechanical experiments, the spatial distributions of mechanical parameters were obtained by seismic inversion, and a 3D heterogeneous geomechanical model was established using a finite element stress analysis approach to simulate the in-situ stress 1

ACCEPTED MANUSCRIPT fields. The effects of depth, faults, rock mechanics, and layer variations on the principal stresses, horizontal stress difference (Δσ), horizontal stress difference coefficient (Kh), and stress type coefficient (Sp) were determined. The results show that the direction of the maximum principal stress is ESE 120°. Additionally, the development zones of natural fractures appear to correlate with regions with high principal stress differences. At depths shallower than 375 m, the stress type is

PT

mainly a thrust faulting stress regime. At depths ranging from 375 to 950 m, the stress type is mainly a strike-slip faulting stress regime. When the depth is greater than 950 m, the stress type is

RI

mainly a normal faulting stress regime. Depth, fault orientation, and rock mechanics all affect the

SC

type of stress. The knowledge regarding the Cen'gong block is reliable and can improve borehole stability, casing set point determination, well deployment optimization, and fracturing area

NU

selection.

MA

Keywords: Geomechanical modeling; in-situ stress field; shale reservoir; numerical simulation; lower Cambrian; Cen'gong block

PT E

D

1. Introduction

Rock stresses are mainly formed by gravitational and tectonic forces, and these forces can be determined via drilling, grooving, and coring techniques (Zang and Stephansson, 2010). A

CE

geomechanical model is developed early in the development of an oil and gas field, and it can be

AC

applied to address a wide range of problems that are encountered during the life cycle of a hydrocarbon reservoir (Zoback, 2007; Tutuncu and Bui, 2015; Wang and Samuel, 2016a). These issues may arise: (1) during the exploration and assessment phase of reservoir development, such as evaluating the fault seal (or leakage) potential (Zoback, 2007;Tingay et al, 2009; Çiftçi, 2013), borehole stability (Mount et al., 1988; Haghi et al., 2013; Rajabi et al., 2017), and hydrocarbon migration and accumulation (Tan et al., 2001; Zoback, 2007; Khair et al., 2013); (2) during the development phase, such as determining optimal well trajectories (Zoback, 2007; Tao et al., 2016a), casing set points and mud weights, as well predicting anisotropic permeability in fractured 2

ACCEPTED MANUSCRIPT reservoirs; (3) throughout the production phase of the reservoir, which requires the selection of optimal completion methodologies such as hydraulic fracturing or repeated hydraulic fracturing (Zhao et al., 2005; Zoback, 2007; Wang et al., 2011), as well as the assessment of stress sensitivity in reservoirs and fractures (Smart et al., 2001; Wang et al., 2015); and (4) during the secondary and tertiary recovery phases of reservoir development, which include the optimization of processes such

PT

as water flooding and steam injection (Zoback, 2007). After nearly a hundred years of development, the stress state in a reservoir is mainly determined

RI

by multiple methods. Research techniques involving in-situ stress are mainly divided into five types:

SC

structural trajectory analyses (Mitsakaki et al., 2013; Zhao et al., 2015), core test methods (Lehtonen et al., 2012; Ito et al., 2013; Goodfellow and Young, 2014), computational methods

NU

based on logging data (Nie et al., 2013a; Tao et al., 2016b), mine stress measurement methods

MA

(Zhao et al., 2005; Wang, 2011) and in-situ stress simulation methods (including physical and numerical simulations) (Jiu et al., 2013; Wu et al., 2017; Liu et al., 2017). In the published literature

D

(Hashimoto and Matsu'Ura, 2006; Jiu et al., 2013; Zeng et al., 2013; Guo et al., 2016), most of the

PT E

stress field simulations are 2D or simple 3D geological models, and the geological units are regarded as homogeneous and assigned the same mechanical parameters. However, such assumptions are often inconsistent with the actual geological conditions (King et al., 2012; Liu et al.,

CE

2015; Rajabi et al., 2016b, 2017), leading to bias in stress field simulations.

AC

The Niutitang formation in the TX-1 well generally contains shale gas and has a thickness of 41 m. During the vertical well development phase, the shale gas production in this well was 3,000 m3 per day. In addition, large-angle deviated well and horizontal well technology will be implemented to improve shale gas production; however, the World Stress Map (WSM) database does not contain high-quality information on the maximum horizontal principal stress (σH) (Table 1). Consequently, a systematic understanding of in-situ stress remains lacking. In this paper, the data obtained from hydraulic fracturing, drilling-induced fractures, borehole breakout, GPS, and well deviation statistics were used to determine the orientation of the maximum horizontal principal stress in the 3

ACCEPTED MANUSCRIPT Cen'gong block and its circumjacent regions. Additionally, hydraulic fracturing and multi-pole array acoustic (XMAC) logging were used to determine the vertical variations in the in-situ stress magnitude. Under the constraints of logging interpretation and mechanical experiments, the spatial distribution of the mechanical parameters was obtained using seismic inversion. Based on the 3D heterogeneous geomechanical model, the spatial distribution of principal stresses, the horizontal

PT

stress difference (Δσ), the horizontal stress difference coefficient (Kh), and the stress type

RI

coefficient (Sp) were predicted.

SC

2. Geological setting

NU

2.1. Location and structure

The Cen'gong block is situated in northeastern Guizhou province, South China. The block is

MA

located in a trough-like fold belt in the southeastern Upper Yangtze plate (Fig. 1A) and has undergone multiple phases of tectonic activity, resulting in an extremely complex structure. The

D

main structural characteristics of the Cen'gong block are as follows: (1) complex structures, faults

PT E

and folds mainly trend NE and NNE; (2) most faults are reverse faults with high dip angles of 60–85°; and (3) the deformation of the strata and late fault activity were mainly dominated by

CE

Yanshanian activity in this area (Zeng et al., 2013; Zhao, 2013) (Fig. 1A and B). GPS data show that the orientation of displacement and strain in the study area is mainly related to the subduction of the

AC

Indian plate, followed by the subduction of the Pacific plate (Fig. 2).

2.2. Stratigraphy

As shown in Fig. 3A, the main strata that developed in the study area are Sinian, Cambrian, Lower Ordovician, and Quaternary (Nie et al., 2013b; Wang et al., 2016b). The Niutitang shale, with a thickness of 50 to 70 m, in the lower Cambrian developed stably and was mainly deposited on a deep shelf. The shale can be divided into three members (Fig. 3B): the Lower Niutitang (Є1n1), which is mainly composed of siliceous shale with a TOC content that ranges from 1.04% to 13.13% 4

ACCEPTED MANUSCRIPT (mean, 4.45%); the Middle Niutitang (Є1n2), which is mainly composed of black shale with a TOC content that ranges from 0.70% to 10.49% (mean, 4.71%) and includes the high-quality source rocks of the Cen'gong block; and the Upper Niutitang (Є1n3), which is mainly composed of gray shale with a TOC content that ranges from 0.36% to 2.43% (mean, 1.33%). Generally, high silica, feldspar, and carbonate contents in the Niutitang formation correspond to high Young's modulus

PT

values, low Poisson's ratios, and high brittleness (Wang et al., 2016b).

RI

2.3. Natural fractures

SC

In this study, core and thin section descriptions, outcrop observations, imaging, and logging analysis (Fig. 4) were the main methods used to investigate the fracture development characteristics

NU

of the shale in the Cen'gong block, and the results indicate that fracture abundance can change

MA

abruptly from bed to bed for the variety of the mechanical parameters (e.g., Young's modulus, shear modulus); sub-vertical fractures are common in the outcrop (Figs. 4B and E); and the dip angles of fractures give priority to vertical fractures (80.1%), followed by high angle fractures (11.2%) and a

PT E

D

small amount of low angle tectonic fractures (2.8%). The Lower and Middle Niutitang shales are relatively well-developed fracture zones compared with the Upper Niutitang shale; the dominant strikes of fractures in the Cen'gong block are NE-SW, NW-SE, and N-S (Figs. 4D and F). Core

CE

fractures mainly include vertical (sub-vertical), high dip angle structural, and horizontal bedding or

and C).

AC

interlayer fractures, and these fractures are mainly either filled or half-filled with calcite (Figs. 4A

3. Samples and testing A mechanical analysis was performed with cores from the TX-1 well, CY-1 well, and TM-1 well. Sixty samples were collected from different depths in the Niutitang formation, and experimental data were used to establish fine-scale mechanical models. All the samples were collected in cylinders with diameters of 25.00 mm, heights of 50.00 mm, and surface roughness 5

ACCEPTED MANUSCRIPT values of greater than class 5 (Ra<2.5 μm; GB1031-83). The non-perpendicular and non-parallel errors were <0.11° and <5 μm, respectively. All mechanical tests were performed at the Civil Engineering Laboratory of Beijing University of Science and Technology, China. Uniaxial compressional experiments were carried out to obtain the Young's modulus and Poisson's ratio values of the rock. The experimental steps were based on the national standard of part 7 of GB/T

PT

23561.7-2009 and part 8 of GB/T 23561.8-2009: Methods for determining the uniaxial compressive strength and deformation parameters of coal and rock. The loading rate of the stress was 0.5

RI

MPa/s-1.0 MPa/s, and the accuracy of the deformation measuring device was 0.001 mm. A uniaxial

SC

compressive test was performed using an YTD-200 electronic pressure testing machine with a test force measurement accuracy of 1.0%. Tensile strength values were determined through Brazilian

NU

splitting tests using a WEP-600 microcomputer screen display testing machine with a test force

MA

measurement range of 12-600 kN. The experimental steps were based on the national standard of the GB/T 23561.10-2010. The loading rate of the stress was 0.03 MPa/s-0.05 MPa/s, the accuracy

D

of the length measuring device was 0.02 mm, and the machine accuracy was 1% of the test force

PT E

value. Triaxial compression tests were performed to collect shear strength data for the different types of rocks (Table 2). Generally, shear strength parameters were obtained according to the Mohr-Coulomb theory using 3-4 cylindrical specimens under different conning pressures. Cohesion

CE

values and the internal friction angles of the rock (Table 2) were determined using a TYS-500

AC

triaxial testing machine with a precision of 1.0%. The confining pressures were 0, 5, and 10 MPa. The applicable standard was part 9 of JTG E41-2005: Methods for determining the triaxial strength and deformation parameters of shale. The loading rate of the stress was 0.5 MPa/s-1.0 MPa/s, the accuracy of the deformation measuring device was 0.001 mm, and the machine accuracy was 1% of the test force value.

4. Test results of in-situ stress tensors The state of in-situ stress is described by a stress tensor, which comprises three orthogonal 6

ACCEPTED MANUSCRIPT principal stresses, each with an orientation and magnitude (Engelder, 1993; Houska, 2014; Rajabi et al., 2016a). The in-situ stress tensor is typically described using just four components: the minimum horizontal stress magnitude, maximum horizontal stress magnitude, vertical stress magnitude, and maximum horizontal stress orientation (Bell, 1996).

PT

4.1. Vertical stress magnitude A total of 25 hydraulic fracturing tests in two boreholes were conducted successfully at depths

RI

ranging from -100 to -1300 m. The vertical stress magnitude in Fig. 5 was commonly assumed to be

SC

the pressure exerted by the overlying rocks and could be exactly calculated by integrating the density log from the surface to the depth of interest (Zoback, 2007; Naimi et al., 2015) (Fig. 4D and

NU

E). The horizontal principal stress values at different depths were derived from hydraulic fracturing tests, and the computed vertical stresses are shown in Fig. 5. In Fig. 5, the linear regression

MA

represents the variations in σH and σh with depth. Note that two wells were drilled into the faults at -725 m and -1300 m, respectively, and low differential stresses in fault zones (Figs. 5A and B) can

PT E

D

be observed at depths of -700 to 750 m and -1240 to -1340 m. XMAC logging is an effective and practical method for determining stress fields using borehole cross-dipole acoustic logging data (Wen et al., 2005; Zhao, 2008; Huang et al., 2012). The direction of the fast shear wave and the

CE

maximum principle is consistent, and the value of the horizontal principal stress can be determined

AC

by a stress-velocity equation (Patterson et al., 1999; Wang et al., 2008; Zheng et al., 2009). The horizontal principal stress values at different depths derived from XMAC at depths ranging from -1400 to -1880 m are shown in Fig. 5C and Fig. 6. These figures show that linear regression provides a good representation of variations in σH and σh with depth (Figs. 5A and B).

σh = -0.0098D + 10.149

(1)

σH = -0.0159D + 11.105

(2)

7

ACCEPTED MANUSCRIPT 4.2. Orientation of σH The σH orientation in the Cen'gong block was determined from borehole breakouts (BO), GPS data, deviated well information, hydraulic fracturing, and drilling-induced fractures (DIF). The quality ranking system of the WSM project (Heidbach et al., 2010; Rajabi et al., 2017) was used to assess the reliability of the results in each well (Tables 1 and 3). Breakouts are stress-induced

PT

enlargements of the wellbore cross-section (Bell and Gough, 1979; Reinecker et al., 2010). When a

RI

wellbore is drilled, the cores removed from the subsurface are no longer supporting the surrounding rock; as a result, the stresses become concentrated in the surrounding rock. Borehole breakout

SC

occurs when stresses around the borehole exceed the compressive strength or shear strength of the

NU

borehole wall. Hence, the long axes of borehole breakouts are oriented approximately perpendicular to the σH orientation (Tingay et al., 2010; Zeng et al., 2015; Rajabi et al., 2017). For vertical well

MA

drilling, the ideal trajectory is perfectly vertical, but under the influence of in-situ stresses, the bit often deviates from the vertical direction, and the direction of deviation is oriented approximately

D

parallel to the σH direction. DIFs are created when the stresses concentrated around a borehole

PT E

exceed the tensile strength of the wellbore wall (Aadnoy and Bell, 1998; Tingay et al., 2010). DIFs typically develop as narrow, sharply defined features that are subparallel or at a slight angle to the

CE

borehole axis in vertical wells. The stress concentration around a vertical borehole is at a minimum in the σH direction. Hence, DIFs develop approximately parallel to the σH orientation (Aadnoy and

AC

Bell, 1998; Rajabi et al., 2017). Hydraulic fractures can be documented by the impression packers in open hole wells. The impression orientation system is placed at a specified depth. Then, a high-pressure water pump is pressurized, and the applied pressure value generally meets or exceeds the rock fracturing pressure. When the high pressure inside the impression is removed, the impression retracts away from the wall, but its shape is retained in the fracture. Between the fracture, the pointer of the magnetic indexer, and the "base line" marked on the impression, the direction of the maximum horizontal principal stress can be determined. We catalogued the deviated well information ranging from -50 to -1875 m in the TX-1 well. 8

ACCEPTED MANUSCRIPT The results, based on drilling data, indicate that the deviated well orientations (quality A, B and C) in the Cen'gong block range from 104° to 124° (Table 3 and Fig. 7A). According to the well deviation data from the vertical well, the σH orientation is ESE 114°. A well diameter analysis was performed to determine the spatial pattern of the mean σH orientation around the borehole according to the six-arm caliper log data in the TM-1 well. Based on this logging data, the results indicate that

PT

the main borehole breakout orientation is NNE 21°. Therefore, according to the principle of determining the in-situ stress based on borehole breakout, the σH orientation is ESE 111° (Table 3

RI

and Fig. 7C). After a statistical analysis of the drilling-induced fracture trends in the HY-1 well, we

SC

found that strikes of DIFs were mainly NNE-SSW, and the dip angles were mainly vertical fractures. Thus, the mean σH azimuth was generally ESE 117° (Table 3 and Fig. 8). Based on the impresser

NU

orientation of hydraulic fractures in the ZK4220 well at different depths, the σH orientation at depths

generally ESE 123° (Table 3 and Fig. 9).

MA

ranging from 248.90 m to 653.00 m ranged from 93.5° to 145.5°, and the mean σH azimuth was

D

In 1997, China launched the crustal movement observation network (CMONWC), which is

PT E

based on GPS observation technology. GPS data provide information on the current strain rate and orientation, which may or may not be consistent with the orientation of the total strain or stress field at a point in the Earth. Using the numerical simulation of the present crustal stress of mainland

CE

China as a viscoelastic spherical shell, Fan et al. (2012) found that the direction of the crustal strain

AC

detected via GPS measurements is consistent with the present maximum horizontal principal stress direction in the interior of the plate. However, the current strain orientation derived from GPS data may not be consistent with the orientation of the total stress field at the boundary of the plate. The Cen'gong block is located inside the Yangtze plate; thus, it is reasonable to use GPS data to determine the direction of maximum principal stress. The GPS data indicate the mean σH orientation in North Guizhou is ESE 118°. Regional stress directions (Table 1 and Figs. 2 and 7D), obtained from FM, HF, BO and OC, are consistent with the principal stress directions determined using other methods (Table 3). The results regarding the local and regional stress directions of the in-situ stress 9

ACCEPTED MANUSCRIPT obtained using the methods described above indicated that the direction of the horizontal maximum principal stress in the Cen'gong block was primarily ESE 120° (±10°).

5. Theory and method The finite element method (FEM) and ANSYS computer software were employed to produce a

PT

three-dimensional simulation of the stress field. The basic concept behind the FEM is that a geological body can be discretized into finite continuous elements that are connected by nodes.

RI

Each element is associated with appropriate mechanical parameters, and these elements are

SC

displaced by the loads placed on the model. The material parameters then convert the strain to a stress value for the element. The simulation of the tectonic stress field involves the following steps:

NU

(1) determine the tectonic units (e.g., faults, anticlines, and synclines) and the geometric shapes and

MA

characteristics of strata (e.g., elevation, thickness, and pinch out); (2) determine the spatial distributions of the mechanical parameters (e.g., Poisson's ratio and Young's modulus) of the rock at

D

different locations; and (3) determine the loading regimes for geological bodies and the associated

5.1. Geomechanical model

PT E

boundary conditions.

CE

5.1.1 Geological model and mechanical parameters

AC

The simulation of the 3D tectonic stress field in the Cen'gong block considers spatial variations in the mechanical parameters and changes in the morphology and occurrence of 65 faults. NMR and XMAC logging were employed to determine the mechanical parameters (e.g., Young's modulus, Poisson's ratio, and shear modulus), brittleness index, and rock elemental weight fraction in the profile. The rock mechanical parameters can change abruptly from bed to bed for a variety of lithologies. The siliceous shales at the bottom of the Niutitang formation have high elastic modulus, low Poisson's ratio and high horizontal principal stress (Fig. 6). The relationship between dynamic and static mechanical parameters has been determined based on logging interpretation and 10

ACCEPTED MANUSCRIPT mechanical experiments (Table 2; Figs. 5F, 5G, 6 and 10). The P-wave and S-wave velocities can be obtained according to the time difference between acoustic waves, and the elastic modulus of dynamic Young's modulus (Ed) and dynamic Poisson's ratio (μd) can then be obtained directly according to the following equations:

3Vp 2 - 4Vs 2

(3)

Vp 2 - Vs 2

Vp 2 - 2Vs 2 2(Vp 2 - Vs 2 )

(4)

RI

d =

2

PT

Ed = 10 ρVs -3

SC

where ρ is the rock density, g/cm3; Vp is the P-wave velocity, m/s; and Vs is the S-wave velocity, m/s.

NU

The seismic velocity field data contain information on the wave velocities in the different target layers throughout the study area. Based on the P-wave velocity and S-wave velocity of the

MA

logging interpretation, the 3D distribution of the P-wave velocity at different depths in the Niutitang shale can be separated on the basis of the correction of the seismic layer velocity data (Figs. 11A

D

and B), the time-depth conversion and the sub-layer correlation. Using the mathematical models of

PT E

the dynamic and static mechanical parameters (Fig. 10) and Vp and Vs (Fig. 11C), the 3D mechanics parameters of the whole study area can be determined. The processing method used to define the

CE

mechanical parameters of the fault zones was fundamental to the outcome of the modeling because

AC

it directly affected the results of the numerical simulation. The orientation of σH is expected to swing perpendicular to mechanically stiff or hard structures and parallel to mechanically weak or soft structures (Bell, 1996; Faulkner et al., 2006; Heap et al., 2010; Rajabi et al., 2017). The Young's modulus of a fault zone is typically smaller than the elastic modulus of a normal stratum in North Guizhou (Jiu et al., 2013; Zeng et al., 2013; Wu et al., 2017). In general, this modulus is generally 65–85% of the elastic modulus of a corresponding normal stratum (Hashimoto et al., 2006; Guo et al., 2016; Liu et al., 2017). Different mechanical parameters (Young's modulus and Poisson's ratio) were associated with different faults based on the fault size, namely the larger of the fault length, the greater of the fault throw, and the greater Young's modulus of the rock, the smaller the Poisson's 11

ACCEPTED MANUSCRIPT ratio, and vice versa. Finally, the spatial distribution of the Young's modulus of the Niutitang shale and faults was obtained (Fig. 12). As shown in Fig. 12, the Young's modulus of the Niutitang shale generally decreased from northwest to southeast in the Cen'gong block, which may be related to the paleo-sedimentary environment where the water gradually became shallower from the southeast to the northwest

PT

during the Niutitang shale sedimentary period (Zeng et al., 2013; Liu et al., 2016; Wu et al., 2016). In addition, the Young's modulus of the Niutitang shale is high in the northeast, northwest and

RI

southwest regions of the study area, which may be related to the rock flexure. Based on the

SC

establishment of accurate geological geometrical models (including the faults, folds and buried depths and thickness of the Niutitang shale), rational meshing is the key to stress field simulation.

NU

5.1.2 FEM meshing

MA

As shown in Fig. 13A-1, because of the existence of pinch out and the significant lateral variation in sediment thickness of the Niutitang shale, when the model is meshed with the same

D

length, the quadrilateral elements are of poor quality (Fig. 13B-2) and the triangular elements are

PT E

evenly distributed and uniform in shape (Fig. 13B-1). To reduce the calculation complexity and improve the calculation accuracy, the geological model was meshed with the triangular elements and subdivided into a series of nodes and elements. According to the characteristics of the element

CE

type in ANSYS software, the Solid 45 element type was chosen for the meshing of the shale

AC

reservoirs, which is suitable for the element division of the layered structure. In general, finer elemental divisions produce more accurately computed results, but the computational complexity is also greater. Therefore, we divided the target stratum and fault zone into a number of triangular mesh elements with a side length of 3-10 m (Fig. 13A-1). The surrounding rocks were meshed with a coarser grid with a side length of 10-300 m. The study area around the model was divided into a fine grid (Fig. 13A-2), and the area away from the study area was divided into a rough grid (Fig. 13A-3). The model was subdivided into 1,091,415 elements and 196,452 nodes, and then the elements and nodes were associated with the different mechanical 12

ACCEPTED MANUSCRIPT attributes that corresponded to Figs. 10 and 12, namely from a homogeneous model (Fig. 13B) to a heterogeneous model (Fig. 13C).

5.2. Boundary conditions As outlined above, the good correlation between the regional stress directions and the observed

PT

stress information in the Cen'gong block (Fig. 2-L1, L2, L3 and L4; Table 1; Table 3) indicated that gravity and plate boundary stresses are an important source of stress in the study area. The regional

RI

pattern of the stress field over the Cen'gong block can be predicted by considering the plate-scale

SC

boundary conditions. The boundary conditions within the model include the boundary stresses and displacement conditions. The study area is primarily affected by the Indian plate, with

NU

approximately twice the force as that from the Pacific plate (Wang et al., 2003b). As shown in Fig. 2,

MA

GPS data also indicate that the movement rate of the land decreased from the northwest to the southeast in South China; therefore, it is reasonable to assume that there is a relatively stationary surface (zero displacement) when the ESE boundary is far enough away from the study area. As

PT E

D

shown in Fig. 13A, when l is large enough, the model's ESE 120° and NNE 30° boundaries (ideal infinity) are relatively stationary in the X direction and Y direction, respectively. Therefore, the Y direction constraint (zero displacement in the Y direction) was applied to the NNE boundary of the

CE

model, the X direction constraint (zero displacement in the X direction) was applied to the ESE

AC

boundary (Fig. 13), and the Z direction constraint was applied to the bottom boundary to prevent it from undergoing rotation and rigid displacement. An important model input is the boundary conditions. Compared with the force measured from the borehole, stresses are more easily obtained by hydraulic fracturing and logging data. After the finite element meshing, the boundary stresses at different depths (Figs. 5A, B and C) is more easily applied to the elements with different depths and sizes; therefore, the stress and displacement constraints are used as boundary conditions for the stress field simulation in this study. Minimum and maximum horizontal stress magnitudes measured in the ZK4220, CD-1, and TX-1 wells are 13

ACCEPTED MANUSCRIPT used as boundary conditions. We used point measurements to represent boundary conditions for two main reasons: firstly, local heterogeneities are "smoothed" out, leaving only regional stress values when using a linear model with depth; secondly, because the range of the stress field simulation is much larger than that of the study area, unlike the simple linear model applied to the boundary of the model (Fig. 13A), applying the actual boundary conditions to the study area (Fig. 13B) takes

PT

into account the local and regional stresses, leading to results that are closer to the actual stress conditions (Fig. 13C-EFGH). The variation of the minimum and maximum horizontal principal

RI

stresses with depth (see Figs. 5D and E for variations of the principal stress with depth) is applied to

SC

the SSW and WNW boundaries of the model, respectively (Fig. 13A). The distributions of the maximum, intermediate, and minimum principal stresses within the Niutitang formation were

NU

determined using a three-dimensional FEM numerical simulation.

MA

6. Results and discussion

D

6.1. In-situ stress field

PT E

In this paper, compressional stress was positive and tensile stress was negative. Fig. 14 shows that the distributions of the maximum, intermediate, and minimum principal stresses are similar,

CE

and that the gravity and plate boundary stresses play the major role within the Cen'gong block. The minimum principal stress (σ3) values are generally between 1-42 MPa, and the direction of σ3 is

AC

mainly NNE, except for in three denudation areas, in which σ1 is a vertical stress. The intermediate principal stress (σ2) values are generally between 3-55 MPa, and the direction of σ2 is mainly ESE. Near the denudation zone, the direction of σ2 is NNE. Slightly away from the denudation zone, the direction of σ2 is vertical. The maximum principal stress (σ1) values are generally between 7-81 MPa, and the direction of σ1 is mainly vertical, except in the three denudation areas, in which the direction of σ1 is ESE. The Cen'gong block is currently in the exploration stage, and widespread fracturing and horizontal wells have not yet been deployed. This paper verifies the reliability of the simulation 14

ACCEPTED MANUSCRIPT results by comparing the in-situ stress results calculated for the TM-1, TX-1 and CY-1 wells. As shown in Fig. 15, the relative error between the maximum horizontal principal stress and the minimum horizontal principal stress is not more than 10%, and the relative error of the local layers of the TX-1 well is greater than 10%. Therefore, the heterogeneous geomechanical modeling method proposed in this paper can be used to accurately predict the 3D in-situ stress field.

PT

In the study area, the principal stresses increase as the buried depth increases, but each tectonic zone has different characteristics. As shown in Fig. 16A-F, the gradient of σh is 1.16-1.76 MPa/100

RI

m, and the gradient of σH is 1.52-3.35 MPa/100 m. The gradients of σh and σH are the largest in the

SC

gentle structure areas, and the smallest in the fault areas. The gradient of σh is almost the same as that of the anticline zone, and the gradient of σH (1.56 MPa/100 m) in the syncline area is obviously

NU

smaller than that of the anticline (1.85 MPa/100 m). The simulation results show that, compared

MA

with the NS faults, the orientation of σH is more affected by the E-W faults. The orientation of σH rotates in the direction perpendicular to the fault strike in the E-W fault zone, and the rotation

D

angles are 8-24°. The orientation of σH rotates in the direction parallel to the fault strike in the

PT E

vicinity of faults, and the rotation angles are 0-20° (Figs. 16G and H). 6.2. Stress difference Δσ and stress difference coefficient Kh

CE

In fractured reservoirs, the propagation direction of hydraulic fractures is closely related to the

AC

state of stress, the occurrence of fractures, and the development degree of natural fractures. Various stress parameters, the horizontal stress difference Δσ, and horizontal stress difference coefficient Kh (Zhou et al., 2010; Wen et al., 2016; Taleghani et al., 2016; Fatahihassan et al., 2016) are widely used to determine whether a hydraulic fracture can cross natural fractures in 3D.

σ = σH  σh

(5)

σH  σh σh

(6)

Kh =

Combining Eqs. (5) and (6), the spatial distributions of the two parameters, Δσ and Kh, can be 15

ACCEPTED MANUSCRIPT obtained. The main developmental characteristics of Δσ are listed as follows: (1) the Δσ values are generally between 1.5 and 18 MPa (Fig. 17A); (2) Δσ is mainly controlled by depth, and in structurally high areas, Δσ values are low, and vice versa; and (3) the effects of different strike faults on Δσ are different, and the Δσ values of ESE-trending faults are low (Fig. 17A-F1). Additionally, NNE-trending faults have high Δσ values (Fig. 17A-F3) or little to no effect on Δσ values (Fig.

PT

17A-F2). The concentration of the principal stress is related to the trend of the fault and the mechanical parameters of the surrounding rock. The horizontal principal stress is high near the

RI

NNE-trending faults and low near the ESE-trending faults; thus, the Δσ values of the ESE-trending

SC

faults are low. For the same fault, the mechanical properties of the surrounding rock also affect the concentration of stress. When the Young's modulus of the surrounding rock is larger near the

NU

NNE-trending faults, the increase in the maximum horizontal principal stress is greater than the

MA

increase in the minimum horizontal principal stress. Hence, some of the NNE-trending faults have high Δσ values. In areas where the Young's modulus of the surrounding rock is smaller, the increase in the maximum horizontal principal stress is similar to the increase in the minimum horizontal

PT E

D

principal stress. Hence, some of the NNE-trending faults have little to no effect on the Δσ values. The distribution of Kh can be described as follows: (1) Kh values are generally between 0 and 1.2 (Fig. 17B); (2) Kh values are mainly controlled by depth, and Kh values are high in structurally

CE

high areas, and vice versa; and (3) the effects of different strike faults on the differential stress are

AC

different. In ESE-trending faults, the Kh values are low (Fig. 17B-F4). NNE-trending faults have low Δσ values (Fig. 17B-F5) or little to no effect on Kh values (Fig. 17B-F6). The distributions of Δσ and Kh do not follow simple linear relationships, as observed in Fig. 18. The relationship between Δσ and Kh can be described as follows: (1) when the depth difference between the different layers is small (<500 m), Δσ and Kh exhibit an approximately linear relationship (Figs. 18A, B, C, and D); (2) as the depth increases, the slopes of Δσ and Kh gradually increase (Figs. 18A, B, C, and D), and the blue line in Fig. 18E shows the possible conversion relationship between Δσ and Kh at different depths; (3) as the depth increases further, the correlation 16

ACCEPTED MANUSCRIPT (R2) between Kh and Δσ increases (Figs. 18A, B, C, and D); and (4) Kh is negatively exponentially correlated with depth, and Δσ is linearly positively correlated with depth (Fig. 18F). Thus, in the Cen'gong area, when the depth difference between different layers is small (<500 m), the propagation directions of hydraulic fractures in fractured reservoirs determined by Δσ or Kh are approximately equal. Conversely, when the depth difference between different layers is large (>800

PT

m), the propagation directions of hydraulic fractures predicted by Kh and Δσ may vary. The in-situ stress profile can provide a reliable basis for the stability evaluation of boreholes,

RI

the determination of casing procedures, strength design, the prediction and diagnosis of layers of

SC

deformed casing, perforation scheme optimization and water injection, and reservoir rebuilding. As illustrated in Fig. 19, the variations in horizontal principal stresses, Δσ, and Kh in each well were

NU

similar to variations in the vertical principal stress, and the principal stress value in the siliceous

MA

shale in the bottom section is high. The stress distribution, which is affected by the mechanics of rock layers, varies throughout the profile, and the development zones of natural fractures obtained

PT E

6.3. Types of in-situ stress

D

from cores are consistent with areas of high principal stresses.

Stress values can be used to infer the faulting type (Table 4) (Anderson, 1972; Wang, 1992;

CE

Reinecker et al., 2010; Wasantha et al., 2017), and the stress type coefficient Sp is defined to

AC

determine whether the stress type changes in a specific area:

Sp =

Np NT

(7)

where NT is the total number of stress types at a specified depth and N p is the amount of a certain type of stress at a specified depth. As illustrated in Fig. 20, σ1, σ2, and σ3 (Fig. 20A) in the profile can be described as follows: at depths shallower than 375 m, the stress type is mainly type II, σv<σh<σH, and the horizontal stress is the main stress. At depths ranging from 375 to 950 m, the stress type is mainly type III, and σh<σv<σH.. Finally, when the depth is greater than 950 m, the stress type is mainly type Ia, σh<σH 17

ACCEPTED MANUSCRIPT <σv, and gravity plays the major role. There are two stress transfer bands in the vertical direction (Fig. 20A): transition zone I at 375 m (σv≈σh≈σH) and transition zone II at 950 m (σh<σv≈σH). The predicted results of stress transfer bands are consistent with the statistical results of hydraulic fracturing tests (Fig. 5A-I and Fig. 5B-II). At depths ranging from -740 to -810 m (near transition zone II) in Fig. 20B, the distribution

PT

of mechanical rock parameters affects the stress type transformation. The layers with large elastic modulus values comprise the development sections of stress type Ia, and the layers with small

RI

elastic modulus values comprise the development sections of stress type Ⅲ. At depths ranging from

SC

-380 to -450 m (near transition zone I) in Fig. 20C, the layers with large elastic modulus values comprise the development area of stress type II, and the layers with small elastic modulus values

NU

comprise the development area of stress type Ⅲ. Faults also affect the distribution of stress types,

MA

which can be determined based on the stress transition interpreted from hydraulic fracturing tests (Fig. 5A-Ⅰ and Fig. 5B-Ⅱ). The simulation results show that ESE-trending faults have the largest

D

change of being characterized as stress type I. In the fault zone, the type of in-situ stress changes

PT E

from type III to type Ia (Fig. 20D-F1), and in the vicinity of a fault, the type of in-situ stress generally changes from type III to type Ia (Fig. 20D-F2).

CE

7. Conclusions

AC

(1) The 3D heterogeneous geomechanical model of the Niutitang formation in the Cen'gong block is established by considering the spatial variations in mechanical parameters and changes in the morphology and occurrence of 65 faults. Analyses of hydraulic fracturing, DIFs, borehole breakout, land movement, and well deviation show that the direction of maximum horizontal principal stress is mainly ESE 120° in the Cen'gong block. The simulation results show that the minimum principal stress values in the lower Niutitang formation are generally between 1 and 42 MPa, the intermediate principal stress values are generally between 3 and 55 MPa, and the maximum principal stress values are generally between 7 and 81 MPa. The orientation of σH rotates 18

ACCEPTED MANUSCRIPT in the direction perpendicular to the fault strike in the E-W-trending fault zone, with rotation angles of 8-24°. However, the orientation of σH rotates in the direction parallel to the fault strike in the vicinity of faults, with rotation angles of 0-20° (2) Stress difference Δσ values are generally between 1.5 and 18 MP and positively correlated with depth. The Δσ values of ESE-trending faults are low, while NNE-trending faults have high Δσ

PT

values or little to no effect on Δσ values. The stress difference coefficient Kh values are generally between 0 and 1.2 and are negatively exponentially correlated with depth. Additionally,

RI

ESE-trending faults have high Kh values. The propagation directions of hydraulic fractures in

SC

fractured reservoirs determined by Δσ or Kh are approximately equal when the depth difference between different layers is small (<500 m); however, the propagation directions of hydraulic

NU

fractures predicted by Kh and Δσ may be inconsistent when the depth difference between different

MA

layers is large (>800 m).

(3) At depths shallower than 375 m, the stress type is mainly type Ⅱ, and σv<σh<σH. At depths ranging from 375 to 950 m, the stress type is mainly type Ⅲ, and σh<σv<σH. When the depth is

PT E

D

greater than 950 m, the stress type is mainly type Ia, σh<σH <σv, and vertical stress becomes the most critical factor. There are two stress transfer bands in the vertical direction in the Cen'gong block: transition zone I at 375 m (σv≈σh≈σH) and transition zone II at 950 m (σh<σv≈σH). The

CE

principal stresses increase as the burial depth increases, but each tectonic zone has different

AC

characteristics. The gradients of σh and σH are the largest in the gentle structure areas and the smallest in the fault areas. The gradient of σh in the syncline area is almost the same as that of the anticline zone, and the gradient of σH in the syncline area is smaller than that in the anticline area. (4) Depth, fault orientation, and rock mechanics all affect the stress type. In the vicinity of transition zone Ⅰ, the layers with large elastic modulus values comprise the development sections of stress type Ⅱ, while the layers with small elastic modulus values comprise the development sections of stress type Ⅲ. In the vicinity of transition zone Ⅱ, the layers with large elastic modulus values comprise the development sections of stress type Ia, and the layers with small elastic 19

ACCEPTED MANUSCRIPT modulus values comprise the development area of stress type Ⅲ. The simulation results show that ESE-trending faults exhibit the largest change in stress type. In the ESE-trending fault zone, the type of in-situ stress changes from type III to type Ia, and in the vicinity of the faults, the type of in-situ stress may change from type III to type Ia. Acknowledgments

PT

This research was supported by the National Natural Science Foundation of China (Grant Nos. 41372139 and 41072098) and the National Science and Technology Major Project of China

RI

(2016ZX05046-003-001, 2016ZX05034-004-003). The authors would like to thank the China

SC

Earthquake Data Center for providing the data of present-day stress pattern. We are also grateful to the anonymous reviewers, whose comments improved the quality of this manuscript.

NU

References

MA

Aadnoy, B.S., Bell, J.S., 1998. Classification of drilling-induced fractures and their relationship to in-situ stress directions. Society of Petrophysicists and Well-Log Analysts, 39(6).

D

Anderson E.M., 1972.The dynamics of faulting. New York: Hafner.

PT E

Bell, J. S., Gough, D. I., 1979. Northeast-Southwest compressive stress in Alberta evidence from oil wells. Earth Planet. Sci. Lett., 45, 475–482. Bell, J.S., 1996. In-situ stresses in sedimentary rocks (part 2): applications of stress measurements.

CE

Geosci. Can. 23, 135-153.

AC

Çiftçi, N.B., 2013. In-situ stress field and mechanics of fault reactivation in the gediz graben, western turkey. Journal of Geodynamics 65(1), 136-147. Engelder, T., 1993. Stress Regimes in the Lithosphere. Princeton University Press, New Jersey. Fan, T. Y., Long, C. X., & Yang, Z. Y., 2012. Comprehensive modeling on the present crustal stress of china mainland with the viscoelastic spherical shell. Chinese Journal of Geophysics, 55(4), 1249-1260. Fatahihassan, H., Hossain, M.M., Sarmadivaleh, M., 2016. Numerical and experimental investigation of the interaction of natural and propagated hydraulic fracture. Journal of Natural Gas 20

ACCEPTED MANUSCRIPT Science and Engineering 37: 409-424. Faulkner, D., Mitchell, T., Healy, D., Heap, M., 2006. Slip on weak fault by the rotation of regional stress in the fracture damage zone. Nature, 444(7121), 922-5. Goodfellow, S.D., Young, R.P., 2014. A laboratory acoustic emission experiment under in-situ conditions. Geophysical Research Letters, 41(10), 3422–3430.

PT

Guo, P., Yao, L., Ren, D., 2016. Simulation of three-dimensional tectonic stress fields and quantitative prediction of tectonic fracture within the Damintun Depression, Liaohe Basin, northeast

RI

China. J. of Structural Geology, 86, 211-223.

SC

Haghi, A.H., Kharrat, R., Asef, M.R., Rezazadegan, H., 2013. Present-day stress of the central persian gulf: implications for drilling and well performance. Tectonophysics 608(6), 1429-1441.

NU

Haimson, B.C., Cornet, F.H., 2003. ISRM suggested methods for rock stress estimation-part 3:

MA

hydraulic fracturing (HF) and/or hydraulic testing of pre-existing fractures (HTPF).Int.J.Rock Mech.Min. 40, 1011-1020.

D

Hashimoto, C., Matsu'Ura, M., 2006. 3-D simulation of tectonic loading at convergent plate

PT E

boundary zones: Internal stress fields in northeast Japan. Pure and Applied Geophys, 163(9), 1803-1817.

Heap, M. J., Faulkner, D. R., Meredith, P. G., Vinciguerra, S., 2010. Elastic moduli evolution and

CE

accompanying stress changes with increasing crack damage: implications for stress changes around

AC

fault zones and volcanoes during deformation. Geophysical Journal International, 183(1), 225–236. Heidbach, O., Tingay, M., Barth, A., Reinecker, J., Kurfeß, D., Müller, B., 2010. Global crustal stress pattern based on the world stress map database release 2008. Tectonophysics, 482, 3–15. Houska, J., 2014. Fundamentals of rock mechanics. Geofluids, 9(3), 251–252. Huang, B., Chen, H., Han, J., 2012. Decreasing casing deformation by using cross-dipole acoustic logging data. Journal of the Acoustical Society of America, 131(4), 3370. Ito, T., Funato, A., Lin, W., Doan, M.L., Boutt, D.F., Kano, Y., 2013. Determination of stress state in deep subsea formation by combination of hydraulic fracturing in-situ test and core analysis: A case 21

ACCEPTED MANUSCRIPT study in the IODP Expedition 319. Jiu, K., Ding, W., Huang, W., You, S., Zhang, Y., Zeng, W., 2013. Simulation of paleotectonic stress fields within Paleogene shale reservoirs and prediction of favorable zones for fracture development within the Zhanhua Depression, Bohai Bay Basin, east China. J. Pet. Sci. Eng. 110, 119–131. Khair, H.A., Cooke, D., Hand, M., 2013. The effect of present day in-situ stresses and paleo-stresses

PT

on locating sweet spots in unconventional reservoirs, a case study from Moomba-big lake fields, cooper basin, south Australia. Journal of Petroleum Exploration & Production Technology 3(4),

RI

1-15.

SC

King, R., Backe, G., Tingay, M., Hillis, R., Mildren, S., 2012. Stress deflections around salt diapirs in the gulf of mexico. Geological Society London Special Publications, 367(1), 141-153.

NU

Lehtonen, A., Cosgrove, J.W., Hudson, J.A., Johansson, E., 2012. An examination of in-situ, rock

MA

stress estimation using the Kaiser effect. Engineering Geology, 124(1), 24-37. Liu, J.S. , Dai, J.S., Zou, J., Yang, H.M., Wang, B.F., Zhou, J.B., 2015. Quantitative prediction of

D

permeability tensor of fractured reservoirs. Oil Gas Geology. 36(6): 1022-1029 (In Chinese with

PT E

English abstract).

Liu, J., Ding, W., Wang, R., Yin, S., Yang, H., Gu, Y., 2017. Simulation of paleotectonic stress fields and quantitative prediction of multi-period fractures in shale reservoirs: a case study of the

CE

Niutitang Formation in the Lower Cambrian in the Cen'gong block, South China. Marine and

AC

Petroleum Geology 84, 289-310. Liu, J., Yao, Y., Elsworth, D., Pan, Z., Sun, X., Ao, W., 2016. Sedimentary characteristics of the lower cambrian Niutitang shale in the southeast margin of sichuan basin, China. J. Nat. Gas Sci. Eng., 36, 1140-1150. Mitsakaki, C., Sakellariou, M.G., Tsinas, D., 2013. A study of the crust stress field for the Aegean region (Greece). Tectonophysics, 597(4), 50–72. Mount, V.S., Suppe, J., Mount, V.S., Suppe, J., 1988. Present-day stress directions in california determined through borehole breakout analysis. AAPG Bulletin.72(3), 390-390. 22

ACCEPTED MANUSCRIPT Naimi-Ghassabian, N., Khatib, M.M., Nazari, H., Heyhat, M.R., 2015. Present-day tectonic regime and stress patterns from the formal inversion of focal mechanism data, in the north of central–east Iran blocks. Journal of African Earth Sciences, 111, 113-126. Nie, X., Zou, C., Pan, L., Huang, Z., Liu, D., 2013a. Fracture analysis and determination of in-situ stress direction from resistivity and acoustic image logs and core data in the Wenchuan earthquake

PT

fault scientific drilling borehole-2 (50–1370m). Tectonophysics, 593(3), 161-171. Nie, Y.S., Leng, J.G., Han, J.H., Sun, L., Shen, G.C., 2013b. Exploration potential of shale gas in

RI

Cen'gong block. Southeast. Guizhou Prov. Oil Gas Geol.34 (2), 274–280 (in Chinese with English

SC

abstract).

Niu Z. J, Wang M., Sun H. R., Sun J. Z., YOU X. Z., Gan W. J., Xue G. J., Hao J. X., Xin S. H.,

NU

Wang Y. Q., Wang Y. X., Li B., 2005. Contemporary velocity field of crustal movement of Chinese

MA

mainland from Global Positioning System measurements. Science Bulletin, 50 (8), 839-840. Patterson, D., Skjong, G., Wade, J.M., 1999. Horizontal Stress Orientation Analysis Using Tile

D

Cross-Dipole Acoustic Log In The Ekofisk Field. In SPWLA 40th Annual Logging Symposium.

PT E

Society of Petrophysicists and Well-Log Analysts. Rajabi, M., Tingay, M., Heidbach, O., 2016a. The present-day state of tectonic stress in the darling

776-790.

CE

basin, Australia: implications for exploration and production. Marine & Petroleum Geology 77,

AC

Rajabi, M., Tingay, M., Heidbach, O., 2016b. The present-day stress field of new south wales, australia. Australian Journal of Earth Sciences,63(1), 1-21. Rajabi, M., Tingay, M., King, R. and Heidbach, O., 2017. Present-day stress orientation in the Clarence-Moreton Basin of New South Wales, Australia: a new high density dataset reveals local stress rotations. Basin Res, 29: 622-640. Reinecker, J., Tingay, M., Müller, B., Heidbach, O., 2010. Present-day stress orientation in the molasse basin. Tectonophysics 482(1), 129-138. Smart, B.G.D., Somerville, J.M., Edlman, K., Jones, C., 2001. Stress sensitivity of fractured 23

ACCEPTED MANUSCRIPT reservoirs. Journal of Petroleum Science & Engineering, 29(1), 29-37. Taleghani A.D., Gonzalez M., Shojaei A., 2016. Overview of numerical models for interactions between hydraulic fractures and natural fractures: Challenges and limitations. Computers & Geotechnics 71:361-368. Tan, C., Jin, Z., Zhang, M., Tang, L., Jia, C., Chen, S., 2001. An approach to the present-day

PT

three-dimensional (3d) stress field and its application in hydrocarbon migration and accumulation in the Zhangqiang depression, Liaohe field, China. Marine & Petroleum Geology 18(9), 983-994.

RI

Tao, N., Wang, G., Xiao, C., Zhou, L., Deng, L., Li, R., 2016a.The in-situ stress determination from

SC

borehole image logs in the Kuqa depression. Journal of Natural Gas Science & Engineering 34, 1077-1084.

NU

Tao, N., Wang, G., Xiao, C., Zhou, L., Sun, Y., Song, H., 2016b. Determination of in-situ stress

MA

orientation and subsurface fracture analysis from image-core integration: an example from ultra-deep tight sandstone (Bsjqk formation) in the Kelasu belt, Tarim basin. Journal of Petroleum

D

Science & Engineering 147, 495-503.

PT E

Tingay, M. R. P., Morley, C. K., Hillis, R. R., Meyer, J., 2010. Present-day stress orientation in thailand's basins. Journal of Structural Geology, 32(2), 235-248. Tingay, M.R.P., Morley, C.K., Swarbrick, R.E., Damit, A.R., Hillis, R.R., King, R.C., 2009.

CE

Present-day stress and neotectonics of Brunei: implications for petroleum exploration and

AC

production. AAPG Bulletin 93(1), 75-100. Tutuncu, A.N., Bui, B.T., 2015. A coupled geomechanics and fluid flow model for induced seismicity prediction in oil and gas operations and geothermal applications. Journal of Natural Gas Science & Engineering 29, 110-124. Wang, H., Samuel, R., 2016a. 3d geomechanical modeling of salt-creep behavior on wellbore casing for presalt reservoirs (spe-166144-pa). Wang, M., Shen, Z., Niu, Z., Zhang, Z., Sun, H., Gan, W., 2003a. Contemporary crustal deformation of the Chinese continent and tectonic block model. Science China Earth Sciences 46(2), 25-40. 24

ACCEPTED MANUSCRIPT Wang, P., 1992. A geomechanical technique—types and distribution of geostress under various tectonic forces. Acta Petrolei Sinica 13(1), 1-12(in Chinese with English abstract). Wang, R., Ding, W., Zhang, Y., Wang, Z., Wang, X., He, J., 2016b. Analysis of developmental characteristics and dominant factors of fractures in lower Cambrian marine shale reservoirs: a case study of Niutitang formation in Cen'gong block, southern China. Journal of Petroleum Science &

PT

Engineering 138, 31-49. Wang, S., Ma, M., Ding, W., Lin, M., Chen, S., 2015. Approximate analytical-pressure studies on

RI

dual-porosity reservoirs with stress-sensitive permeability. Spe Reservoir Evaluation &

SC

Engineering.

Wang, S., Xu, Z., Pei, S., 2003b. Pn velocity variation beneath china continent and deep structure

NU

background for strong earthquake preparation. Chinese Journal of Geophysics, 46(6), 1114–1124.

MA

Wang, X.C., 2011. Geological conditions and key rock mechanics issues in the Western Route of South-to-North Water Transfer Project. J. Rock Mech. Geotech. Eng 3 (3), 234-243

D

Wang, X.J., Peng, S.M., Lu, B.X., Ma, J.Y., 2008. Researching earth stress field using cross-dipole

PT E

acoustic logging technology. Journal of China University of Petroleum 32(4), 42-46(in Chinese with English abstract).

Wasantha, P.L.P., Konietzky, H., Weber, F., 2017. Geometric nature of hydraulic fracture

CE

propagation in naturally-fractured reservoirs .Computers & Geotechnics, 83, 209–220.

AC

Wen, Q., Wang, S., Duan, X., Li, Y., Wang, F., Jin, X., 2016. Experimental investigation of proppant settling in complex hydraulic-natural fracture system in shale reservoirs. Journal of Natural Gas Science and Engineering 33, 70-80. Wen, S., Xiangzhi, Z., Song, J., 2005. The application of XMAC ground stress logging data in the water injection development of a oil field. World Well Logging Technology. Wu, C., Tuo, J., Zhang, M., Sun, L., Qian, Y., Liu, Y., 2016. Sedimentary and residual gas geochemical characteristics of the lower Cambrian organic-rich shales in southeastern Chongqing, China. Mar. Petrol. Geol., 75, 140-150. 25

ACCEPTED MANUSCRIPT Wu, Z., Zuo, Y., Wang, S., Chen, J., Wang, A., Liu, L., 2017. Numerical study of multi-period palaeotectonic stress fields in lower cambrian shale reservoirs and the prediction of fractures distribution: a case study of the niutitang formation in feng'gang no. 3 block, south china. Marine & Petroleum Geology, 80, 369-381. Zang, P.D.A., Stephansson, O., 2010. Stress Field of the Earth's Crust. Springer Netherlands.

PT

Zeng, L., Tang, X., Jiang, J., Peng, Y., Yang, Y., Lyu, W., 2015. Unreliable determination of in-situ stress orientation by borehole breakouts in fractured tight reservoirs: a case study of the upper

RI

Eocene Hetaoyuan formation in the Anpeng field, Nanxiang basin, China. AAPG Bulletin 99(11),

SC

1991-2003.

Zeng, W., Ding, W., Zhang, J., Zhang, Y., Guo, L., Kai, J., 2013. Fracture development in Paleozoic

NU

shale of Chongqing area (south China).part two: numerical simulation of tectonic stress field and

MA

prediction of fractures distribution. Journal of Asian Earth Sciences 75(8), 267-279. Zhao, J., Hefny, A.M., Zhou, Y.X., 2005. Hydrofracturing in-situ stress measurements in Singapore

D

granite. Int. J. Rock Mech. Min 42, 577-583

PT E

Zhao, Q., 2008. Analysis of rock mechanics characteristics in the casing damage area of the oilfield with XMAC well logging data. World Well Logging Technology. Zhao, S., 2013. The structural features and Lower Paleozoic of black shale fracture characteristics

CE

and distribution in Qianbei area Beijing, China University of Geosciences. Beijing. pp. 7–29 (in

AC

Chinese with English abstract).

Zhao, X.G., Wang, J., Qin, X.H., Cai, M., Su, R., He, J.G., 2015. In-situ stress measurements and regional stress field assessment in the Xinjiang candidate area for China's HLW disposal. Engineering Geology 197(4), 42-56. Zheng, Y., Moos, D., Tang, X.M., Dubinsky, V., Patterson, D.J., 2009. Identification of stress in formations using angles of fast and slow dipole waves in borehole acoustic logging. Acoustical Society of America Journal CA 2712271 A1. Zhou, J., Jin, Y., Chen, M., 2010. Experimental investigation of hydraulic fracturing in random 26

ACCEPTED MANUSCRIPT naturally fractured blocks. International Journal of Rock Mechanics & Mining Sciences 47(7), 1193-1199.

AC

CE

PT E

D

MA

NU

SC

RI

PT

Zoback, M.D., 2007.Reservoir geomechanics. Cambridge University Press.

27

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 1. (A) Location of the Cen'gong block. (B) Cross-section through the Cen'gong block. The location is

AC

CE

shown in (A-3).

28

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 2. The present-day stress pattern in South China. The GPS data were provided by Wang et al. (2003a) and Niu et al. (2005); the other data were provided by the China Earthquake Data Center. The different symbols

D

show the orientation of the maximum horizontal stress, and the different colors show the different types of stress

AC

CE

PT E

regimes (normal faulting; strike–slip faulting; thrust faulting). The location is shown in (Fig. 1A-1).

29

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 3. (A) Stratigraphic and lithological systems in the study area. (B) Lithological systems and depositional

AC

CE

PT E

D

MA

environments of different members of the Niutitang formation.

30

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 4. Natural fractures in the Niutitang formation of the Cen'gong block. (A) The vertical shear fractures sealed with calcite are well developed in the siliceous shale in the TX-1 well (depth 1784.50 m). (B) Structural fractures with high dip angles in outcrop. (C) Microfractures in a thin section (the sample was observed with a polariscope), fractures sealed with calcite, and unsealed fractures that developed in the calcite in the CY-1 well, where blue (color version) is the epoxy resin (depth 1458.25 m). (D) Strike of the fractures in the outcrop, dominated by NE-SW, N-S, and NW-SE fractures. (E) Dip angles of the fractures in the outcrop showing the dominance of vertical fractures. (F) Strike of the fractures from the FMI log interpretation, which are similar to 31

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

the outcrop fractures.

32

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 5. Variations in the principal stresses with depth as interpreted from the hydraulic fracturing tests conducted in (A) the ZK4220 well, (B) the CD-1 well, and (C) the TX-1 well. (D) Relationship between the minimum horizontal principal stress and the burial depth based on data from the ZK4220, CD-1, and TX-1 wells. (E) Relationship between the maximum horizontal principal stress and the burial depth based on data from the ZK4220, CD-1, and TX-1 wells. (F) Relationship between rock density and depth in the TM-1 well. (G) Relationship between rock density and depth in the TX-1 well. STZ=stress transition zone. 33

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Fig. 6. Comprehensive interpretation of the in situ stress and mechanical parameters based on nuclear

AC

CE

PT E

D

SM=shear modulus; and BM=bulk modulus.

MA

magnetic resonance (NMR) and XMAC logging. DP=density porosity; NMP=nuclear magnetic porosity;

34

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

D

Fig. 7. (A) Rose diagram of the well deviations in the TX-1 well. (B) Frequency distribution of the deviation

PT E

angles in the TX-1 well. (C) Rose diagram of the borehole breakout azimuths derived from the six-arm caliper

AC

CE

logs in the TM-1 well. (D) Rose diagram of the land movement azimuths derived from the GPS data.

35

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 8. (A) DIFs in the resistivity image log of the HY-1 well. (B) Rose diagram of the strikes of the DIFs. (C)

AC

CE

Frequency distribution of the dip angles of the DIFs.

36

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 9. The hydro-fractures based on the impresser orientation in the ZK4220 well. (A) Hydro-fractures at a

D

depth of 248.90 m. (B) Hydro-fractures at a depth of 347.20 m. (C) Hydro-fractures at a depth of 425.40 m. (D)

AC

CE

PT E

Hydro-fractures at a depth of 653.00 m (the blue dotted lines are "base lines").

37

RI

PT

ACCEPTED MANUSCRIPT

Fig. 10. Mathematical relationships between the mechanical parameters in the Niutitang formation of the

SC

Cen'gong block. (A) Relationship between the dynamic Young’s modulus (Ed) and Poisson's ratio (μd) in the TX-1

NU

well. (B) Relationship between the dynamic Young’s modulus (Ed) and static Young’s modulus (E) in the TX-1

AC

CE

PT E

D

MA

well.

38

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 11. (A) Average P-wave velocity in different strata. (B) P-wave velocity distribution of the lower part of

AC

the Niutitang formation. (C) Relationship between P-waves and S-waves in XMAC logging.

39

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

Fig. 12. Planar distribution of the Young’s modulus of the lower part of the Niutitang formation.

40

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 13. (A) Mechanical model and boundary conditions for the stress simulation in the Cen'gong block. (B) Homogeneous geomechanical model of the lower part of the Niutitang formation. In adjacent elements, the uniform colors represent the same mechanical attributes. (C) Heterogeneous geomechanical model of the lower part of the Niutitang formation. In adjacent elements, the colors represent different mechanical attributes; in this model, each unit is given a different mechanical property that corresponds to Fig. 12.

41

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 14. In-situ stress field of the lower Niutitang formation within the Cen'gong block. (A) Minimum principal stress. (B) Direction of the minimum principal stress. (C) Intermediate principal stress. (D) Direction of intermediate principal stress. (E) Maximum principal stress. (F) Direction of the maximum principal stress. 42

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 15. (A) Comparison of the horizontal minimum principal stress results based on the simulation and the

CE

logging calculations. (B) Comparison of the horizontal maximum principal stress results based on the simulation

calculations.

AC

and the logging calculations. (C) Relative error for the principal stresses for the simulation and logging

43

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 16. Relationship between the horizontal principal stresses and burial depth of the Niutitang shale. (A) Anticline area (A1). (B) Anticline area (A2). (C) Anticline area (A3). (D) Syncline area (A4). (E) Gentle area (A5). (F) Fault area (A6). See Figs. 14A and C for the locations of A1, A2, A3, A4, A5 and A6. (G) The rotation of the maximum horizontal stress in the vicinity of the E-W-trending faults in the CD area. (H) Schematic plan view for the rotation of the maximum horizontal stress in the vicinity of the fault zone in the Cen'gong block (modified from Bell, 1996; Rajabi et al., 2017); see Fig. 1A for the location of the CD area. 44

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

45

RI

PT

ACCEPTED MANUSCRIPT

SC

Fig. 17. (A) Planar distribution of the horizontal stress difference Δσ. (B) Planar distribution of the horizontal

AC

CE

PT E

D

MA

NU

stress difference coefficient Kh.

46

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 18. Panels (A)-(E) show the relationships between ∆σ and Kh at depths of 0 to 500 m, -900 to 1400 m,

AC

-1800 to -2300 m, -2700 to -3200 m, and 0 to -3200 m, respectively. (F) The relationships between ∆σ and depth and between Kh and depth.

47

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 19. Distributions of principal stresses, the horizontal stress difference ∆σ, and the horizontal stress difference coefficient Kh in the TM-1, TX-1, and CY-1 wells. 48

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

49

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 20. (A) Vertical distributions of σ1, σ2, σ3, and stress type. (B) Vertical correlation between the Young's modulus and stress type in the P1 area. (C) Vertical correlation between the Young's modulus and stress type in the P2 area. (D) The influence of faults on the stress type distribution in the P3 area. See Fig. 1A for the locations of P1, P2, and P3.

50

ACCEPTED MANUSCRIPT Table 1 Orientation of the maximum horizontal stress around the Cen'gong block, based on the World Stress Map project.

Stress indicators

Quality

1.0

-

U

E

FM

15.0

87

SS

D

119

HF

1.0

-

U

E

116.000

98

FM

15.0

40

SS

D

26.566

115.798

93

OC

0.4

-

U

E

CHUN

29.660

115.710

95

FM

17.8

-

SS

C

CH483

25.850

115.573

135

OC

0.4

-

U

E

CH484

25.659

115.525

123

OC

0.5

-

U

E

CH 43

23.200

115.500

125

FM

15.0

60

SS

D

CH487

30.085

114.958

120

OC

0.0

-

U

E

CH741

30.133

114.917

129

OC

0.0

-

U

E

CH105

23.700

114.700

107

FM

5.0

-

SS

C

CH421

23.583

113.583

133

HF

1.0

-

U

E

CH423

23.583

113.583

124

HF

1.0

-

U

E

CH 44

23.500

112.500

126

FM

15.0

59

SS

D

CH 40

28.000

112.000

80

FM

15.0

27

SS

D

CH136

22.300

111.800

114

FM

5.0

-

SS

D

CH496

27.267

111.575

145

OC

0.0

-

U

E

CH652

27.267

111.563

145

OC

0.0

-

U

E

CH485

25.023

107.167

130

OC

0.0

-

U

E

CH655

25.017

107.167

130

OC

0.0

-

U

E

CH 46

23.500

107.000

129

FM

15.0

44

SS

D

30.500

119.367

140

HF

CH 42

25.000

117.500

134

CH419

24.450

116.667

CH 39

27.500

CH482

Depth (m)

D

MA

NU

CH415

CE

Longitude

PT

Regime

Latitude

SC

Number

RI

SHmax Orientation (°)

PT E

Well name

Well Location

AC

In Table 1, FM=focal mechanism; HF=hydraulic fracturing; BO=borehole breakout; OC is a stress relief method by overscoring; SS=strike–slip faulting; NF=normal faulting; U=unknown. The quality ranking system of the WSM project (Heidbach et al., 2010) was used to assess the reliability of the results in each well (Table 1 and 3).

51

ACCEPTED MANUSCRIPT Table 2 Mechanical properties of the Niutitang shale in TX-1 well.

1785.91 ~1791.80

TX-1N-4

1795.81 ~1796.33

TX-1N-5

TX-1N-6

1801.54 ~1802.37 1808.19 ~1812.70

24.1 32.3

SV 0.18

27.8

0.21

27.1

0.19

23.2

0.23

26.6

23.6

20.9 29.2 25.6

0.19

27.4

0.20 0.18

21.7

0.19 0.21

24.5

0.21 25.6

0.2

26.8

0.24

30.2

0.20

32.2

0.21

0.19

0.17

20.8 25.4

0.19

0.20

20.5 23.8

Mean

31.3

31.5

PT

TX-1N-3

Mean

RI

TX-1N-2

1779.80 ~1785.66

SV

0.19

SC

TX-1N-1

1770.00 ~1776.38

μ

E/GPa

NU

Depth /m

MA

Samples

0.21

0.22

0.22

0.24

AC

CE

PT E

D

In the Table 2, E is the rock Young's modulus; μ is the Poisson's ratio; and SV = single value.

52

ACCEPTED MANUSCRIPT Table 3 Orientation and quality of the maximum horizontal stress in the Cen'gong block and its circumjacent regions. Stress Depth (m)

Combined

Mean σH

lengths

Azimuth

Region

SD

Quality

Local

24°

C

PT

Well/place

Number

Local

46°

E

124°

Local

21°

C

115°

Local

17°

B

111°

Local

38°

D

-

123°

Local

19°

B

-

117°

Local

14°

B

118°

Regional

10°

A

indicators

TX-1

WD

50-550

20

500 m

104°

TX-1

WD

550-1050

20

500 m

141°

TX-1

WD

1050-1550

20

500 m

TX-1

WD

1550-1875

14

325 m

TM-1

BO

480-1521

8353

1044 m

ZK4220

HF

248-653

7

HY-1

DIF

2985-2303

21

NG

GPS

-

MA

NU

SC

RI

name

24

-

D

In Table 3, WD=well deviation; BO=borehole breakout; DIF=drilling-induced fracture;

PT E

HF=hydraulic fracturing; SD=standard deviation of the stress orientation; A, B, C, D and E are the WSM quality ranking; NG=North Guizhou. See Fig. 1A for the well and place locations. The GPS

CE

data were obtained from Wang et al. (2003a) and Niu et al. (2005). Note that due to changes in the stress regimes, the in situ stress orientation determination using the well deviation data is unreliable

AC

at depths ranging from 550 to 1050 m.

53

ACCEPTED MANUSCRIPT Table 4 The stress type classification proposed by Anderson (1972) and Wang (1992). Subtypes

Characteristics

Fault styles

Ⅰ

Ⅰa

σv>σH>σh; σh>0

Normal faulting

Ⅰ

Ⅰb

σv>σH>σh; σh<0

Normal faulting

Ⅱ

Ⅱ

σH>σh>σv

Thrust faulting

Ⅲ

Ⅲ

σH >σv>σh

Strike-slip faulting

AC

CE

PT E

D

MA

NU

SC

RI

PT

Types

54

ACCEPTED MANUSCRIPT Highlights  The direction of maximum horizontal principal stress in the Cen'gong block was determined. 3D geomechanical parameters were established based on seismic inversion. 3D heterogeneous geological models of the in situ stress fields were developed.

PT

 3D distributions of the horizontal stress difference, horizontal stress difference coefficient and stress type coefficient.

AC

CE

PT E

D

MA

NU

SC

RI

 Depth, fault orientation, and rock mechanical parameters all affect the stress type.

55