In the research process of this study, we utilized a variety of sources for our research (Figure 1D). We first used core paleomagnetism and velocity anisotropy experiments, combined with imaging logging and dipole shear wave interpretation, to evaluate the direction of in-situ stress. Using core differential strain experiments and hydraulic fracturing test data, completed the interpretation of in-situ stress magnitude. Secondly, we constructed a transverse isotropic in-situ stress calculation model and assigned values to all parameters within the model. Finally, we performed a detailed evaluation of in-situ stress for wells. The following sections provide a detailed explanation of the core experimental methods.

The samples used in this experiment were primarily obtained from core samples taken from ten wells in the Yingshan-Pingchang area of northeastern Sichuan. These samples underwent paleomagnetic testing, velocity anisotropy testing, differential strain testing, and rock mechanics testing.

Combining paleomagnetic and velocity anisotropy experimenting, the results show that the current in-situ stress orientation of the Jurassic formation in the Pingchang area ranges from N82.7° to N116.65°E, with the main direction concentrated at N95° ± 10°E. In the Yingshan area, the stress orientation range is from N76° to N111.6°E. The maximum principal stress orientation at well YY1H ranges from 71.77° to 106.45°, while at well PY1, it ranges from 59.1° to 85.6° (Table 1). Due to the influence of local NW-trending structures, the in-situ stress orientation is deflected clockwise, and the three-dimensional stress regime exhibits a strike-slip stress regime. The stress orientation at well YY1H shows a slight clockwise deflection due to its proximity to the fault.

The resistivity differences in wellbore imaging logging can explain rock fractures caused by stress. By utilizing these characteristics, the orientation of the maximum horizontal stress can be further determined (Andrew et al., 2016; Ju et al., 2018). During the drilling process, due to the coupling of wellbore stress concentration and mud pressure, induced fractures or wellbore breakouts can form under specific conditions. The former is manifested in wellbore imaging logging as approximately vertical or feather-like dark bands or a single feather-like dark band (Figures 4A, B), with their orientation indicating the orientation of the maximum horizontal principal stress (Zoback et al., 2003; Colleen et al., 2009; Schmitt et al., 2012). The latter appears in wellbore imaging logging as an approximately symmetric distribution of dark bands or patches, which are discontinuous along the longitudinal orientation with blurred boundaries (Figures 4C, D). The long axis of the breakout is generally parallel to the orientation of the minimum horizontal principal stress (Zoback et al., 2003; Ju et al., 2018). The wellbore imaging logging at depths of 1044–1046 m in well YS8 characteristics is represented by feather-like dark bands arranged in a fan-shaped pattern, while at depths of 1581–1583 m in well YS5, a single feather-like dark band is observed. Based on the identified orientation of the induced fractures, the orientation of the maximum horizontal principal stress for the corresponding intervals in the two wells is 185° and 176°, respectively. The wellbore imaging logging at depths of 3,424–3,426 m in well DY1 and 1581–1583 m in well YS8 show significant wellbore breakouts, with breakout orientations of 6° and 80°, respectively. The interpreted orientations of the maximum horizontal principal stress are 96° and 170°, respectively. In intervals with more developed fractures, false breakout phenomena can easily be observed in borehole imaging imaging. Caution should be exercised during the identification process to avoid misinterpretation.

Electrical imaging logging interpretation software was utilized to identify induced fractures and wellbore breakouts in the DY1 well. The interpretation results indicate that the current stress orientation is primarily EW (Figure 4E).

For anisotropic formations, by considering the variation in shear wave velocities under different array configurations, can reveal the anisotropic characteristics of the subsurface rock layers (Heng et al., 2015). In anisotropic strata, in-situ stress primarily acts along different directions of the layers. Since the shear wave velocity is closely related to the anisotropy, particularly when the direction of shear wave propagation aligns with the stress direction, changes in velocity can reflect the direction of the in-situ stress. The polarization orientation of the fast shear wave indicates the maximum horizontal stress orientation. The time difference between the fast and slow waves is commonly used to measure the anisotropy of the formation. The stronger the anisotropy, the higher the reliability of the extracted data. The Formula 2 for calculating the anisotropy coefficient is as follows: I D T = ∆ t s s − ∆ t s f ∆ t s s + ∆ t s f / 2 × 100 % (2)

