Geofluids
Volume 2021 (2021), Article ID 5552287, 20 pages
https://doi.org/10.1155/2021/5552287

Research Article

Numerical Analysis of Multiple Factors Affecting Hydraulic Fracturing in Heterogeneous Reservoirs Using a Coupled Hydraulic-Mechanical-Damage Model

Tao Xu,¹ Mingyang Zhai,² Bo Huang,¹ Liaoyuan Zhang,¹ Aishan Li,¹ Quansheng Zhang,¹ Zilin Zhang,¹ Lei Wang,¹ and Lianchong Li²

¹Shengli Oilfield Branch Company, SINOPEC, Dongying 257000, China
²Center of Rock Instability and Seismicity Research, School of Resources and Civil Engineering, Northeastern University, Shenyang 110819, China

Correspondence should be addressed to Mingyang Zhai; zhaimingyang_2011@163.com

Received 17 February 2021; Revised 6 March 2021; Accepted 19 March 2021; Published 14 April 2021

Academic Editor: Chun Zhu

Copyright © 2021 Tao Xu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Hydraulic fracturing performance, affected by multiple factors, was essential to the economic exploitation of oil and gas in heterogeneous unconventional reservoirs. Multifactor analysis can gain insight into the fracturing response of reservoirs and in turn optimize the treatment design. Based on characterizations of the geological setting of a heterogeneous glutenite reservoir, the hydraulic fracture (HF) initiation and propagation process, as well as the stimulated reservoir volume (SRV), were simulated and analyzed using a coupled hydraulic-mechanical-damage model. The Weibull distribution was employed to describe rock heterogeneity. The numerical model was verified with microseism (MS) interpretation results of HF geometry. A multifactor analysis and optimization workflow integrating response surface methodology, central composite design (CCD), and numerical simulations was proposed to investigate the coupling effects of multiple geomechanical and hydrofracturing factors on SRV and identify the optimum design of fracturing treatment. The results showed that the horizontal stress difference and injection rate were the most significant factors to control the SRV. Increasing the injection rate and reducing fluid viscosity may contribute to improving the SRV. It is more difficult to increase the SRV at higher horizontal stress difference than at lower horizontal stress difference. The multifactor analysis and optimization workflow introduced in this work was a practical and effective method to control the HF geometry and improve the SRV. This study provided a deep understanding of the hydraulic fracturing mechanism and possessed theoretical significance for treatment design.

1. Introduction

Glutenite reservoirs are widely distributed throughout the world, including the Los Angeles Basin in North America, the Argentine basin in South America, and the Bohai Bay Basin in eastern China [1, 2]. The tight heterogeneous glutenite reservoirs are typically characterized by high variable lithology, low/ultralow permeability, poor porosity connectivity, and significant heterogeneity [3–5]. This type of reservoir is generally not suitable for commercial production without hydraulic fracturing [6]. With the popularization of hydraulic fracturing technology, the development of tight glutenite reservoirs has attracted widespread attention in recent years [7–11]. The tight heterogeneous glutenites of Shengli Oilfield are located in Jiyang depression, Bohai Bay Basin, eastern China. The exploration area is nearly 1500 km², and the potential oil and gas reserves were found to exceed 450 million tons [12]. The oil and gas development of deep tight heterogeneous glutenites has become the focus of fossil energy exploitation in Shengli Oilfield.

Both microseismic interpretation and production analysis in Shengli Oilfield have shown that oil production and ultimate recovery in tight glutenite reservoirs are mainly associated with the hydraulic fracture (HF) complexity and stimulated reservoir volume (SRV) [13]. Field operation experiences indicate that simple biwing HFs are typically generated after fracturing in tight heterogeneous glutenites, which may lead to a rapid decline in production and low oil/gas recovery. The geological characteristics of tight glutenite reservoirs are different from other sedimentary reservoirs, such as shale and coal, due to the existence of gravels and the absence of a developed natural fracture system [14, 15]. It is relatively difficult to create complex HF or fracture networks in tight glutenite reservoirs. It is of great significance for the economic development of tight heterogeneous glutenites to optimize hydraulic fracture (HF) geometry and improving stimulated reservoir volume (SRV). With the development of fracturing technology, Shengli Oilfield is no longer satisfied with the traditional double-wing long fractures and has attempted to optimize injection strategies to create complex HFs and improve SRV. To achieve this, it is necessary to understand the mechanism of hydraulic fracturing and optimize the treatment parameters. The hydraulic fracturing performance is jointly affected by geomechanical properties, geological environment, and treatment parameters. These factors can be classified into two categories: uncontrollable factors (i.e., geomechanical properties and in situ stress environment) and controllable factors (i.e., treatment parameters) [16]. The geomechanical conditions, such as in situ stress and mechanical properties of glutenites, may be difficult to be changed, but we can improve fracture complexity and maximize SRV by optimizing the treatment parameters, such as injection rate and fluid viscosity. An optimized hydraulic fracture design requires a comprehensive understanding of the coupling effects of multiple controlling factors on HF propagation and SRV and finds the optimal combinations of geomechanical conditions and fracturing operations. Figure 1 shows the schematic diagram of HF geometries and their corresponding SRVs in tight heterogeneous glutenites. Complex HF usually possesses a larger SRV than simple HF [5, 17–21]. Multiple secondary fractures, even tertiary fractures, generate on the main fractures and eventually form complex HFs. Whether tight glutenites can form complex HFs or not depends on geological factors and fracturing operational factors [22–26].

Figure 1: Schematic diagram of HF geometries and their corresponding SRVs in tight heterogeneous glutenites.

Extensive studies have investigated the propagation mechanism of HFs in tight heterogeneous glutenites based on numerical simulations and laboratory tests. Li et al. [27] simulated the hydraulic fracturing process using a 2D numerical model and summarized five fracturing behaviors when a fracture approaches a gravel particle, namely, termination, penetration, deflection, branching, and attraction. Rui et al. [28] conducted 2D numerical simulations to investigate the influence of rock mechanical parameters and gravel property on the hydraulic fracture propagation and indicated that the HF geometry is controlled by reservoir physical and mechanical parameters, in situ stress, and operational parameters. Ju et al. [29] studied the hydraulic fracturing behaviors of heterogeneous glutenites and their influencing mechanisms using 3D numerical simulation based on bonded particle models and suggested that gravels in heterogeneous glutenites inhibit crack propagation. Ma et al. [30] investigated the growth law of HFs in glutenite reservoirs based on laboratory experiments and examined the effects of gravel size, horizontal differential stress, fluid viscosity, and injection rate on the HF growth behaviors. However, in these studies, the hydraulic fracture propagation under the coupling effects of multiple factors was still not fully understood. Moreover, most of the present researches studied the hydraulic fracturing mechanisms in tight heterogeneous glutenites at laboratory scales, while little researches studied them in field scales. The conclusions drawn from experiments and simulations in the laboratory scale cannot be directly used in field operations due to the scale effects [8, 31–34]. Therefore, it is of great significance to carry out field-scale numerical simulations to optimize the operational parameters and provide guidance for fracturing design and field operation.

