Analysis of Hydraulic Fracturing Efﬁciency Considering the Principal Stress in Brushy Canyon Formation of the Permian Basin

: The purpose of this study is to investigate the effect of principal stress direction on the efﬁciency of hydraulic fracturing treatment. There are two different drilling scenarios: 1. Four horizontal wells drilled in four orthogonal directions regardless of in-situ stress condition (“Actual”). 2. Three horizontal wells drilled equivalent to “Actual” case by considering the direction of principal stress (“Proposed”). The hydraulic fracturing modeling was carried out based on well logging data and completion reports of Brushy Canyon formation, Permian Basin. In the results of “Actual” case, transverse fractures were generated in two horizontal wells drilled parallel to σ hmin -dir (direction of σ hmin ), similar to “Proposed” case. Meanwhile, for two other wells drilled perpendicular to σ hmin -dir, longitudinal fractures were generated. These obliquely deviated fractures signiﬁcantly decreased the fracture spacing between the stages up to 26%. This induced great stress shadow, however, the fractures propagated straight due to the large stress anisotropy of 2000 psi ( σ Hmax / σ hmin = 1.4). Therefore, it was found that due to the different direction of fracture propagation in “Actual” case, “Proposed” case was 14.6% of stimulated reservoir volume (SRV) higher. In conclusion, for successful hydraulic fracturing treatment, the direction of horizontal well must be determined in consideration of the principal stress direction as well as stress anisotropy.


Introduction
In the United States, wells drilled horizontally into tight oil and shale gas formations account for an increasing share of crude oil and natural gas production. By the end of 2018, about 96% of U.S. crude oil production in tight oil formations came from horizontal wells. In the Permian Basin, oil well density reaches up to 56 wells/km 2 [1]. In a typical design of section development of horizontal wells, the laterals have regular spacings and are aligned in consideration of the formation characteristics such as the in-situ principal stresses, permeability, and poroelasticity, etc. [2]. In the normal fault regime, hydraulic fractures (HFs) open in the direction of the minimum horizontal principal stress (σ hmin -dir) and propagate in the direction of the maximum horizontal principal stress (σ Hmax -dir). Therefore, an ideal horizontal well is designed to be parallel to the orientation of σ hmin to induce "transverse fractures", i.e., perpendicular to the wellbore [3]. In general, transverse fractures provide a higher stimulated reservoir volume (SRV) than longitudinal fractures [4]. However, in the target field of this study, the Forty Niner Ridge field located in the Permian Basin in Eddy County, four horizontal wells were drilled orthogonally to each other from 2006 to 2014.
In 2011, studies on the Barnett and Marcellus formations that analyzed the efficiency of hydraulic fracturing according to the drilling azimuth of a horizontal well were published [5,6]. In fact, it was shown that in terms of the efficiency of the hydraulic fracturing treatment (hereinafter described as HFT) method, the more the horizontal wellbore deviated toward the σ Hmax direction (σ Hmax -dir), the lower the estimated ultimate recovery (EUR). Then, Yang et al. (2016) conducted a performance comparison study of transversely and longitudinally fractured horizontal wells over varied reservoir permeability [7]. They found that horizontal wells with transverse fractures are preferable over the entire permeability range of their study (50 nD-5 mD) in oil reservoirs. However, the productivity was predicted assuming that the fractures have the same physical properties and were placed either transversely or longitudinally. A hydraulic fracture induces a stress shadow (SS) in the surrounding rock; in turn, the increased stress causes a reduction in width for an adjacent HF that propagates within the stress shadow effect (SSE) region or causes the fracture path to change [8][9][10][11]. Dohmen et al. (2014) used downhole geophones to monitor microseismic events associated with hydraulic fracturing in horizontal wells with different completion stage spacings in the Middle Bakken Play and the Utica Play. They found that hydraulic fractures in closely spaced stages interfere with each other by cyclically bouncing out-of-zone. In addition, the authors attributed the bouncing phenomena to the SSE that accumulated stress forced subsequent fractures preferentially propagate upward, which could have a significant impact on the production-efficient contact area as fractures could be created outside the target production area [12]. Tan et al. (2019) deployed dual microseismic arrays to monitor fracture geometry in Vaca Muerta shale of North America and observed that due to the SSE among neighbor stages, the regional tectonic stress direction before hydraulic fracturing changed and the stress regimes varied between strike-slip and reverse faulting after the stimulation [13].
Several modeling approaches have been developed to simulate complex fracture networks [14][15][16][17][18]. Aldakheel et al. (2020) and Noii et al. (2020) applied a phase-field fracture method, evaluated as robust and efficient, that predicted the propagation behavior of fractures in anisotropic heterogeneous material and the corresponding pressure changes [19,20]. Meanwhile, Kresse et al. (2012) examined the SSE from the precedent treatment stage on the HF propagation and shape at next stages by using the Unconventional Fracture Model (UFM) model. According to the detailed simulation study, the authors indicated that the HF propagation is more dominated by the SSE when the principal stress anisotropy is smaller. In a formation with small stress anisotropy, fracture interaction can lead to dramatic divergence of the fractures as they tend to repel each other. However, they also pointed out that even when the stress anisotropy is large enough to limit the fracture turning magnitude, the stress shadowing still has a strong effect on the fracture geometry and proppant placement. [21].
The influence of the SSE generated by the fractures may vary as the hydraulic fracture is generated transversely or longitudinally with respect to the horizontal well. However, little work has been performed to understand the relationship between the principal stresses, the lateral direction of a horizontal well, and HFT efficiency. Therefore, we performed this study to achieve the following,

