Back Analysis of Rock Hydraulic Fracturing by Coupling Numerical Model and Computational Intelligent Technology

Hydraulic fracturing is widely used to determine in situ stress of rock engineering. In this paper we propose a new method for simultaneously determining the in situ stress and elastic parameters of rock. The method utilizing the hydraulic fracturing numerical model and a computational intelligent method is proposed and verified. The hydraulic fracturing numerical model provides the samples which include borehole pressure, in situ stress, and elastic parameters. A computational intelligent method is applied in back analysis. A multioutput support vector machine is used to map the complex, nonlinear relationship between the in situ stress, elastic parameters, and borehole pressure. The artificial bee colony algorithm is applied in back analysis to find the optimal in situ stress and elastic parameters. The in situ stress is determined using the proposed method and the results are compared with those of the classic breakdown formula.The proposedmethod provides a good estimate of the relationship between the in situ stress and borehole pressure and predicts the maximum horizontal in situ stress with high precision while considering the influence of pore pressure without the need to estimate Biot’s coefficient and other parameters.


Introduction
Hydraulic fracturing is widely used in the recovery of oil, gas, geothermal, and mineral resources [1].In petroleum engineering it is important to determine the in situ stresses and elastic parameters of the rock mass when using hydraulic fracturing in fracturing operations, wellbore stability analysis, and reservoir simulation [2].While high accuracy is required for the values of the in situ stress and mechanical parameters of the rock mass, determination of these parameters is still one of the most challenging tasks in hydraulic fracturing.
Hydraulic fracturing tests are considered the most effective method for determining the in situ stress and mechanical parameters of rock mass [3][4][5][6][7][8][9].The Hubbert and Willis hydraulic fracturing criterion and Haimson and Fairhurst's hydraulic fracturing criterion are the two classic formulae for hydrofracture breakdown ( [10]; Hubbert et al., 1953).However, the pore pressure term, which is a significant factor in deep boreholes, is ignored in Hubbert and Willis's hydraulic fracturing criterion.Modifications of the original equations were proposed to account for the pore pressure (Detournay et al., 1988; [1,[11][12][13]), but they have not been used in practice because it involves the Biot poroelastic parameters and Poisson's ratio which are difficult to obtain.Schmitt and Zoback built a more useful generalized form of the hydrofracture breakdown equation by considering the poroelastic effects [1].It can be used to provide upper and lower bound to the maximum horizontal in situ stress because it depends on the specific pore and microcrack structure.However, this method requires the poroelastic coefficients which are difficult to determine in practice.
Owing to the limitations of the classic breakdown formulae and the complexity of hydraulic fracturing tests, laboratory and field tests have been commonly used to determine the in situ stress and mechanical parameters of rock mass (Algorithm 1).However, these tests may not always produce the poroelastic parameters or may provide inaccurate results  because of low-quality core samples [14].Alternatively, back analysis, associated with the "in situ" approach, has been widely used to determine the mechanical parameters of rock mass in rock engineering [15][16][17][18][19][20].Zhang and Yin proposed a back analysis method which combined a neural network and a genetic algorithm to simultaneously identify the in situ stresses and elastic parameters [2]; however, this method did not consider the poroelastic effect.To overcome this difficulty, in this paper we extend our proposed displacement back analysis method to determine the in situ stress and mechanical parameters of a rock mass based on measured borehole pressure.The borehole pressure can be easily measured in the field with a pressure gauge installed inside the borehole [21].Back analysis is implemented following an optimization strategy based on the multioutput support vector machine (MSVM) and artificial bee colony algorithm (ABC) model, which is effective in multiple parameter identification [22].
The rest of this paper is organized as follows.The classic breakdown formula is presented in detail in Section 2. The formulation and procedure of back analysis based on borehole pressure are presented in detail in Section 3. In Section 4, a numerical example is used to verify the proposed method, and our conclusions are presented in Section 5.

Hydraulic Breakdown Equations
Hydraulic fracturing is a widely accepted technology used for determining in situ stress magnitude and direction.The principal stress  V has a magnitude equal to the overburden pressure in the vertical direction.The smallest horizontal principal stress  ℎmin is usually determined directly in the experiment from the shut-in pressure.The greatest horizontal principle stress  max must be calculated using a breakdown formula derived from an appropriate hydraulic fracturing model.Hubbert and Willis proposed a classic breakdown formula (1) to calculate  max for hydraulic fracturing in nonporous impermeable rocks [23], ignoring the pore pressure term.
where  is the rock tensile strength.Equations ( 2) and ( 3) are the breakdown formulae of porous impermeable rocks and porous permeable rocks, respectively, including the pore pressure [24]: where   is the breakdown pressure,   is the pore pressure,  is the Biot poroelastic parameter, and  is Poisson's ratio.Although (3) may best describe the conditions under which hydraulic fracturing is conducted from an open borehole, ( 2) is used in practice because of the difficulty of determining  and .Schmitt and Zoback proposed a more generalized form for the equations of hydrofracture breakdown for porous impermeable rocks and porous permeable rocks [1]: where  is the poroelastic effect parameter.

Back Analysis Model Based on Borehole Pressure
The in situ stress can be estimated based on borehole pressure using the hydraulic breakdown equations (Section 2).However, these equations present some limitations in practice.Therefore, we propose a back analysis method that combines a numerical method and an intelligent computational method.A multioutput support vector machine (MSVM) is used to map the complex, nonlinear relationship between the in situ stress, elastic parameters, and borehole pressure.
The numerical method provides the training samples for the MSVM.It is important to use an optimization method in back analysis.Here, we use the ABC algorithm to find the best-fit in situ stress and elastic parameters by comparing the measured pressure data and the MSVM predicted pressure.

Nonlinear Relationships between Pressure and Geomechanical Parameters.
The relationship between the borehole pressure and geomechanical parameters can be derived by the MSVM.The basic idea of MSVM is to extend the single-output support vector machine to a multidimensional output case.Given a set of training samples {( 1 ,  1 ), ( 2 ,  2 ), . . ., (  ,   )},   ∈   ,   ∈   , the MSVM model can be built by solving the following optimization problem based on an iterative reweighted leastsquare algorithm [25].
where  is the number of input,  is the number of output,  is the weight,  and  are constants,  is the sum of constant terms that do not depend on either  or ,  is the tolerant error, and  denotes the tth iteration.A brief description and the MSVM algorithm can be found in the literature [22].The MSVM model can be expressed as

Geofluids
Based on the above MSVM model, the nonlinear relationship between the borehole pressure and geomechanical parameters can be described as where X = ( 1 ,  2 , . . .,   ) is the -dimensional vector of the identified parameter, for example, the in situ stress, Young's modulus, or Poisson's ratio.Y = ( 1 ,  2 , . . .,   ) is the dimensional vector of the borehole pressure.
To build the MSVM model, the necessary training or learning samples are constructed and the MSVM parameters are determined.The samples are constructed by numerical analysis which computes the corresponding borehole pressure for a given set of tentative determined parameters.The MSVM parameters have a strong influence on the performance of the MSVM.In this study, we determined these parameters using the formulation presented by Meza et al. [26].

Optimization Method.
The back analysis ABC algorithm, developed by Karaboga [27], was adopted to search for the optimal geomechanical parameters of the rock mass.In the algorithm, the colony of artificial bees consists of three groups: employed bees, onlookers, and scouts.The ABC algorithm involves a cycle of four phases: the initialization phase, employed bees phase, onlooker bee phase, and scout bee phase.
In the initialization phase, the ABC generates a randomly distributed initial population of SN solutions and calculates the fitness of each solution.
(, ) = where (, ) is the candidate solution of the problem;  = 1, 2, . . ., SN/2 and SN/2 denotes the size of the population;  = 1, 2, . . .,  and  is the dimension number of each solution; rand(0, 1) is a random number between [0, 1];   min and   max are the upper and lower bounds of each solution.Once initialization is completed, the employed bees search for a solution and calculate the fitness value (see Section 3.3) in the employed bees phase.A candidate solution is produced according to the following equation: V (, ) =  (, ) +   ( (, ) −  (, )) , where  is different from  and is a randomly chosen index from {1, 2, . . ., SN/2},  is also an index randomly chosen from {1, 2, . . ., }, and   is a random number in the range [−1, 1] that controls the generation of food sources around (, ) and represents the comparison of two food positions seen by a bee.
In the onlooker bee phase, the onlooker bees choose a solution based on the fitness value, determine which solution will be abandoned, and allocate its employed bees as scout bees.The probability of being selected for each fitness value can be expressed as where fitness  is the fitness value of the solution.
Finally, in the scout bee phase the scout bees randomly search for a new solution in the determined ranges.A solution that cannot be improved further through a predetermined number of cycles is assumed to be abandoned by the onlookers.

The Fitness Function.
In order to find the optimal solution, it is necessary to build the fitness function for the ABC algorithm; that is, where MSVM() is the predicted pressure using the MSVM model, Y is the vector of the monitored pressure, and  is the number of monitored points.

Procedure of the MSVM-ABC Based Method.
If the MSVM model can establish the nonlinear relation between the borehole pressures and determined parameters, the model can be used to predict the borehole pressures.Then, the ABC algorithm is utilized to find the optimal parameters through error minimization between the pressures predicted by the MSVM model and the measured pressures.The back analysis flowchart is shown in Figure 1 and the algorithm is described as follows.
Step 1. Determine the general information and data such as the unknown (determined by back analysis) and known parameters of the numerical model, the MSVM and ABC algorithm parameters, and the range of parameters to be determined.
Step 2. Generate the combination of determined parameters, calculate the borehole pressure for each combination, and then build the learning samples for MSVM.
Step 3. Based on the learning samples of Step 2, construct the MSVM model using the MSVM algorithm and activate the ABC algorithm.
Step 4. Search for the optimal determined parameters based on the monitored pressure.

Validation and Application
To verify the proposed method, a numerical example is adopted to determine the in situ stress and elastic mechanical parameters of the elastic rock.The numerical experiment is conducted based on 2D hydraulic fracturing model of water injection into a hypothetical deep formation.Further details of the physical and numerical model can be found in the literature published by Schmitt and Zoback [1].The parameters to be determined are the maximum and minimum horizontal in situ stress  max and  ℎmin , respectively, Young's modulus E, and Poisson's ratio .The rock mass in all the zones is considered to be elastic.The mechanical parameters of the joints and the permeability Step 1: general information and data Step 2: numerical model Step 3: MSVM and ABC model Step 4: back analysis model Figure 1: Flowchart of the back analysis process to obtain the rock parameters.
of the rock mass are known; the parameter values can be seen in the literature published by Schmitt and Zoback [1].Thirty sets of training samples and ten testing samples derived in previous studies [1,2] were selected.Based on the MSVM algorithm, the MSVM code was written in Excel and VBA.The MSVM parameters and some of the weight   and constant   values and samples are shown in Figure 2. Good agreement between the measured data and the pressures estimated by the MSVM is shown in Figure 3, indicating the good performance of the MSVM model.Thus, the proposed model can accurately estimate the borehole pressures, replacing the existing numerical analysis method for calculating borehole pressures.The results also confirm that the MSVM model provides an accurate representation of the nonlinear relationship between the pressures and the determined parameters.The ABC code is also written in Excel and VBA.The parameters of the ABC algorithm and the calculation results are shown in Figure 4. Based on the proposed method for determining the in situ stress and mechanical parameters of rock mass, the results are shown in Table 1.we obtained the values of  max ,  ℎmin , , and  as 24.46 MPa, 14.33 MPa, 44.02 GPa, and 0.25, respectively.The elastic mechanical parameters of the rock agree with the results calculated by the Genetic Algorithm-Neural Network (ANN-GA) [2].A comparison of in situ stress values calculated using four different formulations is shown in Figure 5.  ℎmin agrees well with the value estimated by ANN-GA, (1), and (2).The relative error is only 1.07%. max agrees well with the value estimated by (2), which considers the pore pressure, but differs considerably from the values calculated by (1) and ANN-GA which do not consider the poroelastic effects.The relative error is up to 31.8%.Using (4), we obtain the upper and lower limits of the maximum horizontal in situ stress  max (14.95-34.55MPa).The value 24.46 MPa is within this range.Thus, the proposed method can be used in back analysis as an alternative numerical analysis method, which considers the poroelastic effects and provides rational, high-precision results.Note that the proposed method can determine the maximum horizontal in situ stress without estimating the poroelastic coefficient, which is a difficult parameter to obtain.Moreover, there are four borehole pressures, namely, formation breakdown pressure (FBP) 1, fracture propagation pressure (FPP) 2, instantaneous shut-in pressure (ISIP) 3, and leak-off pressure (LOP) 4.A comparison of back analysis on borehole pressures obtained by three different methods is presented in Figure 6.The borehole pressure calculated by the proposed MSVM method is very close to the measured pressure.The relative error is less than 3%.On the other hand, the convergence processes of the algorithm and fitness variations are shown in Figures 7 and 8. Initially, the data was distributed randomly in the searching space (Figure 8) and then converged to the solution of the problem in the 500th generation.This indicates that the proposed method can determine both the in situ stress and elastic mechanical parameters of the elastic rock with excellent converging performance.

Conclusions
In this paper, a new borehole pressure-based back analysis approach to determine the stress and mechanical parameters of rock mass is proposed.The method combines a coupling numerical model of hydraulic fracturing and a computational intelligent method.The method is applied to a numerical example to successfully determine both the in situ stress and mechanical parameters of a rock mass.In this approach, the MSVM is adopted to represent the nonlinear relationship between the borehole pressure and mechanical parameters of the rock mass, proving more efficient than existing numerical models.The ABC algorithm is used to search for the optimal parameters in the search space.The proposed approach is implemented in Excel with VBA.
In the classic breakdown formula, it is difficult in practice to determine the maximum horizontal in situ stress while considering the poroelastic coefficient.The proposed back analysis method can predict the maximum horizontal in situ stress based on the borehole pressures without the need to obtain the poroelastic coefficient.Thus, it is a more practical method for determining the in situ stress from hydraulic fracturing.The proposed method is practical and accurate and can be conveniently applied to simultaneously determine the in situ stress and mechanical parameters of rock from hydraulic fracturing.

Symbols 𝜎 V :
Principal stress in  direction  max : The greatest horizontal principle stress  ℎmin : The smallest horizontal principal stress : Poisson's ratio : Poroelastic effect parameter : Th ew e i g h tv e c t o r MSVM(): The predicted pressure using the MSVM model : Hyper parameter that determines trade-off between the regularization and the error reduction term : T o l e r a n te r r o r : N u m b e ro fo u t p u t   : A random number in the range [−1, 1] : A constant for classification threshold   : Breakdown pressure   : Pore pressure : Y o u n g ' sm o d u l u s : Biot poroelastic parameter : R o c kt e n s i l es t r e n g t h : N u m b e ro fi n p u t : Sum of constant terms that do not depend on either  or  rand(0, 1): A random number between [0, 1] fitness: Fitness value of the solution.no.2016YFC0600702), and Innovative Research Team (in Science and Technology) in University of Henan Province (no.15IRTSTHN029).

