Uncertainty Analysis of Factors Influencing Stimulated Fracture Volume in Layered Formation

Hydraulic fracture dimension is one of the key parameters affecting stimulated porous media. In actual fracturing, plentiful uncertain parameters increase the difficulty of fracture dimension prediction, resulting in the difficulty in the monitoring of reservoir productivity. In this paper, we established a three-dimensional model to analyze the key factors on the stimulated reservoir volume (SRV), with the response surface method (RSM). Considering the rock properties and fracturing parameters, we established a multivariate quadratic prediction equation. Simulation results show that the interactions of injection rate (Q), Young’s modulus (E) and permeability coefficient (K), and Poisson’s ratio (μ) play a relatively significant role on SRV. The reservoir with a high Young’s modulus typically generates high pressure, leading to longer fractures and larger SRV. SRV reaches the maximum value when E1 and E2 are high. SRV is negatively correlated with K1. Moreover, maintaining a high injection rate in this layered formation with high E1 and E2, relatively low K1, and μ1 at about 0.25 would be beneficial to form a larger SRV. These results offer new perceptions on the optimization of SRV, helping to improve the productivity in hydraulic fracturing.

Fracture parameters are closely related to reservoir production in a low permeable formation. Micro-seismic [17] and numerical simulation [18] of hydraulic fracturing indicate that geometric and hydromechanical characteristics of fractures possess the ability to change the well pressure and transient rate performance [19][20][21]. Successful fracturing (large stimulated reservoir volume) leads to greater cumulative production within a short production time [22,23].

Theoretical Background
Basic rock deformation-seepage coupling runs through the whole process of hydraulic fracturing. The fluid flow in the rock causes the change of pore pressure, which is acting on the rock matrix deforming the rock; thus, the subsequent structural deformation, in turn, affects the seepage flow. In this paper, it is assumed that the plastic evolution of reservoir rock as porous media follows the Mohr-Coulomb criterion, and the pores of the rock are completely saturated and filled with incompressible fluid, so the mechanical equilibrium equation of the rock deformation is: Fluid flow continuity equation:

