Finite Element Simulation of Multi-Scale Bedding Fractures in Tight Sandstone Oil Reservoir

: Multi-scale bedding fractures, i.e., km-scale regional bedding fractures and cm-scale lamina-induced fractures, have been the focus of unconventional oil and gas exploration and play an important role in resource exploration and drilling practice for tight oil and gas. It is challenging to conduct numerical simulations of bedding fractures due to the strong heterogeneity without a proper mechanical criterion to predict failure behaviors. This research modiﬁed the Tien–Kuo (T–K) criterion by using four critical parameters (i.e., the maximum principal stress ( σ 1 ), minimum principal stress ( σ 3 ), lamina angle ( θ ), and lamina friction coe ﬃ cient ( µ lamina )). The modiﬁed criterion was compared to other bedding failure criteria to make a rational ﬁnite element simulation constrained by the four variables. This work conducted triaxial compression tests of 18 column samples with di ﬀ erent lamina angles to verify the modiﬁed rock failure criterion, which contributes to the simulation work on the multi-scale bedding fractures in the statics module of the ANSYS workbench. The cm-scale laminated rock samples and the km-scale Yanchang Formation in the Ordos Basin were included in the multi-scale geo-models. The simulated results indicate that stress is prone to concentrate on lamina when the lamina angle is in an e ﬀ ective range. The low-angle lamina always induces fractures in an open state with bigger failure apertures, while the medium-angle lamina tends to induce fractures in a shear sliding trend. In addition, the regional bedding fractures of the Yanchang Formation in the Himalayan tectonic period tend to propagate under the conditions of lower maximum principal stress, higher minimum principal stress, and larger stratigraphic dip. distribution of tight sandstone reservoirs of Yanchang Formation in Ordos Basin during the Himalayan episode.