In this paper, a coupled hydraulic-mechanical-damage model was presented to simulate the hydraulic fracturing process in tight heterogeneous glutenite reservoirs. The MS interpretation results of HF in tight glutenites were used to verify the validity and accuracy of the numerical model. The central composite design (CCD) was adopted to design the numerical experiments. The response surface methodology was employed to establish the response surface in terms of SRV with various hydrofracturing-related parameters, based on which the coupling effects of five factors, including brittleness index, tensile strength, horizontal stress difference, injection rate, and fluid viscosity on SRV, were analyzed. Finally, the SRV of tight heterogeneous glutenites was optimized, and the optimal combinations of these influential parameters were obtained based on the acquired response surface model. This study provided a practical approach for multifactor analysis and SRV optimization in unconventional reservoirs by integrating numerical modeling and optimization algorithm. The results had a theoretical significance for understanding the mechanism of hydraulic fracturing in tight heterogeneous glutenites.

2. Geologic Setting

The tight heterogeneous glutenites of Shengli Oilfield are located in Jiyang depression, Bohai Bay Basin, eastern China. The exploration area is nearly 1500 km², and the potential oil reserves were found to exceed 450 million tons. In recent years, the oil development of deep tight heterogeneous glutenites has become the focus of fossil energy exploitation in Shengli Oilfield. The glutenite reservoir has a low porosity with an average porosity of 7.69%. The reservoir is dominated by gravel-bearing coarse sand and gravel-bearing medium sand, with a particle size distribution range of 0.25~20 mm and a particle size concentration distribution range of 0.5 ~ 4.5 mm, with an average value of about 1 mm. The glutenite contains mainly fine gravels, and the gravel particle size is mostly concentrated between 2 and 4 mm, with an average content of about 15.1% and a maximum of 40%. The overall sorting property is poor to medium, the roundness is subangular, particle support, linear contact is the main one, and occasional point contact is seen. The cementation type is mainly porous cementation, occasionally conjunctival-porous cementation. Lithologically, the glutenite mainly consists of fine-grained lithic arkose. The mean composition is 32% quartz, 16% feldspar, 35% plagioclase, 3% calcite, 6% dolomite, 5% clay, and minor anhydrite. The gravel is mainly composed of quartz, feldspar, etc. The pores between gravels are filled with unequal grain sand, which is mainly mixed unevenly with medium sand, coarse sand, and fine sand.

The tight heterogeneous glutenites in Shengli Oilfield are characterized by poor physical properties, low/ultralow permeability, and significant heterogeneity, which result in various challenges in reservoir development. Furthermore, the burial depth of glutenite reservoirs is approximately between 3300 m and 4050 m, and in situ stress presents a high horizontal stress difference, ranging from 7 MPa to 11 MPa, which is not conducive to the formation of complex HFs. To improve the glutenite oil development in the Bohai Bay Basin, eastern China, it is important to understand the mechanism of hydraulic fracturing and optimize the treatment parameters. How to maximize the SRV and improve the range of the stimulated zone plays an important role in successful economic oil production for tight glutenites.

3. Methodology

3.1. Numerical Simulation

3.1.1. Description of RFPA3D-Petrol Model

Numerical modeling can overcome the limitations of the physical model in simulation scales and costs and provide valuable insight into the hydraulic fracturing process. To study the fracturing response to multiple factors involved in hydraulic stimulation in tight heterogeneous glutenites, a series of field-scale numerical simulations are performed based on RFPA3D-petrol software [35–37]. The reservoir heterogeneity, in addition to the stress-seepage-damage coupling, is considered to simulate the evolution process and complex propagation behavior of HF. An eight-node element was employed as the fundamental element. The mechanical parameters of elements (such as Young’s modulus and strength properties) are assumed to follow a Weibull distribution:

where is the mechanical parameters of rock obtained from laboratory tests, is the mean value of the variable , and is the degree of homogeneity, which is called the homogeneity index. A larger indicates a more homogeneous material and vice versa. is the distribution density of the variable .

The stress and fluid-flow fields are obtained using the Finite Element Method (FEM), and the coupling between stress and fluid flow is governed by Biot’s constitutive theory. The basic governing equations in RFPA3D are as follows:

where and are the stress tensor and strain tensor of the solid, respectively; is the volume strain; and are the body force and displacement components, respectively; is the effective stress; is the pore pressure coefficient; is the pore pressure; is the Kronecker delta; is the Lame constant; is the shear modulus of the solid; is Biot’s coefficient; is the permeability coefficient that varies with the stress-damage evolution; is the initial permeability coefficient; is an abrupt coefficient of the permeability, which is induced by the rock element damage; and is a coupling coefficient that reflects the influence of the normal stress on the permeability.

As the damage evolves, Young’s modulus of the element gradually decreases, as defined below:

where is the damage variable and and are Young’s modulus of the undamaged and damaged elements, respectively.

When the tensile stress in the element reaches its tensile strength, , the variable of the element under uniaxial tension can be expressed as

where is the tensile strain at the elastic limit, as shown in Figure 2(a). is the ultimate tensile strain, which is defined as . is the ultimate strain coefficient. is the residual tensile strength coefficient.

(a) The case under uniaxial tensile stress

(b) The case under uniaxial compressive stress

Figure 2: Elastic-brittle damage constitutive law of element subject to uniaxial stress.

When the equivalent principal tensile strain reaches the threshold strain,, under multiaxial stress states, the element is assumed to be damaged in the tensile mode. is defined as follows:

where , , and are the three principal strains. is a function defined as

The variable is then expressed as

When the element is under uniaxial compression, the constitutive law is shown in Figure 2(b). The Mohr-Coulomb criterion is selected as the damage threshold to describe the element damage under a compressive or shear stress condition:

where and are the principal stresses. is the internal friction angle. is the uniaxial compressive strength. The variable under uniaxial compression can be expressed as

where is the compressive strain at the elastic limit.

When the Mohr-Coulomb criterion is satisfied for the strength of the element subjected to multiaxial stress, the greatest compressive principal strain can be obtained as follows:

where is Poisson’s ratio. Thus, the shear damage under triaxial stress conditions can be extended from Equation (10) as follows:

Young’s modulus of the damaged element under distinct stress levels can then be obtained as per the previously derived equations related to the damage variable .

The fracture aperture is represented as a linear function of the net pressure in the RFPA3D-petrol model:

where is fracture aperture, is net pressure in the HF, and is fracture height.

3.1.2. Model Setup

A field-scale numerical model is established based on RFPA3D-petrol software as shown in Figure 3. It has 3 layers including 1 oil layer and 2 barrier layers. The numerical model was established based on a vertical well named Y-11 in Jiyang depression, Bohai Bay Basin, eastern China, by Shengli Oilfield. The mechanical parameters of different layers are shown in Table 1. In this simulation, 2 barrier layers share the same mechanical parameters. The dimension of the reservoir model is , and the model is divided into 5040000 FEM units. Its mesh size, fracturing parameters, and boundary conditions are listed in Table 2. The thickness of the oil layer is 20 m. To reduce the influence of boundary effect on fracture propagation, the barrier layers with a thickness of 20 m are set up above and below the oil layer, respectively. The vertical stress, maximum horizontal stress, and minimum horizontal stress of the reservoir are set to 72 MPa, 64 MPa, and 57 MPa, respectively, according to logging interpretation. The perforation section with a height of 2 m is prefabricated in the center of the model, which is used to inject fracturing fluid. The viscosity of fracturing fluid is 0.023 Pa·s, the injection rate is 8.0 m³/min, and the total fluid volume is 480 m³. The mechanical parameters of different layers are derived from experiment tests of 121 core samples drilled from well Y-11 [38]. The homogeneity index is obtained by statistical analysis of the uniaxial compressive strengths and elastic moduli of glutenite cores. The hydraulic parameters and the in situ stresses are interpreted from the well logging data of well Y-11 [39].

