Evaluation of Favorable Fracture Area of Deep Coal Reservoirs Using a Combination of Field Joint Observation and Paleostress Numerical Simulation: A Case Study in the Linxing Area

Abstract

The development of fractures under multiple geological tectonic movements affects the occurrence and efficient production of free gas in deep coal reservoirs. Taking the No.8 deep coal seam of the Benxi formation in the Linxing area as the object, a method for evaluating favorable fracture areas is established based on the combination of field joint staging, paleogeological model reconstruction under structural leveling, finite element numerical simulation, and fracture development criteria. The results show that a large number of shear fractures and fewer tensile joints are developed in the Benxi formation in the field and mainly formed in the Yanshanian and Himalayan periods. The dominant strikes of conjugate joints in the Yanshanian period are NWW (100°~140°) and NNW (150°~175°), with the maximum principal stress magnitude being 160 MPa along the NW orientation. Those in the Himalayan period are in the NNE direction (0°~40°) and the EW direction (80°~110°), with the maximum principal stress magnitude being 100 MPa along the NE orientation. The magnitudes of the maximum principal stress of the No. 8 deep coal seam in the Yanshanian period are between −55 and −82 MPa, indicative of compression; those in the Himalayan period are from −34 to −70 MPa in the compressive stress form. Areas with high shear stress values are mainly distributed in the central magmatic rock uplift, indicating the influence of magmatic rock uplift on in situ stress distribution and fracture development. Based on the comprehensive evaluation factors of fractures, the reservoir is divided into five classes and 24 favorable fracture areas. Fractures in Class I areas and Class II areas are relatively well developed and were formed under two periods of tectonic movements. The method for evaluating favorable fracture areas is not only significant for the prediction of fractures and free gas contents in this deep coal reservoir but also has certain reference value for other reservoirs.

1. Introduction

Deep coalbed methane (deeper than 1500 m) has enormous resource potential and has achieved multiple breakthroughs in basins including Piceance basin in the U.S., the Alberta basin in Canada, the Cooper basin in Australia, the Junggar basin, the Ordos basin, and the Sichuan basin in China [1,2,3]. Exploration practice has shown that most deep coal reservoirs belong to dry coal systems with low water content and relatively high free gas content [4,5,6,7,8]. As the main spaces for the occurrence of free gas, the formation and development of fractures are directly influenced by the stress field of different geological periods, and clarification of paleostress field distribution characteristics is of great significance for the prediction of favorable fracture areas [9,10,11]. Previous studies have mainly focused on the formation mechanism, controlling factors, characterization methods and prediction methods of fractures [12,13,14,15,16], and the numerical simulation method for 2D or 3D models combined with field fracture data, acoustic emission and rock mechanics experiments, imaging logging, and fracture criteria is widely used to predict paleostress [17,18,19]. However, geological models used in the numerical simulation are based on the current burial depth of reservoirs, which is a superimposed product of multiple tectonic movements. Additionally, the influence of subsequent tectonic movements on the stress numerical simulation of the early period cannot be eliminated. Therefore, two key problems still need to be examined in depth: (1) the formation and distribution characteristics of fractures under the superposition of multiple structure movements of deep coal reservoirs; and (2) techniques for evaluating favorable fracture areas constrained by multi-stage structural movements of deep coal reservoirs.

In this study, two periods of tectonic fractures are identified based on field joint observations; the maximum principal stress direction and magnitude are obtained using a joint staging and conjugate joint angle estimation method; 3D heterogeneous geological models of deep coal reservoirs in two periods are established using the tectonic trace recovery method; the distributions of stress fields of deep coal reservoirs during different tectonic periods are obtained based on the finite element method; and the fracture development characteristics of different periods in deep coal reservoirs are predicted based on Mohr–Coulomb and Griffith criteria. The results will be beneficial for the effective exploration and development of deep coalbed methane.

2. Geological Setting

The Linxing area is located in the northern part of Jinxi Fold at the eastern margin of the Ordos Basin, with a generally southwest oriented monocline structure that dips westward at 1–5° (Figure 1a,b). It has Cenozoic, Mesozoic, and Paleozoic strata from top to bottom, and the Upper Carboniferous and Lower Permian strata are exposed in the area (Figure 1c). The No. 8 + 9 coal seam in the Benxi formation is the main production layer of coalbed methane, with a thickness varying between 2.5 and 15.5 m and a depth varying between 1093 and 2114 m. The No. 8 + 9 coal seam has undergone four uplifts under multi-stage tectonic movements, resulting in the development of multi-directional folds and faults, and uplift amplitudes in the latter two periods are relatively larger (Figure 1d), reflecting stronger tectonic movements during the Yanshanian and Himalayan periods [20,21,22]. In addition, due to the magmatic activity during Yanshanian period, the Zijinshan magmatic pluton formed an uplift (magmatic rock uplift) in the central part of the study area [23,24].

3. Methodology

3.1. Field Joint Observations

Field joint observation is performed to obtain the fracture characteristics (strike, dip direction, and dip angle) of Carboniferous–Permian rocks in the study area and its surrounding areas, including the Baode–Palougou Section, Fugu–Sunjiagou Section, Xingxian–Guanjiaya Section, and Liulin–Chengjiazhuang Section. The observed fractures mainly include shear and tension joints from sandstone, mudstone, limestone, and the coal seam. The GPS is used for positioning, and the compass is used to measure the strike, dip direction, and dip angle of joints.

3.2. Staging of Joints