Theoretical Background
Basic rock deformation-seepage coupling runs through the whole process of hydraulic fracturing. The fluid flow in the rock causes the change of pore pressure, which is acting on the rock matrix deforming the rock; thus, the subsequent structural deformation, in turn, affects the seepage flow. In this paper, it is assumed that the plastic evolution of reservoir rock as porous media follows the Mohr-Coulomb criterion, and the pores of the rock are completely saturated and filled with incompressible fluid, so the mechanical equilibrium equation of the rock deformation is: (1) Fluid flow continuity equation: where V is the integral space; S is the integral space surface; σ is the effective stress matrix; p w is the pore pressure; I is the matrix vector for the unit; δε is the strain matrix; T is the region outside the surface force; δv is the virtual velocity matrix; f is the unit volume force without regard to fluid gravity; ϕ is the porosity; ρ w is the pore fluid density; n represents the direction parallel to the normal of the outer surface of the integral; v w is the velocity of fluid flow between rock pores; t is the computing time.
Coupling of fracture propagation and fluid flow is a typical nonlinear behavior; thus, reservoir stimulation is a complex process. The cohesive zone model (CZM), which overcomes the stress singularity before the crack tip, is widely used in the nucleation and growth of hydro-fractures in porous elastic formations [37,38]. The constitutive relation in the bond zone reflects the separation and slip of traction and fault surface, and CZM can be adopted using the finite element method (FEM) to investigate hydraulic fracture propagation [39][40][41], including the fracture tip material softening.
whereσ n is the normal stress component at the specified point; σ 0 n is critical tensile stress in the normal direction of the element, which represents the tensile strength of reservoir rock; τ 1 and τ 2 are the shear stress at specified point; τ 0 1 and τ 0 2 are the critical shear stress. Cohesion is a function of the relative displacement of the upper and lower surface of the crack and is controlled by the interface damage factor. When D is 0, the unit exhibits no damage; when D is 1, the crack opens completely. The cohesive element method adopts the damage factor D to describe the crack extension: where d f m is the displacement value corresponding to the complete failure of the element; d max m is the maximum displacement of the element; d 0 m is the displacement value corresponding to the occurrence of damage at the beginning of the element.
As the propagation, the process of hydraulic fracturing is usually accompanied by the shear slip effect, and the crack mode is a composite crack. When the element is completely damaged, it means that the cohesion disappears completely, and at the same time, it marks the generation of a new macroscopic crack surface, and the energy consumed in this damage process is the fracture energy. When describing the propagation behavior of composite fractures, the Benzeggagh and Kenane [42] is applied to give the criteria of the critical energy release rate of fracture propagation: where the B-K criterion assumes that the first and second tangential equals the fracture energy release rate, namely G C s = G C t ; G C n is the critical strain energy release rate of normal fracture; η is the constants relating to the properties of the material itself; G C is critical fracture energy release rate for composite fractures.
Fracturing is hydraulically driven, in which the high-pressure fluid in the wellbore is injected into the fracture. The fracturing fluid pressure acting on the fracture surface is the driving force of the fracture expansion. The flow within the cohesive pore pressure element includes a tangential flow along the fracture surface and normal flow along the vertical fracture surface (shown in Figure 2). The mathematic fracture tip corresponds to the point that is about to separate; the cohesive fracture tip is the damage initiation point where the separation reaches the critical value (δ 0 ) or the traction reaches the cohesive strength; the material fracture tip refers to the complete failure point where the separation reaches the critical value (δ f ), and the traction acting across the unbroken cohesive zone vanishes [43]. The fracturing fluid is regarded as a Newtonian fluid, and the tangential flow on the fracture surface can be written as: where is the injection rate; is the flow coefficient; ∇ is the fluid pressure gradient along the cohesive zone.
According to the Reynolds number equation, the flow coefficient can be expressed as: where is the fracture width, the fluid viscosity coefficient. The filtration coefficient is set to form a permeable layer on the fracture surface, and the normal percolation of the fracturing fluid on the fracture surface can be expressed as: where is the interfacial fluid pressure of cohesive fracture element; and are, respectively, the seepage flow rate of the fluid on the upper and lower surfaces of the crack; and are, respectively, the fluid loss coefficients on the upper and lower surfaces of the crack; and are, respectively, the fluid pore pressure above and below the fracture surface.
The cohesive fracture zone is filled with fluid, and the fracture of the rock does not produce a traction effect, but the fluid pressure acts on the open fracture surface.
( ) is the injection rate. Substituting Equations (7)-(10) into Equation (11) results in the Reynold's lubrication equation: The fracturing fluid is regarded as a Newtonian fluid, and the tangential flow on the fracture surface can be written as: where q is the injection rate; k t is the flow coefficient; ∇p is the fluid pressure gradient along the cohesive zone. According to the Reynolds number equation, the flow coefficient k t can be expressed as: where d is the fracture width, µ the fluid viscosity coefficient. The filtration coefficient is set to form a permeable layer on the fracture surface, and the normal percolation of the fracturing fluid on the fracture surface can be expressed as: where p i is the interfacial fluid pressure of cohesive fracture element; q t and q b are, respectively, the seepage flow rate of the fluid on the upper and lower surfaces of the crack; c t and c b are, respectively, the fluid loss coefficients on the upper and lower surfaces of the crack; p t and p b are, respectively, the fluid pore pressure above and below the fracture surface.
The cohesive fracture zone is filled with fluid, and the fracture of the rock does not produce a traction effect, but the fluid pressure acts on the open fracture surface.
Q(t) is the injection rate. Substituting Equations (7)-(10) into Equation (11) results in the Reynold's lubrication equation: Equation (12) is used to describe the fluid flow in the fracture plane. It can be seen from the equation that there is a pressure-pull-separation relationship between the cohesive fracture zone and the compressed fracture.