'
The value of the numbers of parameters of training samples N = Range(''N'').Cells.Value Dim x = Range(''Dim x'').Cells.Value Dim y = Range(''Dim y'').Cells.Value 'The parameters of SVM C = Range('' C'').Cells.Value epsilon = Range(''epsilon'').Cells.Value sigma = Range(''sigma'').Cells.Value ReDim X input(1 To N, 1 To Dim x) As Double ReDim Y output(1 To N, 1 To Dim y) As Double ReDim xi(1 To Dim x) As Double ReDim xj(1 To Dim x) As Double ReDim W(1 To Dim y, 1 To N) As Double ReDim b(1 To Dim y) As Double ReDim W k(1 To Dim y, 1 To N) As Double ReDim b k(1 To Dim y) As Double ReDim W s(1 To Dim y, 1 To N) As Double ReDim b s(1 To Dim y) As Double ReDim u(1 To N) As Double ReDim ai(1 To N) As Double ReDim ui(1 To Dim y) As Double ReDim u new(1 To N) As Double ReDim A(1 To N + 1, 1 To N + 1) As Double 'ReDim X(1 To N + 1, 1 To Dim y) As Double ReDim BB(1 To N + 1, 1 To Dim y) As Double ReDim A last row(1 To N) As Double ReDim B last row(1 To Dim y) As Double ReDim P W(1 To Dim y, 1 To N) As Double ReDim P b(1 To Dim y) As Double ReDim k(1 To N, 1 To N) As Double ReDim D a(1 To N, 1 To N) As Double ReDim D a 1(1 To N, 1 To N) As Double ReDim kf(1 To N) As Double 'Read the input of training samples For I = 1 To N For j = 1 To Dim x X input(I, j) = Range(''xi'').Cells(I, j) Next Next 'Read the output of training samples For I = 1 To N For j = 1 To Dim y Y output(I, j) = Range(''yi'').Cells(I, j).Value Next Next ' The initial value of Wk and bk For I = 1 To Dim y b(I) = 0 For j = 1 To N W(I, j) = 0 Next Next delta u = 1 Algorithm 1: Continued.
and determine the parameters: Specify the determined parameter Build the numerical model Specify the parameters of MSVM algorithm Specify the parameters of ABC algorithm

Figure 2 :
Figure 2: Parameters and model of MSVM in the Excel VBA platform.

Figure 3 :
Figure 3: Comparison of pressure data estimated by MSVM and borehole measured data.

Figure 4 :
Figure 4: MSVM-ABC-based back analysis, its parameters, and results in the Excel VBA platform.

Figure 5 :Figure 6 :
Figure 5: Comparison of recognized in situ stress using different models (  ,  ℎ are the maximum and minimum in situ stress).

Figure 7 :
Figure 7: Fitness variation with increased cycles in the ABC analysis.

Figure 8 :
Figure 8: Fitness distribution in the searching space at different cyclic stages.

Table 1 :
Comparison of in situ stresses and mechanical parameters between MSVM model and other preexisted models.