The tectonic movements of different stages form different types of joint combinations and the staging of joints is essential for the determination of the paleostress direction. In the article, the regional geological setting, fracture occurrence, and the intersection relationship (staggering, limiting, intercutting, tracking, utilizing, transforming) of joints are combined to determine the formation sequence of joints and clarify the joints combination in the same tectonic period.

3.3. Numerical Simulation of Paleostress Fields

The paleostress field is the stress field of the paleogeological period, which is influenced by the geological body, mechanical properties, and boundary conditions (direction and magnitude of principal stress) [25,26,27]. Based on the restoration of the geological model, mechanical parameters and boundary conditions of the paleogeological period, the finite element method, and ANSYS software (18.0 version) are used to produce 3D simulations of the paleostress fields.

3.3.1. Paleogeological Model

The geological body went through a process of sedimentation, compaction, folding, and fracture during the sedimentary evolution. To study the paleostress state of a geological body, the restoration of the paleogeological model is necessary, including defaultization, defolding, and decompaction. The construction-flattening method [20], based on the superposition theory of waves, is used to flatten the folds of geological bodies and restore the geological bodies of different geological periods, which is effective in paleostress studies.

In this study, the structural flattening method is used to restore paleogeological models based on the contour lines of the 8 + 9 # coal seam floor of current period. Additionally, a combination strata of roof–coal seam–floor is adopted in the model to consider the influence of the roof and floor on the in situ stress of the coal seam. Additionally, the geological model is discretized into hexahedron elements, with 70,529 elements and 24,430 nodes.

3.3.2. Mechanical Parameters

The mechanical parameters of the 8 + 9 # coal seam are calculated by logging interpretation under the constraints of triaxial compression tests, including Young’s modulus, Poisson’s ratio, compressive strength, tensile strength, shear strength, cohesion, and internal friction angle of the coal seam (Formulas (1)–(6)) [28,29,30].

$$E_d=\frac{{10}^3\rho v_s^2\left[{3\left(v_p/v_s\right)}^2-4\right]}{\left[{\left(v_p/v_s\right)}^2-1\right]}$$

(1)

$${\mu }_d=\frac{{\left(v_p/v_s\right)}^2-2}{2\left[{\left(v_p/v_s\right)}^2-1\right]}$$

(2)

$${\sigma }_c=12{\sigma }_t=0.0045E_d\left(1-V_{sh}\right)+0.008E_d{V}_{sh}$$

(3)

$$k=0.026{\sigma }_c/\left[\frac{3\left(1-2{\mu }_d\right)}{E_d}\times \frac{\psi }{1-\psi }\times {10}^6\right]$$

(4)

$$C=5.44\times {10}^{-3}{\rho }^2{v}_p^4{\left(\frac{1+{\mu }_d}{1-{\mu }_d}\right)}^2\left(1-2{\mu }_d\right)\left(1+0.78V_{sh}\right)$$

(5)

$$\phi =90-\frac{360}{\pi }arctan\left(1/\sqrt{4.73-0.098\psi }\right)$$

(6)

where E_d represents Young’s modulus, MPa; ρ represents logging density, g/cm³; V_p represents primary wave velocity of logging, km/s; vs. represents secondary wave velocity of logging, km/s; μ_d represents Poisson’s ratio; σ_c represents uniaxial compressive strength, MPa; σ_t represents uniaxial tensile strength, MPa; k represents shear strength, MPa; V_sh represents volume percent of shale, %; ψ represents porosity, %; C represents cohesion, MPa; and φ represents internal friction angle, °.

3.3.3. Boundary Conditions

Boundary conditions include the determination of stress direction, stress magnitude, and displacement constraints on the boundary of geological model. The stress direction is determined by conjugate shear joint strike of different tectonic stages. The stress magnitude is calculated based on the relationship expression between the conjugate shear joint angle and the stress magnitude (Formulas (7) and (8)) [31,32].

$${\sigma }_1={\sigma }_t-\frac{k^2}{4{\sigma }_t}+\frac{k^2}{2{\sigma }_t}\left(\frac{1}{cos\theta }-\frac{1}{2{cos}^2\theta }\right)$$

(7)

$${\sigma }_3={\sigma }_t-\frac{k^2}{4{\sigma }_t}-\frac{k^2}{2{\sigma }_t}\left(\frac{1}{cos\theta }+\frac{1}{2{cos}^2\theta }\right)$$

(8)

where σ₁ represents the maximum principal stress, MPa; σ₃ represents the minimum principal stress, MPa; and θ represents the conjugate angle of conjugate shear joint, °.

3.4. Fracture Development Criteria

The fracture is formed in coal reservoirs when the paleostress field reaches the initiate threshold [33,34,35]. Based on the Mohr–Coulomb and Griffith criteria, the shear fracture coefficient ($C_S$, Formulas (9)–(11)) and the tensile fracture coefficient ($C_T$, Formula (2)) are established to evaluate the development degree of shear fracture and tensile fracture, respectively.

$$C_S=\left\{{\sigma }_1-\left[{\sigma }_3{tan}^2\left({45}^{°}+\frac{\phi }{2}\right)+2C·tan\left({45}^{°}+\frac{\phi }{2}\right)\right]\right\}/{\sigma }_1$$

(9)

$$\mathrm{When} {\sigma }_1+3{\sigma }_3\le 0 C_T=\left({\sigma }_3-{\sigma }_t\right)/{\sigma }_3$$

(10)

$$\mathrm{When} {\sigma }_1+3{\sigma }_3>0 C_T=\left(\frac{{\left({\sigma }_1-{\sigma }_3\right)}^2}{8\left({\sigma }_1+{\sigma }_3\right)}+{\sigma }_t\right)/{\sigma }_3$$