•
construction of the reservoir model in Brushy Canyon formation and reflect the in-situ principal stress, • investigation of the effects of σ Hmax -dir in the Permian Basin and drilling direction of horizontal wells on the hydraulic fracturing treatment efficiency, • analysis of the stress shadow influence according to the transverse or longitudinal fracture direction.

Geology and Geophysical Model of Forty Niner Ridge Field in Brushy Canyon Formation
The target field in this study is the Forty Niner Ridge field in the Permian Basin in New Mexico, U.S. The Brushy Canyon formation consists of fine-grained sandstone layers of 150 m which represent a significant reservoir interval. However, the low permeability and petrophysical heterogeneity of the formation rock limits primary recovery to as little as 10-15% of the original oil in place (OOIP), suggesting that horizontal wells should be drilled and hydraulically fractured [22].
In the target field, four horizontal wells were drilled in four different directions from 2006 to 2014. To perform this study, the G&G (Geology and Geophysics) reservoir model was first constructed based on well logging data from 17 wells, completion reports, and production reports. First, logging data interpretation was performed by using Schlumberger's Petrel and Techlog simulators to construct the reservoir and to calculate petrophysical properties. As shown in Figure 1, the formation layering was conducted by interpreting the gamma ray patterns from measured physical log data. After determining the existence of bad holes from the caliper log and bit size using Techlog, shale volume, effective porosity, and total water saturation were calculated by synthesizing bad hole, gamma ray, neutron and bulk density, and resistivity logs. Then, based on the petrophysical properties, geomechanical properties of the reservoir were calculated by applying the following correlation equations. From the core analysis results of the same Brushy Canyon formation 6 km west from the target area, a regression analysis was performed to determine the porosity/ permeability relationship as Equation (1) [22].
where ∅ is porosity. From the field technical progress document, the pore pressure gradient was reported as 9.856 MPa/km [22]. Due to the absence of the field measurement data for the rock strength and elastic properties, correlations based on the similar rock mineralogy were adopted to calculate geomechanical properties such as the uniaxial compressive strength (UCS) and the Young's modulus [23]. UCS correlation (2) derived from multivari-  [24,25].
where ρ is density (g/cm 3 ) and S w is water saturation. The σ Hmax /σ hmin acting on our reservoir model was reflected in the model based on the research results published by Snee and Zoback (2018) [26]. They provided regional faulting regime value (A ∅ ) proposed by Simpson (1997) to represent the relative magnitudes of the three principal stresses (σ v , σ Hmax , and σ hmin ) [27]. The value ranges from 0 (the most extensional possible condition of radial normal faulting) to 3 (the most compressive possible condition of radial reverse faulting). The value is defined mathematically by Equations (4) and (5) and where, σ 1 , σ 2 , and σ 3 are the magnitudes of the maximum, intermediate, and minimum principal stresses, respectively, and n is 0 for normal faulting, 1 for strike-slip faulting, and 2 for reverse faulting. According to the report, in the region to which our target field belongs, normal or strike-slip faults mainly occur, in which the faulting regime value is 0.6, and the σ Hmax -dir is N075 • E. Thus, the relationship between the σ v /σ Hmax /σ hmin was summarized as Equations (6) and (7).
where ϑ represents the Poisson's ratio, α is the Biot's constant which is the dimensionless value, and P p is pore pressure (MPa). The reservoir G&G model constructed by reflecting the above properties is shown in Figure 2. From the horizontal well trajectory information, well logging interpretation, and production reports, the true vertical depth (TVD) ranges from 2316-2377 m in depth, with a σ v gradient of 24.9 MPa/km, and the σ Hmax /σ hmin anisotropy of the field is calculated to be about 13.8 MPa (Table 1).