Model Setup
The target formation, in central-western China, is composed of shale, mudstone, and sandstone formed by sedimentation. The upper lithology is mainly mudstone with the inclusion of siltstone, and the lithology of the middle lower part is mainly mudstone, sandstone, and shale. Based on logging data and mechanical rock tests, to study the hydraulic response of SRV to vertical fracture propagation and clarify its key influencing factors, a fully three-dimensional model was established (Figure 3). The overburden of the model was composed of mudstone with a depth range of 1000 to 1020.06 m, while the reservoir was composed of sandstone with a depth range of 1020.06 to 1029.94 m, and the vertical well fracturing was performed in a reservoir. An interval with a depth of 1029.94 to 1050 m was mudstone.
The hydraulic fracturing numerical model follows the following assumptions: (1) The formation is isotropic pore material.
(2) Initial stress gradient is evenly distributed in each layer.
(3) Interface slippage between reservoir and barrier is not considered.
(4) Fluid in the pores is incompressible.
The unopened vertical fracture surface is described by a nested zero-thickness cohesive zone with the same stiffness as the surrounding rocks, and no initial separation occurred. A total of 19,500 pore pressure structural elements and 23,409 nodes were the basis of the whole numerical model analysis, and the mesh was encrypted around the perforation. The geo-stress field acted on the whole formation model, fixing the lower boundary and constraining the normal displacement of the remaining boundary. The pore pressure boundary condition was set to a certain value of 20 MPa. The minimum horizontal principal stress and the maximum horizontal principal stress were 8 and 12 MPa, respectively. The vertical stress was equal to 10 MPa. The fracturing fluid viscosity was 0.001 Pa·s. Considering the cohesion and internal friction angle of the mudstone layer and sandstone layer, we adopted Mohr-Coulomb failure criterion to analyze the plastic failure of the rock and calculate the percolation-stress failure process of the rock. The basic parameters used for the numerical simulation are shown in Table 1.
The whole hydraulic fracturing analysis consists of four steps. Firstly, in order to eliminate the displacement caused by gravity and initial conditions, the boundary and initial stress conditions are applied to achieve the geo-stress balance. The results are used in the next step of the hydraulic fracturing calculation. Then, fluid is injected into the central perforation in the reservoir of the model for 100 s: the injection rate increases linearly to 1 × 10 −2 m 3 /s within 30 s. In the next step, the injection is stopped, and the hydraulic fracture keeps open to simulate proppant behavior, with this stage lasting for 10 s. The final stage totals 1200 s and aims to simulate the fluid back through the fracture, and a production pressure difference of 20 kPa is applied to the perforation node. After stopping the injection, the accumulative high pore pressure in the fracture bleed off into the formation.

Experimental Design
SRV increases with the optimal combination of various parameters. To study the influence rules of several factors on SRV, the response surface method (RSM) and Box-Behnken Design (BBD) were used to establish multiple quadratic regression models with SRV as a response value.

Response Surface Method
Response surface analysis is a mathematical method based on statistics, which is commonly used in modeling and the analysis of engineering problems [44,45]. Based on the acquired response surface, RSM can quantify the significance of each index and its interaction [46]. In addition, each index and its interaction relationship can be directly displayed by a 3D response surface diagram [33].
SRV was taken as an optimization objective to search for corresponding parameter conditions, and the reasonable experimental design method was adopted to obtain certain data through experiments. The multivariate quadratic regression equation was adopted to fit the functional relationship between factors and response values: whereŶ is the predicted response value; x i is the independent variable; K, K i , K j , and K ij are the parameters to be determined in the second-order polynomial mathematical model; ε is the random error. In order to discuss whether the model based on RSM approximates the real function, an F statistic was constructed to test significance as follows: where S SR is the regression sum of squares; S SE is the residual sum of squares; Y i ,Ŷ i , Y are the test value, predicted value, and test mean of the response at point i, respectively; n is the sample size; p is the number of independent variables.