(11)

where C_S represents the shear fracture coefficient, and C_T represents the tensile fracture coefficient.

It is obvious that when C_S is less than 0, the rock has not undergone shear fracture. When C_S is greater than or equal to 0, the rock has undergone shear fracture, and the degree of fracture development increases with the increase in the C_S value; Similarly, when C_T is less than 0, the rock has not undergone tensile fracture. When C_T is greater than or equal to 0, the rock has undergone tensile fracture, and the fracture development degree increases with the increase in the C_T value.

4. Results and Discussions

4.1. Field Joints’ Characteristics

Field joint observations are conducted on 81 geological points, and a total of 230 sets of joint orientation observation data are obtained (Appendix A 

Table A1. Field joint characteristics in the Linxing area and surroundings.

Number	Latitude	Longitude	Dip Direction and Angle	Number	Latitude	Longitude	Dip Direction and Angle
1	N 38°45′32.96″	E 111°08′13.97″	45°∠78°	116	N 39°03′31.40″	E 111°06′11.67″	325°∠69°
2	N 38°45′32.96″	E 111°08′13.97″	75°∠84°	117	N 39°02′56.54″	E 111°06′38.80″	235°∠70°
3	N 38°45′30.10″	E 111°08′9.74″	255°∠77°	118	N 39°02′56.54″	E 111°06′38.80″	151°∠86°
4	N 38°45′30.10″	E 111°08′9.74″	102°∠55°	119	N 39°02′56.54″	E 111°06′38.80″	316°∠73°
5	N 38°45′27.78″	E 111°08′6.06″	75°∠78°	120	N 39°02′56.54″	E 111°06′38.80″	225°∠69°
6	N 38°45′27.78″	E 111°08′6.06″	168°∠82°	121	N 39°02′25.24″	E 111°05′50.44″	187°∠84°
7	N 38°45′28.63″	E 111°08′3.04″	65°∠57°	122	N 39°02′25.24″	E 111°05′50.44″	86°∠81°
8	N 38°45′28.63″	E 111°08′3.04″	178°∠64°	123	N 39°02′21.39″	E 111°05′45.30″	184°∠76°
9	N 38°45′27.89″	E 111°07′51.18″	78°∠65°	124	N 39°02′21.39″	E 111°05′45.30″	77°∠54°
10	N 38°45′27.89″	E 111°07′51.18″	342°∠76°	125	N 39°02′21.39″	E 111°05′45.30″	189°∠71°
11	N 38°45′28.89″	E 111°07′47.44″	260°∠85°	126	N 39°02′21.39″	E 111°05′45.30″	76°∠67°
12	N 38°45′28.89″	E 111°07′47.44″	50°∠71°	127	N 39°02′57.69″	E 111°05′17.15″	102°∠81°
13	N 38°45′26.42″	E 111°07′43.21″	47°∠73°	128	N 39°02′57.69″	E 111°05′17.15″	182°∠85°
14	N 38°45′26.42″	E 111°07′43.21″	268°∠68°	129	N 39°02′57.69″	E 111°05′17.15″	98°∠80°
15	N 38°45′26.94″	E 111°07′41.98″	50°∠53°	130	N 39°02′57.69″	E 111°05′17.15″	65°∠81°
16	N 38°45′26.94″	E 111°07′41.98″	260°∠79°	131	N 39°02′49.33″	E 111°04′32.01″	324°∠82°
17	N 38°45′25.98″	E 111°07′40.60″	48°∠78°	132	N 39°02′49.33″	E 111°04′32.01″	345°∠78°
18	N 38°45′25.98″	E 111°07′40.60″	242°∠69°	133	N 39°02′49.33″	E 111°04′32.01″	60°∠79°
19	N 38°45′24.68″	E 111°07′34.20″	76°∠70°	134	N 39°02′49.33″	E 111°04′32.01″	68°∠75°
20	N 38°45′24.68″	E 111°07′34.20″	340°∠86°	135	N 39°02′49.33″	E 111°04′32.01″	347°∠76°
21	N 38°45′25.90″	E 111°07′29.62″	55°∠82°	136	N 39°02′21.47″	E 111°03′17.23″	335°∠83°
22	N 38°45′25.90″	E 111°07′29.62″	159°∠88°	137	N 39°02′21.47″	E 111°03′17.23″	62°∠78°
23	N 38°45′22.51″	E 111°07′24.62″	37°∠77°	138	N 39°02′21.47″	E 111°03′17.23″	325°∠87°
24	N 38°45′22.51″	E 111°07′24.62″	292°∠54°	139	N 39°02′21.47″	E 111°03′17.23″	84°∠77°
25	N 38°45′21.51″	E 111°07′22.61″	22°∠82°	140	N 39°02′11.29″	E 111°02′35.42″	335°∠82°
26	N 38°45′21.51″	E 111°07′22.61″	258°∠66°	141	N 39°02′11.29″	E 111°02′35.42″	54°∠51°
27	N 38°46′1.61″	E 111°04′9.33″	315°∠87°	142	N 39°02′11.29″	E 111°02′35.42″	63°∠54°
28	N 38°46′1.61″	E 111°04′9.33″	43°∠82°	143	N 39°02′11.29″	E 111°02′35.42″	358°∠85°
29	N 38°46′2.58″	E 111°04′11.09″	12°∠73°	144	N 38°30′14.85″	E 111°10′14.05″	97°∠56°
30	N 38°46′2.58″	E 111°04′11.09″	115°∠61°	145	N 38°30′14.85″	E 111°10′14.05″	23°∠82°
31	N 38°46′0.10″	E 111°04′11.78″	149°∠72°	146	N 38°30′14.85″	E 111°10′14.05″	352°∠87°
32	N 38°46′0.10″	E 111°04′11.78″	52°∠74°	147	N 38°30′14.85″	E 111°10′14.05″	94°∠73°
33	N 38°45′59.16″	E 111°04′12.44″	56°∠71°	148	N 38°30′15.85″	E 111°10′14.41″	22°∠85°
34	N 38°45′59.16″	E 111°04′12.44″	313°∠79°	149	N 38°30′15.85″	E 111°10′14.41″	87°∠81°
35	N 38°45′58.49″	E 111°04′20.51″	343°∠76°	150	N 38°30′12.00″	E 111°10′11.72″	70°∠82°
36	N 38°45′58.49″	E 111°04′20.51″	65°∠78°	151	N 38°30′12.00″	E 111°10′11.72″	142°∠83°
37	N 38°45′58.10″	E 111°04′22.05″	337°∠89°	152	N 38°30′1.18″	E 111°10′11.06″	78°∠75°
38	N 38°45′58.10″	E 111°04′22.05″	75°∠83°	153	N 38°30′1.18″	E 111°10′11.06″	20°∠82°
39	N 38°45′58.41″	E 111°04′23.18″	52°∠63°	154	N 38°30′1.02″	E 111°10′8.31″	84°∠82°
40	N 38°45′58.41″	E 111°04′23.18″	73°∠76°	155	N 38°30′1.02″	E 111°10′8.31″	341°∠81°
41	N 38°45′58.41″	E 111°04′23.18″	352°∠51°	156	N 38°30′0.25″	E 111°10′7.68″	57°∠78°
42	N 38°45′58.41″	E 111°04′23.18″	70°∠81°	157	N 38°30′0.25″	E 111°10′7.68″	120°∠71°
43	N 38°45′58.41″	E 111°04′23.18″	353°∠76°	158	N 38°29′57.27″	E 111°10′4.44″	43°∠75°
44	N 38°45′58.41″	E 111°04′23.18″	74°∠85°	159	N 38°29′57.27″	E 111°10′4.44″	86°∠86°
45	N 38°45′58.41″	E 111°04′23.18″	355°∠87°	160	N 38°29′57.27″	E 111°10′4.44″	43°∠83°
46	N 38°45′58.41″	E 111°04′23.18″	75°∠76°	161	N 38°29′57.27″	E 111°10′4.44″	335°∠81°
47	N 38°45′58.97″	E 111°04′40.51″	321°∠88°	162	N 38°29′55.65″	E 111°10′2.87″	64°∠86°
48	N 38°45′58.97″	E 111°04′40.51″	46°∠67°	163	N 38°29′55.65″	E 111°10′2.87″	26°∠82°
49	N 38°45′58.97″	E 111°04′40.51″	330°∠79°	164	N 37°26′27.77″	E 110°54′12.61″	47°∠79°
50	N 38°45′58.97″	E 111°04′40.51″	53°∠77°	165	N 37°26′27.77″	E 110°54′12.61″	125°∠64°
51	N 38°46′2.86″	E 111°05′14.02″	201°∠88°	166	N 37°26′39.75″	E 110°53′51.22″	145°∠84°
52	N 38°46′2.86″	E 111°05′14.02″	76°∠68°	167	N 37°26′39.75″	E 110°53′51.22″	62°∠89°
53	N 38°46′2.86″	E 111°05′14.02″	152°∠88°	168	N 37°26′41.60″	E 110°53′49.24″	16°∠55°
54	N 38°46′2.86″	E 111°05′14.02″	77°∠72°	169	N 37°26′41.60″	E 110°53′49.24″	107°∠80°
55	N 38°46′2.86″	E 111°05′14.02″	206°∠86°	170	N 37°26′41.60″	E 110°53′49.24″	117°∠83°
56	N 38°46′2.86″	E 111°05′14.02″	106°∠59°	171	N 37°26′41.60″	E 110°53′49.24″	44°∠75°
57	N 38°45′27.48″	E 111°06′44.02″	321°∠84°	172	N 37°26′41.60″	E 110°53′49.24″	86°∠88°
58	N 38°45′27.48″	E 111°06′44.02″	73°∠66°	173	N 37°33′30.36″	E 110°53′51.46″	298°∠80°
59	N 38°45′26.41″	E 111°06′46.27″	144°∠81°	174	N 37°33′30.36″	E 110°53′51.46″	195°∠84°
60	N 38°45′26.41″	E 111°06′46.27″	81°∠69°	175	N 37°33′35.10″	E 110°53′33.25″	290°∠87°
61	N 38°45′26.41″	E 111°06′46.27″	346°∠88°	176	N 37°33′35.10″	E 110°53′33.25″	75°∠75°
62	N 38°45′26.41″	E 111°06′46.27″	74°∠74°	177	N 37°33′35.07″	E 110°53′33.50″	297°∠86°
63	N 38°45′26.41″	E 111°06′46.27″	155°∠76°	178	N 37°33′35.07″	E 110°53′33.50″	194°∠89°
64	N 38°45′26.41″	E 111°06′46.27″	74°∠65°	179	N 37°33′8.90″	E 110°51′58.80″	52°∠83°
65	N 38°45′27.10″	E 111°06′47.15″	4°∠65°	180	N 37°33′8.90″	E 110°51′58.80″	141°∠68°
66	N 38°45′27.10″	E 111°06′47.15″	97°∠73°	181	N 37°33′9.01″	E 110°51′58.85″	183°∠89°
67	N 38°45′27.10″	E 111°06′47.15″	341°∠82°	182	N 37°33′9.01″	E 110°51′58.85″	81°∠66°
68	N 38°45′27.10″	E 111°06′47.15″	74°∠78°	183	N 37°33′10.65″	E 110°51′55.53″	75°∠70°
69	N 38°45′24.64″	E 111°06′56.57″	322°∠82°	184	N 37°33′10.65″	E 110°51′55.53″	155°∠78°
70	N 38°45′24.64″	E 111°06′56.57″	74°∠67°	185	N 37°33′9.98″	E 110°51′54.10″	196°∠89°
71	N 38°45′22.63″	E 111°07′0.53″	123°∠68°	186	N 37°33′9.98″	E 110°51′54.10″	81°∠75°
72	N 38°45′22.63″	E 111°07′0.53″	52°∠74°	187	N 37°33′10.09″	E 110°51′51.38″	85°∠79°
73	N 38°45′21.75″	E 111°07′3.08″	18°∠61°	188	N 37°33′10.09″	E 110°51′51.38″	184°∠83°
74	N 38°45′21.75″	E 111°07′3.08″	86°∠74°	189	N 37°33′9.76″	E 110°51′50.04″	183°∠84°
75	N 38°45′21.75″	E 111°07′3.08″	150°∠78°	190	N 37°33′9.76″	E 110°51′50.04″	285°∠81°
76	N 38°45′21.75″	E 111°07′3.08″	76°∠72°	191	N 37°33′9.76″	E 110°51′50.04″	78°∠89°
77	N 38°45′20.86″	E 111°07′14.95″	146°∠75°	192	N 37°33′9.76″	E 110°51′50.04″	152°∠73°
78	N 38°45′20.86″	E 111°07′14.95″	85°∠89°	193	N 37°32′55.11″	E 110°49′30.59″	57°∠81°
79	N 39°03′23.56″	E 111°07′6.74″	188°∠56°	194	N 37°32′55.11″	E 110°49′30.59″	129°∠84°
80	N 39°03′23.56″	E 111°07′6.74″	100°∠76°	195	N 37°32′55.11″	E 110°49′30.59″	51°∠79°
81	N 39°03′19.84″	E 111°07′4.13″	219°∠82°	196	N 37°32′55.11″	E 110°49′30.59″	142°∠80°
82	N 39°03′19.84″	E 111°07′4.13″	127°∠84°	197	N 37°35′33.50″	E 110°53′10.92″	290°∠74°
83	N 39°03′20.39″	E 111°07′3.14″	133°∠76°	198	N 37°35′33.50″	E 110°53′10.92″	185°∠86°
84	N 39°03′20.39″	E 111°07′3.14″	221°∠84°	199	N 37°35′33.50″	E 110°53′10.92″	190°∠87°
85	N 39°03′20.73″	E 111°07′2.48″	212°∠79°	200	N 37°35′33.50″	E 110°53′10.92″	290°∠82°
86	N 39°03′20.73″	E 111°07′2.48″	139°∠87°	201	N 37°35′30.47″	E 110°53′0.46″	294°∠82°
87	N 39°03′20.32″	E 111°07′2.01″	222°∠79°	202	N 37°35′30.47″	E 110°53′0.46″	193°∠82°
88	N 39°03′20.32″	E 111°07′2.01″	143°∠84°	203	N 37°35′30.47″	E 110°53′0.46″	190°∠89°
89	N 39°03′19.69″	E 111°07′1.11″	56°∠81°	204	N 37°35′30.47″	E 110°53′0.46″	290°∠84°
90	N 39°03′19.69″	E 111°07′1.11″	127°∠79°	205	N 37°35′30.47″	E 110°53′0.46″	191°∠87°
91	N 39°03′18.83″	E 111°07′0.56″	53°∠74°	206	N 37°35′30.47″	E 110°53′0.46″	275°∠88°
92	N 39°03′18.83″	E 111°07′0.56″	131°∠87°	207	N 37°35′26.86″	E 110°52′50.82″	190°∠89°
93	N 39°03′21.35″	E 111°07′4.26″	42°∠81°	208	N 37°35′26.86″	E 110°52′50.82″	285°∠87°
94	N 39°03′21.35″	E 111°07′4.26″	320°∠83°	209	N 37°35′26.86″	E 110°52′50.82″	198°∠88°
95	N 39°03′21.35″	E 111°07′4.26″	45°∠71°	210	N 37°35′26.86″	E 110°52′50.82″	281°∠71°
96	N 39°03′21.35″	E 111°07′4.26″	324°∠86°	211	N 37°35′25.39″	E 110°52′48.02″	80°∠73°
97	N 39°03′35.25″	E 111°06′19.22″	182°∠67°	212	N 37°35′25.39″	E 110°52′48.02″	193°∠88°
98	N 39°03′35.25″	E 111°06′19.22″	277°∠90°	213	N 37°35′25.28″	E 110°52′41.15″	291°∠84°
99	N 39°03′35.25″	E 111°06′19.22″	359°∠71°	214	N 37°35′25.28″	E 110°52′41.15″	189°∠87°
100	N 39°03′35.25″	E 111°06′19.22″	272°∠77°	215	N 37°35′22.00″	E 110°52′35.99″	192°∠89°
101	N 39°03′35.25″	E 111°06′19.22″	272°∠57°	216	N 37°35′22.00″	E 110°52′35.99″	289°∠68°
102	N 39°03′35.25″	E 111°06′19.22″	183°∠88°	217	N 37°35′20.02″	E 110°52′30.63″	189°∠89°
103	N 39°03′34.15″	E 111°06′20.15″	267°∠76°	218	N 37°35′20.02″	E 110°52′30.63″	282°∠79°
104	N 39°03′34.15″	E 111°06′20.15″	183°∠74°	219	N 37°35′20.02″	E 110°52′30.63″	186°∠85°
105	N 39°03′30.99″	E 111°06′19.14″	355°∠83°	220	N 37°35′20.02″	E 110°52′30.63″	82°∠66°
106	N 39°03′30.99″	E 111°06′19.14″	87°∠81°	221	N 37°35′19.18″	E 110°52′26.46″	82°∠81°
107	N 39°03′30.99″	E 111°06′19.14″	230°∠82°	222	N 37°35′19.18″	E 110°52′26.46″	193°∠86°
108	N 39°03′30.99″	E 111°06′19.14″	325°∠46°	223	N 37°35′13.14″	E 110°51′57.89″	179°∠86°
109	N 39°03′30.26″	E 111°06′18.75″	359°∠82°	224	N 37°35′13.14″	E 110°51′57.89″	270°∠79°
110	N 39°03′30.26″	E 111°06′18.75″	273°∠79°	225	N 37°35′13.14″	E 110°51′57.89″	175°∠84°
111	N 39°03′30.15″	E 111°06′13.40″	224°∠75°	226	N 37°35′13.14″	E 110°51′57.89″	283°∠84°
112	N 39°03′30.15″	E 111°06′13.40″	82°∠87°	227	N 37°35′13.70″	E 110°51′57.07″	184°∠82°
113	N 39°03′31.40″	E 111°06′11.67″	337°∠66°	228	N 37°35′13.70″	E 110°51′57.07″	272°∠85°
114	N 39°03′31.40″	E 111°06′11.67″	54°∠75°	229	N 37°35′14.15″	E 110°51′51.27″	177°∠72°
115	N 39°03′31.40″	E 111°06′11.67″	209°∠73°	230	N 37°35′14.15″	E 110°51′51.27″	280°∠69°

), showing that sub-vertical joints are developed in different geological points. Additionally, based on the strike rose diagram of field joints (Figure 2), the dominant strikes of joints are NNE (10°~25°), NEE (50°~80°), near-EW (85°~110°), NW (130°~150°), and NNW (160°~175°) trending. Most of the triangular rock blocks sandwiched between the joint surfaces have fallen off.