Hydraulic Fracturing Application Case Study
The horizontal wells in the Brushy Canyon formation of our study were drilled in four orthogonal directions regardless of in-situ stress conditions. While two horizontal wells, SANDY01 and Roadrunner Fed 2H, were drilled parallel to σ hmin -dir, the other two, SANDY03 and Roadrunner Fed 1H, were drilled perpendicular to σ hmin -dir. We implemented this field as "Actual Well Pattern" model as depicted in Figure 3a. On the other hand, the "Proposed Well Pattern" model drilling horizontal wells in consideration of the σ Hmax -dir was implemented as shown in Figure 3b. In both cases, the total number of fracturing stages was 45 and the total length of lateral parts was set to be 6400 m. The completion of the well was performed based on the completion report at the site, and Figure 4 is the completion design of SANDY03. There are 10 fracturing stages and the spacing between stages is designed to range from 76.5 to 213.4 m. In the case of the proposed model, each horizontal well was designed to have 15 fracturing stages with the same 137.2 m spacing. Then, the pumping design to generate HFs is as shown in Table 2. First, in stage #1, the breakdown linear gel was injected to break the formation, followed by an acid treatment, and the crosslinked gel and proppants were mixed and injected. Then, flushing was followed. This process has become simpler, as the acid treatment and flushing process are omitted in the later stages as shown in Table 2 stage #3.

Hydraulic Fracturing Model
In this study, Schlumberger's Kinetix simulator was used for hydraulic fracturing modeling. The "Unconventional Fracture Model," which is a complex fracture network model that creates planar fractures with complex stress profiles, was applied as the hydraulic fracture model.

Unconventional Fracture Model (UFM)
The fracturing model applied in this study is UFM, which simulates vertical fracture propagation, rock deformation, and fluid flow in the complex fracture network created during a fracture treatment. The basic governing equation describes the fluid flow in the fracture network, the fracture deformation, and the fracture propagation/interaction criterion [21]. The Equations (8) and (9) are the continuity equations assuming that fluid flow propagates along fracture network with mass conservation, ∂q ∂s where q is the local flow rate inside the hydraulic fracture along the length, w is an average width, H f l is the local height of the fracture occupied by fluid, and q L is the leak-off volume rate through the wall of the hydraulic fracture into the rock matrix per unit length. Through the Carter's leak-off model, q L is yielded by h L , leak-off height, times velocity, u L , at which fracturing fluid infiltrates into a surrounding permeable medium. For a vertical fracture ( Figure 5) in a layered formation with variable σ hmin and fluid pressure p, the width profile can be determined from an analytical solution given as and the varying parameter, height of the fracture h, is calculated based on the approach described in Kresse et al. (2011). The fracture height is calculated based on the pressure at each location of the fracture by matching "stress intensity factors, K Iu and K Il " as in Equations (11) and (12).
where K Iu and K Il are the stress intensity factors at the top and bottom tips of the fracture, σ i and h i are the minimum stress and distance from top of the ith layer to fracture bottom tip, p cp is the fluid pressure at a reference (perforation) depth, h cp , measured from the bottom tip. Then fracture width w(z) at any position z measured from the bottom tip is given by Equation (13) [21,28,29].