Box-Behnken Design
The Box-Behnken experimental design is a cyclic second-order, three-level high (high, medium, and low levels), partial factor experimental design method. This cubic design method is characterized by selecting the midpoints on each side of the multidimensional cube and the midpoints of the cube as design experiments to avoid data loss [47]. As shown in Figure 4, it can be seen that the Box-Behnken design is more efficient than the full three-level factorial design. The total number of tests (N) is defined as: N = 2k(k − 1) + C 0 , where k is the number of parameters, and C 0 is the number of center points. This method can not only overcome the limitations of orthogonal experiments but also reduce the number of experiments and save computing time. Parameters of the model are divided into high, medium, and low levels (140%, 100%, and 60%), and the statistical experimental design is shown in Table 2. These parameters include the injection rate (Q), Young's modulus (E1, E2), Poisson's ratio (µ1, µ2), the friction angle (ϕ1, ϕ2), the cohesion (C1, C2), the permeability coefficient (K1, K2), the void ratio (ν1, ν2), the dilatancy angle (ψ1, ψ2), and the viscosity (vis1, vis2), and the corresponding meanings of each parameter are shown in Table 2. The Box-Behnken experimental design is a cyclic second-order, three-level high (high, medium, and low levels), partial factor experimental design method. This cubic design method is characterized by selecting the midpoints on each side of the multidimensional cube and the midpoints of the cube as design experiments to avoid data loss [47]. As shown in Figure 4, it can be seen that the Box-Behnken design is more efficient than the full three-level factorial design. The total number of tests (N) is defined as: = 2 ( − 1) + , where is the number of parameters, and is the number of center points. This method can not only overcome the limitations of orthogonal experiments but also reduce the number of experiments and save computing time. Parameters of the model are divided into high, medium, and low levels (140%, 100%, and 60%), and the statistical experimental design is shown in Table 2. These parameters include the injection rate (Q), Young's modulus (E1, E2), Poisson's ratio (μ1, μ2), the friction angle (φ1, φ2), the cohesion (C1, C2), the permeability coefficient (K1, K2), the void ratio (ν1, ν2), the dilatancy angle (ψ1, ψ2), and the viscosity (vis1, vis2), and the corresponding meanings of each parameter are shown in Table 2.

Results
At the beginning of fracturing in the benchmark model, the injection rate of the fracturing fluid at the perforation rapidly reaches a high level; thus, rock element damage occurs after the pore water pressure rises to the fracture pressure of 39.8 MPa ( Figure 5). Pressure waves travel through the reservoir, causing pressure and rate responses [48], which affect the surrounding rock stress and pore pressure in the coupling equation. As the fracture extends upwards and forwards, the fracture height reaches stability outside the reservoir after the fracture crosses the interface of the thin reservoir, and the range of pore water pressure increases. at the perforation rapidly reaches a high level; thus, rock element damage occurs after the pore water pressure rises to the fracture pressure of 39.8 MPa ( Figure 5). Pressure waves travel through the reservoir, causing pressure and rate responses [48], which affect the surrounding rock stress and pore pressure in the coupling equation. As the fracture extends upwards and forwards, the fracture height reaches stability outside the reservoir after the fracture crosses the interface of the thin reservoir, and the range of pore water pressure increases.
When t = 16.8 s, the fracture height reaches stability outside the reservoir after the fracture crosses the interface of the thin reservoir ( Figure 6). Injection pressure at the fracturing point decreases sharply with the passing of time, and fluctuations appear in the injection pressure curve, but injection pressure gradually stabilizes at about 33.5 MPa. At this time, under the condition of a constant pump injection rate, the fracture begins to expand steadily in the formation, which can be considered as the stable expansion stage of hydraulic fractures [49]. At the end of the fracturing, the resulting fracture size is 32.6 m in length, 13.6 m in height, and 5.8 mm in width. Figure 5. The change in pore pressure during the fluid injection process in the benchmark model. ① Fluid is injected into the central perforation in the reservoir of the model for 100 s, and pore pressure rapidly rises to the breakdown pressure and gradually stabilizes at a stable level. ② Injection is Figure 5. The change in pore pressure during the fluid injection process in the benchmark model. 1 Fluid is injected into the central perforation in the reservoir of the model for 100 s, and pore pressure rapidly rises to the breakdown pressure and gradually stabilizes at a stable level. 2 Injection is stopped, and the hydraulic fracture keeps open, proppant acts at this stage, and then dissipated pore pressure activates the production pressure difference to make the fluid move back through the fracture. The research focus is on the fluid injection stage when SRV is formed. When t = 16.8 s, the fracture height reaches stability outside the reservoir after the fracture crosses the interface of the thin reservoir ( Figure 6). Injection pressure at the fracturing point decreases sharply with the passing of time, and fluctuations appear in the injection pressure curve, but injection pressure gradually stabilizes at about 33.5 MPa. At this time, under the condition of a constant pump injection rate, the fracture begins to expand steadily in the formation, which can be considered as the stable expansion stage of hydraulic fractures [49]. At the end of the fracturing, the resulting fracture size is 32.6 m in length, 13.6 m in height, and 5.8 mm in width. stopped, and the hydraulic fracture keeps open, proppant acts at this stage, and then dissipated pore pressure activates the production pressure difference to make the fluid move back through the fracture. The research focus is on the fluid injection stage when SRV is formed.