In this formula: ∆ t s s is the time difference for the fast shear wave, in μs/ft; ∆ t s f is the time difference for the slow shear wave, in μs/ft.

The interpretation of dipole shear wave logging indicates that the current in-situ stress orientation is primarily EW (Figure 5). The anisotropy analysis results show that as the difference coefficient between the two principal stresses increases, shear wave anisotropy becomes stronger, reducing the number of interpreted maximum principal stress orientations and increasing the reliability of the interpretation.

The current stress orientation in the Yingshan-Pingchang area of the Jurassic Formation, determined by comprehensive analysis of differential strain, electrical imaging logging, and dipole shear wave logging, is primarily in the EW orientation (Figure 6A). The orientation ranges from N 87°–120° E, with the main direction maintaining a consistent range of NE 95° ± 10°.

The overall in-situ stress orientation of the Jurassic formation in the study area is consistent with that of the neighboring Yuanba and Langzhong areas (Figure 6B). In the Yingshan area, the local tectonic influence causes a clockwise deviation of the in-situ stress orientation, primarily ranging from N60° to 80°E. Taking YS6 Well as an example (Figure 6C), the in-situ stress orientation at the target layer between 1,650 and 1,900 m is mainly influenced by the NWW-trending fault, resulting in a deviation of the stress orientation. Figure 6D presents a simulation showing the area’s maximum principal stress deviating due to its 25° angle with a local NW-trending fault. The direction of the maximum horizontal principal stress clearly deviates along the fault strike and propagates in a localized manner along it. Fault activity leads to the displacement of rock layers. When a fault forms a specific angle with the orientation of the maximum horizontal principal stress, the distribution of the original stress field undergoes significant alteration. In such cases, the fault exerts a strong guiding effect on the in-situ stress orientation. In areas distant from the fault (depths less than 1,650 m and greater than 1,900 m), the in-situ stress orientation at Yingshan Well six remains generally stable, exhibiting a characteristic EW orientation without noticeable deviation. This stability is because, away from the fault zone, the rock layers experience relatively uniform compression or tension and are not significantly affected by fault activity. Regional tectonic stress plays a dominant role during this period.

During the differential strain experiment, strain data were directly collected from nine directions on three orthogonal surfaces of the rock sample. The ratio of the three-dimensional stress was then calculated using elasticity theory. The pressure at the microfracture closure point was used to estimate the maximum principal stress, or alternatively, the vertical stress was inferred from the depth of the core sample. Finally, the three-dimensional stress was estimated by combining the principal stress ratios (Figure 7).

The results of the three-dimensional stress from the differential strain experiment indicate that the maximum horizontal principal stress ranges from 42.33 MPa to 102.56 MPa, averaging 74.89 MPa; the minimum horizontal principal stress ranges from 39.20 MPa to 84.04 MPa, averaging 67.20 MPa; the vertical principal stress ranges from 31.91 MPa to 91.39 MPa, averaging 60.23 MPa (Table 2).

Based on the stress-strain relationship along the three principal orientations, the ratio of the magnitudes of the three-dimensional stress can be obtained. The stress ratio of the tested samples indicates that the maximum horizontal principal stress value > vertical principal stress value > minimum horizontal principal stress value, suggesting a strike-slip stress regime. Statistical analysis of the minimum principal stress gradient for different lithologies shows that, the minimum principal stress gradient of shale is higher than that of limestone and sandstone (Figure 8A). Due to the high clay mineral content in shale, the minimum horizontal principal stress gradient tends to be relatively high. According to the formula for calculating the minimum horizontal principal stress, a higher content of plastic minerals results in a larger tectonic strain coefficient, which in turn leads to an increase in the minimum horizontal stress gradient. Overall, the Jurassic stress values are higher in the Pingchang area than in the Longgang area and higher than in the Yingshan area due to the burial depth (Figure 8B).