Stress Shadow Effect (SSE)
Stress interaction occurs between HFs when they propagate in proximity, known as the "stress shadow effect." An HF induces a SS in the surrounding rock, proportional to the net pressure inside the fracture, in a region of approximately the size of the fracture height. Then, the increased stress causes a reduction in width for an adjacent HF that propagates within the SSE region or causes the fracture path to change. The stress field around a 2D fracture with internal pressure p can be calculated from Sneddon (1946) and Sneddon and Elliott (1946) solutions. The stress normal to the fracture is σ x (Figure 6) and can be calculated from Equation (14) [30,31]. And x, y are the coordinates and x, y, L, L 1 , L 2 are the distances in Figure 6.

Analysis of HFT Results
To analyze the effect of the in-situ principal stress condition and horizontal well direction on the hydraulic fracturing efficiency, we investigated the hydraulic fracturing results in the "Actual Well Pattern" model and the "Proposed Well Pattern" model. We analyzed the fracture characteristics such as the growth direction of HFs, the size, width, and conductivity of individual fractures, as well as the stress shadows generated in the formation adjacent to these HFs. In addition, by calculating the SRV according to the HF generation, the influence of the direction of the formation principal stress and the azimuth of horizontal wells on the HFT efficiency was analyzed.

Hydraulic Fracture Geometry
First, as a result of performing HFT on the actual and proposed models, all of the HFs were facing a certain direction, as shown in Figure 7a,b, although the horizontal wells were drilled in four different directions. The σ Hmax -dir of the reservoir model is N075 • E, and all HFs generated in both the actual and the proposed models have grown toward σ Hmax -dir. In other words, in actual case, transverse (T-dir) fractures were generated in the horizontal wells drilled parallel to σ hmin -dir (SANDY01, ROADRUNNER FED 2H), and in the horizontal wells drilled perpendicular to σ hmin -dir (SANDY03, ROADRUNNER FED 1H), the HFs formed obliquely and almost longitudinally (L-dir) in the well. Meanwhile, in "Proposed" case, all three horizontal wells were drilled in the σ hmin -dir, so that all HFs were formed in the T-dir.

Fracturing Efficiency
This time, the loss of the injected fracturing fluid was examined to inspect the fracturing efficiency by fracturing fluid injection. The fracturing fluid injected to generate HFs in the formation can be lost by being pushed out through the fracture surface into the contacting formation due to high internal pressure. Since the fracturing fluid as much as the lost amount cannot be used for fracturing, it is considered a factor affecting the fracturing efficiency. The continuity equation assuming that fluid flow propagates long fracture of the UFM model includes q L term, which is the leak-off volume rate as described in Equation (8). The q L is yielded by h L , leak-off height, times velocity, u L , at which fracturing fluid infiltrates into surrounding permeable medium as in the Equation (9). The permeability of formation may affect this velocity, u L . For this reason, the average porosity and permeability of the adjacent well formation and the leak-off ratio determined in each model were depicted in Figure 8. In the actual model, the average porosity and permeability were 7.7% and 3.9 × 10 −15 m 2 , which are 5.5% and 88.7% higher in comparison to that of the proposed model. The leak-off ratio, which indicate the ratio of lost fluid to the amount of injected fluid, was calculated to be 11.1% higher. Although the porosity of the two models is similar, the difference in permeability of 88.7% might have affected the calculation of leak-off ratio, which represents the loss of fracturing fluid to the surrounding formation. Under this influence, the HFs generation results of both models were examined. By investigating the results of the length, height, average width, and conductivity of the HFs created per unit injection amount of the fracturing fluid, a detailed analysis was performed to identify how the difference in the horizontal well orientation of the two models had influenced the fracturing efficiency. The HF volume of the actual model was about 3% larger, the HF width of the proposed model was slightly thicker, and the conductivity was almost the same, so the difference between the results was insignificant as listed in the Table 3. In other words, the difference between the two models is insignificant in terms of the HF length, height, width, and conductivity, so it was found that the azimuth of the horizontal well had little effect on the fracturing efficiency.