Figure 3: The elastic modulus map of the numerical model.

Table 1: Mechanical parameters of different layers.

Parameter
Symbol
Unit
Oil layer
Barrier layers

Macroscopic elastic modulus

GPa
33
28

Uniaxial compressive strength

MPa
72
93

Homogeneity index

—
4
4

Compression-tension ratio
—
—
10
10

Friction angle

°
30
35

Poisson ratio

1
0.21
0.23

Permeability

m²

Porosity

%
12.1
0.1

Storativity
Ss
m⁻¹

Residual strength coefficient

—
0.1
0.1

Ultimate tensile strain coefficient

—
5
5

Ultimate compressive strain coefficient

—
200
200

Table 2: Mesh size, fracturing parameters, and boundary conditions of the numerical model.

Parameter
Symbol
Unit
Value

Volume

m³

Element number
nele
—
5040000

Injection rate

m³·min⁻¹
8

Fluid density

kg m⁻¹
1000

Fluid viscosity

mPa·s
53

Coupling coefficient

—
1.8

Maximum horizontal stress

MPa
64

Minimum horizontal stress

MPa
57

Vertical stress

MPa
72

3.1.3. HF Propagation and Model Validation

In RFPA3D-HF, the hydraulic fracture geometry is reflected by the aperture of elements. In the numerical simulation process, different elements possess different apertures. The hydraulic fractures are represented by full failure elements (damage variable ) and damage elements (damage variable ). The matrix is represented by elastic elements (). The fracture geometry and aperture of elements at different fracturing times are illustrated in Figure 4. It is shown that the HF exhibits an X-shaped complex fracture geometry. The fracture length is 218 m, and the height is 32 m. The main fracture begins to bifurcate at a fracturing time of 10 minutes due to the strong heterogeneity of the glutenite formation. The secondary fractures at both ends of the main fracture are characterized by highly asymmetric propagation originally. The secondary fractures continue to propagate over time. At first, the two secondary fractures extend in the same direction at a close range; the stress interference between the fractures is quite severe. Thus, they are mutually exclusive due to the stress shadow effects, which makes the other one propagate far away from its path. As they continue to propagate away from each other, the stress interference between them is weakened, and their extension paths tend to be parallel to the maximum principal stress.

Figure 4: The X-shaped HF propagation process obtained from numerical simulation.

According to the definition of SRV [40, 41], the value of SRV is estimated by building up a three-dimensional geometric body to completely envelope the monitored MS events and then calculating the volume of the geometric envelope. Due to the geological uncertainty of tight glutenite reservoirs, it is challenging to simulate accurate SRV in numerical modeling. In our numerical model, the simulated SRV was calibrated with the MS interpretation results with consideration of the MS data uncertainty. The SRV for stage 2 of well Y-11 is 229000 m³. In our numerical model, the hydraulic fracture is represented by full failure elements (damage variable ) and damage elements (damage variable ). When calculating SRV, the pressurized elements around the hydraulic fracture are automatically searched and accumulated in the model until the total number of these elements equals the SRV of microseismic interpretation. Thus, the SRV for a certain HF in the numerical model can be calculated, which will be used to evaluate the hydraulic fracturing performance subsequently. Figure 5 illustrates the 3D view of the fracture aperture and simulated SRV in the numerical model.

(a)

(b)

Figure 5: The 3D view of the (a) HF geometry and (b) SRV obtained from numerical simulation.

To verify the validity and accuracy of the numerical model, the MS monitoring data of the HF of stage 2 in well Y-11 in Shengli Oilfield were interpreted and analyzed. Figure 6 illustrates the HFs and their corresponding SRVs obtained from MS interpretation and numerical simulation. For stage 2 of well Y-11, a total of 146 microseismic events were detected during fracturing operations, and the HF presents an X-shaped morphology, as shown in Figure 6(a). Comparing Figures 6(a) and 6(b), it is obvious that the HF obtained from numerical simulation is consistent with the results of MS interpretation. The length and height of HF interpreted from MS events are 223 m and 36 m, respectively. The geometric parameters of HFs obtained from MS monitoring and numerical simulation are quite close, which proves the validity and accuracy of the numerical model.

(a)

(b)

(c)

Figure 6: The comparison of the HF geometry and SRV obtained from MS interpretation and numerical simulation.

3.2. Response Surface Methodology

3.2.1. CCD Experimental Design

The CCD experimental design [42] was employed to obtain the mathematical relationship between various hydrofracturing-related factors and SRV responses. The CCD experimental design is a common experimental design method for response surface analysis. It is suitable for multifactor and multilevel experiments, which can effectively reduce the workload and ensure the calculation accuracy simultaneously. A five-factor five-level experimental design was employed in response surface analysis. The response surface was calculated from the results of numerical experiments. The hydraulic fracturing performance is controlled by geological conditions and treatment parameters. The reservoir tensile strength, brittleness, and horizontal stress difference are the main geological factors affecting hydraulic fracture geometry. As for treatment parameters, the injection rate and viscosity of fracturing fluid are the key parameters controlling fracture geometry. Therefore, these five parameters were selected to perform the experimental design and SRV optimization. For CCD experimental design, 5 levels of each parameter are usually needed. The geomechanical and hydraulic parameters of the numerical model established in this work are consistent with the field fracturing operation. The upper and lower limits of each parameter are determined according to the field fracturing operation. Based on this, five levels are selected between the upper and lower limits of each parameter for CCD experimental design. According to the geological and hydrofracturing parameters of the tight glutenite in well Y-11 of Shengli Oilfield, each uncertain parameter is given a reasonable range with five value levels, respectively, as listed in Table 3.

Table 3: Value levels of different parameters for CCD experimental design.

Parameter
Coded symbol
Unit
Level -2
Level -1
Level 0
Level 1
Level 2

Brittleness index

—
0.15
0.30
0.45
0.60
0.75

Tensile strength

MPa
3
5
7
9
11

Horizontal stress difference

MPa
5
7
9
11
13

Injection rate

m³/min
4
6
8
10
12

Fluid viscosity

mPa·s
3
28
53
78
103

The definition of the brittleness index proposed by Rickman et al. [43] was used to quantify the glutenite brittleness in numerical simulations. It was generalized as the following equation:

where and are the maximum and minimum values of Young’s modulus , respectively. and are the maximum and minimum values of Poisson’s ratio , respectively. According to well logging data and core laboratory test results of well Y-11, the following values were employed in brittleness index calculation: , , , and .

According to the principle of CCD optimal design, 50 groups of model tests were required for five-factor five-level experimental design. Table 4 shows the 50 combinations of these parameters generated by the CCD experimental design. The parameters of each case are input into the numerical model for hydraulic fracturing established in Section 3.1.2; then, the HF geometry and corresponding SRV can be obtained. Table 4 also lists the response value of SRV for each case. The data listed in Table 4 will be used for sensitivity analysis, multifactor analysis, and SRV optimization subsequently.