Hurricane Analysis of Multiple Factors
Before the response surface analysis, independent factor analyses were carried out for multiple parameters of reservoir and barrier. Geo-mechanical parameters (E, μ, C, φ, ψ), hydromechanical parameters (K, v and vis), and operational parameters (Q) were taken into account to study the influence of each factor. Thirty-four simulation results of SRV indicated that SRV formed by changing diverse variables is dissimilar, as shown in Figure 7. From the hurricane analysis result, the key

Hurricane Analysis of Multiple Factors
Before the response surface analysis, independent factor analyses were carried out for multiple parameters of reservoir and barrier. Geo-mechanical parameters (E, µ, C, ϕ, ψ), hydromechanical parameters (K, v and vis), and operational parameters (Q) were taken into account to study the influence of each factor. Thirty-four simulation results of SRV indicated that SRV formed by changing diverse variables is dissimilar, as shown in Figure 7. From the hurricane analysis result, the key influencing factors are Q, E2, E1, K1, µ2, µ1, and K2. The BBD experiment was designed for these seven key parameters. influencing factors are Q, E2, E1, K1, μ2, μ1, and K2. The BBD experiment was designed for these seven key parameters.

Response Surface Analysis of Variables to SRV
The interaction between the factors shall not be underestimated [50]. Due to the limitation of the single factor experiments, the continuous point cannot be analyzed, and the interaction among various factors is not considered. Thus, we used the response surface method to comprehensively analyze the results.
SRV is regarded as the function of several influencing factors. The numbers of the simulation experiments were designed to obtain response values, which were established for the mathematical model of the response to each factor level (60%, 100%, 140%). From the above-mentioned hurricane analysis results, seven significant factors (Q, E2, E1, K1, μ2, μ1, K2) were considered in the Box-Behnken experiment design. Sixty-two groups of experiments with different combinations of variables are shown in Table 3, along with corresponding SRV simulation results.

Response Surface Analysis of Variables to SRV
The interaction between the factors shall not be underestimated [50]. Due to the limitation of the single factor experiments, the continuous point cannot be analyzed, and the interaction among various factors is not considered. Thus, we used the response surface method to comprehensively analyze the results.
SRV is regarded as the function of several influencing factors. The numbers of the simulation experiments were designed to obtain response values, which were established for the mathematical model of the response to each factor level (60%, 100%, 140%). From the above-mentioned hurricane analysis results, seven significant factors (Q, E2, E1, K1, µ2, µ1, K2) were considered in the Box-Behnken experiment design. Sixty-two groups of experiments with different combinations of variables are shown in Table 3, along with corresponding SRV simulation results.  based on the actual values of the seven parameters, but in the response surface calculation, the equation is expressed as a combination of coded factors with different coefficients. The high levels of the factors are coded as +1, and the low levels are coded as −1. Parameter coefficients and their standard errors are shown in Figure 8. Comparing parameter coefficients is a supplementary means to identify the relative impact of the seven parameters. In order to discuss the significance of all the model terms, the F-value and p-value were introduced to statistical analysis, as shown in Figure 9. A p-value less than 0.05 implies that the model terms are significant. The statistical results of the F-value and p-value indicate that A, B, C, AB, AC, C 2 are the most significant model terms.