Stress Shadow Influenced by Fracture Propagating Direction
This time, the characteristics of the HFs generated according to the drilling direction of the horizontal well were examined. As shown in Figure 7, in the horizontal well drilled toward σ hmin , HFs grew toward σ Hmax ; thus, fractures placed transverse (T-dir) to the well were created. When a fracture grows and the width increases, a SS is generated in the adjacent rock. These stresses locally increase the total stress and, as they are aggregated to σ hmin , the change of σ Hmax /σ hmin anisotropy may induce the rotation of the fracture propagating direction. In this model, the SS, which was less than 0.3 MPa generated by T-dir fractures, was weighted on the minimum horizontal principal stress. However, the maximum principal stress in the formation was on average 13.8 MPa higher than the minimum, so it did not affect the propagating direction of the fractures. As a result, the facing distance between HFs remained the same as the distance between the initially designed fracture stages. On the other hand, the HFs of the horizontal well drilled toward σ Hmax in Figure 9 grew almost longitudinally as the fracture deviated obliquely from the horizontal well. Accordingly, the initially designed spacing between the stages was 76.5 m to 213.4 m, but the distance that the actual HFs faced decreased to 26%, resulting in closer spacing of 19.8 m to 55.2 m ( Table 4). As shown in the right side of Figure 9, the 2nd stage HF region, which locally faces the maximum fracture width (~2.2 cm) of the 1st stage HF, was stressed by a maximum of 6.4 MPa due to the SS. Consequently, growth of the 2nd stage HF was biased toward the non-SS area as shown in Figure 10, resulting in a relatively narrow proppant distribution in the fracture. Moreover, the width of the 2nd stage fracture facing the 1st stage fracture of yellow circle of SS area in Figure 10 was 0.03 cm or less, and the conductivity converged to 0 m 3 , such that the weighted SS eventually caused the HFs to be generated only on one side of the horizontal well.

Stimulated Reservoir Volume (SRV)
We have calculated the difference between the above T-dir and L-dir HFs induced stimulated reservoir volume (SRV) for each "Well Pattern" model. The concept of SRV provides a link between well productivity, the observed distribution of microseismic events, and HF networks that are created. The SRV dimensions are affected by both reservoir geomechanical properties and hydraulic fracturing parameters [32]. In this study, using the "Kinetix Stimulation" simulator, as shown in Figure 11, the SRV is calculated by creating unstructured grids centered on the HFs and summing the volume of the grids at a certain distance from the HFs [33]. The unstructured gridding created a Voronoi (PEBI; perpendicular bisection) grid which can conform to complex HF geometries generated by UFM [34]. The fracture cell width was upscaled to 3.28 ft. The minimum zone height was 10 ft and layers smaller than 10 ft were merged. As shown in Figure 12a,b and Table 5, the total SRV result of the proposed case, which was designed by considering the in-situ principal stress, was 99,926,737 m 3 , which was 16.4% higher than the actual case (116,291,809 m 3 ). This is due to the L-dir HFs, which were very closely spaced, formed in SANDY03 and ROADRUNNER 1H of the actual model drilled toward σ Hmin -dir. Additionally, it is because the fractures created almost parallel to the horizontal well were limited to a narrow range close to the wellbore, so that the contact area to the reservoir could not be formed widely. Figure 13 shows the monthly oil production and cumulative oil production of SANDY03 from 2012 to 2019. The monthly oil production for the 1st year indicates that it decreased to 16.2%, from the initial 1848 m 3 to 300 m 3 . In addition, the cumulative production is only 2.7% (8777.2 m 3 ) of the calculated OOIP. If the T-dir HFs were formed in SANDY03 and the SRV was 16.4% higher, it is expected that 1437.7 m 3 of oil can be additionally produced when the oil recovery of 2.7% is incorporated. Figure 11. Concept diagram of unstructured gridding for stimulated reservoir volume (SRV) calculation [33]. Figure 12. SRV of (a) "Actual Well Pattern" case and (b) "Proposed Well Pattern" case.  Figure 13. Monthly oil production and cumulative oil production history of SANDY03 well for 7 years (October 2012-September 2019).