Table 4: CCD experimental design for the numerical simulation cases.

Run
(BI)
(T) (MPa)
(HSD) (MPa)
(IR) (m³/min)
(FV) (mPa·s)
SRV (m³)

1
0.6
5
11
6
78
173751

2
0.3
5
11
10
78
186716

3
0.3
5
7
6
78
210852

4
0.45
7
9
8
53
223992

5
0.3
5
7
6
28
252678

6
0.6
5
11
10
78
199153

7
0.6
9
11
10
28
214317

8
0.3
5
7
10
28
307224

9
0.6
9
11
10
78
194971

10
0.6
5
7
10
78
271889

11
0.45
7
9
8
3
250824

12
0.3
9
11
10
78
208920

13
0.6
9
7
6
28
247717

14
0.6
5
11
10
28
237263

15
0.45
7
9
12
53
318294

16
0.45
3
9
8
53
220928

17
0.45
7
13
8
53
173216

18
0.3
9
7
10
28
291938

19
0.3
5
11
6
78
161220

20
0.45
7
9
4
53
192264

21
0.45
7
9
8
53
238992

22
0.3
9
7
10
78
251072

23
0.45
7
9
8
53
223992

24
0.6
9
7
10
78
271496

25
0.3
5
7
10
78
265012

26
0.6
9
11
6
28
198490

27
0.45
7
9
8
103
180620

28
0.3
9
7
6
28
236669

29
0.6
5
7
6
28
255033

30
0.3
9
7
6
78
211866

31
0.6
9
7
10
28
315737

32
0.45
7
9
8
53
223992

33
0.45
7
9
8
53
223992

34
0.45
11
9
8
53
204572

35
0.3
5
11
10
28
211851

36
0.3
9
11
10
28
213919

37
0.6
9
11
6
78
179018

38
0.6
9
7
6
78
240270

39
0.3
5
7
10
78
279644

40
0.45
7
9
8
53
223992

41
0.15
7
9
8
53
207008

42
0.3
5
11
6
28
179522

43
0.45
7
9
8
53
223992

44
0.3
9
11
6
78
145813

45
0.6
5
7
10
28
318080

46
0.45
7
9
8
53
223992

47
0.3
9
11
6
28
169840

48
0.45
7
5
8
53
299884

49
0.6
5
11
6
28
185501

50
0.75
7
9
8
53
233064

Figure 7 shows the numerical simulation results of the partially selected 9 run cases in Table 4. The HFs exhibit different morphologies under different geological and fracturing conditions. Multiple secondary and tertiary fractures can be observed in the HFs of Case 9, Case 12, Case 15, and Case 31. The common feature of these four cases is that the large injection rate is employed. However, for Case 20, a simple biwing HF is created due to the small injection rate. Therefore, it can be preliminarily concluded that a large injection rate is beneficial to the formation of complex HFs. From Case 9, Case 12, Case 14, Case 17, and Case 42, it can be seen that even if the horizontal stress difference is high (11 MPa~13 MPa), it is still possible to create complex HFs. The HF geometry and corresponding SRV are affected by the coupling effects of various factors. Understanding the mechanism of hydraulic fracturing under the influence of multiple factors is of great significance to the optimization of fracturing design.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 7: The numerical simulation results of the partial 9 run cases in Table 4.

3.2.2. RSM Setup

To enhance the productivity and ultimate recovery of tight heterogeneous glutenites, the response surface methodology is employed to optimize the fracturing treatment parameters and improve the HF complexity and SRV. Response surface methodology is a statistical-based mathematical method, which is commonly used in modeling and analysis of engineering problems. The basic principle of response surface methodology is to establish an approximate explicit function between the output response (SRV) and the input variables (influencing factors) and analyze the sensitivity of output response (SRV) to multiple input variables (influencing factors). Figure 8 shows the schematic diagram of the response surface. It can effectively reduce the numerical calculation workload and ensure the result accuracy simultaneously. Moreover, once an accurate response surface model (RSM) is established, the significance of factors and their interactions can be quantified. Meanwhile, the optimal SRV and its corresponding parameters can be easily obtained, as shown in Figure 8. In this paper, the SRV was taken as the optimization objective to search the corresponding parameters and optimize the design of treatment parameters.

Figure 8: The schematic diagram of the response surface.

Figure 9 shows the workflow of the RSM-based multifactor analysis and SRV optimization for tight heterogeneous glutenite reservoirs. Five factors, including mechanical properties (brittleness index and tensile strength), geological environment (in situ stress), and treatment parameter (injection rate and fluid viscosity), are considered in the response surface analysis. The workflow is aimed at sensitivity analysis of SRV, multifactor analysis, and SRV optimization by integrating RSM, CCD experimental design, and numerical modeling.

Figure 9: The workflow of the RSM-based multifactor analysis and SRV optimization for glutenites.

Generally, the polynomial functions are often used as the response surface. The statistical approach was employed to select an appropriate model as the response surface, among the linear function, two-factor interaction function, quadratic function, and cubic function. The results are listed in Table 5. It can be seen that the quadratic function is suggested, and the cubic function is aliased. Therefore, the quadratic polynomials with crossing terms are used to describe the mathematical relationships between the SRV and the uncertain factors:

where is the response value, and are the uncertain variables, and ,,, and are the coefficients to be determined of the RSM.

Table 5: Statistical approach to select the RSM model for SRV.

Source
Sequential value
Lack of fit value
Adjusted
Predicted

Linear
<0.0001
0.0136
0.9156
0.9001

2FI
0.4498
0.0124
0.9159
0.8839

Quadratic
0.0010
0.0613
0.9497
0.8948
Suggested

Cubic
0.2232
0.0726
0.9602
0.6949
Aliased

According to the numerical results listed in Table 4, the multivariate equations are established, and the coefficients of the response surface function are solved by the least square method. The equation of the RSM for SRV with coded symbol is

where is the brittleness index, is tensile strength (MPa), is horizontal stress difference (MPa), is the injection rate (m³/min), and is fluid viscosity (mPa·s).

3.2.3. RSM Validation

To test and verify the reliability of the acquired RSM, the normal probability plot of the externally Studentized residuals is analyzed to check for normality of residuals. The results are illustrated in Figure 10(a). It can be seen that all the points fall on the straight line, meaning the residuals are normally distributed.

(a)

(b)

Figure 10: Validation of the RSM model. (a) The normal plot of residuals and (b) predicted SRV versus actual SRV.

The plot of “predicted SRV versus actual SRV” is shown in Figure 10(b). The correlation coefficient of the mathematical model is 0.954, and the fitting degree is good. The acquired RSM is capable of providing accurately predicted values of SRV, as compared with the actual values of SRV. This means that the established RSM for SRV analysis is reliable.

4. Results and Analysis

4.1. Sensitivity Analysis

