A method for predicting the probability of formation of complex hydraulic fracture networks in shale reservoirs: development and application

Abstract

Shales can form a complex fracture network during hydraulic fracturing, which greatly increases the stimulated reservoir volume (SRV) and thus significantly increases oil or gas production. It is therefore important to accurately predict the probability of formation of the hydraulic fracture network for shale gas exploration and exploitation. Conventional discriminant criteria are presented as the relationship curves of stress difference vs. intersection angle. However, these methods are inadequate for application in the field. In this study, an effective and quantitative prediction method relating to the probability of complex fracture network formation is proposed. First, a discriminant criterion of fracture network was derived. Secondly, Monte Carlo simulation was applied to calculate the probability of the formation of the complex fracture network. Then, the method was validated by applying it to individual wells of two active shale gas blocks in the Sichuan Basin, China. Results show that the probabilities of fracture network are 0.98 for well JY1 and 0.26 for well W204, which is consistent with the micro-seismic hydraulic fracturing monitoring and actual gas production. Finally, the method was further extended to apply for the regional scale of the Sichuan Basin, where the general probabilities of fracture network formation are 0.32−1 and 0.74−1 for Weiyuan and Jiaoshiba blocks, respectively. The Jiaoshiba block has, therefore, an overall higher probability for formation of fracture network than the Weiyuan block. The proposed method has the potential in further application to evaluation and prediction of hydraulic fracturing operations in shale reservoirs.

Graphical abstract

Keywords

shale gas / complex fracture network / shale reservoir / Monte Carlo simulation / Sichuan Basin

1 Introduction

In recent years, the exploitation of shale gas has attracted widespread interest among the international energy community (Ross and Bustin, 2007; Curtis et al., 2012; Gasparik et al., 2014; Zhao et al., 2018). Compared with conventional gas reservoirs such as sandstones, shale reservoirs have ultra-low permeability, and thus hydraulic fracturing is necessary to make economic gas production (Dahi-Taleghani and Olson., 2011; Xu et al., 2015; Yuan et al., 2015; Zhao et al., 2016; Hou et al., 2018). Unlike the bi-wing and symmetric fractures found in conventional reservoirs, hydraulic fractures required in shale reservoirs commonly form complex networks (i.e., connected fractures of irregular orientation and distribution, see Fig.1), with large stimulated reservoir volumes (SRV). These networks can greatly increase the matrix permeability and resulting gas production (Page and Miskimins, 2009; Shapiro, 2009; Hou et al., 2019). Therefore, it is important to understand how complex hydraulic fracture networks form in shale gas reservoirs.

At present, it is widely accepted that the complex hydraulic fracture network in the naturally fractured formations is mainly caused by the interaction between hydraulic and natural fractures (Renshaw and Pollard, 1995; Jeffrey et al., 2009; Weng et al., 2011; Gu et al., 2012; Liu et al., 2018c). In other words, when a hydraulic fracture approaches a pre-existing natural fracture, the propagation of the hydraulic fracture is not always in the direction of the maximum principal stress, because of the interference of natural fractures. If a hydraulic fracture crosses a natural fracture, the hydraulic fracture may simply continue to propagate along its original direction. However, if a hydraulic fracture does not cross a natural fracture but propagates along it, the natural fracture may be reactivated and a complex fracture network may form. There is considerable research on whether hydraulic fractures cross, or propagate along natural fractures. For example, Zhou et al. (2008, 2010) concluded that the factors controlling hydraulic fractures propagating along natural fractures are the horizontal stress difference, the intersection angle between the hydraulic fracture and the natural fracture, the internal friction coefficient of the natural fracture, and the net pressure of the fracture. Weng et al. (2011) suggested that for the same initial natural fractures, the stress anisotropy can result in a change in the induced-fracture geometry from a wing-fracture to a complex-fracture network. Dehghan et al. (2015) concluded that the dip and strike of pre-existing fractures have a substantial influence on the hydraulic fracture propagation behavior and geometry. Fatahi et al. (2017) found that an increase in the angle between the natural fracture plane and the direction of maximum horizontal principal stress can improve the possibility of a hydraulic fracture crossing a natural fracture. Liu et al. (2018a) found that orthogonal intersection of natural fractures is favorable for creating fracture networks. These studies show that multiple factors influence the interaction between hydraulic fractures and natural fractures, which makes it difficult to understand the formation of complex fracture networks in shale reservoirs.