The hydraulic fracturing method for in-situ stress measurement relies primarily on the pressure-time curve recorded during the fracturing process. By integrating relevant rock mechanics parameters and fracturing theory, the magnitudes of in-situ stresses can be calculated. Compared to other measurement methods, this approach does not require prior knowledge of the rock’s mechanical parameters, yet it enables the determination of multiple stress parameters within the formation. Additionally, it offers several advantages, including simple equipment, convenient operation, significant measurement depth, the ability for continuous or repeated testing at any depth, rapid measurement speed, intuitive measurement values, and comprehensive representation of the stress field (Schmitt et al., 2012; Wang et al., 2020). The use of fracturing data to determine in-situ stress is one of the most direct and reliable methods currently available.

Through the analysis and calculation of fracturing data (Table 3), it is concluded that the shut-in pressures for wells PA101, PY1-1H, PY1-2H, and PY1-3H mainly range from 35 to 42 MPa, averaging 37.01 MPa. Using the fracturing data, the in-situ stress magnitudes calculated by the hydraulic fracturing method are obtained (Table 3). Taking well PA101H as an example, the calculated maximum horizontal principal stress is 75.52 MPa, and the minimum horizontal principal stress is 66.187 MPa. A comparison of the results from the differential strain experiment with the hydraulic fracturing stress calculation results at the same depth (Table 3) shows that the in-situ stress values measured in the differential strain experiment are close to those calculated by hydraulic fracturing, demonstrating considerable accuracy.

Based on the calculated magnitudes of in-situ stresses derived from construction parameter data, it was determined that the relationship among the stresses follows the pattern of maximum principal stress > vertical principal stress > minimum principal stress. Comprehensive analysis suggests that the stress regime in the Jurassic formation in this area is characterized by a strike-slip faulting stress type, which is consistent with the results obtained from the differential strain experiment.

The magnitude of horizontal in-situ stress is influenced by factors such as vertical in-situ stress, formation pressure, and tectonic stress. Therefore, accurate calculation of vertical in-situ stress and pore pressure is essential. Typically, vertical stress is determined by integrating the density log curve. However, not all well sections possess a density log curve. Consequently, vertical stress is initially calculated using the average density of sections lacking a density log curve and subsequently augmented by the vertical stress obtained from integrating the density log curves of sections with density logs. The Formula 3 for calculating the vertical stress at a given depth in the formation is as follows: p o = 0.00981 ρ av H 1 + ∫ H 1 H 2 ρ h d h (3)

In this formula: P _o is the vertical stress, in MPa; ρ av is the average overburden density of the intervals without density log curves, in g/cm²; p is the density of each individual rock layer, in g/cm²; H ₁ is the vertical depth of the deepest point of the interval without density log curves, in meters; H ₂ is the vertical depth of the overburden pressure calculation point (where density log curves are available), in meters; h is the formation thickness corresponding to the sampling interval depth of the log curve, in meters (commonly taken as 0.125 m).

The overburden pressure gradient in the Ying Shan area mainly ranges from 2.46 to 2.5 MPa/100 m, while in the Longgang area, it ranges from 2.38 to 2.40 MPa/100 m, and in the Pingchang area, it ranges from 2.46 to 2.48 MPa/100 m.

Pore pressure is a critical parameter in drilling processes, fracturing design, and various other applications, and it is essential for evaluating in-situ stresses. Based on an analysis of the target formation in the study area, the method adheres to Teraghi’s effective stress theory. According to this theory, the overburden pressure within a given layer is partitioned into pore fluid pressure, supported by the pore fluids, and rock stress, supported by the solid grains of the bedrock. The relationship between overburden pressure, rock stress, and pore fluid pressure is calculated by Formula 4: P p = P 0 − σ (4)

