Numerical Evaluation and Optimization of Multiple Hydraulically Fractured Parameters Using a Flow-Stress-Damage Coupled Approach

Multiple-factor analysis and optimization play a critical role in the the ability to maximizethe stimulated reservoir volume (SRV) and the success of economic shale gas production. In this paper, taking the typical continental naturally fractured silty laminae shale in China as anexample, response surface methodology (RSM) was employed to optimize multiple hydraulic fracturing parameters to maximize the stimulated area in combination with numerical modeling based on the coupled flow-stress-damage (FSD) approach. This paper demonstrates hydraulic fracturing effectiveness by defining two indicesnamelythe stimulated reservoir area (SRA) and stimulated silty laminae area (SLA). Seven uncertain parameters, such as laminae thickness, spacing, dip angle, cohesion, internal friction angle (IFA), in situ stress difference (SD), and an operational parameter-injection rate (IR) with a reasonable range based on silty Laminae Shale, Southeastern Ordos Basin, are used to fit a response of SRA and SLA as the objective function, and finally identity the optimum design under the parameters based on simultaneously maximizingSRA and SLA. In addition, asensitivity analysis of the influential factors is conducted for SRA and SLA. The aim of the study is to improve the artificial ability to control the fracturing network by means of multi-parameteroptimization. This work promises to provide insights into the effective exploitation of unconventional shale gas reservoirs via optimization of the fracturing design for continental shale, Southeastern Ordos Basin, China.


Introduction
Marine shale gas has recently been explored and has achieved significant success in the USA, which has triggered a worldwide fever for shale gas exploration and development.Recently, more and more attention has been paid to the continental shale hydrocarbons in China.The typical continental shale plays in China are mainly distributed at the southeast Ordos Basin.The Ordos Basin covers 26 ˆ104 km 2 and is a huge, hydrocarbon prolific basin located in the middle of northern China, possessing giant gas fields in the Upper Palaeozoic and oil fields in the Ordovician, Triassic, and Jurassic strata [1][2][3].The main stratigraphy of the basin is the Yanchang formation, which consists of the upper Triassic strata, can be subdivided into ten members, named Member Chang1 through Chang10 [4].All members are composed of sandstones, shales, and mudstone, among which Member Chang7 of the Yanchang Formation is composed of shales predominantly interbedded with silt sandstones, characterized by a typical laminae structure [5][6][7][8][9][10].Chang 7 is an important exploration target for conventional shale gas in the basin.Continental shales in Ordos Basin China are generally characterized by high clay content and low thermal maturity characteristics that are significantly different from those of marine shales documented in the USA and elsewhere [4,11].It is really a question that whether such continental shales have similar commercial hydrocarbon potential as the marine shale in the USA or not.However, in the Yanchang formation, the lacustrine shales have produced an average of 1000-3000 m 3 of gas and two to five tons of oil per day by hydraulic fracturing treatment.This proves that the continental shale formations have a potential exploration prospect.As a result, it is in urgent need to study the continental shales and combination of horizontal drilling and multistage hydraulic fracturing technology to gas production from the gas shale Yanchang formation to flourish.
Current studies are mainly focused on the occurrence of silty laminae and the depositional environments [12], hydrodynamic conditions [13], origin of silty laminae [14][15][16], the shale gas generation, and distribution [17,18].For the Shale of the Upper Triassic Yanchang Formation, Ordos Basin, China, Lei et al. [10] have conducted a large study on the microscopic sedimentary characteristics, mineralogy, pore structure, porosity, and gas and oil content of silty laminae.Tang et al. [4] have studied the kerogen type, thermal maturity, heterogeneity and lithological, porosity, and clay mineral content of Chang 7 and Chang 9. Li et al. [19] have researched the geology characteristics of the Chang 7 reservoir in the Xiasiwan oil field.Zhang et al. [8] proposed a method to determine the oil/gas phase in the Chang 7 shale of the Yanchang formation in the Ordos Basin.Zeng et al. [20] analyzed the gas content of continental Yanchang Shale and the other main controlling factors.Liu et al. [21] used X-ray diffraction analysis to investigate the geochemistry, pore structures and fractal characteristics of 16 shale samples from fourwells in the Ordos Basin, andthe relationships between total organic carbon (TOC), mineralogical compositions and pore structure parameters were alsoinvestigated.From the exploration and experimental results of these scholars, it is obvious that silty laminae forms significant migration conduits and contribute significantly to reservoir storage [10,22,23].Thus, it is urgent to exploit the continental laminated gas shale formation.A decisive factor in the successful economic shale gas production is to maximize the total stimulated reservoir volume.It has brought a new perspective for the unprecedented growth of shale reservoirs, which focus on the optimization of multiple hydraulically fractured parameters in silty laminaed shales.
Numerical modeling techniques can provide us a feasible alternative solution to understand the fracture initiation and propagation processes [24][25][26][27][28][29][30][31][32][33].In this paper, two-dimensional numerical modeling based on the FSD coupled approach (realistic failure process analysis (RFPA)-Flow) was integrated to simulate the hydraulic fracturing process for each case.D-Optimal design was adopted to design the experiment.The response surface methodology (RSM) is employed to build the response surface in terms of stimulated reservoir area (SRA) and stimulated silty laminae area (SLA) with seven parameters including three silty laminae structure parameters (thickness, spacing, dip angle), two mechanical parameters (internal friction angle (IFA), cohesion), one in situ stress parameter (stress difference (SD)) and one operational parameter (injection rate (IR)) to obtain the best economic scenario for a given range of these influential parameters.To the authors' knowledge, so far no studies have been published about multi-parameter optimization for hydraulic fracturing.The aim of this paper is to provide new insights into the effective development of continental shale gas formations by optimization of the hydraulic fracturing design.The results show that IR plays an important role in the interaction between silty laminae and shale matrix and stimulated area.The laminae structural thickness parameter and the mechanical parameters of cohesion and friction angle, present a reversed trend fromSRA and SLA.