Discussion
The established response surface equation consists of seven variables, which are geo-mechanical parameters (E and µ), a geo-physical parameter (K), and an operational parameter (Q). During hydraulic fracturing, the damage and failure mechanisms of the rock are controlled by multiple co-effective properties [51]. Fracture height, width, and injection pressure are sensitive to the variation of material parameters of the reservoir and the barrier (E, µ, and K) [52][53][54]. The effect of variables and their interaction on SRV are emphatically discussed in this section.

Effect of the Injection Rate
Variations of fracture height, width, and SRV with Q are studied, as shown in Figure 11. When the pump rate increases, the net pressure in hydraulic fractures will increase, resulting in more evident perturbation of formation stress and more fracturing fluid leak-off [55][56][57], and the stimulated reservoir volume will increase. Our statistical results highlight that the injection rate is the most significant factor towards SRV, which is consistent with other researches [34,58,59]. A high pumping rate provides enough energy to overcome resistance, resulting in a larger fracture height and width [11]. Thus, when the pumping rate is low, the hydraulic fracture is incapable of expanding with a wide range. As the injection rate increased to 0.012 m 3 /s, the created fracture height reaches the maximum level ( Figure 11a). However, due to the fracture length at this time being less than the length of the condition of the injection rate of 0.014 m 3 /s, the value of SRV is not at its maximum. Moreover, a fracture occurs when the local pore pressure at the perforation point reaches the breakdown pressure. When the injection rate is high, the volume of fluid pumped into the fracture per unit time increases, and the pore pressure and skeleton stress at the tip and sides of the hydraulic fracture becomes larger, which causes the breakdown pressure to become higher (Figure 11b).

Discussion
The established response surface equation consists of seven variables, which are geo-mechanical parameters (E and μ), a geo-physical parameter (K), and an operational parameter (Q). During hydraulic fracturing, the damage and failure mechanisms of the rock are controlled by multiple coeffective properties [51]. Fracture height, width, and injection pressure are sensitive to the variation of material parameters of the reservoir and the barrier (E, μ, and K) [52][53][54]. The effect of variables and their interaction on SRV are emphatically discussed in this section.

Effect of the Injection Rate
Variations of fracture height, width, and SRV with Q are studied, as shown in Figure 11. When the pump rate increases, the net pressure in hydraulic fractures will increase, resulting in more evident perturbation of formation stress and more fracturing fluid leak-off [55][56][57], and the stimulated reservoir volume will increase. Our statistical results highlight that the injection rate is the most significant factor towards SRV, which is consistent with other researches [34,58,59]. A high pumping rate provides enough energy to overcome resistance, resulting in a larger fracture height and width [11]. Thus, when the pumping rate is low, the hydraulic fracture is incapable of expanding with a wide range. As the injection rate increased to 0.012 m 3 /s, the created fracture height reaches the maximum level ( Figure 11a). However, due to the fracture length at this time being less than the length of the condition of the injection rate of 0.014 m 3 /s, the value of SRV is not at its maximum. Moreover, a fracture occurs when the local pore pressure at the perforation point reaches the breakdown pressure. When the injection rate is high, the volume of fluid pumped into the fracture per unit time increases, and the pore pressure and skeleton stress at the tip and sides of the hydraulic fracture becomes larger, which causes the breakdown pressure to become higher (Figure 11b).

Effect of the Young's Modulus
Sensitivity analyses for the fracture profile and SRV of Young's modulus are shown in Figures  12 and 13. The interface between strata can produce a barrier effect on the fracture height [60]. When Young's modulus of the reservoir (E1) is high, the layered formation model tends to be hard, which limits the width and height of the fracture (Figure 12a). The main reason for this control is the

Effect of the Young's Modulus
Sensitivity analyses for the fracture profile and SRV of Young's modulus are shown in Figures 12  and 13. The interface between strata can produce a barrier effect on the fracture height [60]. When Young's modulus of the reservoir (E1) is high, the layered formation model tends to be hard, which limits the width and height of the fracture (Figure 12a). The main reason for this control is the coupling effect of the change in crack width and net fracture pressure [61]. Under the condition of a constant injection rate, the net pressure at the crack tip is determined by the fracture width. When Young's modulus is low, the fracture width in the reservoir and the fluid flow resistance is small. Thus, the corresponding net treatment pressure is low, the fracture length is short, and the resulting fracture volume is small (Figure 12b). The formations with high Young's modulus typically generate high pressure, leading to longer fractures.