4.2. Paleo Tectonic Stress Characteristics

Based on the development characteristics of conjugate joints and the tectonic evolution history in the study area, four combinations of different types and periods are classified using the stereographic projection method [36] (Figure 3): the type I combination has a set of conjugated joints with NWW (50°~80°) and NNW (150°~175°) strikes, with joint surfaces approximately perpendicular to the geological strata formed in the Yanshanian period and the maximum stress direction in the near-NW orientation (Figure 3a,b).

The conjugated joint strikes of the type II combination are in the NNE direction (0°~40°) and the near-EW direction (80°~110°), corresponding to the Himalayan period, with the maximum stress direction being in the near-NE orientation (Figure 3c,d).

It is obvious that the joints in the study area are mainly formed in the Yanshanian period and the Himalayan period. Additionally, the paleostress magnitude is calculated using the conjugate shear angle estimation method, which shows the maximum horizontal principal stress in the Yanshanian period is 160 MPa, and the minimum principal stress is 10 MPa. The maximum horizontal principal stress in the Himalayan period is 100 MPa, and the minimum principal stress is 20 MPa.

4.3. Paleogeological Models

The compression with NW orientation in the Yanshanian period causes the coal reservoir to produce NW oriented folds. After the compression with NE orientation in the Himalayan period, which is nearly perpendicular to the Yanshanian period, the NW oriented folds are superimposed on the NE oriented folds, forming a superimposed fold structure. The 8 + 9 # coal reservoir in the study area forms a large number of synclines, anticlines, and saddle-shaped structures. For the restoration of the ancient geological model in the study area, the core is the products of these superimposed structures. Firstly, taking the Benxi formation 8 + 9 # coal reservoir in the Linxing area as the research object, the structural traces of the current coal seam floor contour lines are analyzed, and structures including anticlines, synclines, and folds are categorized. The structural superposition method [20] is applied to flatten the superimposed products of anticlines, synclines, and saddle structures and obtain the contour map of the coal seam floor in paleo periods (Figure 4).