Hydraulic Fracturing Model Setup
In shale fracturing, one of the challenges that influence the stimulation and gas productivity predication is the heterogeneity of reservoirs.Firstly, the rock matrix is composed of many different Energies 2016, 9, 325 3 of 19 kinds of minerals, crystal grains, pores and micro-cracks etc. from the mesoscopic perspective [34,35].Secondly, from a macro perspectivenatural fractures are greatly developed in natural reservoirs.Some parts of very anisotropic rocks fracture at a much lower pressure than other parts, thus the frac will have a preferential path through the shale fabric [36,37].The impact of anisotropic rock fabric and the flow paths that they create is one of the major unknowns for any shale.Even fracs in shale with extensive natural fracture systems have shown preferential fracture development that is traceable from well to well.The naturally fractured formations under hydraulic loading exhibit a unique feature: the flow and transport behavior within sound developed fractures are very different from those in rock mass with existing fractures under the same loading.When the rocks are not damaged, the permeability of rocks with existing fractures does not change, but it can vary dramatically due to the evolution of damage in fracturing rocks.The influence of damage on the variation of permeability as well as the original nature of the existing fractures in reservoirs is critical to shear stimulation of natural fractures.Another major issue in the description of the hydraulic fracturing behaviors of shale reservoirs concerns the irregular flow paths that depend on the heterogeneity of mechanical properties.In working with heterogeneous reserviors, an important aspect is to determine the specific data that are used to judge the influence of heterogeneity on the complicated flow paths in naturally fractured reserviors.To solve the coupled flow-damage problems, Tang et al. [38] proposed a flow-stress-damage (FSD) coupling model by taking into account the growth of existing fractures and the formation of new fractures.A flow-stress-damage coupled model for heterogeneous rocks that considering the growth of existing fractures and the propagation of new fractures is implemented in RFPA-Flow, this FSD model is used to investigate the behaviour of fluid flow and damage evolution, and their coupling behaviors, in specimens that are subjected to external loadings.The flow-stress-damaged coupled model has a unique ability to reappear the evolution of fluid flow; also, the fracturing phenomena from the local properties to the global properties can be realized.In the model, the initiation, propagation of natural fractures are dynamic, the interaction behaviors (cross, offset and capture) between hydraulic fracture and natural fracture are automatic realized according to the mechanical characteristics.