Effect of the Poisson's Ratio
T-stress at the crack tip and Poisson's ratio affect the three typical types of fractures in fracture mechanics. Fracture parameters, including stress intensity factor and t-stress, are related to Poisson's ratio.
T-stress can be defined as: where is the thickness strain component, I is the interaction integral, and f is the point force [63]. is positively correlated with μ.
The fracture morphology and SRV with the change of μ are shown in Figure 14. When the reservoir Poisson's ratio (μ1) is low, the fracture is wide (Figure 14a), but the SRV reaches the maximum and is negatively correlated with μ1 ( Figure 14b). That is because the material stress intensity factor and T-stress decrease with the decrease of the Poisson's ratio. Once the tip of the fracture crosses the interface, the high barrier Poisson's ratio (μ2) has a relatively weak effect on inhibiting the fracture growth (Figure 14c), and consequently, increasing the fracture extension range (Figure 14d). Lower breakdown pressure leads to easier initiation of hydraulic fractures. Combined with Figures 14 and 15, when μ1 is 0.25, the breakdown pressure is relatively low, and as a result, SRV is optimized. Fracture profiles with the change in Young's modulus of the barrier (E2) are shown in Figure 12c. The effective Young's modulus decreased after the fracture penetrated the material interface [62], and the net fracture pressure positively correlated with the effective Young's modulus.
The stress intensity factor is related to the net crack pressure: where P net is the net fracture pressure, and H is the fracture height. The decrease of P net results in a lower stress intensity factor. Thus, after the penetration of the fracture through the barrier, P net decreases with the decrease in effective Young's modulus. The fracture length in the reservoir increases with the increase in E2, resulting in larger SRV (Figure 12d), but there is no significant change in breakdown pressure.
In general, the fracture width depends on the joint action of Young's modulus of reservoir and barrier. Our statistical results show that SRV reaches the maximum value when E1 and E2 are high, which indicates that E1 and E2 are positively correlated with SRV ( Figure 13).

Effect of the Poisson's Ratio
T-stress at the crack tip and Poisson's ratio affect the three typical types of fractures in fracture mechanics. Fracture parameters, including stress intensity factor and t-stress, are related to Poisson's ratio.
T-stress can be defined as: where ε zz is the thickness strain component, I is the interaction integral, and f is the point force [63].
T is positively correlated with µ.
The fracture morphology and SRV with the change of µ are shown in Figure 14. When the reservoir Poisson's ratio (µ1) is low, the fracture is wide (Figure 14a), but the SRV reaches the maximum and is negatively correlated with µ1 (Figure 14b). That is because the material stress intensity factor and T-stress decrease with the decrease of the Poisson's ratio. Once the tip of the fracture crosses the interface, the high barrier Poisson's ratio (µ2) has a relatively weak effect on inhibiting the fracture growth (Figure 14c), and consequently, increasing the fracture extension range (Figure 14d). Lower breakdown pressure leads to easier initiation of hydraulic fractures. Combined with Figures 14 and 15, when µ1 is 0.25, the breakdown pressure is relatively low, and as a result, SRV is optimized.

Effect of the Permeability Coefficient
Matrix porosity and permeability are sensitive to effective stress and pore pressure [64]. Figure  16 shows the fracture morphology and changes of SRV with the permeability coefficient. When the