In this formula: P ₀ is the overburden pressure, in MPa; P _P is the pore fluid pressure, in MPa; σ is the effective stress, in Mpa.

Based on the measured formation pressure data, well logging interpretation results of lithology, physical properties, and electrical response analysis, and considering the relationships among clay content, porosity, P-wave velocity, and effective stress, a multi-parameter regression algorithm was employed. To establish a formation pressure logging interpretation model that integrates data from different regions. The model is as follows in Formulas 5, 6:

Pingchang σ = 304.8 A C + 23.1 − 1.31 φ + 0.9 V s h (5)

Yingshan σ = 304.8 A C + 18.33 − 0.13 φ − 21.45 V s h (6)

In this formula: AC is the acoustic time delay from the sonic log, in μs/ft; φ is the porosity, in %; V _sh is the shale content, in %.

The relationship between stress and strain in a transversely isotropic medium (VTI medium) satisfies the generalized Hook’s law 4, and can be expressed as follows in Formula 7: C 11 C 11 − 2 C 66 C 13 0 0 0 C 11 − 2 C 66 C 11 C 13 0 0 0 C 13 C 13 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 66 (7)

In this formula: C is the stiffness matrix, in GPa; ε is the strain; C ₁₁ is the longitudinal modulus along the parallel to the bedding plane, in GPa; C ₃₃ is the longitudinal modulus perpendicular to the bedding plane, in GPa; C ₄₄ is the shear modulus perpendicular to the bedding plane, in GPa; C ₆₆ is the shear modulus along the parallel to the bedding plane, in GPa; and C ₁₃ is the stiffness modulus, in GPa.

The symmetry type of formation rocks can be determined by the number of independent elastic parameters. An isotropic formation has two independent elastic parameters, while a transversely isotropic formation has five independent elastic parameters. To describe the stress-strain relationship in a VTI formation, the stiffness coefficients C ₁₁ , C ₃₃ , C ₄₄ , C ₆₆ , and C ₁₃ need to be determined. These coefficients are calculated in Formulas 8-12, respectively, as follows: C 11 = ρ V p h   2 (8) C 33 = ρ V p v   2 (9) C 44 = ρ V s v   2 (10) C 66 = ρ V s h   2 (11) 4 ρ V p 45     2 − C 11 − C 33 − 2 C 44 2 − C 11 − C 33 2 4 1 2 − C 44 (12)

In this formula: V p h is the P-wave velocity parallel to the bedding plane, in km/s; V p v is the P-wave velocity perpendicular to the bedding plane, in km/s; V s h is the S-wave velocity parallel to the bedding plane, in km/s; V s v is the S-wave velocity perpendicular to the bedding plane, in km/s; V p 45 is the P-wave velocity at a 45° angle to the bedding plane, in km/s.

Using an ultrasonic transducer system, the velocity of sound is measured on core plug samples extracted from large natural shale samples with vertical stratification. The core plug samples are taken at 0°, 45°, and 90° orientations. The corresponding five independent elastic parameters at the specified depth points can be calculated (Table 4).

C ₃₃ is calculated using the compressional wave velocity obtained from density and sonic data from well logging data. C ₄₄ is calculated from the shear wave velocity derived from density and sonic data. C ₆₆ can be determined by inversing the Stoneley wave velocity to compute the horizontal shear wave velocity and density. The parameters C ₁₃ and C ₁₁ need to be obtained from core experiments or empirical formulas. Based on core sample experimental data, Suarez-Rivera and Bratton proposed an improved ANNIE model called MANNE. By combining experimental data and using multiple linear regression to determine the coefficients of the MANNE model, resulting in the relationship Formulas 13-16 for C ₁₃ and C ₁₁ in vertical wells, sandstone, and shale as follows:

Vertical Well (Sandstone) c 13 = 0.795947 c 33 − 1.24059 c 44 (13) c 11 = 1.008131 c 13 (14)

Vertical Well (Shale) c 13 = 0.928115 c 33 − 1.66159 c 44 (15) c 11 = 0.916093 c 13 (16)