To quantitatively analyze the factors affecting HF geometry and SRV, the analysis of variance (ANOVA) method is adopted for sensitivity analysis using the numerical experimental results in CCD experimental design. The results of ANOVA are listed in Table 6. The model value of 47.3 indicates that the model is significant. In statistics, value was usually used to evaluate the significance of parameters. Taking 5% as the threshold before testing, when the value was more than 5%, it was considered that the observed data is inconsistent with the zero hypothesis; the zero hypothesis should be rejected, and its opposite hypothesis should be accepted. When the value is less than 0.001, the parameters are extremely significant. When the value was more than 0.05, the parameters are not significant. In this case, factors of , , , and are significant terms. The influential order of these five factors is . The horizontal stress difference is the most sensitive factor, followed by the injection rate. This implies that horizontal stress difference and injection rate are the most dominant factors to hydraulic fracturing performance. This conclusion is consistent with King [44], Nagel et al. [45], and Wang et al. [46]. The horizontal stress difference determines the HF propagation path and final morphology. Under high horizontal stress difference, HF tends to propagate along with the maximum horizontal stress; deflection and bifurcation behavior are significantly limited. Besides, HF is prone to pass through natural discontinuities, thus significantly decreasing its complexity. On the contrary, low horizontal stress difference is beneficial to the formation of complex HFs. With a low injection rate, fluid is prone to leak off into the natural discontinuities, and the pressure can rise above the confining stress without inducing new hydraulic fractures. The net pressure in HFs increases with the increase of the injection rate, resulting in that the HF tends to connect more natural discontinuities due to the increase of the pressurization rate. In this case, the HF complexity may increase.

Table 6: ANOVA of response surface for SRV.

Source
Sum of squares
df
Mean square
value
value

Model
85782693635.2
20
4289134681.8
47.3
0.0001
Significant

(BI)
1699352244.5
1
1699352244.5
18.7
0.0003

(T)
204709990.1
1
204709990.1
12.3
0.0340

(HSD)
45947592466.4
1
45947592466.4
506.6
0.0000

(IR)
21866665859.3
1
21866665859.3
241.1
0.0001

(FV)
7848762991.7
1
7848762991.7
86.5
0.0002

80383456.3
1
80383456.3
0.9
0.3543

273343.0
1
273343.0
0.0
0.9571

110249754.0
1
110249754.0
1.2
0.2794

654923.6
1
654923.6
0.0
0.9332

47595197.2
1
47595197.2
0.5
0.4756

29324248.9
1
29324248.9
0.3
0.5747

101004299.1
1
101004299.1
1.1
0.3009

666473821.5
1
666473821.5
7.3
0.0114

338548304.8
1
338548304.8
3.7
0.0635

161631189.8
1
161631189.8
1.8
0.1923

116686135.2
1
116686135.2
1.3
0.2661

445138750.0
1
445138750.0
4.9
0.0354

157114793.6
1
157114793.6
1.7
0.1987

1521126591.6
1
1521126591.6
16.8
0.0001

285557410.8
1
285557410.8
3.1
0.0872

Residual
2630402476.7
29
90703533.7

Lack of fit
2326476838.3
21
110784611.3
2.9
0.061
Not significant

Pure error
303925638.4
8
37990704.8

Cor total
88413096111.9
49

4.2. Multiple Factor Analysis

To reveal the coupling effects of multiple geomechanical and hydraulic parameters on the HF propagation and resulting SRV, multiple factor analysis is conducted based on the acquired RSM. The influence of interaction between two factors on SRV is analyzed by fixing the values of other factors to zero levels. The 3D response surface and 2D contour map are drawn according to the regression equation, as illustrated in Figures 11–14. The response surfaces all exhibit curved shapes, indicating the nonlinear influential mechanism of these factors on SRV.

(a) 3D response surface

(b) 2D contour map

Figure 11: Response surface curve of SRV showing the interaction between injection rate and horizontal stress difference.

(a) 3D response surface

(b) 2D contour map

Figure 12: Response surface curve of SRV showing the interaction between fluid viscosity and horizontal stress difference.

(a) 3D response surface

(b) 2D contour map

Figure 13: Response surface curve of SRV showing the interaction between brittleness index and horizontal stress difference.

(a) 3D response surface

(b) 2D contour map

Figure 14: Response surface curve of SRV showing the interaction between fluid viscosity and injection.

Figure 11 shows that the SRV increases with the increase of the injection rate and decreases with the increase of horizontal stress difference. The highest red color is at an injection rate of 12 m³/min and a horizontal stress difference of 5 MPa. The results are inconsistent with Li et al. [15] and Beugelsdijk et al. [47]. Additionally, at a lower injection rate and a higher horizontal stress difference, the values tend to be blue, indicating smaller SRV. The effects of the injection rate and fluid viscosity on SRV are more significant under lower horizontal stress difference than that under higher horizontal stress difference. For example, when the injection rate increases from 4 m³/min to 12 m³/min, the SRV increases from 285000 m³ to 421000 m³ at a horizontal stress difference of 5 MPa, while the SRV increases from 160000 m³ to 234000 m³ at a horizontal stress difference of 13 MPa. It is more difficult to increase the SRV at higher horizontal stress differences than at lower horizontal stress differences. Nonetheless, the SRV can still be elevated within a certain range by increasing the injection rate.

Figure 12 illustrates that the SRV decreases with the increase of fluid viscosity and decreases with the increase of horizontal stress difference. The highest red color is in the fluid viscosity range of 3 mPa·s ~53 mPa·s and at a horizontal stress difference of 5 MPa. At a higher fluid viscosity and horizontal stress difference, the values tend to be blue, indicating smaller SRV. Comparing the value of SRV in different conditions, it can be concluded in Figure 11 that under lower horizontal stress difference, the effect of fluid viscosity on SRV is more significant than that under higher horizontal stress difference. Although low-viscosity fluid has a strong filtration capacity and is helpful to create complex HFs, its sand-carrying capacity is limited, which is not conducive to supporting fractures. Therefore, fracturing fluid with a viscosity of 28 mPa·s ~53 mPa·s is recommended under low horizontal stress difference, and fracturing fluid with a viscosity of 3 mPa·s ~28 mPa·s is recommended under high horizontal stress difference.

Figure 13 shows that the SRV increases with the increase of the brittleness index and decreases with the increase of horizontal stress difference. The highest red color is at a brittleness index of 0.75 and a horizontal stress difference of 5 MPa. Additionally, at a lower brittleness and a higher horizontal stress difference, the values tend to be blue, indicating smaller SRV. Many previous studies have shown that reservoir brittleness has a significant impact on the HF complexity, and reservoirs with strong brittleness are more likely to generate complex HFs. With the increase of horizontal stress difference, it is difficult to create complex HFs even in brittle glutenite reservoirs. Therefore, it is necessary to optimize the treatment parameters, reasonably elevating the injection rate, and lowering fluid viscosity might help.

Figure 14 shows that the SRV increases with the increase of the injection rate and decreases with the increase of fluid viscosity. The highest red color is at an injection rate of 12 m³/min and a fluid viscosity of 3 mPa·s. The results are consistent with Yao [48] and Zhou et al. [49]. Additionally, at a lower injection rate and a higher fluid viscosity, the values tend to be blue, indicating smaller SRV. It is also shown in Figure 10 that the effect of the injection rate on SRV is more significant than that of fluid viscosity, which is inconsistent with the result of sensitivity analysis in Section 4.1.

4.3. SRV Optimization