Effect of the Permeability Coefficient
Matrix porosity and permeability are sensitive to effective stress and pore pressure [64]. Figure 16 shows the fracture morphology and changes of SRV with the permeability coefficient. When the reservoir permeability coefficient (K1) is low, the formed fracture is relatively narrow (Figure 16a). This is because the lower the permeability, the greater the fluid flow resistance in the formation, and the lower openness the fracture will develop. However, SRV is negatively correlated with K1 (Figure 16b). The mechanism is that the lower the permeability, the easier it is for fluid to flow in the fracture, and the more fluid remains in the fracture rather than being lost in the formation. Thus, a larger fracture volume can be formed.  The barrier permeability coefficient (K2) is an order of magnitude smaller than K1; thus, the fluid is mainly in the fracture rather than in the barrier. Increasing K2 can reduce fluid resistance in the barrier, leading to an appropriate increase in the fracture length and width (Figure 16c). So, a larger SRV is formed with the increase of K2 (Figure 16d). The simulation result, as shown in Figure 17, indicates that when K2 is greater than 1.1 × 10 −8 m/s, SRV increases with the increase of K1. That is because K2 and K1 almost remain at the same order of magnitude, the interlayer barrier effect in layered formation is no longer significant, and the fracture length in the reservoir increases.

Conclusions
In this paper, the evolution of the morphological characteristics of hydraulic fractures during the whole process of hydraulic fracturing were studied. Using the cohesive zone method, we investigated the effect of each characteristic parameter on SRV and established a quadratic model. The main conclusions are as follows: (1) To screen the significant influencing factors and determine the optimal level combination, the quadratic response surface model was established. The statistical analysis of the F value and pvalue of regression equations implies that the model has a high degree of fitting, which is beneficial for analyzing and predicting the hydraulic behavior of fracturing in a layered formation.
(2) The injection rate and Young's modulus are the most significant factors to SRV among all the independent variables. A high pumping rate provides enough energy to form a hydraulic fracture with a wide range. The reservoir with high Young's modulus typically generates high pressure, leading to longer fractures and larger SRV. SRV reaches the maximum value when E1 and E2 are high, which indicates that E1 and E2 are positively correlated with SRV.
(3) When the reservoir Poisson's ratio (μ1) is low, the fracture is wide, but the SRV reaches the maximum and is negatively correlated with μ1. Once the tip of the fracture crosses the interface, the high barrier Poisson's ratio (μ2) has a relatively weak effect on inhibiting the fracture growth, and consequently, increasing the fracture extension range. Moreover, the lower the reservoir permeability coefficient (K1), the more fluid remains in the fracture rather than being lost in the formation; thus, SRV is negatively correlated with K1.
(4) Different interactions affect the SRV to varying degrees. In general, maintaining a high injection rate in this layered formation with high E1 and E2, and relatively low K1 and μ1 at about 0.25, would be beneficial to form a larger SRV.

Conclusions
In this paper, the evolution of the morphological characteristics of hydraulic fractures during the whole process of hydraulic fracturing were studied. Using the cohesive zone method, we investigated the effect of each characteristic parameter on SRV and established a quadratic model. The main conclusions are as follows: (1) To screen the significant influencing factors and determine the optimal level combination, the quadratic response surface model was established. The statistical analysis of the F value and p-value of regression equations implies that the model has a high degree of fitting, which is beneficial for analyzing and predicting the hydraulic behavior of fracturing in a layered formation. (2) The injection rate and Young's modulus are the most significant factors to SRV among all the independent variables. A high pumping rate provides enough energy to form a hydraulic fracture with a wide range. The reservoir with high Young's modulus typically generates high pressure, leading to longer fractures and larger SRV. SRV reaches the maximum value when E1 and E2 are high, which indicates that E1 and E2 are positively correlated with SRV. (3) When the reservoir Poisson's ratio (µ1) is low, the fracture is wide, but the SRV reaches the maximum and is negatively correlated with µ1. Once the tip of the fracture crosses the interface, the high barrier Poisson's ratio (µ2) has a relatively weak effect on inhibiting the fracture growth, and consequently, increasing the fracture extension range. Moreover, the lower the reservoir permeability coefficient (K1), the more fluid remains in the fracture rather than being lost in the formation; thus, SRV is negatively correlated with K1. (4) Different interactions affect the SRV to varying degrees. In general, maintaining a high injection rate in this layered formation with high E1 and E2, and relatively low K1 and µ1 at about 0.25, would be beneficial to form a larger SRV.