Based on the contour map of coal thickness and ancient coal seam floor, an isotropic idealized geological model is established using triangular meshes in ANSYS software (18.0 version). In addition, the stress of the reservoir is greatly affected by the roof and floor rock layers, so the geological model is constructed as a roof–coal–floor combination type (Figure 5), and a cube is considered to surround the roof–coal–floor combination model to simulate the stress condition of the surrounding rock.

4.4. Paleostress Field

Based on the tectonic stress direction and magnitude, overburden stress, Young’s modulus, Poisson’s ratio, and density of geological model (

Table 1. Mechanical parameters for numerical simulation of the paleostress field in the Linxing area.

Geologic Bodies	Poisson’s Ratio	Young’s Modulus (GPa)	Density (g/cm3)
Roof	0.22	21.33	2.730
Coal seam	0.36	6.2	1.480
Floor	0.21	21.55	2.750
Others	0.23	20	1.655

), geomechanical models with boundary conditions for different paleogeological periods are applied, and the finite element method and ANSYS software (18.0 version) are used to obtain the stress distribution of the 8 + 9 # coal reservoir in different geological periods.

4.4.1. Yanshanian Period

The stress distribution of the 8 + 9 # coal reservoir in the Yanshanian period is shown in Figure 6, where the values of maximum principal stress are between −55 and −82 MPa, with an average of −67 MPa, which is indicative of compression (Figure 6a). Excluding the influence of model boundaries on the results, the overall maximum principal stress gradually increases from north to south, and high-value zones are mainly distributed in the southern syncline (mostly −76 to −82 MPa) and central magmatic rock uplift (mainly between −73 and −80 MPa). The minimum principal stress values are distributed between −8 and −35 MPa, with an average of −20 MPa, which is also indicative of compression (Figure 6b). The value of minimum principal stress gradually decreases from east to west, and the central magmatic rock uplift has a low minimum principal stress value, with values below −10 MPa.