Introduction
Multi-scale bedding fractures, including km-scale regional bedding fractures and cm-scale lamina-induced fractures, are caused by lamina dissolution or induced by regional tectonic stress [1][2][3][4][5][6][7][8][9][10]. The formation of sedimentary lamina and bedding underground exerts a strong control on the following fracture propagation [11][12][13], so the term lamina-induced fractures is used to indicate the cm-scale fractures forming along the core lamina under the influence of external forces and internal rock mechanical properties [14]. Moreover, the lamina-induced fractures also correspond to the km-scale regional bedding fractures, which are induced by the formation beddings. These fractures have been the focus of conventional and unconventional oil and gas reservoir characterization including shale, tight oil [14][15][16][17][18][19][20][21][22][23], and Carboniferous rocks [24]. The fractures have also played an important role in methane gas emissions from coal seams [25][26][27]. It has been found that the open state of bedding fractures in complicated tectonic zones could play an essential role in oil and gas diffusion, emission, migration, and accumulation [14,17,[28][29][30][31][32][33]. In addition, the lamina or bedding in the tight reservoirs coupled with the hydraulic fractures could induce a more complicated in-situ fracture network [34]. However, the quantitative simulation work of bedding fractures, despite the qualitative description, as mentioned above, is still lacking due to the complicated stress distribution caused by the lamina heterogeneity. It is challenging to distinguish the lamina or bedding from surrounding rocks and predict the failure behaviors under the effects of heterogeneity.
Multiple methods and softwares have been used to conduct numerical simulations and make predictions on the energy resources, e.g., artificial neural networks, ant tracking algorithms, petrophysical logging, and microseisms [35][36][37][38]. The ant tracking algorithm can be applied to detect small faults, but it is difficult to extract detailed formation of the fractures due to the extremely low resolution and low coherence of data [39]. Even though the microseism is often used for artificial fracture detection and modeling in the hydraulic fracturing process of the horizontal wells, its roles are confined to the oil-gas exploration and development period [40]. In addition, the accuracy of logging interpretation strongly depends on the data amount, and the resolution of seismic data is too low to conduct a precise fracture simulation. Therefore, it is necessary to find out an efficient method to make an accurate prediction for the multi-scale fractures without the need for a great amount of data.
In this study, a finite element simulation of bedding fractures was conducted based on the dynamic propagation condition for the tectonic bedding fractures. The finite element simulation is widely used to predict the present geological stress and fracture index, because it can provide a platform for researchers to focus on the geological or mechanical model, with the adjustability of the rock failure criterion and boundary stress conditions [41][42][43]. In addition, this study focused on the modified failure criterion of tectonic bedding fractures in the dynamic simulation environment.
The most challenging work for bedding fracture simulation is to build rational failure criteria in the finite element simulation environment considering four factors, i.e., maximum stress, minimum stress, lamina angle (or the km-scale regional stratigraphic dip), and differences between the lamina and surrounding rock in different lamina lithofacies. Although some failure criteria have been proposed for the bedding fractures based on the mechanical tests of laminated rock, as shown in Table 1 [6,[44][45][46][47], a coefficient should be additionally proposed to indicate the different lamina lithofacies and the difference between the lamina and the surrounding rock in the laminated cores or the rock formations with bedding fractures.
The µ lamina values were tested based on the triaxial compression tests of different lamina lithofacies.
(2) A finite element simulation was conducted based on the modified T-K criterion to study the stress distribution in the laminated rock model composed of the independent lamina and the surrounding rock bodies. (3) The regional bedding fractures distribution of the Upper Triassic Yanchang Formation in the Ordos Basin was clarified based on the regional stratigraphic dip distribution, as well as the regional laminated rock lithofacies distribution proposed in this work, and the regional stress field, which was simulated in our prior simulation work [14].

Materials
The six laminated samples studied in this work were collected from the Chang 8 to Chang 6 member, the principal hydrocarbon-accumulating intervals in the Upper Triassic Yanchang Formation of the southern Ordos Basin. These samples could represent the different petrophysical facies in the delta-lacustrine environment. Core-1 and Core-2 are the tight siltstone-and sandstone-deposited distributary channel facies in tractive current hydrodynamic conditions. Core-1 represents the tight channel siltstone with the inconspicuous lamina, i.e., micro-cross-bedding lamina. Core-2 represents the tight channel sandstone with the macro-cross-bedding lamina. Core-3 represents the tight argillaceous siltstone with the horizontal lamina. Core-4 and Core-5 represent the tight sheet siltstone with the obvious lamina, i.e., wavy or lenticular lamina. Core-6 represents the tight siltstone deposited in the lacustrine turbidite fan facies with slumping-induced deformation lamina. Generally, the six laminated cores are typical representatives for all tight sandstone oil reservoirs, especially for the tight sandstone oil reservoir of the Ordos Basin in China.

Triaxial Compression Tests
To indicate the heterogeneity of laminated rocks and lamina lithofacies in the geo-models, we tested 18 column samples with different lamina angles from the six different laminated cores in triaxial compression experiments. The triaxial compression tests were performed with an RTR apparatus under a 20 MPa confining pressure. The 20 MPa approximates the underground pressure at a 2000 m depth, so it can be used to simulate real geological conditions. Sample treatment was conducted in strict accordance with the International Society of Rock Mechanics (ISRM) rock triaxial test requirements and the China National Standard (GB/T 50266-99).

Finite Element Simulation
Finite element simulation was performed in the statics module of the ANSYS workbench. In this method, the geological model is divided into different specific elements linked to nodes. In addition, the approximate value of the node displacement could be calculated using the element functions in the equilibrium state. Based on the computed result, we can obtain the stress and strain of these elements [48][49][50][51][52]. For the geological model building work, it is generally assumed that the stress and strain of the shallow crust rocks are linearly dependent during the period of rock elastic deformation [53,54]. In this study, the principle of the finite element approach is composed of the following three relationships: (1) The strain-displacement relationship; (2) strain-stress relationship; and (3) stress-external force relationship. The principle and its detailed explanations can be found in our previous study [17].
In this work, to calculate the boundary stress conditions for the finite element simulation, we modified the Newberry formula to make it more suitable for the tensional and compressive stress fields. Although the rock spatial stress computed by the Newberry formula is appropriate for the tight reservoir simulation [55], the equation is only defined in compressive stress fields and leaves the tensional stress field out of consideration [56][57][58]. Hence, the positive sign of the compressive stress field in the Newberry formula is changed to a negative sign for the tensional stress field. In addition, this work constrained the calculative process of the boundary stress state, based on the assumption that the 3D geo-model can be regarded as a point mass with the corresponding well burial history area. In this way, we can calculate the boundary stress state by combining the regional tectonic stress, tectonic direction, and rock density. The detailed modification of the Newberry formula and the calculative process of the boundary stress state can also be found in our previous study [17].

Stress Balance along the Lamina
A finite element simulation of lamina-induced fractures is conducted in the statics module of the ANSYS workbench based on two geo-bodies, i.e., the lamina and surrounding rock. The contact relationship was settled to the frictional model characterized by a friction coefficient representing the difference between the two separated geo-bodies. Here, this study defined the coefficient as the lamina friction coefficient (µ lamina ) in the laminated rock, where µ lamina could be used to indicate different lamina lithofacies and modify the T-K failure criterion in the finite element simulation. The laminated rock model is shown in Figure 1, which was constructed with the laminated body and the surrounding rock body. Two parameters were introduced into the geo-model to characterize the laminated rock, i.e., the lamina angle (θ) and µ lamina . It should be emphasized that µ lamina of the laminated rock differs from the internal friction coefficient of intact rocks. The µ lamina is not only a physical factor representing the physical sliding behavior in simulation but also a geological factor indicating different lamina lithofacies. Therefore, the range of µ lamina may vary from the traditional internal friction coefficient. Its detailed values for different laminated rock types were tested using the triaxial compression technologies in Section 2.2.
Energies 2020, 13, x FOR PEER REVIEW 4 of 20 this way, we can calculate the boundary stress state by combining the regional tectonic stress, tectonic direction, and rock density. The detailed modification of the Newberry formula and the calculative process of the boundary stress state can also be found in our previous study [17].

Stress Balance along the Lamina
A finite element simulation of lamina-induced fractures is conducted in the statics module of the ANSYS workbench based on two geo-bodies, i.e., the lamina and surrounding rock. The contact relationship was settled to the frictional model characterized by a friction coefficient representing the difference between the two separated geo-bodies. Here, this study defined the coefficient as the lamina friction coefficient (μlamina) in the laminated rock, where μlamina could be used to indicate different lamina lithofacies and modify the T-K failure criterion in the finite element simulation. The laminated rock model is shown in Figure 1, which was constructed with the laminated body and the surrounding rock body. Two parameters were introduced into the geo-model to characterize the laminated rock, i.e., the lamina angle (θ) and μlamina. It should be emphasized that μlamina of the laminated rock differs from the internal friction coefficient of intact rocks. The μlamina is not only a physical factor representing the physical sliding behavior in simulation but also a geological factor indicating different lamina lithofacies. Therefore, the range of μlamina may vary from the traditional internal friction coefficient. Its detailed values for different laminated rock types were tested using the triaxial compression technologies in Section 2.2. Another critical factor for the failure criterion is the lamina angle (θ), indicating the intersection angle between the axial stress direction and the lamina dip direction, which is the complementary angle to the β angle (the acute angle between the direction of maximum principal stress and the discontinuity) in previous works [6,47]. A stress equilibrium equation was built based on the decomposition of stress along the lamina surface (Equation (1), Figure 1). In this equation, θ indicates the lamina angle, σ1 indicates the axial stress, σ3 indicates the confining pressure, and Clamina indicates the lamina cohesion, which is different from the cohesion of homogeneous rock. Another critical factor for the failure criterion is the lamina angle (θ), indicating the intersection angle between the axial stress direction and the lamina dip direction, which is the complementary angle to the β angle (the acute angle between the direction of maximum principal stress and the discontinuity) in previous works [6,47]. A stress equilibrium equation was built based on the decomposition of stress along the lamina surface (Equation (1), Figure 1). In this equation, θ indicates the lamina angle, σ 1 indicates the axial stress, σ 3 indicates the confining pressure, and C lamina indicates the lamina cohesion, which is different from the cohesion of homogeneous rock.
The Mohr-Coulomb rock failure criterion (Equation (2)) is usually used to explain the propagation condition for shear fractures. Here, an assumption was made that the propagation condition for bedding fractures meets the Mohr-Coulomb rock failure criterion. In other words, the lamina-induced fractures could be regarded as specially compressed fractures along the lamina surface in the intact rock. Then, Equation (3) was achieved when we solved the simultaneous Equations (1) and (2). It describes the relationship among the internal friction angle (ϕ) of intact rock, the lamina friction coefficient (µ lamina ), and the lamina angle (θ) of laminated rock. Equation (4) describes the relationship among the cohesion (S 0 ), the internal friction angle (ϕ) of intact rock, and the lamina angle (θ) of laminated rock. The final relationship (Equation (5)) was built among the cohesion (S 0 ), the internal friction angle (ϕ) of intact rock, the lamina angle (θ), and the lamina friction coefficient (µ lamina ) of laminated rock, when Equation (3) was put into Equation (4). Thus, an important function (Equation (5)) was used to modify the T-K criterion in the following section. Furthermore, the detailed derivation and physical interpretations are shown in Appendix A.

The Failure Criterion of Bedding Fractures
Tien and Kuo (2001) proposed a common failure criterion (Equation (6)) for the intact bedded rocks, which demonstrates the characteristics of strength anisotropy revealed by the laboratory experiments, where k and n indicate the elastic constants of laminated rocks [6]. In addition, the T-K criterion is based on the nonlinear Hoek-Brown criterion for the homogeneous rock with no lamina.
Zhou et al. (2017) proposed a more efficient modified T-K criterion (Equation (7)) to replace the nonlinear Hoek-Brown criterion with the linear Mohr-Coulomb criterion, where ϕ indicates the cohesion of surrounding rock beyond the lamina and S 0 indicates the cohesion of laminated rocks [47]. The main problem for Zhou et al.'s modified criterion is that the cohesion of laminated rocks (S 0 ) was a constant in Equation (7), but it actually varied with the lamina angle (θ) and the lamina friction coefficient (µ lamina ) when the compressed fractures primarily propagated along the lamina surface in the laminated rocks. Thus, θ and µ lamina could be used together, as two modified factors of S 0 in the T-K criterion, to predict the simulated failure behaviors in cm-scale laminated cores to the km-scale bedding formation.
Here, Equation (5) was put into Equation (7) to modify the T-K criterion. Thus, a modified T-K criterion is proposed by Equation (8), where ϕ indicates the surrounding rock properties and the other parameters indicate the lamina properties. The modified criterion is composed of four critical variables, i.e., the maximum principal stress (σ 1 ), minimum principal stress (σ 3 ), lamina angle (θ), and lamina friction coefficient (µ lamina ), which are advantageous and could also be used in the ANSYS finite element simulation.

Index of Bedding Fractures
Li et al. (2018) simulated the distribution of the structural fracture of the Upper Triassic Yanchang Formation in the Ordos Basin by the ANSYS software with a simplified fracture index set. The simulation was proposed based on the relationship between the maximum principal stress and the failure strength, and represents the total fracture possibility based on the distance between the envelope line and the stress state point (Equation (9)) [14]. In the bedding fracture simulation, σ 1 (the maximum normal stress satisfied the critical rupture condition) is replaced by the failure criterion for the bedding fractures (i.e., Equation (8)) to indicate the index of bedding fractures, as shown in Equation (10). (10) where T 0 is the tensile strength, T 0-lamina is the tensile strength of laminated sandstone; µ is the coefficient of internal friction; f LF is the fracture index of lamina induced fractures.

Comparisons of Different Failure Criteria For Bedding Fractures
To verify the modified T-K criterion, the testing data were compared to the calculated data using three criteria, i.e., the modified T-K criterion in this study, Zhou et al.'s modified T-K criterion, and the Jaeger criterion (Table 1). Figure Table 1, and the detailed values of each criterion's parameters are listed in Table 2. Detailed peak strength values are listed in Table 3       As shown in Figure 2, there was a much lower peak strength for those samples with the lamina angle of 30 • . Notably, black square points in Figure 3  Generally, the modified T-K failure criterion has almost the same trend as the T-K criterion and even agrees better with the testing data in , as shown in Figure 3, especially for the medium-low-angle range (nearly 20 • -50 • ) [59]. Moreover, the four variables in this criterion are in accordance with the ANSYS parameter environment and could provide us with a better failure criterion in the finite element simulation. Generally, the modified T-K failure criterion has almost the same trend as the T-K criterion and even agrees better with the testing data in , as shown in Figure 3, especially for the medium-low-angle range (nearly 20°-50°) [59]. Moreover, the four variables in this criterion are in accordance with the ANSYS parameter environment and could provide us with a better failure criterion in the finite element simulation.

Lamina Friction Coefficients of Different Lamina Lithofacies
It is necessary to use the lamina friction coefficient (μlamina) to indicate the regional heterogeneity of lamina lithofacies in the geo-simulation work. In order to test the lamina friction coefficient (μlamina) of the different lamina rock types, we tested 18 column samples with different lamina angles from six different laminated cores in triaxial compression experiments. Triaxial compression stress-strain curves are shown in

Lamina Friction Coefficients of Different Lamina Lithofacies
It is necessary to use the lamina friction coefficient (µ lamina ) to indicate the regional heterogeneity of lamina lithofacies in the geo-simulation work. In order to test the lamina friction coefficient (µ lamina ) of the different lamina rock types, we tested 18 column samples with different lamina angles from six different laminated cores in triaxial compression experiments. Triaxial compression stress-strain curves are shown in  As indicated in Figure 5, the following results could be drawn. (1) Different lamina friction coefficients (μlamina) indicate different laminated core types. The cores more strongly characterized by the obvious lamina, such as Core-4, Core-5, and Core-6, are typically correlated with a lower μlamina.
(2) The fitting data agreed well with the testing data, except Core-1, which is characterized by the inconspicuous and micro-cross-bedding lamina. Here, the strength values of Core-1 are nearly located in the same range due to the weak influence of these inconspicuous and micro-cross-bedding lamina within Core-1. In other words, the smaller the difference between the lamina and surrounding rock, the weaker the influence of the lamina in different angle ranges. Furthermore, the values of μlamina of Core-4 and Core-5 are around 0.5, which are the laminated siltstone deposited by traction currents in the delta front sedimentary environment. The value of μlamina of Core-6 is around 0.3, which is the laminated siltstone sedimented in the lacustrine turbidite environment. The value of μlamina of Core-3 is around 0.9, which is the laminated mudstone sedimented in a still water environment. From the perspective of the regional finite element simulation, the μlamina, as a critical variable, should be assigned differently according to the laminated rock types deposited in the different sedimentary environments during Yanchang Formation in the Ordos Basin. As indicated in Figure 5, the following results could be drawn. (1) Different lamina friction coefficients (µ lamina ) indicate different laminated core types. The cores more strongly characterized by the obvious lamina, such as Core-4, Core-5, and Core-6, are typically correlated with a lower µ lamina .
(2) The fitting data agreed well with the testing data, except Core-1, which is characterized by the inconspicuous and micro-cross-bedding lamina. Here, the strength values of Core-1 are nearly located in the same range due to the weak influence of these inconspicuous and micro-cross-bedding lamina within Core-1. In other words, the smaller the difference between the lamina and surrounding rock, the weaker the influence of the lamina in different angle ranges. Furthermore, the values of µ lamina of Core-4 and Core-5 are around 0.5, which are the laminated siltstone deposited by traction currents in the delta front sedimentary environment. The value of µ lamina of Core-6 is around 0.3, which is the laminated siltstone sedimented in the lacustrine turbidite environment. The value of µ lamina of Core-3 is around 0.9, which is the laminated mudstone sedimented in a still water environment. From the perspective of the regional finite element simulation, the µ lamina , as a critical variable, should be assigned differently according to the laminated rock types deposited in the different sedimentary environments during Yanchang Formation in the Ordos Basin.

Stress Distribution in the Laminated Rock
The laminated core model (cm-scale) in the finite element simulation is composed of two bodies, as shown in Figure 6, where the contact part between two bodies is settled as frictional with a sliding friction coefficient in accordance with the lamina friction coefficient in the failure criteria of this study. The failure behavior of the laminated rock in the special lamina angle range is determined by the stress distribution constrained by the constant value of μlamina. Therefore, it is necessary to study the stress distribution of the laminated rock using the finite element simulation technology, especially for the different lamina angle ranges, i.e., the low, medium, and high lamina angle ranges.
Researchers have tried to simulate the compression behavior of tight intact rock and the laminated rock using finite element simulation technology [60,61]. To better understand the stress distribution during the compression process, the failure distribution around the compression samples was simulated by using the modified criterion in the ANSYS environment. The simulated samples were divided into two parts: The surrounding rock and the lamina body, with the lamina surface acting as the sliding surface. Additionally, the mechanical parameters of the surrounding rock used in the simulation were different from those of the laminated rock. Here, the surrounding rock was set as siltstone, while the laminated rock was claystone in the ANSYS software. In addition, the contact relationship, i.e., the lamina sliding surface, was set as frictional in the ANSYS software with the fractional coefficient value of 0.2. Another critical parameter, the lamina angle, was designed in the computer model, as shown in Figure 6a,c,e. Accordingly, the failure criterion was defined according to the modified T-K criterion in the ANSYS software. Figure 6a,c,e show the simulated results, in contrast to the CT (computerized tomography) results of fractured samples (Figure 6b,d,f) with the low, medium, and high lamina angle ranges, respectively.
By contrast, for the high lamina angle (Figure 6a,b) outside the effective lamina angle range, stress was prone to concentrate on the surrounding rock body compared to the samples simulated in the medium and low lamina angle ranges (Figure 6c,e). Therefore, the typical X-conjugate fractures

Stress Distribution in the Laminated Rock
The laminated core model (cm-scale) in the finite element simulation is composed of two bodies, as shown in Figure 6, where the contact part between two bodies is settled as frictional with a sliding friction coefficient in accordance with the lamina friction coefficient in the failure criteria of this study. The failure behavior of the laminated rock in the special lamina angle range is determined by the stress distribution constrained by the constant value of µ lamina . Therefore, it is necessary to study the stress distribution of the laminated rock using the finite element simulation technology, especially for the different lamina angle ranges, i.e., the low, medium, and high lamina angle ranges.
Researchers have tried to simulate the compression behavior of tight intact rock and the laminated rock using finite element simulation technology [60,61]. To better understand the stress distribution during the compression process, the failure distribution around the compression samples was simulated by using the modified criterion in the ANSYS environment. The simulated samples were divided into two parts: The surrounding rock and the lamina body, with the lamina surface acting as the sliding surface. Additionally, the mechanical parameters of the surrounding rock used in the simulation were different from those of the laminated rock. Here, the surrounding rock was set as siltstone, while the laminated rock was claystone in the ANSYS software. In addition, the contact relationship, i.e., the lamina sliding surface, was set as frictional in the ANSYS software with the fractional coefficient value of 0.2. Another critical parameter, the lamina angle, was designed in the computer model, as shown in Figure 6a,c,e. Accordingly, the failure criterion was defined according to the modified T-K criterion in the ANSYS software. Figure 6a,c,e show the simulated results, in contrast to the CT (computerized tomography) results of fractured samples (Figure 6b,d,f) with the low, medium, and high lamina angle ranges, respectively.
By contrast, for the high lamina angle (Figure 6a,b) outside the effective lamina angle range, stress was prone to concentrate on the surrounding rock body compared to the samples simulated in the medium and low lamina angle ranges (Figure 6c,e). Therefore, the typical X-conjugate fractures were prone to propagate without associated bedding fractures. Notably, the effective lamina angle indicates that those compressed fractures were generally induced by the lamina body on which the stress was primarily concentrated.
were prone to propagate without associated bedding fractures. Notably, the effective lamina angle indicates that those compressed fractures were generally induced by the lamina body on which the stress was primarily concentrated. For the medium lamina angle in the effective lamina angle range, a shear fracture propagated between the two bedding fractures, as shown in the CT photo (Figure 6c,d), and was accompanied by some shear fractures between the laminated bodies. The fracture aperture was smaller than that in the core with a small lamina angle and was characterized by the shear sliding distance, as shown in Figure 6d. In addition, stress was prone to concentrate around the rock body rather than the sample simulated at a low lamina angle, which could be indicated by the CT photos (Figure 6a,c,d).
For the low lamina angle (Figure 6e,f) in the effective lamina angle range, the CT photo showed a bigger fracture aperture in the darker color. In addition, the simulated fracture index of the lamina was much higher than that of the surrounding rock, indicating a stronger propagation tendency for the bedding fractures caused by the stress concentration on the lamina.
Overall, stress in these compression samples was prone to concentrate on the lamina when the lamina angle was in the effective lamina angle (low-and medium-angle range). The low-angle lamina always induces fractures in an extensively open state with a bigger failure aperture, and the medium-angle lamina always induces fractures in a shear sliding trend. On the contrary, stress was prone to concentrate on the surrounding rock body when the lamina angle was outside of the effective range, and contributed less to the compressed bedding fractures in the laminated rock. For the medium lamina angle in the effective lamina angle range, a shear fracture propagated between the two bedding fractures, as shown in the CT photo (Figure 6c,d), and was accompanied by some shear fractures between the laminated bodies. The fracture aperture was smaller than that in the core with a small lamina angle and was characterized by the shear sliding distance, as shown in Figure 6d. In addition, stress was prone to concentrate around the rock body rather than the sample simulated at a low lamina angle, which could be indicated by the CT photos (Figure 6a,c,d).
For the low lamina angle (Figure 6e,f) in the effective lamina angle range, the CT photo showed a bigger fracture aperture in the darker color. In addition, the simulated fracture index of the lamina was much higher than that of the surrounding rock, indicating a stronger propagation tendency for the bedding fractures caused by the stress concentration on the lamina.
Overall, stress in these compression samples was prone to concentrate on the lamina when the lamina angle was in the effective lamina angle (low-and medium-angle range). The low-angle lamina always induces fractures in an extensively open state with a bigger failure aperture, and the medium-angle lamina always induces fractures in a shear sliding trend. On the contrary, stress was prone to concentrate on the surrounding rock body when the lamina angle was outside of the effective range, and contributed less to the compressed bedding fractures in the laminated rock.

Distribution of the Regional Bedding Fractures
In order to simulate the regional bedding fractures generated in a particular tectonic period using the modified criteria, previous simulation studies of the conventional structural fractures at the regional scale are used as the tectonic stress field [14,19,20,62]. Thus, the bedding fractures of the Yanchang Formation in the Ordos Basin generated during the latest orogenic episode (i.e., Himalayan episode) have been simulated by ameliorating   [14]. Specifically, as a result, the bedding fractures were calculated (Figure 7d) based on the improved rock failure criteria for bedding fractures (Equation (10)), the stress field simulated in   (Figure 7a,b) [14], the stratigraphic dip distribution map (Figure 7c), and the lamina fraction coefficient in different lithofacies (Figure 7d). Here, the stratigraphic dip at the regional scale (i.e., the intersection angle between the horizontal direction and the dip direction at the regional scale) is equal to the lamina angle at the cm-scale, because the maximum principal stress is in the horizontal direction during the orogenic periods, where the stratigraphic dips of the Yanchang Formation are valued in the relatively effective range for the bedding fracture failure criteria, as mentioned before (Figure 7c). In addition, µ lamina is assigned different values according to the different fitting results in Figure 5. The laminated siltstone sedimented in a delta front environment ( Figures 5 and 7d) is assigned 0.5, the laminated turbidite sandstone sedimented in a lacustrine turbidite environment ( Figures 5 and 7d) is assigned 0.3, and the mud sedimented in a still water environment (Figures 5 and 7d) is assigned 0.9.
The finite element simulation environment, including the boundary condition, Himalayan tensile stress field, geological model, and petrophysical model, remained consistent with those in   [14]. The bedding fracture index in the Yanchang Formation was calculated based on four variables, i.e., the maximum principal stress (σ 1 ), minimum principal stress (σ 3 ), lamina angle (θ), and lamina friction coefficient (µ lamina ). The calculated bedding fracture distribution, as shown in Figure 8, is broadly consistent with the structural fractured regions simulated by   [14], such as the typical oilfield regions of "Dingbian," "Xin'anbian," and "Zhiluo" (Figure 7c). The difference is that the simulated bedding fractures spread over an even larger area, which reveals that bedding fractures are more easily induced than the conventional structural fractures (i.e., the tensile, shear, and hybrid fractures) under the same stress condition. For the Himalayan tensile stress field of the Yanchang Formation in the Ordos Basin, it is obvious that the bedding fractures tend to propagate under the conditions of lower maximum principal stress, higher minimum principal stress, and higher stratigraphic dip values (Figure 7a-d). Generally, the bedding fractures distribution simulated in this work (Figure 7d) would be a supplement for the structural fractures simulated by   [14], and both of the two research works contribute to the further study of the natural fracture network underground (including the bedding fractures and the conventional structural fractures) for the Yanchang Formation of the Ordos Basin.

Conclusions
1. Two critical parameters were introduced to modify the T-K criterion, i.e., the lamina angle (θ) based on the stress equilibrium equation along the lamina surface, and the lamina friction coefficient (μlamina) related to the different laminated rock type.
2. Stress is prone to concentrate on the lamina when the lamina angle is in the effective low and medium range. The low-angle lamina always induces fractures in the expansively open state with the bigger failure aperture, and the medium-angle lamina always induces fractures in the shear sliding trend. On the contrary, stress is prone to concentrate on the surrounding rock body, when the lamina angle is outside of the effective range, and contributes less to the compressed fractures in the laminated rock.
3. The bedding fractures of the Yanchang Formation in the Himalayan stress field were simulated based on the new rock failure criteria proposed in this work. It is concluded that these bedding fractures tend to propagate under the conditions of lower maximum principal stress, higher minimum principal stress, and larger stratigraphic dip.

Conclusions
1. Two critical parameters were introduced to modify the T-K criterion, i.e., the lamina angle (θ) based on the stress equilibrium equation along the lamina surface, and the lamina friction coefficient (µ lamina ) related to the different laminated rock type.
2. Stress is prone to concentrate on the lamina when the lamina angle is in the effective low and medium range. The low-angle lamina always induces fractures in the expansively open state with the bigger failure aperture, and the medium-angle lamina always induces fractures in the shear sliding trend. On the contrary, stress is prone to concentrate on the surrounding rock body, when the lamina angle is outside of the effective range, and contributes less to the compressed fractures in the laminated rock.
3. The bedding fractures of the Yanchang Formation in the Himalayan stress field were simulated based on the new rock failure criteria proposed in this work. It is concluded that these bedding fractures tend to propagate under the conditions of lower maximum principal stress, higher minimum principal stress, and larger stratigraphic dip. Then, ϕ and S 0 can be expressed as Equations (A6) and (A7): ϕ = arcsin µ lamina cos θ + sin θ (µ lamina + 1) cos θ + (1 − µ lamina ) sin θ (A6) S 0 = (1 − sin ϕ)C lamina 2 cos ϕ(cos θ − µ lamina sin θ) Thus, S 0 can be further expressed as Equation (A8) when Equation (A6) is put into Equation (A7).