The Kuqa Depression is located in the north of Tarim Basin. It is connected with the South Tianshan fault fold belt by the thrust fault in the north, Tabei Uplift in the south, Yangxia sag in the east, and Wushi sag in the West. It is a superimposed foreland basin dominated by Mesozoic and Cenozoic sedimentation and can be divided into seven secondary structural units (Figure 1A). Among them, the Kelasu tectonic belt is a typical deep and ultra-deep natural gas resource enrichment area. From west to east, it can be divided into Awate segment, Bozi segment, Dabei segment, and Keshen–Kela segment. The Kelasu tectonic belt develops four first-order faults from north to south, which are divided into Bozi–Kela fault tectonic belt, Keshen fault tectonic belt, Baicheng fault tectonic belt, and Baicheng South fault tectonic belt. Multiple secondary thrust faults are developed between the faults.

The Bozi gas field is located in the west of Kelasu tectonic belt and is characterized by rows of faulted anticlines. The strike of anticline is generally the same, which is in the NE–SW direction. The northern Bozi gas field is a pop-up structure sandwiched by two thrust faults F1 and F2, and the southern part is a series of imbricate structures, with gradually increasing burial depth (Figures 1B,C). From east to west, the overlapping degree of the imbricate structure increases gradually.

The drilling of the Bozi gas field shows that the drilling strata from the top to bottom are: Quaternary Xiyu Formation (Q₁ x), Neogene Kuqa Formation (N₂ k), Kangcun Formation (N_1-2 k), Jidike Formation (N₁ j), Paleogene Suweiyi Formation (E_2-3 s), Kumugeliemu Group (E_1-2 km), Cretaceous Bashijiqike Formation (K₁ bs), and Baxigai Formation (K₁ bx). Among them, The E_1-2 km is a set of gypsolith rock formation with great thickness distribution and plastic flow characteristics. In the process of compression deformation, layered shrinkage deformation occurs, and the strata above this plastic stratum folds, thrusts, and rises sharply, forming a “roof like” structure. Under the control of “roof” structure, salt rock flow results in the characteristics of abovementioned salt Kelasu tectonic belt, salt body, and subsalt differential structural deformation.

The target strata of the study are the Cretaceous K₁ bs and K₁ bx. The K₁ bs can be divided into three lithologic segments. The first lithologic segment is missing due to the influence of paleosedimentary environment. The second lithologic segment mainly comprises medium thick ∼ super thick layered brown and grayish brown medium sandstone and fine sandstone, which is deposited at the front of braided river delta. The third lithologic segment mainly comprises medium thick ∼ thick layered brown fine sandstone, thin ∼ medium thick layered siltstone, and argillaceous siltstone, which are deposited at the front of fan delta. The K₁ bx comprises thin medium thick brown mudstone and medium thick layered brownish gray, grayish brown fine sandstone, siltstone, and argillaceous siltstone, which are interbedded in an isopachous-slight isopachous manner and are deposited at the front of braided river delta (Figure 2).

Reservoir pore types in the Bozi gas field are mainly primary intergranular pores, with a small amount of microfractures and micropores. The porosity of Bashijiqike formation mainly ranges between 4% and 10.0%, with an average of 6.30%, and the permeability mainly ranges between 0.1 and 0.5 mD, with an average of 0.256 mD. The porosity of Baxigai formation sandstone ranges between 4 and 9%, with an average porosity of 6.06%, and the permeability mainly ranges between 0.035 and 0.1 mD, with an average permeability of 0.067 mD. The reservoir fractures are relatively developed, mostly being partially filled—unfilled, high angle, and shear fractures. The fracture strike is mainly near E–W, and the fracture density ranges between 0.2/m and 0.4/m.

At present, there are many methods to measure and test in situ stress, including the drill-based borehole collapse method, relaxation strain measurement, stress relief technique, acoustic emission, and differential strain method (Zang and Stephansson, 2010). The Bozi gas field is deep buried and has high bottom hole temperature and high pressure. Hence, the measuring instrument has poor applicability, and there is no constraint of confining pressure after core removal, causing rapid stress relief. In mild cases, it may lead to inaccurate measurement, and in serious cases, it may lead to possible core fracture. Therefore, it is not suitable to use downhole measurement and core testing to obtain the in situ stress state.

This article utilizes the drilling and logging data to obtain the in situ stress state. Different basins are in different states, and scholars have put forward a variety of in situ stress calculation models considering more comprehensive factors. As a result, the accuracy of the calculation has gradually improved, and the applicability has become more extensive. At present, there is little dispute about the calculation of vertical stress (S_v). It is considered that its magnitude is equal to the gravity of the overlying rock mass, which can be obtained by integrating the density curve from the ground to the target depth such as Eq. 1: S V = ∫ 0 z ρ ( z ) g d z , (1) where g is the gravitational acceleration (m/s²), ρ(z) is the density as a function of burial depth (kg/m³), and z is the burial depth (m).