Figure 6c indicates that shear stress values in the Yanshanian period are generally between 11 MPa and 27.5 MPa. The stress value in the northern zone of the study area gradually increases from northeast to southwest, exhibiting a strip-shaped distribution; the central magmatic rock uplift and southwest zones show high shear stress values, which means shear fractures are prone to occur in those zones.

4.4.2. Himalayan Period

Figure 7 shows the stress distribution of the 8 + 9 # coal reservoir in the Himalayan period, in which the maximum and minimum principal stress values are from −34 to −70 MPa and from −2 to −27 MPa, respectively, both in the form of compressive stress. The overall trends of high in south zones and low in the north zones of maximum and minimum principal stress are displayed in the study area (Figure 7a,b), and local low value zones of maximum and minimum principal stress are distributed in the central magmatic rock uplift.

As is shown in Figure 7c, the values of shear stress range from −12 to −18 MPa, with an average of −15 MPa. The value of shear stress gradually increases from southeast to northwest zones, and high shear stress value occurs in the central magmatic rock uplift, indicating the influence of magmatic rock uplift on in situ stress distribution and fracture development.

4.5. Evaluation of Favorable Fracture Area

Based on the numerical simulation results, combined with the tensile strength, cohesion, and internal friction angle calculated from logging, the tensile and shear fracture coefficients are calculated. The results are shown in