The Biot coefficient is a parameter used to describe the coupling relationship between the solid matrix and the pore pressure in a porous medium. It reflects the contributions of both the solid matrix and the pore space to the overall properties of the porous material. To obtain a more accurate estimate of the Biot coefficient, this study utilizes the effective stress theory. By employing experimentally measured parameters such as porosity, bulk modulus, and elastic modulus of the rock (Table 4), the Biot coefficient is estimated through theoretical derivation. The Biot coefficient is typically denoted by the symbol α and is defined as Formula 17: α = 1 − K d r y K s o l i d (17)

In this formula: K _dry is the rock’s bulk modulus; K _solid is the bulk modulus of the skeleton.

After obtaining all the stiffness coefficients, the anisotropic rock mechanical parameters, such as dynamic elastic modulus and Poisson’s ratio, are calculated using the corresponding formulas. Based on the rock mechanical parameters derived from both mechanical experiments and acoustic logging data, a correction relationship model is established for the Pingchang and Yingshan areas, Formulas 18-25 are shown below:

Pingchang μ v = 1.6588 μ d v − 0.2339 (18) μ h = − 0.9972 μ d h + 0.4638 (19) E v = 0.8981 E d v − 14.502 (20) E h = 0.9003 E d h − 15.632 (21)

Yingshan μ v = 1.848 μ d v − 0.3127 (22) μ h = − 1.9646 μ d h + 0.6507 (23) E v = 1.7827 E d v − 26.312 (24) E v = 1.7827 E d v − 26.31 (25)

In this formula: E _dv and E _dh is the dynamic Young’s modulus in the vertical and horizontal directions, in GPa; μ d v and μ d h is the dynamic Poisson’s ratio in the vertical and horizontal directions, with no units.

Based on experimental measurements of the skeleton compression coefficient and the rock volume compression coefficient for the rock samples, a statistical analysis is conducted on the Biot coefficient in relation to parameters such as the acoustic travel time, rock bulk density, and other factors. A regression Formula 26 for the Biot coefficient with the GR and triple porosity logging curves is established as follows: B i o t = 2.791863 − 0.006 G R + 0.002568 A C − 0.00245 C N L − 0.65217 D E N (26)

Using core experimental data combined with previously obtained parameters, including vertical stress, formation pore pressure, Young’s modulus, Poisson’s ratio, and Biot’s coefficient, the maximum and minimum tectonic stress coefficients are subsequently inverted based on the transversely isotropic (Sn) model Equations 27, 28.

Higgins et al. analyzed the anisotropic characteristics of VTI formations and developed the Sn model for interpreting the horizontal principal stress in VTI formations. The formula for calculating the horizontal principal stress in a transversely isotropic medium is as follows: σ h = E h o r z   E v e r t   v v e r t   1 − v h o r z   σ z − a σ P + E h o r z   1 − v h o r z   2 ε h + E h o r z   v h o r z   1 − v h o r z   2 ε H + a σ P (27) σ H = E h o r z   E v e r t   v v e r t   1 − v h o r z   σ z − a σ P + E h o r z   1 − v h o r z   2 ε H + E h o r z   v h o r z   1 − v h o r z   2 ε h + a σ P (28)

In this formula: σ h is the minimum horizontal principal stress, in MPa; σ H is s the maximum horizontal principal stress, in MPa; ε h is the minimum tectonic stress coefficient, with no units; ε H is the maximum tectonic stress coefficient, with no units; σ z is the overburden pressure, in MPa; σ P is the pore pressure (MPa); V _vert is the vertical static Poisson’s ratio under the transversely isotropic model, with no units; V _horz is the horizontal static Poisson’s ratio under the transversely isotropic model, with no units; E horz   is the horizontal static Young’s modulus under the transversely isotropic model in MPa; E vert is the vertical static Young’s modulus under the transversely isotropic model in MPa; a is the Biot coefficien, with no units.