Many methods based on numerical models and laboratory experiments have been proposed for predicting whether a hydraulic fracture will propagate across a natural fracture. For example, Blanton (1982, 1986) proposed a method based on the criterion considering stress difference and intersection angle. It was found that hydraulic fractures tend to cross pre-existing fractures when the intersection angle is 60−90 degrees and the stress difference is high enough. Renshaw and Pollard (1995) proposed a criterion considering crossing stress ratio and friction coefficient to determine whether a hydraulic fracture will cross a pre-existing fracture orthogonal to it, and they concluded that the crossing will occur when the tip stress of the hydraulic fracture is sufficiently high. However, natural fractures and hydraulic fractures are not always orthogonal. Gu et al. (2012) extended the criterion for non-orthogonal intersection angles by theoretical derivation, and they suggested that the hydraulic fracture is more likely to propagate along the natural fracture as the intersection angle between the hydraulic fracture and natural fracture decreases from 90°. Sarmadivaleh and Rasouli (2014) further expanded Renshaw and Pollard’s criterion to a case for a non-cohesive interface with a non-orthogonal intersection angle. In general, the above criteria were derived by analyzing the forces in the fracture planes, when several factors such as stress difference, intersection angle, and friction coefficient were considered. As far as the application is concerned, these criteria or methods are generally presented as a relationship between two parameters such as stress difference and intersection angle. This relationship can be used to judge whether hydraulic fractures will propagate along natural fractures to form fracture networks. However, due to the complexity of the relationship, these criteria or methods are hard to apply in the field. Thus, it is worthwhile to explore an effective and applicable method for predicting the probability of the formation of complex fracture networks during hydraulic fracturing of shale reservoirs.

In this study, we propose an improved method to predict the probability of formation of complex fracture networks, which utilizes the development of a discriminant criterion and Monte Carlo simulation. First, the discriminant criterion was derived by introducing a normal stress calculation formula into the derivation process, where the triaxial principal stress, the natural fracture plane strike and dip angle are considered. Secondly, the Monte Carlo simulation is applied. Through a large number of simulations, the probability of formation of complex fracture networks could be calculated, which makes it possible to apply the method to actual wells and gas fields. Then, the method was validated by applying it to individual wells of two active shale gas blocks in China. Finally, the application of this method was further extended to the Weiyuan and Jiaoshiba Blocks in the Sichuan Basin, China to explore the formation conditions of complex fracture networks.

2 Methodology

The methodology for predicting the probability of formation of a complex hydraulic fracture network consists of two main parts: 1) establishing the discriminant criterion for the formation of the complex fracture network, and 2) conducting a Monte-Carlo simulation.

2.1 Discriminant criterion for the formation of the fracture network

The formation of the fracture network is mainly controlled by geological and engineering factors. In this study, we assume that all dynamic engineering parameters, such as fluid viscosity, flow rate, number of perforation clusters, and staged fracturing are the same. Fig.2 shows the interaction between a hydraulic fracture (HF) and a natural fracture (NF) viewed from above. The whole process of interaction between the HF and the NF can be divided into multiple modes, such as approaching, slippage, being captured, crossing, and deflecting (Gu et al., 2012). This study mainly focuses on the behavior at the moment when the HF and NF meet, and the results of the instantaneous interaction can be divided into crossing and deflecting. In other words, the HF is parallel to the maximum principal stress (σ_max) and perpendicular to the minimum principal stress (σ_min) before it encounters the NF, and when the HF encounters the NF, the HF may cross the NF with its original direction or deflect along the NF (Fig.2). The effective stress on the NF can be decomposed into the normal stress (σ_n) perpendicular to the fracture plane and the shear stress (τ_n) parallel to the fracture plane (Fig.2).

The HF is the fracture of tensile nature and it always develops in the direction perpendicular to the minimum compressive stress. When the net pressure (p_net) in the fracture at the time of interaction exceeds the minimum principal stress (σ_min) and the tensile strength (S_t) of the rock (Mokhtari et al., 2014), the HF will occur at the direction perpendicular to the minimum principal stress (σ_min), as given in Eq. (1):

(1)p_{\mathrm{n}\mathrm{e}\mathrm{t}}>{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}+S_{\mathrm{t}}.

The NF is a weak plane and its S_t is generally considered to be 0. When the p_net in the fracture at the time of interaction is greater than the σ_n acting on the NF, the NF will be open and reactivate (Blanton, 1982; Blanton, 1986), as given in Eq. (2):

(2)p_{\mathrm{n}\mathrm{e}\mathrm{t}}>{\sigma }_{\mathrm{n}}.

If the HF encounters the NF, the criterion for choosing the direction at the moment of contact is the principle of least resistance (Wang, 2019). Therefore, when the σ_n on the NF is less than the sum of the σ_min and the S_t, the HF deflects along the NF, and a complex fracture network is formed, as given in Eq. (3):

(3)

where σ_n is the normal stress on the NF, MPa; σ_min is the minimum principal stress, MPa; S_t is the rock tensile strength, MPa.

On the contrary, when the σ_n on the NF is greater than the sum of the σ_min and the S_t, the HF crosses the NF and propagates perpendicular to the direction of σ_min.

For a natural fracture plane with a certain triaxial principal stress, its σ_n can be calculated by Eq. (4) (Tong et al., 2010):