Table 2. Fracture coefficients in two paleogeological periods.

Nodes	CT1	CT2	CS1	CS2	F	Nodes	CT1	CT2	CS1	CS2	F
1149	−0.06	−0.60	−0.06	0.04	0.55	11,651	−0.05	−0.51	−0.09	−0.12	0.45
1409	−0.07	−0.83	−0.05	0.02	0.54	11,916	−0.04	−0.53	−0.09	−0.30	0.41
1565	−0.07	−0.95	−0.06	0.04	0.53	12,124	−0.04	−0.44	−0.09	−0.30	0.41
1669	−0.07	−1.01	−0.06	0.05	0.53	12,332	−0.04	−0.33	−0.09	−0.14	0.48
1773	−0.08	−1.09	−0.04	0.19	0.59	12,701	−0.02	−0.32	−0.09	−0.32	0.41
1877	−0.07	−1.14	−0.06	0.06	0.52	12,805	−0.03	−0.42	−0.09	−0.40	0.39
1981	−0.07	−1.19	−0.06	−0.01	0.49	12,909	−0.02	−0.35	−0.08	0.00	0.54
2033	−0.08	−1.24	−0.05	0.15	0.56	13,387	−0.02	−0.33	−0.10	−0.38	0.39
2298	−0.06	−0.69	−0.06	0.01	0.53	13,491	0.01	−0.21	−0.06	0.15	0.65
2402	−0.07	−0.76	−0.05	0.14	0.59	14,021	−0.01	−0.30	−0.09	−0.48	0.37
3291	−0.06	−0.62	−0.06	0.01	0.54	14,125	0.07	0.07	0.02	0.37	0.92
3499	−0.07	−0.80	−0.07	0.00	0.52	14,230	0.01	0.00	−0.12	0.28	0.59
3707	−0.07	−0.95	−0.07	−0.09	0.47	14,707	0.03	−0.08	−0.04	0.05	0.66
3915	−0.07	−1.08	−0.07	−0.14	0.44	14,811	0.05	0.12	−0.07	0.07	0.65
4123	−0.07	−1.15	−0.06	0.11	0.53	15,237	−0.01	−0.22	−0.10	−0.63	0.32
5225	−0.05	−0.52	−0.07	−0.05	0.52	15,341	0.02	−0.09	−0.07	−0.33	0.48
5433	−0.06	−0.67	−0.07	−0.06	0.50	15,445	0.05	0.11	−0.08	−0.26	0.51
5641	−0.07	−0.84	−0.08	−0.09	0.47	15,923	−0.01	−0.24	−0.09	−0.75	0.29
5849	−0.07	−0.98	−0.08	−0.08	0.45	16,027	0.00	−0.15	−0.07	−0.56	0.40
6057	−0.07	−1.10	−0.08	−0.02	0.46	16,132	0.00	−0.18	−0.10	−0.48	0.37
7367	−0.05	−0.53	−0.08	−0.12	0.48	16,505	−0.01	−0.25	−0.10	−0.76	0.27
7575	−0.06	−0.71	−0.08	−0.13	0.46	16,609	−0.01	−0.30	−0.10	−0.82	0.25
7783	−0.06	−0.86	−0.08	−0.18	0.42	16,713	−0.01	−0.34	−0.10	−0.85	0.24
7991	−0.07	−0.93	−0.08	−0.25	0.39	17,295	−0.02	−0.39	−0.10	−0.89	0.22
8199	−0.07	−1.09	−0.09	−0.16	0.40	17,400	−0.02	−0.41	−0.10	−0.92	0.21
9301	−0.04	−0.46	−0.08	−0.21	0.45	17,959	−0.01	−0.34	−0.10	−0.87	0.23
9509	−0.05	−0.59	−0.08	−0.11	0.47	18,115	−0.02	−0.47	−0.09	−1.00	0.18
9717	−0.06	−0.80	−0.09	−0.25	0.40	18,271	−0.03	−0.62	−0.09	−1.05	0.16
9925	−0.05	−0.70	−0.05	0.15	0.60	18,428	−0.04	−0.67	−0.09	−0.84	0.23
10,133	−0.07	−0.97	−0.10	−0.17	0.39	18,848	−0.03	−0.46	−0.11	−0.88	0.20
11,443	−0.04	−0.49	−0.09	−0.24	0.43	18,952	−0.03	−0.51	−0.11	−0.91	0.19

