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

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.


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][4][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][8][9][10][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 2 , 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][18][19][20][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][23][24][25][26].
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][32][33][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

Geofluids
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.

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 2 , 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. 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 RFPA3Dpetrol software [35][36][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 u is the mechanical parameters of rock obtained from laboratory tests, u 0 is the mean value of the variable u, and m is the degree of homogeneity, which is called the homogeneity index. A larger m indicates a more homogeneous material and vice versa. f ðuÞ is the distribution density of the variable u.
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 σ ij and ε ij are the stress tensor and strain tensor of the solid, respectively; ε v is the volume strain; F i and U i are the body force and displacement components, respectively; σ e ij is the effective stress; α is the pore pressure coefficient; p is the pore pressure; δ ij is the Kronecker delta; λ is the Lame constant; G is the shear modulus of the solid; Q is Biot's coefficient; k is the permeability coefficient that varies with the stress-damage evolution; k 0 is the initial permeability coefficient; ξðξ ≥ 1Þ 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.

Geofluids
As the damage evolves, Young's modulus of the element gradually decreases, as defined below: where D is the damage variable and E 0 and E are Young's modulus of the undamaged and damaged elements, respectively. When the tensile stress in the element reaches its tensile strength, σ 3 ≤ −f t0 , the variable D of the element under uniaxial tension can be expressed as where ε t0 is the tensile strain at the elastic limit, as shown in Figure 2(a). ε tu is the ultimate tensile strain, which is defined as ε tu = ηε t0 . η is the ultimate strain coefficient. λ is the residual tensile strength coefficient. When the equivalent principal tensile strain ε reaches the threshold strain,ε t0 , under multiaxial stress states, the element is assumed to be damaged in the tensile mode. ε is defined as follows: where ε 1 , ε 2 , and ε 3 are the three principal strains. < > is a function defined as The variable D 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 σ 1 and σ 3 are the principal stresses. φ is the internal friction angle. f c0 is the uniaxial compressive strength. The variable D under uniaxial compression can be expressed as where ε c0 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 ε c0 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 D.
The fracture aperture is represented as a linear function of the net pressure in the RFPA3D-petrol model: where w is fracture aperture, p net is net pressure in the HF, and h f is fracture height.

Model Setup.
A field-scale numerical model is established based on RFPA3D-petrol software as shown in  Table 1. In this simulation, 2 barrier layers share the same mechanical parameters.
The dimension of the reservoir model is 420 m × 200 m × 60 m, and the model is divided into 5040000 FEM units. Its mesh size, fracturing parameters, and boundary conditions are listed in Table 2    Vertical stress σ z/ σ v MPa 72 5 Geofluids 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 3 /min, and the total fluid volume is 480 m 3 . The mechanical parameters of different layers are derived from experiment tests of 121 core samples drilled from well Y-11 [38]. The homogeneity index m 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].

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 D = 1) and damage elements (damage variable 0 < D < 1). The matrix is represented by elastic elements (D = 0 ). 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.
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 3 . In our numerical model, the hydraulic fracture is represented by full failure elements (damage variable D = 1) and damage elements (damage variable 0 < D < 1). 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.
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  1.111e+002 Figure 4: The X-shaped HF propagation process obtained from numerical simulation.  [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   Table 3.
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 E max and E min are the maximum and minimum values of Young's modulus E, respectively. ν max and ν min 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: E max = 55 MPa, E min = 25 MPa, ν max = 0:30, and ν max = 0:15.
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. 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.

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 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.
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 Y is the response value, x i and x ij are the uncertain variables, and ω,ω i ,ω ii , and ω ij are the coefficients to be determined of the RSM.
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 SRV = 2:26e 5 where A is the brittleness index, B is tensile strength (MPa), C is horizontal stress difference (MPa), D is the injection rate (m 3 /min), and E is fluid viscosity (mPa·s).

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. The plot of "predicted SRV versus actual SRV" is shown in Figure 10 Table 4.

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 F value of 47.3 indicates that the model is significant. In statistics, p value was usually used to evaluate the significance of parameters. Taking 5% as the threshold before testing, when the p 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 p value is less than 0.001, the parameters are extremely significant. When the p value was more than 0.05, the parameters are not significant. In this case, factors of A, C, D, and E are significant terms. The influential order of these five factors is C ðhorizontal stress differenceÞ > D ðinjection rateÞ > E ðfluid viscosityÞ > A ðbrittleness indexÞ > B ðtensile strengthÞ. 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.

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    Figures 11-14. The response surfaces all exhibit curved shapes, indicating the nonlinear influential mechanism of these factors on SRV.   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 3 /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 3 /min to 12 m 3 /min, the SRV increases from 285000 m 3 to 421000 m 3 at a horizontal stress difference of 5 MPa, while the SRV increases from 160000 m 3 to 234000 m 3 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 sandcarrying 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 hor-izontal 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    [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.

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. 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.
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.

Geofluids
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  16 Geofluids 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.

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 3 , 129000 m 3 , 179000 m 3 , 84000 m 3 , and 151000 m 3 , 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.

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 horizontal stress difference > injection rate > fluid viscosity > brittleness index > tensile strength (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.