(4){\sigma }_{\mathrm{n}}={\sigma }_{\mathrm{v}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2\theta +{\sigma }_{\mathrm{H}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2\theta {\mathrm{c}\mathrm{o}\mathrm{s}}^2\alpha +{\sigma }_{\mathrm{h}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2\theta {\mathrm{s}\mathrm{i}\mathrm{n}}^2\alpha ,

where σ_v is the vertical principal stress; σ_H is the maximum horizontal principal stress; σ_h is the minimum horizontal principal stress; θ is the angle between the NF and the σ_v, that is, the residual angle of the dip angle (β) for the NF; α is the angle between the intersection line of the NF on the σ_H–σ_h plane and the σ_h, which is the angle between the strike of the NF and the σ_h (see Fig.3).

It is worth noting that, since the direction of σ_h is constant, the angle α actually reflects the strike of the NF. Therefore, α can be defined as the relative strike to σ_h of the NF. To simplify the expression, α is referred to as the strike of the NF hereafter.

As shown in Fig.3, there is a relationship of θ + β = 90°, then Eq. (4) can be changed to

(5)\begin{array}{rl}{\sigma }_{\mathrm{n}}=& {\sigma }_{\mathrm{v}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2(90-\beta )+{\sigma }_{\mathrm{H}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2(90-\beta ){\mathrm{c}\mathrm{o}\mathrm{s}}^2\alpha \\ & +{\sigma }_{\mathrm{h}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2(90-\beta ){\mathrm{s}\mathrm{i}\mathrm{n}}^2\alpha .\end{array}

Substituting Eq. (5) into Eq. (3), the discriminant criterion for the complex fracture network can be obtained as:

(6)

If Eq. (6) is true, a complex fracture network can be formed. Once we obtain the values of the σ_H, σ_v, σ_h, S_t, α, and β in Eq. (6), then we can judge whether a complex fracture network system can be formed, or not.

2.2 Monte–Carlo simulation

Monte-Carlo simulation is a broad class of computational algorithms that rely on repeated random sampling and statistical analysis to obtain numerical results. In principle, Monte–Carlo simulation can be used to solve any problem having a probabilistic interpretation. The main advantage of this method is that the error of the method is not related to the scale of the problem, and problems with statistical properties can be solved directly. In the field of oil exploration, it is mainly applied to investment risk analysis and reserves calculation (Stoltz et al., 1998; Nam et al., 2013). In this study, we introduce Monte–Carlo simulation for hydraulic fracture network prediction.

The use of Monte-Carlo simulation to predict the probability of formation of a complex hydraulic fracture network involves four main steps.

The first step is to construct a mathematical probability model (Y) whose numerical solution is the probability distribution of the actual problem. Y is defined as

(7)Y=F\left(X_1,X_2,\dots ,X_i,\dots ,X_n\right),

where X₁, X₂,…, X_i,…, X_n represents n independent random variables whose probability distributions are Q₁(X₁), Q₂(X₂),…, Q_i(X_i),…, Q_n(X_n), respectively.

In Eq. (7), X_i is randomly sampled according to the probability distribution Q_i (X_i), and a series of target values can be obtained. When the sampling numbers are large enough, the probability distribution of system values can be obtained. The values reflect the numerical solution characteristics of the problem to be solved.

In actual shales, the σ_H, σ_v, σ_h, S_t, α, and β in Eq. (6) have random characteristics within a certain range. As a result, the prediction results are governed by statistical rules with a certain probability significance. Therefore, the Monte-Carlo method can be used, and the left side of Eq. (6) can be used for a mathematical probability model as follows:

(8){Y=\sigma }_{\mathrm{v}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2\beta +{\sigma }_{\mathrm{H}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2\beta {\mathrm{c}\mathrm{o}\mathrm{s}}^2\alpha +{\sigma }_{\mathrm{h}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2\beta {\mathrm{s}\mathrm{i}\mathrm{n}}^2\alpha -{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}-S_{\mathrm{t}.}

In Eq. (8), σ_v, σ_H, σ_h, S_t, α, and β are independent random variables.

The second step is to set the value range and mathematical distribution type of each variable. The value ranges of variables can be chosen according to the actual data in the study area. There are three commonly-used distribution types, i.e., the normal, uniform, and BetaPERT distributions. The BetaPERT distribution describes a situation where the minimum, maximum, and most likely values are known, and it is similar to the triangular distribution, except the curve of probability vs random variable is smoothed to reduce the peak.

To determine the most appropriate distribution type for each variable, a Chi-square test is performed. The Chi-square test is a widely-used function that tests the deviation degree between the actual observed values and the theoretical inferred values of a sample (Bryant and Satorra, 2012). The Chi-square test statistic formula can be expressed by Eq. (9):

(9){\chi }^2={\sum }_{i=1}^k\frac{{{(f}_{\mathrm{o}i}-f_{\mathrm{e}i})}^2}{f_{\mathrm{e}i}},

where χ² is Chi-square test value; f_oi is the observation frequency of the ith group; f_ei is the expected frequency of the ith group; k is the number of data groups. In the chi-square test, smaller χ² means smaller fitting error and better fitting result.

The third step is to generate random values of the variables based on the distribution types, and to calculate the Y value in Eq. (8). Then, the sampling and calculation steps are repeated until the number of samplings is sufficient. Generally speaking, ≥ 6000 samplings are needed to achieve a statistically valid outcome.

The fourth step is to read the probability that it is less than 0 from the generated cumulative probability distribution curve of the Y value, which is the probability of the fracture network formation.

Based on the above four steps of Monte-Carlo simulation, a prediction method for formation probability of fracture network can be obtained and then applied in the field.

3 Case studies for the Lower Longmaxi Formation in individual wells in the Sichuan Basin

Wells W204 in the Weiyuan block and JY1 in the Jiaoshiba block were selected as the study locations to validate the proposed method.

3.1 Geological background

The Sichuan Basin is located in the northwest of the Upper Yangtze Craton in southwest China (Fig.4), with an area of about 16 × 10⁵ km² (Zheng et al., 2019). As shown in Fig.4, the Sichuan Basin extends from Micang Mountain (fold belt) and Daba Mountain (thrust belt) in the north to Emei-Leshan (fold belt) in the south, and from Longmen Mountain (orogenic belt) in the west to the Hubei-Hunan-Guizhou fold belt in the east. Tectonically, the Sichuan Basin is divided into three tectonic regions by the Longquan Mountain and Huaying Mountain, and it is further subdivided geographically into six sublevel structural domains: west, north, central, south-east, east, and south (Fig.4). The study areas were the Weiyuan block in the south-west domain of the Sichuan Basin, and the Jiaoshiba block in the east domain of the Sichuan Basin (Fig.4). Two individual wells studied, JY1 and W204, are located in the northern Jiaoshiba block, and the central Weiyuan block, respectively (Fig.4).

Pre-Sinian metamorphic rocks and magmatic rocks form the basement in the Sichuan Basin, and a series of Sinian - Cretaceous strata are developed from bottom to top. The Upper Ordovician Wufeng Formation and the Lower Silurian Longmaxi Formation are two major marine sedimentary strata. The depositional settings of the Wufeng and Longmaxi Formations in the two blocks were mainly deep-water shelf to shallow-water shelf (Xu et al., 2020). The Wufeng Formation underlies the Longmaxi Formation, and comprises a set of thin siliceous shales. The Longmaxi Formation, characterized by black shale and abundant graptolites, is widely distributed in the study area. The Lower Member of the Longmaxi Formation with has high contents of organic matter and graptolites is the focus of this article, and it is the target unit for shale gas development in the Weiyuan and Jiaoshiba blocks (Chen et al., 2015).

3.2 Data acquisition

To quantify the probability of formation of a complex fracture network for the Lower Longmaxi Formation in wells W204 and JY1, the variables in Eq. (8) need to be determined. In this study, the σ_H, σ_v, and σ_h of well W204 were obtained by logging interpretation in the ranges of 85−88.3 MPa, 83−85.2 MPa, and 69.6−70 MPa (Zhang, 2015; Ma et al., 2020); those of well JY1 were obtained from the Kaiser rock mechanics experiments, and they were 52.19−55.52 MPa, 49.25−53.68 MPa, and 48.63−49.92 MPa (Guo, 2014). The strike and dip angles of NF were obtained by image logging interpretations. The strike of NF in Well W204 is 40°−150° and that in well JY1 is 50°−130°; the dip angles of NF are 20°−65° and 15°−75° for wells W204 and JY1, respectively (Li, 2014). Note that the strike mentioned here is the true strike of the natural fracture, which needs to be converted to the relative strike (α) to the σ_h for application. In this study, we failed to collect the data for S_t for wells W204 and JY1 due to the lack of systematic study of the tensile failure of shale formations in the Weiyuan and Jiaoshiba blocks, so as an alternative we used ranges of values of S_t for different cores from the Sichuan Basin (see Tab.1) for simulation. All the results of the S_t were obtained from Brazilian tests on shale samples.

3.3 Simulation and results

3.3.1 Determination of the mathematical probability model

The three-direction principal stress of Lower Longmaxi Formation in the two wells W204 and JY1 satisfies the relationship of σ_H > σ_v > σ_h, belonging to the strike-slip stress regime (Guo, 2014; Zhang, 2015). According to Eq. (6), the σ_min equals σ_h, thus, the fracture network discriminant criterion can be further modified as

(10)

The mathematical probability model of fracture network formation (Eq. (8)) can be further modified as

(11){Y=\sigma }_{\mathrm{v}}{\mathrm{c}\mathrm{o}\mathrm{s}}^2\beta +{\sigma }_{\mathrm{H}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2\beta {\mathrm{c}\mathrm{o}\mathrm{s}}^2\alpha +{\sigma }_{\mathrm{h}}{\mathrm{s}\mathrm{i}\mathrm{n}}^2\beta {\mathrm{s}\mathrm{i}\mathrm{n}}^2\alpha -{\sigma }_{\mathrm{h}}-S_{\mathrm{t}}.

3.3.2 Determination of the variables in the probability model

The range of values and mathematical distribution type of each variable in the probability model are presented in this section. The range of values was determined from the original data in Section 3.2, and the distribution type was determined by the Chi-square test. Chi-square tests were carried out for S_t, σ_H, σ_h, σ_v, α, and β, using actual values of each variable for the Lower Longmaxi Formation in the two wells, or the Sichuan Basin for S_t using data in Tab.1. Tab.2 shows an example for a Chi-square test of the values for S_t. The results show that a uniform distribution has the smallest χ² of 13.78 among the various distribution types, and therefore this was selected to represent the distribution type for S_t.

After Chi-square tests for all variables, the final simulation input parameters for wells W204 and JY1 are summarized in Tab.3.

3.3.3 Simulation results

In this study, we conducted 6000 simulations for each variable using the Crystal Ball software. The simulation process in the software consists of generating random numbers and calculating the Y value in Eq. (11). The simulation results are shown in Fig.5. We can read the probability of Y < 0 from Fig.5, which is the probability of formation of a complex fracture network.

As shown in Fig.5, the probability of Y < 0 is 0.26 for well W204, while it is 0.98 for Well JY1. This means that the probability of complex fracture network formation for the Lower Longmaxi Formation in well JY1 is significantly higher than that in Well W204. We compared our simulated results with the hydro-fracturing results that were obtained by 3D micro-seismic monitoring. As shown in Fig.6, well W204 has a simple hydraulic fracture morphology, which propagates along the direction perpendicular to the σ_h, while Well JY1 has complex hydraulic fracture morphology and network. Moreover, the initial gas production from W204 (16.5 × 10⁴ m³/d) was less than that of JY1 (20.3 × 10⁴ m³/d) (Jin et al. 2018). Thus, the simulation results for the two wells were validated by both the micro-seismic monitoring results and gas production rates.

3.4 Sensitivity analysis

Fig.7 shows the effect of varying different parameters on the probability of complex fracture network formation. For wells W204 and JY1, changes in S_t have the greatest impact on the probability values, and reach 76.7% and 95.3%, respectively. The effect of varying triaxial principal stress (σ_H, σ_v, and σ_h) to the results of complex fracture network formation are 0.7% and 2.2% for W204 and JY1, respectively. The impact of varying β and α for well W204 (12.1% and 10.5%, respectively) is one order of magnitude higher than those for Well JY (1.2% and 0.3%, respectively).

These results indicate that varying S_t is the main control on the simulation results. As a possible reason, the S_t is a critical parameter affecting the formation of the fracture network, which controls the initiation and propagation of hydraulic fractures (Gale et al., 2007; Li et al., 2017). The impact of varying β and α is higher for W204 than for JY1 (Fig.7). In the next section, the effect of changing both β and α on the fracture network formation probability is investigated through application of the method to the blocks as a whole.

4 Application at the regional scale

The method was also applied to the regional evaluation of the Weiyuan and Jiaoshiba blocks. The values of α and β can vary in different wells in the same region. In this study, since it is difficult to obtain α and β values for all wells in the entire region, we assumed that both α and β values are uniformly distributed in the range of 0°−90° in steps of 10°. In this case, for α and β in each block, there are in total 100 groups of values, such as (α = 0, β = 0), (α = 0, β = 10), (α = 0, β = 20), ......, (α = 90, β = 70), (α = 90, β = 80), (α = 90, β = 90). For σ_H, σ_v, σ_h, and S_t, the values are taken as random variables in a certain range. The ranges and distribution types of these random variables were determined in the same way as mentioned in the section 3.3.2 and not discussed in detail here. The input data for Monte-Carlo simulation, including the σ_H, σ_v, σ_h, and S_t are summarized in Tab.4.

For each given group of α and β, the other four variables are randomly sampled 6000 times to calculate the Y value. Fig.8 shows an example for the simulation result of α = β = 0° in Weiyuan block. As a result, the probability of fracture network formation (Y < 0) is 0.32 in Weiyuan block (Fig.8).

Other situations where α and β range from 0° to 90° with a step-length of 10° were also simulated, and in total 100 cumulative probability curves were obtained by the simulation software, and these results for the Weiyuan and Jiaoshiba blocks are summarized in Tab.5. As shown in Tab.5, the probability of complex fracture network formation in Weiyuan block is 0.28−1, and that in Jiaoshiba block is 0.74−1. The probability of complex fracture network formation in Jiaoshiba block is obviously higher than that in Weiyuan block.

Note that the conditions to be met for α and β when the probability is greater than 0.9 of two blocks are exhibited in Fig.9. It is clearly seen that in Weiyuan block, only when α is between 70° and 90°, and β is between 60° and 90°, it is possible to form a complex fracture network (Fig.9(a)). However, complex fracture networks may be formed in the range of 0°−90° for both α and β in Jiaoshiba block (Fig.9(b)). This implies that fracture networks can be formed more easily in Jiaoshiba than in Weiyuan during hydraulic fracturing. The gas production in Jiaoshiba block is higher than that in Weiyuan block (Jin et al., 2018), which also supports the simulated results.

It should be emphasized that with an increase in α and β, it is easier to form fracture networks in both blocks (Fig.9). In other words, the hydraulic fracture networks are concentrated in the areas with larger α and β. Therefore, when natural fractures with high dip angles that are almost parallel to the σ_H develop in the study area, they are more likely to be reactivated and connected to form a complex hydraulic fracture network.

Fig.10 shows discriminant curves proposed by different researchers (Blanton,1986; Gu et al., 2012; Liu et al., 2014; Tan et al., 2017; Dou et al., 2019). As shown in Fig.10, the HF is more likely to propagate along the NF to form fracture networks in the case of small intersection angles (< 60°) or small horizontal principal stress differences (Δσ < 6 MPa). This is consistent with the research results obtained by our method. In the Jiaoshiba block, the horizontal principal stress difference is small (Δσ = 2.98−6.89 MPa; data from Guo, 2014), and the HF can deflect along the NF in all directions (Fig.9(b)); in the Weiyuan block, the principal stress difference is relatively large (Δσ = 11−19.4 MPa; data from Zhang, 2015), only when α > 60° (intersection angle < 30°) can a fracture network be formed (Fig.9(a)). This further confirms the reliability of our simulation results.

As discussed, conventional methods use the discriminant curves that are presented as a function of stress difference vs. intersection angle in order to judge whether a fracture work can be formed (as shown in Fig.10). The only advantage of these conventional methods lies in their simplicity. For an ideal case with a certain horizontal principal stress difference and intersection angle, it is easy to judge from Fig.10 whether hydraulic fractures can deflect along natural fractures to form fracture networks. However, the relationships in actual wells and blocks are usually complicated, making these methods hard to apply in the field. Moreover, the conventional methods only consider the factors of stress difference and intersection angle, but not other factors such as S_t and β (Dehghan et al., 2015; Fatahi et al., 2017; Li et al., 2017; Liu et al., 2018b), which may bring uncertainty to the prediction. In contrast, the method proposed here introduces Monte-Carlo simulation to transform the fracture network prediction into a probability outcome, making it possible for application in actual wells and fields. Most importantly, the proposed method considers several additional aspects such as σ_H, σ_v, σ_h, S_t, α, and β. In other words, the conventional methods can be recognized as special cases in our method. Meanwhile, the new method has great potential for application in hydraulic fracturing efficiency evaluation by combining with the stimulated reservoir volume (Mayerhofer et al., 2010) and fracture complexity index (Cipolla et al., 2008) obtained through micro-seismic monitoring. This will be investigated in future research.

5 Conclusions

A method for the prediction of the probability of formation of complex fracture networks in shale reservoirs is proposed based on derivation of the discriminant criterion and Monte-Carlo simulation.

The method has been applied to individual wells (JY1 and W204) and two active shale gas blocks (Jiaoshiba and Weiyuan) of the Sichuan Basin in China. The predictions are in good agreement with data from micro-seismic fracture monitoring and the actual gas production.

The factors σ_H, σ_v, σ_h, S_t, α, and β are considered in this method. A complex fracture network can be formed more easily when α is closer to the direction of the σ_H and β is closer to the vertical direction.

Compared with conventional methods, the proposed method can be directly applied to actual oil/gas wells and fields, which can provide guidance for the design and implementation of hydraulic fracturing protocols.

References

[1] Blanton T L (1982). An experimental study of interaction between hydraulically induced and pre-existing fractures. In: SPE-10847 Presented at the SPE Unconventional Gas Recovery Symposium, Pittsburgh, Pennsylvania, USA

[2] Blanton T L (1986). Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs.In: SPE-15261 Presented at the SPE Unconventional Gas Technology Symposium, Louisville: Kentucky

[3] Bryant F B, Satorra A (2012). Principles and practice of scaled difference chi-square testing.Struct Equ Modeling, 19(3): 372–398

[4] Chen L, Lu Y C, Jiang S, Li J Q, Guo T L, Luo C (2015). Heterogeneity of the lower Silurian Longmaxi marine shale in the southeast Sichuan basin of China.Mar Pet Geol, 65: 232–246

[5] Chen S B, Zhu Y M, Qin Y, Wang H Y, Liu H L, Fang J H (2014). Reservoir evaluation of the Lower Silurian Longmaxi Formation shale gas in the southern Sichuan Basin of China.Mar Pet Geol, 57: 619–630

[6] Cipolla C L, Warpinski N R, Mayerhofer M J, Lolon E, Vincent M C (2008). The relationship between fracture complexity, reservoir properties, and fracture treatment design. In: SPE Annual Technical Conference and Exhibition. Denver, Colorado, USA. Paper number: SPE-115769

[7] Curtis M E, Sondergeld C H, Ambrose R J, Rai C S (2012). Microstructural investigation of gas shales in two and three dimensions using nanometer-scale resolution imaging.AAPG Bull, 96(4): 665–677

[8] Dahi-Taleghani A, Olson J E (2011). Numerical modeling of multistranded-hydraulic-fracture propagation: accounting for the interaction between induced and natural fractures.SPE J, 16(3): 575–581

[9] Dehghan A N, Goshtasbi K, Ahangari K, Jin Y (2015). The effect of natural fracture dip and strike on hydraulic fracture propagation.Int J Rock Mech Min Sci, 75: 210–215

[10] Dou F, Wang J G, Wang H M, Hu B W, Li C X (2019). Discrete element analysis for hydraulic fracture propagations in laminated reservoirs with complex initial joint properties.Geofluids, 2019(1): 1–23

[11] Fatahi H, Hossain M M, Sarmadivaleh M (2017). Numerical and experimental investigation of the interaction of natural and propagated hydraulic fracture.J Nat Gas Sci Eng, 37: 409–424

[12] Gale J F W, Reed R M, Holder J (2007). Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments.AAPG Bull, 91(4): 603–622

[13] Gasparik M, Bertier P, Gensterblum Y, Ghanizadeh A, Krooss B M, Littke R (2014). Geological controls on the methane storage capacity in organic-rich shales.Int J Coal Geol, 123: 34–51

[14] Gu H R, Weng X W, Lund J, Mack M, Ganguly U, Suarez-Rivera R (2012). Hydraulic fracture crossing natural fracture at nonorthogonal angles: a criterion and its validation.SPE Production & Operations, 27(1): 20–26

[15] Guo X S (2014). Shale Gas Formation Mechanism and Exploration Technologies in Jiaoshiba Block, Fuling Shale Gas Field. Beijing: Science Press (in Chinese)

[16] Hou B, Zhang R X, Tan P, Song Y, Fu W N, Chang Z, Kao J W, Muhadasi Y, Chen M (2018). Characteristics of fracture propagation in compact limestone formation by hydraulic fracturing in central Sichuan, China.J Nat Gas Sci Eng, 57: 122–134

[17] Hou B, Zhang R X, Chen M, Kao J, Liu X (2019). Investigation on acid fracturing treatment in limestone formation based on true tri-axial experiment.Fuel, 235: 473–484

[18] Jeffrey R G, Bunger A, Lecampion B, Zhang X, Chen Z R, Van As A, Allison D P, Beer W D, Dudley J W, Siebrits E, Thiercelin M J, Mainguy M (2009). Measuring hydraulic fracture growth in naturally fractured rock. In: SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana. Paper number: SPE-124919

[19] Jin Z J, Nie H K, Liu Q Y, Zhao J H, Jiang T (2018). Source and seal coupling mechanism for shale gas enrichment in upper Ordovician Wufeng Formation-Lower Silurian Longmaxi Formation in Sichuan Basin and its periphery.Mar Pet Geol, 97: 78–93

[20] Li H, Lai B T, Liu H H, Zhang J L, Georgi D (2017). Experimental investigation on Brazilian tensile strength of organic-rich gas shale.SPE J, 22(1): 148–161

[21] Li Y N (2014). Shale gas reservoir logging evaluation and its application. Dissertation for Doctoral Degree. Beijing: China University of Mining and Technology Beijing (in Chinese)

[22] Liu C, Jin X, Shi F, Lu D T, Liu H, Wu H A (2018a). Numerical investigation on the critical factors in successfully creating fracture network in heterogeneous shale reservoirs.J Nat Gas Sci Eng, 59: 427–439

[23] Liu J, Yao Y B, Liu D M, Xu L L, Elsworth D, Huang S P, Luo W J (2018b). Experimental simulation of the hydraulic fracture propagation in an anthracite coal reservoir in the southern Qinshui Basin, China.J Petrol Sci Eng, 168: 400–408

[24] Liu Z Y, Chen M, Zhang G Q (2014). Analysis of the influence of a natural fracture network on hydraulic fracture propagation in carbonate formations.Rock Mech Rock Eng, 47(2): 575–587

[25] Liu Z Y, Wang S J, Zhao H Y, Wang L, Li W, Geng Y D, Tao S, Zhang G Q, Chen M (2018c). Effect of random natural fractures on hydraulic fracture propagation geometry in fractured carbonate rocks.Rock Mech Rock Eng, 51(2): 491–511

[26] Ma X H, Li X Z, Liang F, Wan Y J, Shi Q, Wang Y H, Zhang X W, Che M G, Guo W, Guo W (2020). Dominating factors on well productivity and development strategies optimization in Weiyuan shale gas play, Sichuan Basin, SW China.Pet Explor Dev, 47(3): 594–602

[27] Mayerhofer M J, Lolon E, Warpinski N R, Cipolla C L, Walser D W, Rightmire C M (2010). What is stimulated reservoir volume?.SPE J, 25(1): 89–98

[28] Mokhtari M, Bui B T, Tutuncu A N (2014). Tensile failure of shales: impacts of layering and natural fractures. In: SPE Western North American and Rocky Mountain Joint Meeting. Denver, Colorado.

[29] Nam S G, Park C, Yoo J (2013). Uncertainty quantification of an asset evaluation for an oilfield property incorporating response-surface Monte Carlo simulation with stochastic oil price models.Energy Explor Exploit, 31(5): 783–795

[30] Page J C, Miskimins J L (2009). A Comparison of hydraulic and propellant fracture propagation in a shale gas reservoir.J Can Pet Technol, 48(5): 26–30

[31] Renshaw C E, Pollard D D (1995). An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic materials.Int J Rock Mech Min Sci Geomech Abstr, 32(3): 237–249

[32] Ross D J, Bustin R M (2007). Shale gas potential of the lower Jurassic Gordondale member, northeastern British Columbia, Canada.Bull Can Pet Geol, 55(1): 51–75

[33] Sarmadivaleh M, Rasouli V (2014). Modified Reinshaw and Pollard criteria for a non-orthogonal cohesive natural interface intersected by an induced fracture.Rock Mech Rock Eng, 47(6): 2107–2115

[34] Shapiro S A (2009). Microseismicity of non-linear fluid-rock interactions: from stimulations of geothermic reservoirs to hydraulic fracturing of shales. In: Ling H I, Smyth A, Betti R, eds. Poro-mechanics IV. Lancaster: DES Tech Publications Inc, 575–581

[35] Stoltz L R, Jones M S, Wadsley A W (1998). Probabilistic reserves assessment using a filtered Monte Carlo method in a fractured limestone reservoir. In: Asia Pacific Conference on Integrated Modelling for Asset Management Kuala Lumpur, Malaysia

[36] Tan P, Jin Y, Han K, Zheng X J, Hou B, Gao J, Chen M, Zhang Y Y (2017). Vertical propagation behavior of hydraulic fractures in coal measure strata based on true triaxial experiment.J Petrol Sci Eng, 158: 398–407

[37] Tong H M, Cai D S, Wu Y P, Li X G, Li X S, Meng L J (2010). Activity criterion of pre-existing fabrics in non-homogeneous deformation domain.Sci China Earth Sci, 53(8): 1115–1125

[38] Wang H Y (2019). Hydraulic fracture propagation in naturally fractured reservoirs: complex fracture or fracture networks.J Nat Gas Sci Eng, 68: 102911

[39] Weng X W, Kresse O, Cohen C E, Wu R T, Gu H R (2011). Modeling of hydraulic fracture network propagation in a naturally fractured formation.SPE J, 26(4): 368–380

[40] Xu D, Gao W, Hu R L (2016). A new method for hydraulic fracturing evaluation and its application in shale reservoirs. J Pet Geol Eng, 30(4): 138–141 (in Chinese)

[41] Xu J C, Guo C H, Wei M Z, Jiang R Z (2015). Production performance analysis for composite shale gas reservoir considering multiple transport mechanisms.J Nat Gas Sci Eng, 26: 382–395

[42] Xu S, Gou Q Y, Hao F, Zhang B Q, Shu Z G, Zhang Y Y (2020). Multiscale faults and fractures characterization and their effects on shale gas accumulation in the Jiaoshiba area, Sichuan basin, China.J Petrol Sci Eng, 189: 107026

[43] Yang J, Fu Y Q, Chen H F, Zeng L X, Li J S (2012). Rock mechanical characteristics of shale reservoirs. J Nat Gas Ind, 32(7): 12–14 (in Chinese)

[44] Yao G H, Chen Q, Liu H, Tan Y H, Wang L S, Du H Y, Zhu H L (2015). Experiment study on mechanical properties of bedding shale in Lower Silurian Longmaxi shale Southeast Chongqing. Chinese J Rock Mechan Eng, 34(S1): 3313–3319 (in Chinese)

[45] Yuan B, Su Y L, Moghanloo R G, Rui Z H, Wang W D, Shang Y Y (2015). A new analytical multi-linear solution for gas flow toward fractured horizontal wells with different fracture intensity.J Nat Gas Sci Eng, 23: 227–238

[46] Zhang H (2015). Preferred Shale gas production targets in the Longmaxi Formation in the Weiyuan area of the Sichuan Basin. Dissertation for Master Degree. Chengdu: Southwest Petroleum University (in Chinese)

[47] Zhao Y L, Shan B C, Zhang L H, Liu Q G (2016). Seepage flow behaviors of multi-stage fractured horizontal wells in arbitrary shaped shale gas reservoirs.J Geophys Eng, 13(5): 674–689

[48] Zhao Y L, Zhang L H, Shan B C (2018). Mathematical model of fractured horizontal well in shale gas reservoir with rectangular stimulated reservoir volume.J Nat Gas Sci Eng, 59: 67–79

[49] Zheng D Y, Pang X Q, Ma X H, Li C R, Zheng T Y, Zhou L M (2019). Hydrocarbon generation and expulsion characteristics of the source rocks in the third member of the Upper Triassic Xujiahe Formation and its effect on conventional and unconventional hydrocarbon resource potential in the Sichuan Basin.Mar Pet Geol, 109: 175–192

[50] Zhou J, Chen M, Jin Y, Zhang G Q (2008). Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs.Int J Rock Mech Min Sci, 45(7): 1143–1152

[51] Zhou J, Jin Y, Chen M (2010). Experimental investigation of hydraulic fracturing in random naturally fractured blocks.Int J Rock Mech Min Sci, 47(7): 1193–1199

RIGHTS & PERMISSIONS

Higher Education Press