Notes: C_T1 and C_T2 represent tensile fracture coefficients of Yanshanian and Himalayan periods, respectively; C_S1 and C_S2 represent shear fracture coefficients of Yanshanian and Himalayan periods, respectively; and F represents the comprehensive evaluation factor.

.

It is obvious that tensile fracture coefficients are far lower than shear fracture coefficients in the coal reservoir, which indicates that shear behaviors are more likely to occur than tensile behaviors, and that shear fractures are more developed than tensile fractures. In addition, the shear fracture coefficients of the Himalayan period are higher than those of the Yanshanian period, reflecting that fractures are prone to developing during the Himalayan period. Thus the comprehensive evaluation factor F was established and calculated to show the development degree of fracture under multiple structural movements.

$$F=a{N}_{T1}+{bN}_{T2}+c{N}_{S1}+{dN}_{S2}$$

(12)

$$N=\frac{C_i-C_{min}}{C_{max}-C_{min}}$$

(13)

where F represents the comprehensive evaluation factor; N represents the normalization value of the fracture coefficient; C_i, C_max, and C_min, represent fracture coefficients of node i, the maximum fracture coefficient, and the minimum fracture coefficient, respectively; N_T1 and N_T2 represent normalization values of the tensile fracture coefficients of the Yanshanian and Himalayan periods, respectively; N_S1 and N_S2 represent normalization values of the shear fracture coefficients of the Yanshanian and Himalayan periods, respectively; and a, b, c, and d represent weights of normalization values of fracture coefficients, with values of 0.1, 0.1, 0.3, and 0.5, respectively.

Additionally, the greater the value of F, the greater the development degree of fracture. According to the comprehensive evaluation results, the study area is divided into five classes in the degree of fracture development: Class I (F > 0.50), Class II (0.45 < F ≤ 0.50), Class III 0.40 < F ≤ 0.45), Class IV (0.35 < F ≤ 0.40), and Class V (F ≤ 0.35), thereby performing a quantitative evaluation of the favorable fracture area.

As shown in Figure 8, the study area is divided into 24 fracture development areas, of which Class I areas are mainly distributed in the northwest and surround the magmatic rock uplift. Fractures in Class I areas are well developed and can be formed under two periods of tectonic movements. Class II areas are mainly distributed in the central and western regions of the research area, and fractures are relatively well developed. Class III areas are mainly distributed in the central, western, and southwestern parts of the study area with moderately developed fractures. Class IV and Class V areas have less-developed fractures, and fractures can only be formed through one tectonic movement period in some areas. Generally, fractures in the study area are controlled by tectonic stress, magmatic rock uplift, and buried depth, and favorable fracture areas are located surrounding magmatic rock uplift. Additionally, the division of reservoir fractures area can help to predict the development degree of reservoir fractures and provide a basis for coalbed methane extraction.

5. Conclusions

Areas with high shear stress values are mainly distributed in the central magmatic rock uplift, indicating the influence of magmatic rock uplift on in situ stress distribution and fracture development. Based on the comprehensive evaluation factors of fractures, the reservoir is divided into five classes and 24 favorable fracture areas, and the following conclusions can be made:

(1)

The 8+9 # coal reservoir in the Linxing area has mainly undergone two stages of tectonic movements, which are the compression in the Yanshanian period in the NW direction and the compression in the Himalayan period in the NE direction. The maximum horizontal principal stress during the Yanshanian period is 160 MPa, and the minimum principal stress is 10 MPa. The maximum horizontal principal stress during the Himalayan period is 110 MPa, and the minimum principal stress is 20 MPa.

(2)

The degree of fracture development in deep coal reservoirs in the research area is directly influenced by the paleostress field, with the main fracturing periods being the Yanshanian and Himalayan periods. Based on the distribution of the paleostress field obtained from numerical simulation, the Mohr–Coulomb fracture criterion and Griffith fracture criterion are used to predict shear and tension fractures. It is found that the fracture threshold of shear fracture is smaller than that of tension fracture, and shear fractures are formed earlier than tensile fractures.

(3)

Based on the comprehensive evaluation factors of fractures, the 8 + 9 # coal reservoir is divided into 24 favorable fracture areas from Class V to Class I. Fractures in Class I areas and Class II areas are relatively well developed and were formed under two periods of tectonic movements. Additionally, there are nine favorable zones in Class I and Class II, mainly distributed in the northwest of the study area and the magmatic rock uplift area.