Sensitivity Analysis of the Principal Stress Anisotropy
As mentioned earlier, this reservoir model has a very high σ Hmax /σ hmin anisotropy of about 13.8 MPa, and no rotation occurred in the HFs growth direction. However, in the case of a formation with anisotropy of 2.1 MPa, the fracture propagating direction may be affected depending on the magnitude of SSE. To confirm the effect, HFT was performed again in the reservoir with a lower σ Hmax /σ hmin anisotropy of 2.1 MPa for the SANDY03 well where heavy SSE occurred. In the results, as the fracturing stage continued to 2nd, 3rd, and 4th, the change in the propagating direction was confirmed, as shown in Figure 14. As the total stress in the vicinity of the 2nd stage HF increased due to SSE induced by the 1st stage HF, the propagating direction of the 2nd HF was rotated toward the higher stressed zone (yellow arrow to blue arrow in Figure 14); consequently, the HF was formed closer to the horizontal well. It is reported that for fractures not transverse to a horizontal wellbore, the stress interaction between the adjacent fractures may attract and even connect each other [35][36][37]. In other words, depending on the state of principal stresses, HFT could create HFs grown in a direction that is unfavorable for the well productivity, which is L-dir in this study. In particular, when the principal stress anisotropy is relatively small, the HF propagating direction may change due to the occurrence of SSE, and the disadvantage to production may be more severe. Then, the influence of the horizontal well drilling azimuth determination by considering the in-situ stress condition on the efficiency of the HFT becomes more critical.

Conclusions
The main purpose of this study is to investigate the effect of σ Hmax -dir in the Permian Basin on the efficiency of HFT by using two different drilling cases: "Actual Well Pattern" vs. "Proposed Well Pattern". The following conclusions were made.
In the simulation results of actual case, transverse (T-dir) fractures, which are perpendicular to the σ hmin -dir, were generated in the horizontal wells drilled parallel to σ hmin -dir. In the meantime, two other horizontal wells were drilled perpendicular to σ hmin -dir, which is almost a longitudinal fracture in the well. These obliquely deviated fractures from the well significantly decreased the fracture spacing in between up to 26%, so that the SSE between fractures was greatly increased. With the large stress anisotropy of the principal stresses in Brushy Canyon field, the fracture propagating direction did not change. However, in the case of a formation with moderate anisotropy, SS induced the change in the fracture propagating direction. As total stress in the vicinity of the prior stage fracture increased due to SSE, the propagating direction of the following fracture rotated toward the higher stressed zone even closer to the wellbore. In the results of stimulated reservoir volume (SRV) and original oil in place (OOIP), it was found that the direction of fracture propagation greatly influenced the SRV and OOIP. The proposed case with T-dir fractures was 16.4% and 13.2% higher than the actual case with both T-dir and L-dir fractures, respectively. The fractures created parallel to the well were limited to a narrow range close to the wellbore, resulting in the small contact area to the reservoir. In conclusion, to apply hydraulic fracturing effectively, the horizontal well direction must be determined in consideration of the principal stress direction.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.