There are many calculation models for obtaining horizontal principal stress magnitudes, including the uniaxial strain model, Mohr–Coulomb failure model, Coulomb–Navier failure model, Huang’s model, combined spring model, porous elastic strain model, and biaxial strain model (Li and Zhang, 1997; Liu and Luo, 1999; Rasouli and Sutherland, 2014; Li et al., 2017). In this article, the calculation model is as follows (Li et al., 1997): { S H = μ 1 − μ ( S V − α P p ) + E ξ H 1 − μ 2 + μ E ξ h 1 − μ 2 + α P p S h = μ 1 − μ ( S V − α P p ) + E ξ h 1 − μ 2 + μ E ξ H 1 − μ 2 + α P p , (2) where S_H is the maximum horizontal principal stress, MPa; S_h is the minimum horizontal principal stress, MPa; S_V is the vertical principal stress, MPa; P_P is the pore pressure, MPa; μ is Poisson’s ratio, dimensionless; E is the modulus of elasticity, GPa; α is the Biot coefficient, dimensionless; and ξ H ξ h are the maximum and minimum principal stress coefficient, respectively, dimensionless.

The method of calculating the mechanical parameters with the logging data is as follow (Lu et al., 2015): E d = ρ b Δ t s 2 ⋅ 3 Δ t s 2 − 4 Δ t p 2 Δ t s 2 − Δ t p 2 , (3) μ d = Δ t s 2 − 2 Δ t p 2 2 ( Δ t s 2 − Δ t p 2 ) , (4) where ρ _b is the rock density, kg/m³, and Δt _p and Δt _s are the P-wave time difference and S-wave time difference, respectively, μs/ft.

Generally, we may determine S_h at a specific location through the hydraulic fracturing construction data, which can be used as the constraint and scale basis to determine the value of ξ h . However, the maximum horizontal principal stress (S_H) of ultra-deep wells cannot be determined through hydraulic fracturing (Hickman, 1983; Zoback, 2007). In this article, the S_H is determined with the wellbore fracture information. In the process of drilling, as the borehole is drilled, the borehole wall may produce stress concentration under the action of confining pressure. When the stress concentration exceeds the fracture strength of the rock around the borehole, the borehole wall will collapse (Bell and Gough, 1979; Zoback et al., 2003).

Because wells in different structural locations have different stress distributions (stress concentration and tensile range) around the well in different lithologic sections, and with the increase of radial depth, the stress effect will decrease and the degree of borehole wall collapse will also vary. Therefore, it is a more intuitive method to distinguish the direction of in situ stress by analyzing the changes of formation characteristics around the borehole from the scanning images of formation micro resistivity (FMI). The collapse width of the borehole wall can be judged from FMI, and the collapse width has the following mathematical calculation relationship with rock uniaxial compressive strength and in situ stress state (Zoback, 2007): S H = ( C 0 + 2 P p + Δ p + σ Δ T ) − S h ⋅ [ 1 + 2 ⁡ cos ( π − W B O ) ] 1 − 2 ⁡ cos ( π − W B O ) , (5) where C₀ is the uniaxial compressive strength, P_p is the pore pressure, Δ p is the difference between pore pressure and bottom hole pressure, W_BO is the collapse width, and σ Δ T is the thermal stress, which can be ignored generally.

Therefore, the gradient of the maximum horizontal principal stress and the minimum horizontal principal stress can be inverted according to the collapse width W_BO and the uniaxial compressive strength (C₀) of the borehole wall collapse position.

In this article, a 1D mechanical earth modeling has been carried out in the Bozi gas field. Figure 4 shows the 1D MEM result along the well section of A2–A101-2–A102-4–A105, which indicates the distribution and changes of Cretaceous geomechanical parameters in one dimension of the wellbore.

Taking Well A2 as an example, Young’s modulus is 20 ∼ 30 GPa, the average uniaxial compressive strength is approximately 100 MPa, and Poisson’s ratio changes little, at 0.25 ∼ 0.30. The minimum horizontal, maximum horizontal, and vertical principal stress gradients are 2.11 ∼ 2.18 MPa/100 m, 2.55 ∼ 2.72 MPa/100 m, and 2.45 MPa/100 m, respectively, which are in strike–slip stress regime (S_H>S_V>S_h). For a better comparison, the S_h is filled with colors, with warm color representing low value and cold color for high value. It can be seen that the in situ stress shows obvious stratification vertically. The stress in the upper part of K₁ bs ₂ and K₁ bs ₃, that is, 7456 ∼ 7555 m is low, and the S_h magnitude is approximately 165 MPa. The stress in the lower part of K₁ bs ₃, the upper part of K₁ bx ₁ and K₁ bx ₂, that is, 7555 ∼ 7700 m is high, and the S_h magnitude is approximately 175 MPa, and local areas can reach 180 MPa. The stress in the lower part of K₁ bx ₂ decreases, and the S_h magnitude is approximately 170 MPa. This shows the distribution characteristics of “low–high–low”.

The in situ stress of single wells in the Bozi gas field has certain similar characteristics. The stress has the distribution characteristics of “low–high–low ” in the Cretaceous system. That is, the stress in the upper part of K₁ bs ₂ and K₁ bs ₃ is low, the stress in the lower part of K₁ bs ₃, and the upper part of K₁ bx ₁ and K₁ bx ₂ is high, and the stress in the lower part of K₁ bx ₂ is reduced. In addition, compared with the Baxigai Formation, the stress characteristics of Bashijiqike Formation are lower, and stress difference between the layers is relatively small. The interlaminar stress difference in the Baxigai Formation is large, resulting in a stronger interlaminar property.

The 3D distribution of the in situ stress in the Cretaceous layes of Bozi gas field is shown in Figure 5. The distributions of three principal stresses are similar, showing a trend of low in the north and high in the south. The S_h is mainly 153 ∼ 180 MPa, and S_H is 180 ∼ 200 Pa, and S_V is 175 ∼ 200 MPa, which is a strike slip stress state. There are great differences in stress distribution among the different fault blocks. The S_h in A104 and A102 fault blocks is 153 ∼ 160 MPa and the S_H is approximately 180 ∼ 188 MPa, while the S_h in A8 and A9 fault blocks is as high as 175 MPa and the S_H even exceeds 200 MPa.

Figure 6 shows the in situ stress distribution of five typical sections in the Bozi gas field. It can be seen that the in situ stress distribution is evidently related to the tectonic form. Each fault block in the Bozi gas field is basically a fault anticline structure, and different structural styles vary in terms of stress distribution. Figure 7 summarizes the structural style of the fault anticline in the Bozi gas field and the corresponding in situ stress distribution mode. The structural style of the fault anticline is divided into five categories: north steep and south gentle type (A8 and A10), north gentle and south steep type (A7), symmetrical imbricate type (A104), complete pop-up structure (A7), and semi complete pop-up structure (A1). Generally, the in situ stress distribution of symmetrical imbricate faulted anticline and pop-up structure is regular and symmetrical. The stress is low at the top of the anticline and high at the wings. In north steep and south gentle type, the stress is high in the north wing and low in the south wing. In the north gentle and south steep type, the stress in low in the north wing and high in the south wing.

The in situ stress prediction results of seven wells in different fault blocks of the Bozi gas field are compared with the logging calculation results (Table 1). The calculation method of error (r) is as follows: r = S c - S m , (9) where S_m is the in situ stress value actually calculated, MPa and S_c is the in situ stress predicted by numerical simulation, MPa. It can be seen that the average error of the S_H value is 6.5 MPa, and the average error of the S_h value is 6.3 MPa. Both the errors are within 10%, indicating that the in situ stress simulation results have high reliability.

The in situ stress distribution of the thrust imbricate structure is complex, especially the disturbance effect of the fault on the in situ stress (Meyer et al., 2007; Rasouli and Sutherland, 2014). In this study, the in situ stress shows a special distribution feature in the thrust imbricate structure (Figure 6). In the area around a fault, the in situ stress in the footwall is lower than that of the hanging wall, even if the footwall is buried deeper. Moreover, the greater the fault offset, the greater is the stress difference between the hanging wall and footwall, that is, the in situ stress value within the control range of the footwall is lower. According to the statistics of stress values from the hanging wall and the footwall of 20 faults (Figure 8), when the fault offset is less than 50 m, the stress difference is less than 5 MPa, and while it exceeds 500 m, the stress difference exceeds 15MPa, and there is a good correlation between the fault distance and stress difference between the hanging wall and footwall. In addition, the footwall of the overthrust generally shows low stress value. The higher the degree of overthrust, the larger is the range of low stress area of the footwall.

Previous studies (Liu et al., 2004; Li et al., 2012; Liu JS. et al., 2022) indicate that the wellbore trajectory is the most stable along the direction of the S_H under the strike–slip stress regime. The practice of Tarim Oilfield also confirms that under the ultra-deep strike–slip stress regime, highly deviated wells have fewer complex drilling accidents and shorter drilling cycle compared with the vertical wells (Xu et al., 2020; Cai et al., 2022). Therefore, in the Bozi gas field, highly deviated wells or horizontal wells have better wellbore stability.

The in situ stress field in the Bozi gas field is heterogeneous, and the in situ stress distributes in the belt horizontally and stratified vertically. The directional wells (inclined or horizontal wells) are more likely to drill in favorable low-stress zones than the vertical wells, and the directional trajectory or sidetrack drilling may be adjusted during the drilling based on in situ stress field.

Due to the high degree of overthrust in the Bozi gas field and the favorable low-stress area is often in the footwall, the vertical well method cannot meet the drilling conditions, while the directional well method can avoid the obstacles of nontarget layer and drill into the favorable area of target layer. The existing practice in Tarim Oilfield also proves that the production of directional wells is higher than that of the vertical wells. Therefore, it is considered that the highly deviated wells are more suitable for hydrocarbon exploitation in the Bozi gas field.