Brief Description of the Numerical Model
RFPA-Flow code is a numerical simulation tool using finite element analysis to handle the progressive failure of heterogeneous, permeable rock.In the model, the coupled effects of flow, stress and damage on the extension of existing/new fractures and the permeability change on account of damage evolution of the rocks were addressed.This coupled flow-stress-damage model in RFPA-Flow has been validated in previous publications [24,35,36,[38][39][40].There are two features distinguishing RFPA-Flow from other numerical approaches: (1) the RFPA-Flow code can simulate the non-linear deformation of a quasi-brittle behavior by introducing the heterogeneity of rock properties into the model, with an ideal brittle constitutive law for the local material; (2) by introducing a reduction of material parameters after element failure, the RFPA code can simulate strain-softening and discontinuous mechanics problems in a continuum mechanics mode.
For heterogeneity, the material properties (failure-strength σ c and elastic modulus E c ) for elements are randomly distributed throughout the model by following a WeIibull distribution: where σ is the element strength and σ 0 is the mean strength of the elements for the specimen.For the elastic modulus, E, the same distribution is used.We define m as the homogeneity index of the rock [38].According to the definition, a larger m implies a more homogeneous material and vice versa.In general, it is assumed that the strength properties and elastic modulus follow two separate distributions with the same homogeneity index.Systematic studies of the homogeneity index "m" have been detailed conducted previously [35][36][37][38].
In RFPA2D, the finite element method is employed toobtain the stress and fluid flow fields.The relationship between the fluid flow, stress, and strain in the fractured rock mass is deterimined by the Biot's consolidation theory [6].The main governing equations in RFPA2D are as follows: Equilibrium equation: Strain-displacement equation: Constitutive equation: Seepage equation: Coupling equation: where σ is stress; ρ is density; u is displacement; ε is strain; X is component of body force; α is coefficient of pore water pressure; λ is Lame coefficient; p is pore water pressure; δ is Kronecker constant; G is shear modulus; Q is Biot's constant; k is coefficient of permeability; k 0 is the initial coefficient of permeability; β is a coupling parameter that describe the effect of stress on the permeable property; and ζ (>1) is a mutation coefficient of permeability to describe the increase in permeability of the specimen during the formation of fractures.The value of ζ can be obtained from experimental tests of Thallak et al. [39] and Noghabai [40].Equations ( 2)-( 5) are based on Biot's theory of consolidation, and Equation (6) represents the effect of stress on permeability, which is introduced to describe the dependency of permeability on stress and damage, and the relationship between stress and permeability coefficient is assumed to follow a negative exponential equation.When the stress of the element satisfies a certain strength criterion (e.g., Coulomb criterion), the element starts to damge.In elastic damage mechanics, the elastic modulus of the element could degrade gradually with the incease of the damage degree, and the elastic modulus of the damaged element is defined as follows: E " p1 ´Dq E 0 (7) where D is the damage variable, E 0 and E are elastic modulus of the undamaged and the damaged material, respectively.When the tensile stress in an element reaches its tensile strength f 1 t , the constitutive relationship is adopted.This is: The damage variable can be described by Tang et al. [38] as: where f tr is the residual tensile strength of the element, and ε is equivalent principal strain of the element, ε to is the strain at the elastic limit, or threshold strain, and ε tu is the ultimate tensile strain of the element at which the element would be completely damaged.
In this case the permeability can be described as: where ξ (ξ > 1) is the damage factor of permeability, which reflects the damage-induced permeability increase [38].The value of ξ can be obtained from experimental tests.
In the model, both tensile and shear failure modes are considered.An element is considered to have failed in the tension mode when its minor principal stress exceeds the tensile strength of the element, as described by Equation ( 8), and have failed in shear mode when the compressive or shear stress has satisfied the Mohr-Coulumb failure criteriongiven by Tang et al. [38]: where σ 1  1 is the major effective principal stress, σ 1 3 is the minor effective principal stress, φ 1 is the minor effective angle of friction, f 1 t is the tensile strength of the element, and f 1 c is the compressive failure strength of the element.The damage factor under uniaxial compression is described as: where f 1 cr is the residual compressive strength, ε cu is the ultimate compressive strain of the element at which the element would be completely damaged.In this case, the permeability can be described by:

Realistic Failure Process Analysis-Flow Model Setup
It is believed that the characteristics of silty laminae have a critical effect on the response of naturally fractured reservoirs to fluid injection, and therefore the effectiveness of hydraulic fracturing.Core observations indicate that laminae are found developed in the producing shales.From the core statistical results of the Zhangjiatan Shale, in the southeastern of Ordos Basin, the frequency ranges from 4 to 32 laminae/m and 16 laminae/m on average, and the thickness of silty laminea thickness ranges from 0.2 to 4 mm and is 1.5 mm on average.The thickness proportion of laminae in the measured segments ranges from 6% to 17% [10].By measuring the typical core sample of Chang 7 wells, we found that thickness of silty laminae ranges from 0.5 mm to 6 mm (Figure 1b).Silty laminae was extremely heterogeneous, and composed of quartz, feldspar, mixed-layer montmorillonite, and chlorite, its mechanical properties are controlled by the microheterogeneity characteristics of these minerals.To study the multiple factors affecting optimization for a fracturing network, the flow chart shown in Figure 2 was used.Seven influential factors, such as thickness, spacing, dip angle, friction angle, SD, and IR are considered in the model.
Figure 3 shows the geometry and the set-up of the simulation model.The model represents a 2D horizontal section of a reservoir.In the model, the injection is through a vertical wellbore in the center of the model, the increasing injection pressure was imposed on the wellbore at constant rate.A wellbobre with a diameter of 20 mm is assumed to be located at the centre of a 200 ˆ200 mm model.The model is discretised into 40,000 (200 ˆ200) identical square elements.For all the simulations, the maximum horizontal stress (SHmax) was 20 MPa, the minimum horizontal stress (Shmin) was 15 MPa.The initial pore pressure was set to 10 MPa.The input material mechanical parameters for the numerical models are referred to Wang et al. [36], as shown in Table 1.For all the simulations, the slick-water-frac treatment is selected during the simulations, fluid rheology is 1 centipoise (cp).
Energies 2016, 9, 325 15 of 19 horizontal stress (SHmax) was 20 MPa, the minimum horizontal stress (Shmin) was 15 MPa.The initial pore pressure was set to 10 MPa.The input material mechanical parameters for the numerical models are referred to Wang et al. [36], as shown in Table 1.For all the simulations, the slick-waterfrac treatment is selected during the simulations, fluid rheology is 1 centipoise (cp).horizontal stress (SHmax) was 20 MPa, the minimum horizontal stress (Shmin) was 15 MPa.The initial pore pressure was set to 10 MPa.The input material mechanical parameters for the numerical models are referred to Wang et al. [36], as shown in Table 1.For all the simulations, the slick-waterfrac treatment is selected during the simulations, fluid rheology is 1 centipoise (cp).horizontal stress (SHmax) was 20 MPa, the minimum horizontal stress (Shmin) was 15 MPa.The initial pore pressure was set to 10 MPa.The input material mechanical parameters for the numerical models are referred to Wang et al. [36], as shown in Table 1.For all the simulations, the slick-waterfrac treatment is selected during the simulations, fluid rheology is 1 centipoise (cp).

D-Optimal Design for Response Surface Methodology
Manymethods have been proposed to conduct resolve uncertainty analysis [41][42][43][44][45].In this work, the RSM approach is applied to the optimization of the index SRA and SLA.RSM is utilized to approximate a response, in terms of maximum SRA and SLA, over a region of interest specified by the range of variability of input factors based on the least squares criterion.The RSM model can be in fully quadratic or linear.It can offer a cost-effective and efficient way to manage the uncertainties for shale gas reservoir development.A more detailed statistical and mathematical theories of RSM can be referred in the reports by Myers and Montgomery [46].The shale gas reservoir model is shown in Figure 2. Seven uncertain parameters such as laminae thickness (T), laminae spacing (S), laminae dip angle (θ), laminae IFA, laminae cohesion (C), SD and IR are given a reasonable range with the actual maximum and minimum values or coded symbol of "+1" and "´1" respectively, as shown in Table 2.For the convenience of building up the numerical model, and discuss the influence of thickness on fracturing effectiveness, the minimum thickness is set as 2mm.According to the sevenvariables, 46 cases were required, based on the approach of D-Optimal design, which originated from the optimal design theory [47,48].Table 3 shows the 46 combinations of the seven studied uncertain factors generated by the D-Optimal design.After the numerical simulation of each case, results of hydraulic fracturing are shown in Figure 4 for partial selected cases with various designs in Table 3. From Table 3, in order to qualitatively evaluate the effect of hydraulic fracturing, two of the indices for the model responses were evaluated during the fluid injection.They are: (1) SRA: defined as the interaction area of hydraulic fractures and silty laminae that has experienced a fluid pressure increase due to injection; (2) SLA: defined as the area of silty laminae that has been damage during hydraulic fracturing.SRA is the overall reflection of hydraulic fracturing effectiveness, and SLA reflects the hydraulic fracturing effectiveness of silty laminae.In Table 3, the quantitative maximum stimulated area was Energies 2016, 9, 325 8 of 19 listed in columns 9 and 10; they are the two responses in optimization.As there is no precise criteria for defining the SRA, a criteria based on fracture pressure change was employed in this work.The interaction area of hydraulic fractures and silty laminae having a pore pressure increament below the maximum value of 5 MPa was considered as SRA, it refers to the "red" region in Figure 5a.For SLA, during hydraulic fracturing, because of the lower mechanical properties and lower brittleness of the silty laminae, the silty laminae was easily damaged.From the simulation results, the picture labeled "black" (Figure 5) corresponds to the silty laminae that experienced injection pressure.The total area of the "black" area was considered as SLA.Determination of SRA and SLA are based on digital image process (DIP) method.Taking case 14 for example, after gray segmentation, the area of SRA and SLA are shown in Figure 5. are based on digital image process (DIP) method.Taking case 14 for example, after gray segmentation, the area of SRA and SLA are shown in Figure 5.

3.2.Response Surface MethodologyModel Analysis
According to the research results of [10], the amount of free gas and solution gas in shale increases with the increase of silty laminae, and the amount of adsorbed gas decreases.Silty laminae plays a significant role in free-state shale gas formation, migration, and production.Therefore, in this work, not only SRA but also SLA is considered for the optimization of the fracturing design across multiple parameters to determine the optimal stimulation design.Once the maximum simulated area of 46 run cases was obtained, the RSM method is used to analyze the relationship between the response value and influential factors.To choose the siutable model, the statistical method was adopted to judge which polynomial fit the equation with a linear model, two factor model, interaction model (2FI), quadratic model, or cubic model, as shown in Tables 4 and 5, which are response surface models for SRA and SLA, respectively.The criterion for selecting the appropriate model is choosing the highest polynomial model, where the additional terms are significant and the model is not aliased.Although the cubic model is the highest polynomial model, it is not selected because is it aliased.Aliasing is a result of reducing the number of experimental runs.When it occurs, several groups of  3. are based on digital image process (DIP) method.Taking case 14 for example, after gray segmentation, the area of SRA and SLA are shown in Figure 5.

3.2.Response Surface MethodologyModel Analysis
According to the research results of [10], the amount of free gas and solution gas in shale increases with the increase of silty laminae, and the amount of adsorbed gas decreases.Silty laminae plays a significant role in free-state shale gas formation, migration, and production.Therefore, in this work, not only SRA but also SLA is considered for the optimization of the fracturing design across multiple parameters to determine the optimal stimulation design.Once the maximum simulated area of 46 run cases was obtained, the RSM method is used to analyze the relationship between the response value and influential factors.To choose the siutable model, the statistical method was adopted to judge which polynomial fit the equation with a linear model, two factor model, interaction model (2FI), quadratic model, or cubic model, as shown in Tables 4 and 5, which are response surface models for SRA and SLA, respectively.The criterion for selecting the appropriate model is choosing the highest polynomial model, where the additional terms are significant and the model is not aliased.Although the cubic model is the highest polynomial model, it is not selected because is it aliased.Aliasing is a result of reducing the number of experimental runs.When it occurs, several groups of

Response Surface MethodologyModel Analysis
According to the research results of [10], the amount of free gas and solution gas in shale increases with the increase of silty laminae, and the amount of adsorbed gas decreases.Silty laminae plays a significant role in free-state shale gas formation, migration, and production.Therefore, in this work, not only SRA but also SLA is considered for the optimization of the fracturing design across multiple parameters to determine the optimal stimulation design.Once the maximum simulated area of 46 run cases was obtained, the RSM method is used to analyze the relationship between the response value and influential factors.To choose the siutable model, the statistical method was adopted to judge which polynomial fit the equation with a linear model, two factor model, interaction model (2FI), quadratic model, or cubic model, as shown in Tables 4 and 5 which are response surface models for SRA and SLA, respectively.The criterion for selecting the appropriate model is choosing the highest polynomial model, where the additional terms are significant and the model is not aliased.Although the cubic model is the highest polynomial model, it is not selected because is it aliased.Aliasing is a result of reducing the number of experimental runs.When it occurs, several groups of effects are combined Energies 2016, 9, 325 10 of 19 into one group and the most significant effect in the group is used to represent the effect of the group.Essentially, it is important that the model is not aliased.In addition, other criteria are needed to select the model that has the maximum "Adjusted R-Squared" and "Predicted R-Squared".Thus, the fully quadratic model and linear model are selected to build the SRA and SLA response surface in the subsequent optimization process, respectively.From the results of the anova for response surface quadratic model of SRA, the model F-value of 4.37 implies the model is significant.There is only a 0.86% chance that a "Model F-value" this large could occur due to noise.As shown in Figure 6, values of "Prob > F" less than 0.05 indicates the model terms are significant.In this case, factor of A, D, E, and G, are significant model terms.The influential order of these seven factors is: G-IR > A-thickness > D-cohesion > E-friction angle > C-dip angle > B-spacing > F-SD.effects are combined into one group and the most significant effect in the group is used to represent the effect of the group.Essentially, it is important that the model is not aliased.In addition, other criteria are needed to select the model that has the maximum "Adjusted R-Squared" and "Predicted R-Squared".Thus, the fully quadratic model and linear model are selected to build the SRA and SLA response surface in the subsequent optimization process, respectively.From the results of the anova for response surface quadratic model of SRA, the model F-value of 4.37 implies the model is significant.There is only a 0.86% chance that a "Model F-value" this large could occur due to noise.As shown in Figure 6, values of "Prob > F" less than 0.05 indicates the model terms are significant.In this case, factor of A, D, E, and G, are significant model terms.The influential order of these seven factors is: G-IR > A-thickness > D-cohesion > E-friction angle > C-dip angle > Bspacing > F-SD.From Figures 6 and 7, it can be seen that the operational parameter of the IR is the most sensitive factor in the models.This indicates that IR is crucial to hydraulic fracturing effectives.This conclusion is consistent with Wang et al. [35,36], King [49], Gil et al. [50] and Nagel et al. [51].With low IR, fluid is prone to leak into the natural fractures (bedding faces) despite the effect of fluid pressure and once the natural fractures accepts fluid, the pressure can rise far above the confining pressure with no induced fractures generation.With large IR, the hydraulic fracture is prone to cross the natural fractures resulting from the increase of the pressurization rate.In addition, the factor A (thickness) is From Figures 6 and 7 it can be seen that the operational parameter of the IR is the most sensitive factor in the models.This indicates that IR is crucial to hydraulic fracturing effectives.This conclusion is consistent with Wang et al. [35,36], King [49], Gil et al. [50] and Nagel et al. [51].With low IR, fluid is prone to leak into the natural fractures (bedding faces) despite the effect of fluid pressure and once the natural fractures accepts fluid, the pressure can rise far above the confining pressure with no induced fractures generation.With large IR, the hydraulic fracture is prone to cross the natural fractures resulting from the increase of the pressurization rate.In addition, the factor A (thickness) is also sensitive to both SRA and SLA.Silty laminae contains plenty of free gas, with the increase of thickness, the free gas recovery improves.
Energies 2016, 9, 325 20 of 19 also sensitive to both SRA and SLA.Silty laminae contains plenty of free gas, with the increase of thickness, the free gas recovery improves.
where A is laminae thickness; B is laminae length; C is laminae dip angle; D is laminae cohesion; E is laminae IFA; F is in situ SD, and G is IR.
The normal plot of residuals, reflecting the distribution of the residuals, for SRA and SLA are shown in Figure 8.All the points in the "Normal Plot of Residuals" fall on a straight line, meaning the residuals are normally distributed.Figure 9 shows the plot of "Predicted versus Actual" for SRA and SLA, illustrating whether the generated equation of stimulated area response surface accurately predicts the actual SRA and SLA values.It can be seen that generated SRA and SLA response surface models provide reliable predicted SRA and SLA values, as compared with the actual values of SRA and SLA.This means that the generated SRA and SLA response surface models are reliable.Figure 10 shows the 3D surface of the studied factors for SRA, the response surface represents a run case among the 46 cases.It shows the influential tread of the seven factors to SRA.As shown inFigure 10a, SRA decreases with the increase of thickness, and increases with the increase of spacing.Figure 10b shows the influence of spacing and dip angle on SRA, and SRA increases with increasing spacing.Figure 10c shows the influence of cohesion and friction angle on SRA, and SRA increases where A is laminae thickness; B is laminae length; C is laminae dip angle; D is laminae cohesion; E is laminae IFA; F is in situ SD, and G is IR.
The normal plot of residuals, reflecting the distribution of the residuals, for SRA and SLA are shown in Figure 8.All the points in the "Normal Plot of Residuals" fall on a straight line, meaning the residuals are normally distributed.Figure 9 shows the plot of "Predicted versus Actual" for SRA and SLA, illustrating whether the generated equation of stimulated area response surface accurately predicts the actual SRA and SLA values.It can be seen that generated SRA and SLA response surface models provide reliable predicted SRA and SLA values, as compared with the actual values of SRA and SLA.This means that the generated SRA and SLA response surface models are reliable.
Energies 2016, 9, 325 20 of 19 also sensitive to both SRA and SLA.Silty laminae contains plenty of free gas, with the increase of thickness, the free gas recovery improves.
where A is laminae thickness; B is laminae length; C is laminae dip angle; D is laminae cohesion; E is laminae IFA; F is in situ SD, and G is IR.
The normal plot of residuals, reflecting the distribution of the residuals, for SRA and SLA are shown in Figure 8.All the points in the "Normal Plot of Residuals" fall on a straight line, meaning the residuals are normally distributed.Figure 9 shows the plot of "Predicted versus Actual" for SRA and SLA, illustrating whether the generated equation of stimulated area response surface accurately predicts the actual SRA and SLA values.It can be seen that generated SRA and SLA response surface models provide reliable predicted SRA and SLA values, as compared with the actual values of SRA and SLA.This means that the generated SRA and SLA response surface models are reliable.Figure 10 shows the 3D surface of the studied factors for SRA, the response surface represents a run case among the 46 cases.It shows the influential tread of the seven factors to SRA.As shown inFigure 10a, SRA decreases with the increase of thickness, and increases with the increase of spacing.Figure 10b shows the influence of spacing and dip angle on SRA, and SRA increases with increasing spacing.Figure 10c shows the influence of cohesion and friction angle on SRA, and SRA increases with increasing cohesion and friction angle.Figure 10d shows the relationship of SD, IR and SRA, SRA increases with increasing IR and decreasing SD.It can be seen that large IR, cohesion and friction angle, can lead to the increasing SRA.This can be better interpreted that with larger cohesion and friction angle, the strength of silty laminae is higher; during hydraulic fracturing, leak off is relative lower into laminae, and the injected fluid is mainly used to communicate silty laminae and shale matrix.Figure 11 shows the 3D response surface of the studied factors for SLA. Figure 11a plots the relationship of thickness, spacing, and SLA.SLA increases with increasing thickness and decreasing spacing.Figure 11b presents the 3D surface of spacing and dip angle for SLA, and SLA increases with increasing dip angle.Figure 11c presents the 3D surface of the cohesion and friction angle for Figure 10 shows the 3D surface of the studied factors for SRA, the response surface represents a run case among the 46 cases.It shows the influential tread of the seven factors to SRA.As shown inFigure 10a, SRA decreases with the increase of thickness, and increases with the increase of spacing.Figure 10b shows the influence of spacing and dip angle on SRA, and SRA increases with increasing spacing.Figure 10c shows the influence of cohesion and friction angle on SRA, and SRA increases with increasing cohesion and friction angle.Figure 10d shows the relationship of SD, IR and SRA, SRA increases with increasing IR and decreasing SD.It can be seen that large IR, cohesion and friction angle, can lead to the increasing SRA.This can be better interpreted that with larger cohesion and friction angle, the strength of silty laminae is higher; during hydraulic fracturing, leak off is relative lower into laminae, and the injected fluid is mainly used to communicate silty laminae and shale matrix.with increasing cohesion and friction angle.Figure 10d shows the relationship of SD, IR and SRA, SRA increases with increasing IR and decreasing SD.It can be seen that large IR, cohesion and friction angle, can lead to the increasing SRA.This can be better interpreted that with larger cohesion and friction angle, the strength of silty laminae is higher; during hydraulic fracturing, leak off is relative lower into laminae, and the injected fluid is mainly used to communicate silty laminae and shale matrix.Figure 11 shows the 3D response surface of the studied factors for SLA. Figure 11a plots the relationship of thickness, spacing, and SLA.SLA increases with increasing thickness and decreasing spacing.Figure 11b presents the 3D surface of spacing and dip angle for SLA, and SLA increases with increasing dip angle.Figure 11c presents the 3D surface of the cohesion and friction angle for Figure 11 shows the 3D response surface of the studied factors for SLA. Figure 11a plots the relationship of thickness, spacing, and SLA.SLA increases with increasing thickness and decreasing spacing.Figure 11b presents the 3D surface of spacing and dip angle for SLA, and SLA increases with increasing dip angle.Figure 11c presents the 3D surface of the cohesion and friction angle for SLA.While the cohesion and friction angle increase, SLA decreases gradually.Figure 11d plots the relationship of SD, IR, and SLA, and SLA increases with increasing IR, and decreases with increasing SD.Similarly to SRA, IR is a positive factor resulting in the increase of SLA; however, cohesion and friction angle are negative factors leading to decreasing SLA.In silty laminae, shale gas is mainly in free state and solution state, the amount of free-state and solution gas in shale increases with increasing amount of silty laminae [10].SLA.While the cohesion and friction angle increase, SLA decreases gradually.Figure 11dplots the relationship of SD, IR, and SLA, and SLA increases with increasing IR, and decreases with increasing SD.Similarly to SRA, IR is a positive factor resulting in the increase of SLA; however, cohesion and friction angle are negative factors leading to decreasing SLA.In silty laminae, shale gas is mainly in free state and solution state, the amount of free-state and solution gas in shale increases with increasing amount of silty laminae [10].The relatively large porosity that developed in silty laminae facilitates gas storage, and the mechanical properties are relatively weaker.This may be the reason why the hydraulic fracturing effectiveness is good for weak silty laminae, as studied in this paper.

Stimulated Area Optimization
The index of SRA indicates the total interaction area between hydraulic fractures and silty laminae, the larger the index, the better the fracturing network is.Also, the SLA index represents the production of free-state and solution gas.These two indices are very important to shale gas productivity in silty laminae formations, particularly in southeastern Ordos Basin, China.In this work, the numerical optimization was done to select the set of variables that leads to the maximum SRA and SLA value.A total of 71 optimal projects were generated after numerically optimization.Table 6 shows the top ten optimal solutions maximizing the SRA and SLA simultaneously.The relatively large porosity that developed in silty laminae facilitates gas storage, and the mechanical properties are relatively weaker.This may be the reason why the hydraulic fracturing effectiveness is good for weak silty laminae, as studied in this paper.

Stimulated Area Optimization
The index of SRA indicates the total interaction area between hydraulic fractures and silty laminae, the larger the index, the better the fracturing network is.Also, the SLA index represents the production of free-state and solution gas.These two indices are very important to shale gas productivity in silty laminae formations, particularly in southeastern Ordos Basin, China.In this work, the numerical optimization was done to select the set of variables that leads to the maximum SRA and SLA value.A total of 71 optimal projects were generated after numerically optimization.Table 6 shows the top ten optimal solutions maximizing the SRA and SLA simultaneously.
From Table 6, it can be seen that the IR is a critical factor controlling the stimulated area of SRA and SLA.When SRA and SLA reach a maximum, the IR is 0.02 m 3 /d/m.In addition, the smaller the SD, the maximum stimulated area can be reach for the SRA and SLA.Dip angle of silty laminae is about 60 ˝when SRA and SLA reach the maximum.Mechanical properties control the damage and failure mechanism during hydraulic fracturing, according to Mohr-Coulomb shear failure criterion [52], the fraction strength of the silty laminae affects the fracturing evolution to a large extent.Figure 12a-d shows the 3D response surfaces for the optimal SRA and SLA with desirability of 0.780, respectively.The shape of the response surface displays a curved surface, which is different from the shape of SRA and SLA.This reflects the non-linear fracturing process when SRA and SLA are at a maximum simultaneously.
Energies 2016, 9, 325 17 of 19 From Table 6, it can be seen that the IR is a critical factor controlling the stimulated area of SRA and SLA.When SRA and SLA reach a maximum, the IR is 0.02 m 3 /d/m.In addition, the smaller the SD, the maximum stimulated area can be reach for the SRA and SLA.Dip angle of silty laminae is about 60° when SRA and SLA reach the maximum.Mechanical properties control the damage and failure mechanism during hydraulic fracturing, according to Mohr-Coulomb shear failure criterion [52], the fraction strength of the silty laminae affects the fracturing evolution to a large extent.Figure 12a-d shows the 3D response surfaces for the optimal SRA and SLA with desirability of 0.780, respectively.The shape of the response surface displays a curved surface, which is different from the shape of SRA and SLA.This reflects the non-linear fracturing process when SRA and SLA are at a maximum simultaneously.

Discussion
From the RSM analysis of SRA and SLA, the influence factor of thickness presents a reverse trend.SRA has a negative correlation to laminae thickness, however, SLA has a positive correlation to laminae thickness.Also, laminae spacing is a negative factor for both SRA and SLA.Due to the weak mechanical properties of silty laminae compared to shale matrix, silty laminae has a larger porosity than shale matrix, the injected fluid is mainly used to reactive the silty laminae, and leak-off in silty laminae is most serious,thus, generation of hydraulic fractures is unfavorable.With increased spacing and SD, the degree of communicationbecomes weaker during hydraulic fracturing, which is unfavorable both for SRA and SLA.These two conclusions are consistent with the results of [32].Dip angle of silty laminae has a negative correlation with SRA, but a positive correlation with SLA.Generally speaking, the general opinion on the mechanism leading to the success of waterfrac in shale

Discussion
From the RSM analysis of SRA and SLA, the influence factor of thickness presents a reverse trend.SRA has a negative correlation to laminae thickness, however, SLA has a positive correlation to laminae thickness.Also, laminae spacing is a negative factor for both SRA and SLA.Due to the weak mechanical properties of silty laminae compared to shale matrix, silty laminae has a larger porosity than shale matrix, the injected fluid is mainly used to reactive the silty laminae, and leak-off in silty laminae is most serious,thus, generation of hydraulic fractures is unfavorable.With increased spacing and SD, the degree of communicationbecomes weaker during hydraulic fracturing, which is unfavorable both for SRA and SLA.These two conclusions are consistent with the results of [32].Dip angle of silty laminae has a negative correlation with SRA, but a positive correlation with SLA.Generally speaking, the general opinion on the mechanism leading to the success of waterfrac in shale gas reservoirs is that a complex fracture network is created by the stimulation of pre-existing natural fractures and the silty laminae plays the role of natural fractures, which is actually a kind of cement-filled material.Fracture complexity is thought to be enhanced when pre-existing fractures are oriented at an angle to the maximum stress direction, or when both horizontal stresses and horizontal stress anisotropy are low, because the natural fractures in multiple orientations are prone to be stimulated with this combination of stress condition.The two mechanical parameters of cohesion and fraction angle, are both positively correlated with SRA and negative correlation to SLA.The reason for this result is that with lower mechanical strength, silty laminae are prone to damage and failure; once the laminae elements are damaged, it is different to generate hydraulic fractures, which reduces the interaction between hydraulic fractures and silty laminae.IR is a positive factor both for SRA and SLA, at a low IR, fluid tends to leak into the pre-existing discontinuities despite the influence of fluid pressure and once the fluid leak off into the natural fractures, the pressure can rise far above the confining pressure without inducing new fractures.At a large IR, the hydraulic fracture tends to cross the natural fracture due to the increase in the pressurization rate.

Conclusions
In this work, RSM was used to obtain the optimal design for gas shale hydraulic fracturing by optimizing uncertain parameters (i.e., laminae thickness, spacing, dip angle, cohesion, IFA, insituS, and IR) to maximizethe stimulated area.The following conclusions can be drawn from this study: (1) The proposed approach is practical and efficient for the design and optimization of hydraulic fracturing multi-parameters.SRA and SLA response surface models givereliable predicted values compared with the actual values of SRA and SLA; (2) By the RSM optimization, the optimal design combinations for silty laminae shale are obtained.
Among the optimal projects, IR is a positive factor to increase SRA and SLA; dip angle of laminae is about 60 ˝when SRA and SLA reach the maximum for all the optimal projections; (3) The thickness factor is sensitive to both SRA and SLA.Silty laminae contains plenty of free gas and solution gas, and with the increase of thickness, the free gas recovery improves.However, due to the strong leak-off characteristics of silty laminae, the factor of thickness is a negative factor for the overall stimulated area.The influence of mechanical properties is oppositein SRA compared to SLA, with lower mechanical strength, silty laminae is prone to damage and failure; the injected fluid is difficult to generate hydraulic fractures, which reduces the interaction between hydraulic fractures and silty laminae;

Figure 1 .
Figure 1.The typical shale core with silty laminae from Chang 7-3 member: (a)outcrop of silty laminae in a Chang 7 stratum; and (b)core sample of Chang 7 formation with typical silty laminae structure.

Figure 2 .Figure 3 .
Figure 2. The flow chart for the multiple parameters evalation and optimization to silty laminae shale formation.

Figure 1 .
Figure 1.The typical shale core with silty laminae from Chang 7-3 member: (a) outcrop of silty laminae in a Chang 7 stratum; and (b) core sample of Chang 7 formation with typical silty laminae structure.

Figure 1 .
Figure 1.The typical shale core with silty laminae from Chang 7-3 member: (a)outcrop of silty laminae in a Chang 7 stratum; and (b)core sample of Chang 7 formation with typical silty laminae structure.

Figure 2 .
Figure 2. The flow chart for the multiple parameters evalation and optimization to silty laminae shale formation.

Figure 3 .
Figure 3.The geometry and model setup.

Figure 2 .
Figure 2. The flow chart for the multiple parameters evalation and optimization to silty laminae shale formation.

Figure 1 .
Figure 1.The typical shale core with silty laminae from Chang 7-3 member: (a)outcrop of silty laminae in a Chang 7 stratum; and (b)core sample of Chang 7 formation with typical silty laminae structure.

Figure 2 .
Figure 2. The flow chart for the multiple parameters evalation and optimization to silty laminae shale formation.

Figure 3 .
Figure 3.The geometry and model setup.

Figure 3 .
Figure 3.The geometry and model setup.

Figure 4 .
Figure 4.The hydraulic fracturing simulation results for the partial 8 run cases according to Table3.

Figure 5 .
Figure 5. Determination of SRA and SLA based on gray segmentation for case 14: (a) color shadow indicate relative magnitude of the pore water pressure field, andSRA is the overall reflection of hydraulic fracturing effectiveness, it corresponds to the red region; and (b) SLA reflects the hydraulic fracturing effectiveness of silty laminae, it corresponds to the black region in (a).

Figure 4 .
Figure 4.The hydraulic fracturing simulation results for the partial 8 run cases according to Table3.

Figure 4 .
Figure 4.The hydraulic fracturing simulation results for the partial 8 run cases according to Table3.

Figure 5 .
Figure 5. Determination of SRA and SLA based on gray segmentation for case 14: (a) color shadow indicate relative magnitude of the pore water pressure field, andSRA is the overall reflection of hydraulic fracturing effectiveness, it corresponds to the red region; and (b) SLA reflects the hydraulic fracturing effectiveness of silty laminae, it corresponds to the black region in (a).

Figure 5 .
Figure 5. Determination of SRA and SLA based on gray segmentation for case 14: (a) color shadow indicate relative magnitude of the pore water pressure field, andSRA is the overall reflection of hydraulic fracturing effectiveness, it corresponds to the red region; and (b) SLA reflects the hydraulic fracturing effectiveness of silty laminae, it corresponds to the black region in (a).

Figure 6 .
Figure 6.Anova for SRA response surface with linear model.The anova for the response surface linear model of SLA is shown in Figure7.The model F-value of 3.01 implies the model is significant.There is only a 1.29% chance that a "Model F-value" this large could occur due to noise.Values of "Prob > F" less than 0.05 indicates the model terms are significant.In this case A, B, D, and G, are significant model terms.The influential order of these seven factors is: G-IR > A-thickness > D-cohesion > B-spacing > F-SD > C-dip angle > E-friction angle.From Figures6 and 7, it can be seen that the operational parameter of the IR is the most sensitive factor in the models.This indicates that IR is crucial to hydraulic fracturing effectives.This conclusion is consistent with Wang et al.[35,36], King[49], Gil et al.[50] and Nagel et al.[51].With low IR, fluid is prone to leak into the natural fractures (bedding faces) despite the effect of fluid pressure and once the natural fractures accepts fluid, the pressure can rise far above the confining pressure with no induced fractures generation.With large IR, the hydraulic fracture is prone to cross the natural fractures resulting from the increase of the pressurization rate.In addition, the factor A (thickness) is

Figure 6 .
Figure 6.Anova for SRA response surface with linear model.

Figure 9 .Figure 10 .
Figure 9. Predicted stimulated area vs. the actual stimulated area: (a) The predicted area vs. the actural area for SRA; and (b) The predicted area vs. the actural area for SLA.

Figure 9 .
Figure 9. Predicted stimulated area vs. the actual stimulated area: (a) The predicted area vs. the actural area for SRA; and (b) The predicted area vs. the actural area for SLA.

Figure 9 .Figure 10 .
Figure 9. Predicted stimulated area vs. the actual stimulated area: (a) The predicted area vs. the actural area for SRA; and (b) The predicted area vs. the actural area for SLA.

Figure 10 .
Figure 10.The influential of studied factors on SRA response surface: (a) plots the influence of factor A and B on SRA; (b) plots the influence of factor B and C on SRA; (c) plots the influence of factor D and E on SRA; and (d) plots the influence of factor F and G on SRA.

Figure 11 .
Figure 11.The influential factors studied on SLA response surface: (a) plots the influence of factor A and B on SLA; (b) plots the influence of factor B and C on SLA; (c) plots the influence of factor D and Eon SLA; and (d) plots the influence of factor F and G on SLA.

Figure 11 .
Figure 11.The influential factors studied on SLA response surface: (a) plots the influence of factor A and B on SLA; (b) plots the influence of factor B and C on SLA; (c) plots the influence of factor D and Eon SLA; and (d) plots the influence of factor F and G on SLA.

Figure 12 .
Figure 12.The 3D response surfaces for SRA and SLA with optimal project: (a) plots the influence of factor A and B on the desirability of solutions; (b) plots the influence of factor C and D on the desirability of solutions; (c) plots the influence of factor D and E on the desirability of solutions; and (d) plots the influence of factor F and G on the desirability of solutions;

Figure 12 .
Figure 12.The 3D response surfaces for SRA and SLA with optimal project: (a) plots the influence of factor A and B on the desirability of solutions; (b) plots the influence of factor C and D on the desirability of solutions; (c) plots the influence of factor D and E on the desirability of solutions; and (d) plots the influence of factor F and G on the desirability of solutions;

( 4 )
The influential order of the studied factors to SRA is: G-IR > A-thickness > D-cohesion > E-friction angle > C-dip angle > B-spacing > F-SD.And the influential order to SLA is: G-IR > A-thickness > D-cohesion > B-spacing > F-SD > C-dipangle > E-friction angle.Factors of laminae thickness, cohesion, and IR are the most significant factors for both SRA and SLA.

Table 1 .
Parameters used in multi-parameter optimization.IFA: internal friction angle.

Table 2 .
Uncertainty parameters of multiple parameter optimizations in this study.SD: stress difference; IR: injection rate.

Table 3 .
D-Optimal design table for the hydraulic fracturing simulations.SRA: stimulated reservoir area; SLA: silty laminae area.

Table 4 .
Statistical approach to select the response surface methodology (RSM) model for SRA.

Table 5 .
Statistical approach to select the RSM model for SLA.

Table 4 .
Statistical approach to select the response surface methodology (RSM) model for SRA.

Table 5 .
Statistical approach to select the RSM model for SLA.