In the process of hydraulic fracturing operation, engineerings are usually concerned about how to select the reasonable treatment parameters to obtain a maximum SRV for a tight heterogeneous glutenite under specific geological conditions. As mentioned above, the RSM-based numerical optimization can not only analyze the influence of factors on HF complexity and SRV but also optimize them by selecting the set of variables that leads to the maximum SRV. The index of SRV indicates the overall response of the tight heterogeneous glutenites subjected to hydraulic fracturing. For a specific glutenite reservoir, the best HF configuration corresponds to the largest SRV. Therefore, in this section, the desirable HF geometry is obtained through SRV optimization.

A total of 100 optimal solutions were generated after response surface analysis. Table 7 lists the 25 out of 100 desirable solutions to show the optimal SRVs and corresponding parameter combinations. In general, the tight glutenites with strong brittleness, low tensile strength, and low horizontal stress difference, combined with a large injection rate and low-viscosity fracturing fluid, are beneficial to obtain the maximum SRV.

Table 7: 25 out of 100 desirable solutions for SRV optimization.

Number
(BI)
(T) (MPa)
(HSD) (MPa)
(IR) (m³/min)
(FV) (mPa·s)
SRV (m³)

1
0.57
5.37
5.38
11.42
57.46
369439

2
0.62
3.23
7.4
11.61
42.18
331242

3
0.53
7.55
5.18
11.51
48.53
375204

4
0.51
10.51
5.05
8.97
11.66
325712

5
0.48
4.56
5.21
9.41
34.19
330647

6
0.43
7.81
5.16
11.54
59.06
365943

7
0.6
5.98
5.41
8.03
7.4
330018

8
0.54
7.34
6.08
9.99
52.29
320311

9
0.56
6.31
5.35
9.31
4.41
356766

10
0.57
7.77
6.12
10.28
29.93
341022

11
0.43
7.34
5.36
9.84
20.51
353506

12
0.48
4.24
5.22
10.49
52.21
352337

13
0.38
7.99
5.75
10.07
30.94
337232

14
0.54
7.47
7.6
11.68
53.67
323797

15
0.73
7.7
5.38
7.98
25.03
318901

16
0.57
6.19
5.46
11.84
59.06
377211

17
0.49
5.88
5.19
9.21
6.33
357410

18
0.67
3.58
5.03
8.78
13.51
347520

19
0.39
5.34
5.46
11.84
59.06
376016

20
0.57
7.85
5.57
10.47
63.66
331249

21
0.4
8.49
5.47
10.51
43.99
343790

22
0.55
5.76
8.19
11.58
31.79
326107

23
0.44
6.49
5.4
10.44
74.22
323514

24
0.72
7.76
5.05
7.71
5.41
330651

25
0.37
3.81
7.03
11.62
53.39
333795

Figure 15 shows the effects of various factors on the desirable solutions of SRV according to Table 7. The mechanical properties of glutenites control the damage and failure mechanism during hydraulic fracturing. From Figures 15(a) and 15(b), the value of the brittleness index is set to close to the maximum and the value of tensile strength is set to below medium for the 25 desirable solutions. The strong brittleness and low tensile strength may contribute to the formation of complex HFs, thus increasing the SRV.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15: The effects of various factors on the desirable solutions of SRV.

From Figure 15(c), the value of horizontal stress difference is set to close to the minimum for all the desirable solutions. The horizontal stress difference is a critical factor controlling the SRV. For deep tight glutenite reservoirs with high in situ stress difference, the improvement of SRV by optimizing the injection rate and fluid viscosity is limited. The local in situ stress distribution can be modified by preinjection of low-viscosity fluid (e.g., slickwater) to artificially lower the horizontal stress difference, to improve the HF complexity and SRV [50, 51].

From Figures 15(d) and 15(e), the value of the injection rate is set to close to the maximum, and the value of fluid viscosity is set to below medium for all the desirable solutions. The injection rate and fluid viscosity have a remarkable influence on SRV. The large injection rate and low-viscosity or medium-viscosity fracturing fluid are required to obtain a maximum SRV. As mentioned above, it is difficult to change the mechanical properties and geological environment of tight glutenites, but we can optimize the treatment parameters to obtain a desirable SRV. The RSM-based SRV optimization discussed in this work may help. Once an accurate response surface was established, the treatment parameters could be optimized according to geological conditions to obtain a maximum SRV for tight glutenites. It should be noted that other factors that may significantly affect HF morphology and SRV, such as shear strength and reservoir thickness, are not covered in this work. These matters will be discussed in further study.

The multifactor analysis of SRV provided a deeper understanding of the hydraulic fracturing mechanism and SRV enhancement strategy. Firstly, on the basis of accurate characterization of geological conditions, engineers can conduct SRV prediction and analysis using the response surface method. Secondly, once an accurate response surface is established, the treatment parameters (injection rate and fluid viscosity) can be optimized ahead of the fracturing operation. Thirdly, according to the different sensitivity of SRV to treatment parameters under different geological conditions (e.g., formation brittleness or in situ stress), appropriate adjustments are made to the treatment parameters according to the multifactor analysis results.

5. Field Application

Using the numerical simulation and RSM-based SRV optimization method introduced in this paper, the treatment parameters for stage 3, stage 4, and stage 5 of well Y-11 tight glutenites in Shengli Oilfield are optimized. On the basis of accurate characterization of geological conditions (such as brittleness and in situ stresses) of the target layer, the hydraulic fracturing process and the SRV are numerically simulated. Then, the CCD experimental design is conducted, and the response surface model is established. Within the framework of the response surface model, the maximum SRV and its corresponding treatment parameters can be calculated. Therefore, the treatment parameters (injection rate and fluid viscosity) can be optimized ahead of the fracturing operation.

The MS interpretation results of HFs of different fracturing stages are shown in Figure 16. The SRVs of fracturing stage 1 to fracturing stage 5 are 114000 m³, 129000 m³, 179000 m³, 84000 m³, and 151000 m³, respectively. It can be seen that the HFs of three fracturing stages all display complex morphologies. Compared with stage 2, the fracture widths of stage 3 and stage 5 increased, and the SRV increased by 38.7% and 17.1%, respectively. As for stage 3, the large horizontal stress difference of 10.8 MPa limits the HF complexity; however, the complex HF is generated after optimization of treatment parameters. In conclusion, the field application results show that the RSM-based SRV optimization method is practical and effective in tight glutenite reservoirs and has a good application effect.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 16: The MS interpretation results of HFs of different fracturing stages for well Y-11 tight glutenites.

6. Conclusion

In this paper, the response surface methodology combined with a coupled stress-seepage-damage model was proposed for multifactor analysis and SRV optimization in hydraulic stimulation of tight heterogeneous glutenite reservoirs. The proposed approach provided a valuable insight into the fracturing mechanism and SRV enhancement strategy in tight heterogeneous glutenites, which was of innovative value for fracturing design and treatment parameter optimization. The conclusions are as follows:

(1)A three-dimensional coupled stress-seepage-damage numerical model of hydraulic fracturing was presented to reproduce the complex HF evolution process in tight heterogeneous glutenite reservoirs. The numerical model was validated with the actual HF of well Y-11 interpreted from MS monitoring data. The numerical model was practical and effective for simulating the complex HF propagation in tight glutenites(2)A multifactor analysis and SRV optimization workflow was proposed by integrating the response surface method, CCD experimental design, and numerical simulations. The analysis of variance (ANOVA) method is adopted for sensitivity analysis based on the established response surface. The influential order of the factors affecting the SRV in this work was (3)The coupling effects of five factors on SRV were analyzed. The effects of the injection rate and fluid viscosity on SRV are more significant under lower horizontal stress difference than that under higher horizontal stress difference. It is more difficult to increase the SRV at higher horizontal stress differences than at lower horizontal stress differences. Nonetheless, the SRV can still be elevated within a certain range by increasing the injection rate and decreasing the fluid viscosity(4)The optimal design solutions for tight heterogeneous glutenites were obtained by SRV optimization. Once an accurate response surface was established, the treatment parameters could be optimized according to geological conditions to obtain a maximum SRV for tight glutenites

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was conducted with support from the National Science and Technology Major Project of China (Grant No. 2017ZX05072), the National Natural Science Foundation of China (Grant NOs. 51879041 and 51761135102), and the Fundamental Research Funds for the Central Universities (No. N180105029).

References

Q. Zhang, M. Li, Z. Zhang, Y. Chen, L. Zhang, and A. Li, “Application of volume fracturing technology in tight oil reservoirs of Shengli Oilfield,” China Petroleum Exploration, vol. 24, no. 2, pp. 233–240, 2019.

Q. Wang, H. Gao, B. Jiang, S. Li, M. He, and Q. Qin, “In-situ test and bolt-grouting design evaluation method of underground engineering based on digital drilling,” International Journal of Rock Mechanics and Mining Sciences, vol. 138, article 104575, 2021.

B. Li, R. Bao, Y. Wang, R. Liu, and C. Zhao, “Permeability evolution of two-dimensional fracture networks during shear under constant normal stiffness boundary conditions,” Rock Mechanics and Rock Engineering, vol. 54, no. 1, pp. 409–428, 2021.

B. Chen, B. Xu, B. Li, M. Kong, W. Wang, and H. Chen, “Understanding the performance of hydraulically fractured wells in the laumontite-rich tight glutenite formation,” Journal of Natural Gas Science and Engineering, vol. 185, article 106600, 2020.

C. L. Cipolla, N. R. Warpinski, M. J. Mayerhofer, E. Lolon, and M. C. Vincent, “The relationship between fracture complexity, reservoir properties, and fracture-treatment design,” SPE Production & Operations, vol. 25, no. 4, pp. 438–452, 2010.

G. Zhang, S. Sun, K. Chao et al., “Investigation of the nucleation, propagation and coalescence of hydraulic fractures in glutenite reservoirs using a coupled fluid flow-DEM approach,” Powder Technology, vol. 354, pp. 301–313, 2019.

S. R. Wang, X. G. Wu, Y. H. Zhao, P. Hagan, and C. Cao, “Evolution characteristics of composite pressure-arch in thin bedrock of overlying strata during shallow coal mining,” International Journal of Applied Mechanics, vol. 11, no. 3, article 1950030, 2019.

P. Tan, Y. Jin, K. Han et al., “Analysis of hydraulic fracture initiation and vertical propagation behavior in laminated shale formation,” Fuel, vol. 206, pp. 482–493, 2017.

S. Wang, P. Hagan, D. Xu, B. Hu, Z. Li, and K. Gamage, “Fracture process and energy dissipation analysis of sandstone plates under the concentrated load, Teh,” Tehnički vjesnik, vol. 21, no. 6, pp. 1345–1351, 2014.

Y. Wang, X. Li, Z. H. Zhao, R. Q. Zhou, B. Zhang, and G. F. Li, “Contributions of non-tectonic micro-fractures to hydraulic fracturing—a numerical investigation based on FSD model,” Science China Earth Sciences, vol. 59, no. 4, pp. 851–865, 2016.

Q. Wang, Q. Qin, B. Jiang et al., “Mechanized construction of fabricated arches for large-diameter tunnels,” Automation in Construction, vol. 124, p. 103583, 2021.

Y. Zhao, P. He, Y. Zhang, and C. Wang, “A new criterion for a toughness-dominated hydraulic fracture crossing a natural frictional interface,” Rock Mechanics and Rock Engineering, vol. 52, no. 8, pp. 2617–2629, 2019.

Y. Shao, X. Huang, and Y. Xing, “An integrated study on the sensitivity and uncertainty associated with the evaluation of stimulated reservoir volume (SRV),” Journal of Petroleum Science and Engineering, vol. 159, pp. 903–914, 2017.

S. Lyu, S. Wang, X. Chen et al., “Natural fractures in soft coal seams and their effect on hydraulic fracture propagation: a field study,” Journal of Petroleum Science and Engineering, vol. 192, p. 107255, 2020.

Z. Li, L. Li, M. Li et al., “A numerical investigation on the effects of rock brittleness on the hydraulic fractures in the shale reservoir,” Journal of Natural Gas Science and Engineering, vol. 50, pp. 22–32, 2018.

Z. Tao, C. Zhu, M. He, and M. Karakus, “A physical modeling-based study on the control mechanisms of negative Poisson's ratio anchor cable on the stratified toppling deformation of anti- inclined slopes,” International Journal of Rock Mechanics and Mining Sciences, vol. 138, article 104632, 2021.

Z. Yushi, Z. Shicheng, Z. Tong, Z. Xiang, and G. Tiankui, “Experimental investigation into hydraulic fracture network propagation in gas shales using CT scanning technology,” Rock Mechanics and Rock Engineering, vol. 49, no. 1, pp. 33–45, 2016.

A. N. Dehghan, K. Goshtasbi, K. Ahangari, and Y. Jin, “The effect of natural fracture dip and strike on hydraulic fracture propagation,” International Journal of Rock Mechanics and Mining Sciences, vol. 75, pp. 210–215, 2015.

W. Cheng, Y. Jin, and M. Chen, “Reactivation mechanism of natural fractures by hydraulic fracturing in naturally fractured shale reservoirs,” Journal of Natural Gas Science and Engineering, vol. 23, pp. 431–439, 2015.

T. Guo, S. Zhang, Z. Qu, T. Zhou, Y. Xiao, and J. Gao, “Experimental study of hydraulic fracturing for shale by stimulated reservoir volume,” Fuel, vol. 128, pp. 373–380, 2014.

G. Xiang, W. Zhou, J. Yuan, and X. Chang, “Pore pressure cohesive zone modelling of complex hydraulic fracture propagation in a permeable medium,” European Journal of Environmental and Civil Engineering, vol. 1, pp. 1–17, 2019.

M. Behnia, K. Goshtasbi, G. Zhang, and S. Yazdi, “Numerical modeling of hydraulic fracture propagation and reorientation,” European Journal of Environmental and Civil Engineering, vol. 19, no. 2, pp. 152–167, 2015.

B. Hou, R. Zhang, Y. Zeng, W. Fu, Y. Muhadasi, and M. Chen, “Analysis of hydraulic fracture initiation and propagation in deep shale formation with high horizontal stress difference,” Journal of Petroleum Science and Engineering, vol. 170, pp. 231–243, 2018.

Y. Wang, H. Fu, T. Liang et al., “Large-scale physical simulation experiment research for hydraulic fracturing in shale,” in SPE Middle East oil & gas Show & Conference, Manama, Bahrain, 2015.

N.-Z. Liu, Y.-S. Zou, X.-F. Ma, N. Li, and S. Wu, “Study of hydraulic fracture growth behavior in heterogeneous tight sandstone formations using CT scanning and acoustic emission monitoring,” Petroleum Science, vol. 16, no. 2, pp. 396–408, 2019.

K. K. Chong, W. Grieser, O. Jaripatke, and A. Passman, “A completions roadmap to shale-play development: a review of successful approaches toward shale-play stimulation in the last two decades,” in Proceedings of International Oil and Gas Conference and Exhibition in China, pp. 1–13, Beijing, China, 2010.

L. Li, Q. Meng, S. Wang, G. Li, and C. Tang, “A numerical investigation of the hydraulic fracturing behaviour of conglomerate in glutenite formation,” Acta Geotechnica, vol. 8, no. 6, pp. 597–618, 2013.

Z. Rui, T. Guo, Q. Feng, Z. Qu, N. Qi, and F. Gong, “Influence of gravel on the propagation pattern of hydraulic fracture in the glutenite reservoir,” Journal of Petroleum Science and Engineering, vol. 165, pp. 627–639, 2018.

Y. Ju, Y. Wang, H. Dong, and Y. Yang, “Numerical analysis of the hydrofracturing behaviour of heterogeneous glutenite considering hydro-mechanical coupling effects based on bonded particle models,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 42, no. 13, pp. 1493–1515, 2018.

X. Ma, Y. Zou, N. Li, M. Chen, Y. Zhang, and Z. Liu, “Experimental study on the mechanism of hydraulic fracture growth in a glutenite reservoir,” Journal of Structural Geology, vol. 97, pp. 37–47, 2017.

C. Jiang, B. Niu, G. Yin, D. Zhang, T. Yu, and P. Wang, “CT-based 3D reconstruction of the geometry and propagation of hydraulic fracturing in shale,” Journal of Petroleum Science and Engineering, vol. 179, pp. 899–911, 2019.

J. Liu, Y. Yao, D. Liu et al., “Experimental simulation of the hydraulic fracture propagation in an anthracite coal reservoir in the southern Qinshui basin, China,” Journal of Petroleum Science and Engineering, vol. 168, pp. 400–408, 2018.

A. B. Yost II, W. K. Overby Jr., D. A. Wilkins, and C. D. Locke, “Hydraulic fracturing of a horizontal well in a naturally fractured reservoir: gas study for multiple fracture design,” in SPE 17759 presented at the SPE gas technology symposium, Dallas, TX, USA, 1988.

P. Liu, Y. Ju, P. G. Ranjith, Z. Zheng, and J. Chen, “Experimental investigation of the effects of heterogeneity and geostress difference on the 3D growth and distribution of hydrofracturing cracks in unconventional reservoir rocks,” Journal of Natural Gas Science and Engineering, vol. 35, pp. 541–554, 2016.

C. A. Tang, W. T. Yang, Y. F. Fu, and X. H. Xu, “A new approach to numerical method of modelling geological processes and rock engineering problems -- continuum to discontinuum and linearity to nonlinearity,” Engineering Geology, vol. 49, no. 3-4, pp. 207–214, 1998.

T. Li, L. Li, C.’a. Tang et al., “A coupled hydraulic-mechanical-damage geotechnical model for simulation of fracture propagation in geological media during hydraulic fracturing,” Journal of Petroleum Science and Engineering, vol. 173, pp. 1390–1416, 2019.

L. C. Li, C. A. Tang, and G. Li, “Numerical simulation of 3D hydraulic fracturing based on an improved flow-stress-damage model and a parallel FEM technique,” Rock Mechanics and Rock Engineering, vol. 45, pp. 801–818, 2012.

G. Huang, The Study on the Technologies for Large Scale Fracturing in G892 and F154 Blocks, China University of Petroleum, 2012.

Y. Liu, J. Chen, S. Gao, and S. Peng, Drilling Fluid Technology for the Slim-Hole Horizontal Well of Fan154-Ping, Inner Mongolia Petrochemical Industry, 2012.

M. J. J. Mayerhofer, E. P. P. Lolon, N. R. R. Warpinski, C. L. L. Cipolla, D. Walser, and C. M. M. Rightmire, “What is stimulated reservoir volume?” SPE Production & Operations, vol. 25, no. 1, pp. 89–98, 2010.

G. Yu and R. Aguilera, “3D analytical modeling of hydraulic fracturing stimulated reservoir volume,” in SPE Latin America and Caribbean Petroleum Engineering Conference, pp. 1–17, Mexico City, Mexico, 2012.

G. E. P. Box and D. W. Behnken, “Some new three level designs for the study of quantitative variables,” Technometrics, vol. 2, no. 4, pp. 455–475, 1960.

R. Rickman, M. Mullen, E. Pertre, B. Grieser, and D. Kundert, “A practical use of shale petrophysics for stimulation design optimization: all shale plays are not clones of the Barnett shale,” in Proc. the SPE annual technical conference and exhibition, Denver, CO, USA, 2008.

G. E. King, “Thirty years of gas shale fracturing: what have we learned?” in proceedings of paper SPE 133456 presented at the SPE annual technical conference and exhibition, Florence, Italy, 2010.

N. Nagel, I. Gil, M. Sanchez-Nagel, and B. Damjanac, “Simulating hydraulic fracturing in real fractured rock -overcoming the limits of pseudo 3D models,” in proceedings of SPE 140480, presented at the SPE HFTC, Woodlands. TX, USA, 2011.

Y. Wang, X. Li, R. Hu, C. Ma, Z. Zhao, and B. Zhang, “Numerical evaluation and optimization of multiple hydraulically fractured parameters using a flow-stress-damage coupled approach,” Energies, vol. 9, no. 5, p. 325, 2016.

L. J. L. Beugelsdijk, C. J. de Pater, and K. Sato, “Experimental hydraulic fracture propagation in a multi-fractured medium,” in SPE 59419 presented at the Asia Pacific conference on integrated modeling for asset management, Yokohama, Japan, 2000.

Y. Yao, “Linear elastic and cohesive fracture analysis to model hydraulic fracture in brittle and ductile rocks,” Rock Mechanics and Rock Engineering, vol. 45, no. 3, pp. 375–387, 2012.

J. Zhou, L. Zhang, A. Braun, and Z. Han, “Numerical modeling and investigation of fluid-driven fracture propagation in reservoirs based on a modified fluid-mechanically coupled model in two-dimensional particle flow code,” Energies, vol. 9, no. 9, p. 699, 2016.

L. Zhang, Z. Qu, M. Zhai, A. Li, B. Huang, and Z. Zhang, “Study on effects of low viscosity fluid pre-injection on hydraulic fracture complexity in tight heterogeneous glutenite reservoirs,” Geotechnical and Geological Engineering, vol. 39, no. 2, pp. 839–854, 2021.

W. Mu, L. Li, T. Yang, G. Yu, and Y. Han, “Numerical investigation on a grouting mechanism with slurry-rock coupling and shear displacement in a single rough fracture,” Bulletin of Engineering Geology and the Environment, vol. 78, no. 8, pp. 6159–6177, 2019.