#### Abstract

Using the presented criteria for the breakdown pressure in which the breakdown pressure is related to the horizontal stress, the breakdown pressure in the hydraulic fracturing test is directly used in the estimation of the maximum horizontal in situ stress. However, the classical breakdown pressure criteria do not take into account the effects of near-wellbore stress field, so estimating the horizontal stress using them is accompanied by errors. A plane strain numerical model is presented to obtain the breakdown pressure, in which the effect of initial crack length are considered. In this model, it is assumed that the rock medium is impermeable and the breakdown occurs at the tip of two radially cracks. Thus, in this model, the effect of the initial crack length, which is not present in the classical models, is considered. The results obtained from the numerical model show the significant effect of the initial crack length on the breakdown pressure. Also, the comparison of the numerical model with the classical criteria shows that the condition of using the classical criteria to determine the breakdown pressure is that the ratio of the initial crack length to the borehole radius is very small. Using the numerical model, 1,456 datasets were prepared to train an artificial neural network to predict the maximum horizontal stress. Input parameters include breakdown pressure, minimum horizontal stress, initial crack length, and tensile strength. The evaluation of the results shows that the obtained neural network model has a good ability to predict the maximum horizontal stress.

#### 1. Introduction

Today, however, the applications of hydraulic fracturing in geomechanics and enhancing the geothermal systems are very important [1, 2], but in the petroleum industry, wells are stimulated by hydraulic fracturing technique, to increase the porosity and permeability of the rock. In addition to these applications, today hydraulic fracturing, which is sometimes referred to as hydrofrac or minifrac, is widely used in determining the in situ stresses in the Earth’s crust [3–12]. By using the results of a successful hydraulic fracture test, it is possible to determine the in situ stresses magnitude and their directions in the plane perpendicular to the axis of the borehole. In this method, first a section of the borehole that is a candidate for hydraulic fracturing is sealed off by straddle packer. Then, the fluid which is usually water, is injected at a constant flow rate. As the injection continues, the pressure gradually increases on the borehole wall until the first tensile fracture occurs along the minor principal stress due to the concentration of tensile stress. At this point, the injection is stopped and the borehole pressure is allowed to decay. The injection/shut-in cycle is repeated several times with the same fluid flow rate. Shut-in pressure and breakdown pressure, which are used to determine the in situ stresses, are obtained from the pressure-time record (Figure 1(a)). The breakdown pressure corresponds to the peak pressure, but the shut-in pressure can be determined in the different ways [13, 14].

**(a)**

**(b)**

The use of hydraulic fracturing technique in estimating the in situ stresses has advantages that other methods do not have; first, provided that the borehole is deep, using this method, a continuous stress versus depth diagram can be obtained. Second, there is no theoretical limit to the measurement depth provided that the borehole has access to the target zone and the brittle rock is in the elastic range.

For a vertical borehole the stress component along the borehole axis, which is equivalent to the overburden weight, is one of the principal stresses. Under such conditions, with the hydraulic fracture of the borehole, two vertical plane cracks are formed on both sides of the borehole along and perpendicular to (Figure 1(b)).

In the standard/classical method of determining the in situ stresses from the hydraulic fracturing test, the minimum principal stress is considered equivalent to the shut-in pressure [8]:where is the minimum horizontal stress and is the shut-in pressure. But the maximum horizontal stress can be determined from the minimum horizontal stress, breakdown pressure, and the tensile strength of the rock. The primary breakdown pressure model was presented by Hubbert and Willis [3], assuming linear elastic behavior and impermeable rock. According to this criterion, breakdown occurs when the minimum effective tangential stress in the borehole wall is equal to the rock tensile strength [3]:

In Equation (2), is the breakdown/initiation pressure, are the maximum and minimum horizontal stresses, respectively, is pore pressure in the medium, and is the rock tensile strength. Using this model, as long as there is no initial pore pressure the maximum horizontal stress is given by:

The second model of breakdown pressure is presented in accordance with the same concept of the previous model, with the difference that in this model the theory of poroelasticity is used to analyze the stress around the borehole [15]:where is the poroelastic stress coefficient given by is known as Biot’s poroelastic parameter and is Poisson’s ratio.

Using this model, when the rock is impermeable and there is no initial pore pressure , the maximum horizontal stress is given by:

Comparing Equation (3) with Equation (6), the two equations do not give the same result for the maximum horizontal stress because the breakdown pressure in the Equation (6) has a coefficient of 2; this is due to the effect of effective stress which causes the borehole pressure to be considered twice in the calculation. The tensile strength criterion can also be used in determining the initiation pressure and the location of the fracture initiation in noncircular boreholes [16].

Alternatively, shear failure in the borehole wall can also be used as a criterion for finding the breakdown pressure [17–20]. According to these models, when the shear stress in the borehole wall reaches the rock shear strength, borehole breakdown occurs. According to the shear failure models, the maximum horizontal stress can be obtained based on the rock shear strength and breakdown pressure.

The other category of models proposed for breakdown pressure is based on the theory of fracture mechanics [21–26]. In these models, it is assumed that when at the tip of an initial crack the stress intensity factor (mode I or opening mode) reaches the mode I fracture toughness , the crack starts to extend and breakdown occurs. The stress intensity factor indicates the gradient with which the stress at the crack tip tends to infinity. The stress intensity factor can be obtained by weighted sum of the given tensile stresses perpendicular to the crack surface. Rummel [22] provided an analytical solution for the breakdown pressure of two transverse hydraulic fractures in the borehole wall. In the proposed solution, the stress intensity factor is calculated for each of the stress components and the total stress intensity factor is obtained using the principle of superposition. Using this model, the maximum horizontal stress is given by Rummel [22]:where depends on the distribution of the pressure on the crack surfaces and initial crack length. For zero initial crack length they approach 3, −1, 1, and 0, respectively. By comparing Equation (7) with the standard/classical criteria (Equations (3) and (6)), the term , represents the tensile strength of rock material.

The shortcomings of the standard models made it difficult to estimate the maximum horizontal stress using the breakdown pressure, and this led to the use of other alternatives to determine from the hydraulic fracturing measurements [27]. For example none of the standard breakdown pressure criteria predicts the effect of pressurization rate and size, while, laboratory experiments show that breakdown pressure depends on borehole size and pressurization rate [28–30]. The effect of pressurization rate has also been reported in field tests, in particular Cleary [31] investigated the effect of pressurization rate on initiation pressure of the hydraulic fracture. Bunger et al. [10] used a mathematical–numerical model to investigate the ambiguities and errors associated with estimating in situ stresses using the borehole pressure record in hydraulic fracturing. According to their results, the use of standard breakdown pressure to obtain the major horizontal stress is based on the assumption that the initial crack length is very small compared to the borehole radius and that the crack initiation pressure corresponds to the breakdown pressure.

According to the above, most analytical/numerical models of hydraulic fracturing are so complicated that it is often impractical to use them in predicting the in situ field stresses, on the other hand, the standard models of breakdown pressure are also so simple that the estimation of the in situ stress by them is accompanied by an error. The main weakness of classical models based on rock tensile strength is that they do not consider the preexisting initial crack length in the borehole wall. The preexisting crack in the borehole wall can be caused by the intersection of natural cracks with the borehole or the perforation borehole. Hence, determining the maximum horizontal in situ stress using classical models will be accompanied by errors. These shortcomings are the main motivation of the work presented in this article. In this paper; first, a plane strain numerical–analytical model to obtain the breakdown pressure of hydraulic fracture is presented. Initiation of hydraulic fracturing occurs at the tip of two symmetrical fractures transverse to a borehole. In the numerical model, the medium is assumed to be impermeable with linear elastic behavior. The J-integral technique is used to calculate the stress intensity factor at the crack tip, which is more accurate than the other methods. Furthermore, with this technique, the distribution of pressure on crack surfaces can be considered uniform or nonuniform. This simple model examines the effect of parameters such as initial crack length on breakdown pressure that standard models did not consider.

In the second part of the paper, using the validated numerical model, 1,456 numerical analysis with different values for in situ stresses, dimensionless initial crack length, and fracture toughness are performed. Using 1,456 datasets, a neural network model was developed that accurately predicts maximum horizontal stress by using the input data. Finally, using sensitivity analysis, the effect of each of the input parameters on the maximum horizontal stress is determined.

#### 2. Analytical–Numerical Model, Breakdown Criterion

A borehole with radius and two symmetrical planar cracks of length is considered in an impermeable linear elastic medium assuming plane strain condition (Figure 2(a)). The in situ stress field is nonisotropic and the planar cracks are along the maximum horizontal stress . In this article pressure distribution inside the cracks is assumed to be uniform and equivalent to the borehole pressure. The governing equations are the same as the basic equations in the theory of elasticity and linear elastic fracture mechanics (LEFM), including the relation between strains and displacements:the equilibrium equation:and the stress–strain relation [32]:where and are the stress and strain matrix, respectively, is displacement vector, is body force vector, is Kronecker delta which is a symmetric symbol, and are Lame’s constant and shear modulus, respectively. Finite element numerical technique is used to solve the equations, in summary in this method, the general equation of equilibrium is as follows given by Reddy [33]:where is the global stiffness matrix, is body force vector, and is the loading vector.

**(a)**

**(b)**

According to the fracture toughness criterion from LEFM, fracture initiation occurs when [34]:where is the mode I stress intensity factor at the tip of preexisting crack and is the mode I fracture toughness of the material. In LEFM, the stress intensity factor concept at the crack tip is used to determine the stress, strain, and displacement field. For this reason, accurate evaluation of the stress intensity factor is very important in finite element analysis within the framework of LEFM.

Techniques for calculating the stress intensity factor are different, which are divided into two main categories, direct methods and energy-based methods. Direct methods are more common due to the simplicity of their relationships, but energy methods which are based on calculating the energy release rate, are more accurate. The J-integral technique was first introduced by Rice [35] to calculate the rate of energy release in crack problems. Rice [35] showed that the rate of energy release can be written as a path-independent integral:in this equation, is the strain energy density:and is the unit normal vector to the contour integral (Figure 2(b)), is the traction vector on the contour integral .

It can be shown that for a very small virtual crack, the *J*-integral represents potential energy changes, so for linear elastic materials, the energy release rate is equivalent to the *J*-integral. The energy release rate in two-dimensional problems is also obtained based on the stress intensity factor as follows [35]:where in the plane strain condition . After meshing the problem domain, the *J*-integral can be obtained from the finite element mesh. To increase the accuracy in calculating the stresses, it is better to pass the contour through the Gaussian integration points of the elements. There are also other conditions to consider when choosing an integration domain; first, in complex crack patterns, it is appropriate to be near to the crack tip, second, it should be easy to implement it in an automated procedure of simulation and finally, it must be consistent with the boundary conditions of the problem and geometry. A crack is defined as an edge or face with over-lapped duplicate nodes that can separate during analysis. To calculate the *J*-integral, the contour around the crack tip is considered. Based on the governing equations and using the finite element numerical method and the *J*-integral technique within the framework of LEFM, a program was written in MATLAB software to solve the problem. According to Figure 2(a), for a given value of in situ stresses and the initial crack length, borehole pressure, and pressure on crack surfaces gradually increase until the stress intensity factor equals fracture toughness, i.e., , at this time the borehole pressure is equal to the breakdown pressure.

#### 3. Validation, Comparison of Numerical Model with Integral Formula

In this section, the analytical–numerical model presented in the previous section and the code written in MATLAB are validated. For validation, a comparison is made between the mode I stress intensity factor obtained from the numerical model by the *J*-integral technique and the stress intensity factor obtained from the integral formulas presented by Nilson and Proffer [36]. Based on the weight function method, Nilson and Proffer [36] proposed an integral formula for calculating the mode I stress intensity factor for planar cracks originating from a wellbore in an infinite elastic medium:where is the hole radius, is the fracture length, *x* is distance from the fracture tip, is the pressure distribution, and is the external confining pressure. The weight function is defined as follows:where for plane strain conditions considered in this article and

In this relation, for cylindrical cavity and for spherical cavity.

Figure 3 shows the model geometry, which shows the boundary conditions, loading and two symmetrical exaggerated cracks. The dimensions of the model are with a hole with a radius of in the middle. The meshing is done using isoparametric eight-node elements and the mesh is finer around the boreholes and cracks, to increase accuracy.

Figure 4 shows that a good agreement between the numerical model and the integral formula for the stress intensity factor of opening mode is obtained. In this figure the stress intensity factor is normalized to the stress intensity factor of a Griffith crack . Comparisons are made for three different loading configurations; when the pressure is restricted to the hole, when pressure is applied to crack surfaces and the hole, and in the third case, the model boundaries are subjected to the uniform tensile stress. Because in all three configurations, there is a good agreement between the numerical model and the integral formula, so the numerical model should have a good performance in obtaining the breakdown pressure of hydraulic fracturing.

#### 4. Parametric Study on Breakdown Pressure

In this section the results of the parametric study on breakdown pressure are presented, the parameters given in Table 1 have been used as base case parameters. Various laboratory studies have shown that the mode I fracture toughness of rock materials is directly related to its tensile strength [37, 38]. Zhang [39] showed that in general for all types of rock from soft to hard, in low speed impact loading or quasistatic loading, there is the following relationship between tensile strength and the mode I fracture toughness:

In this paper, using this relationship and for a given tensile strength, the fracture toughness is obtained.

Figure 5 shows the breakdown pressure versus the dimensionless initial crack lengths for different deviatoric stresses . It is observed that at low-deviatoric stresses, the breakdown pressure of hydraulic fracture decreases with increasing . The rate of reduction of breakdown pressure is high up to and then it takes a slower trend. At large deviatoric stresses with increasing , the breakdown pressure does not decrease monotonously. For , all curves in Figure 5 converge to the minimum horizontal stress and the effect of deviatoric stress disappears; this is because at distances away from the borehole, the stress concentration effect disappears and the breakdown pressure is only affected by the minimum horizontal stress. Figure 6 is another representation of Figure 5. In this figure, it is also observed that with increasing deviatoric stress, the breakdown pressure decreases.

Figures 7 and 8 show the breakdown pressure versus dimensionless initial crack length for different fracture toughnesses and for two deviatoric stresses of and , respectively, and Figure 9 illustrates the evolution of breakdown pressure versus fracture toughness. It is seen that the evolution of breakdown pressure versus fracture toughness and deviatoric stress is linear, but the breakdown pressure variation is nonlinear versus dimensionless initial crack length.

#### 5. Comparison with Tensile Strength-Based Breakdown Pressure Models

It is now possible to evaluate the standard breakdown pressure models given in Equations (3) and Equation (6) using the proposed model. Evaluation of standard breakdown pressure models based on the tensile strength are very important because these simple models are usually used in the field to estimate . In the numerical model presented in this article, it is assumed that the effects of poroelasticity are negligible, and the pressure distribution on the crack surfaces is uniform and is equivalent to the borehole pressure, this is almost a condition caused by the injection of a nonviscous fluid in the field [40]. Because the fluid pressure is applied to the all crack surfaces, the model can be compared with the breakdown pressure equation presented by Haimson and Fairhurst [15] (H-F criterion) (Equation (6)).

In breakdown pressure models based on the tensile strength of the rock, it is assumed that the tensile stress is uniformly applied across an infinitesimally small crack of length , therefore, according to the theory of fracture mechanics, considering a finite crack of length under uniform tensile stress , with the assumption where for the edge crack [41]. Provided that , this relation can be substituted in Equation (6) and the other form of the H-F criterion is obtained as follows, in which there is the effect of the initial crack length :

For the base case parameters given in Table 1, the comparison result for is given in Figure 10. In this figure, Curve 2 is the breakown pressure obtained from the numerical method presented in this paper, and Curve 1 is the breakdown pressure obtained from the H-F model in accordance with Equation (6), and Curve 3 is also obtained from H-F model but based on Equation 21.

The H-F criterion according to Equation (6) (Curve 1), because it does not consider the effect, therefore provides a constant value for the breakdown pressure. The H-F criterion based on Equation (18), for intersects with the curve obtained from the numerical model, and for , there is a perfect match between the H-F criterion and the numerical model. However, as increases, the difference between the H-F criterion and the numerical model increases, because, as mentioned, Equation (21) is valid as long as . For example, for , the breakdown pressure obtained from the H-F criterion according to Equation (18) has a error compared to the numerical model, and the breakdown pressure obtained from the H-F criterion according to Equation (6), has an error of compared to the numerical model, and when increases to 0.15, the breakdown pressure obtained from Equation (18) has a 4.65% error and the breakdown pressure obtained from Equation (6) has an error of 8.5%.

#### 6. Neural Network Model for Predicting Major Principal Stress Using Breakdown Pressure

Using the numerical analytical model rendered in this article, maximum horizontal principal stresses can now be predicted by through an artificial neural network (ANN). ANN systems behave similarly to the human brain, this means that they need a learning process using a specific set of input and output data. After neural network training, the network can be used to predict other output parameters using the input data. Using more items for training allows for more reliable network predictions. Another similarity between the neural network and the human brain is related to the internal architecture or structure of a network. There are different types for network architecture, but primarily all of them are composed of perceptrons or nodes that are connected to each other. Nodes in the human brain are called neurons. Each perceptron processes the signals of input parameters and transmits the result to the output neurons. One of the classical network architectures used in the engineering analysis is feed forward backpropagation architecture. In this architecture, there is no connection back between the output neurons and input neurons and the network does not recall its previous output values [42]. This architecture has three layers; input layer, output layer, and hidden layer (Figure 11). Each neuron in the hidden layer is fully connected to all input neurons to receive data, and to transmit processed data must also be connected to output neurons.

BIAS neurons are used to correct bias, they do not receive any input data. Network training is done using an iterative process through synaptic weights, with automatic fitting. In each connection, the data are weighed independently via the reciprocal synaptic weight, thus, each neuron related to the output accepts a distinct value. The advantage of neural networks over analytical fit functions such as the least squares method or regression is the use of a large number of fitting parameters or synaptic weights in the neural networks [43, 44]. Another important advantage of neural networks is that they are not dependent on simplified assumptions such as linear behavior. In the training, synaptic weights are determined using iterative process to create a minimum mean squared error betwixt the network output data and actual target data provided that the inputs are the same. In general, the reply of a neuron is written as follows:where is the result obtained from neuron *j* of current layer, is transfer/activation function, *r* is the number of neurons in the foregoing layer, and is the result of neuron *i* from the previous layer, is the weight of the connection , and is the bias weight. Usually the number of hidden layers and their neurons is obtained by trial and error. The network is optimal when the coefficient of determination between the predicted and the actual is close to 1. on the other hand, its error in predicting is close to zero. The coefficient of determination , the root-mean-square error (RMSPE) and mean square error (MSE) relationships are given below, respectively:where *m* is the number of datasets.

In this paper, after testing several activation functions, the activation function of sigmoid in the interval (0, +1) was chosen for the hidden layer; besides, purelin function which is a linear transfer function was used for the output layer. These functions are defined by Equations (26) and (27), respectively (Figure 11):with these two activation functions, the lowest MSE value and the highest convergence were obtained.

#### 7. Optimal Neural Network Architecture and Results

To obtain a neural network-based model for predicting maximum horizontal stress, 1,456 numerical modeling with different input parameters were built and run and breakdown pressure of the hydraulic fracturing was obtained for each model. The input parameters of the model include minimum horizontal stress , breakdown pressure , dimensionless crack length , and tensile strength of the rock , and the output of the model is the maximum horizontal stress (Figure 11). Input and output data and their range are given in Table 2.

It should be noted that out of 1,456 datasets, 70% of them were randomly used for network training, 15% for network test, and 15% for network validation. In this work, backpropagation algorithm is used to train the network.

Figure 12 shows normal probability distribution and histogram for the input and output variables, from this figure, the distributions of input and output variables is not completely normal and are not similar to each other, indicating that the relationship between output and input variables is nonlinear. Thus the neural network model has an output layer with one neuron, an input layer with four neurons, and a hidden layer. To find the optimal number of hidden layer neurons, several analysis were performed with different neurons and it was observed that with increasing the number of hidden layer neurons, the error decreases and the coefficient of determination increases. In accordance with the results of the analysis given in Table 3, since there is not much difference between topology 4-20-1 and topology 4-50-1, 20 neurons were selected for the hidden layer. Consequently, to predict , a neural network with 25 neurons and 100 synaptic connections is created (Figure 13). The coefficient of determination for this network is 0.9997 and the MSE is equal to 0.16 and the RMSPE is equal to 0.00002 (Table 3). The optimal architectural weights and biases related to connecting the input layer to the hidden layer and the hidden layer to the output layer are given in Tables 4 and 5, respectively.

#### 8. Evaluation of Neural Network Performance

The error histogram for the datasets used in the neural network test is shown in Figure 14, most of the errors are around 0 which is a good indicator of the accuracy of the ANN model in predicting . Figure 15 is also a comparison between the predicted by the neural network model and the obtained from the numerical method, which for the sake of brevity, the comparison is only done for 100 data from the test data. The observed good match indicates a good prediction made by the neural network model.

Figure 16 illustrates the relationship between the actual maximum horizontal stress obtained from the numerical analysis and predicted by the ANN model, in this figure, the line *y* = *x* is also shown for comparison. Using simple linear regression, the following equation shows the relationship between the values obtained from the numerical solution and the predicted values:

#### 9. Strength of Relationship between Neural Network Input Data and Output Data

In the neural network model presented in this article, four groups of input data are involved; minimum horizontal stress , breakdown pressure , dimensionless crack length , and tensile strength of the rock . The important question that arises now is: what is the value and impact of each of the model input parameters on the model output ? A useful way to determine the strength of the relationship among neural network input data and output data is the cosine amplitude method (CAM). The calculation process in the CAM method is such that first an *n*-dimensional space is considered:where *n* is the sum of the input and output parameters of the model, each element in this space is a vector of dimension *m*:where *m* is the total number of datasets, the relationship strength among the model input parameter and the model outputs is given by the participation value expressing that strength:where , and . Equation (31) is similar to the inner product in the cosine function, when the vectors are collinear, their inner product is unity, and when two vectors are perpendicular to each other, their inner product is zero. The results of the sensitivity analysis of the neural network model and as the output of the model are shown in Figure 17. As can be seen in this figure, in the numerical analysis of hydraulic fracturing presented in this article and for the range of values investigated, is the most sensitive parameter affecting the . In the classical relationships presented for estimating based on breakdown pressure (Equations (3) and (6)), it is also observed that has a coefficient of 3. Breakdown pressure has a larger effect on than the tensile strength of rock , in the relation presented by Haimson and Fairhurst [15] (Equation (6)), it is also seen that the breakdown pressure has a coefficient of 2 but, the rock tensile strength has a coefficient of unity. The dimensionless crack length is the least influential parameter on for the values explored. Numerical analysis performed in this article indicated that for small values of , the breakdown pressure obtained from the numerical model is close to the breakdown pressure obtained from the relationship proposed by Haimson and Fairhurst [15]. Therefore, has the least effect on among the other input parameters of the ANN model.

#### 10. Multiple Linear Regression

In this section, multiple linear regression (MLR) analysis is performed to find the relationship between and input parameters of the problem. By normalizing 1,456 datasets from numerical analysis performed in previous sections, using SPSS software, the relationship between and input parameters of the problem was obtained as follows:

The coefficient of determination of this relationship is , which is less than the coefficient of determination obtained by the neural network model . Figure 18 is a comparison between the predicted by the MLR and the obtained from the numerical method, since the coefficient of determination of the MLR is lower than that of the neural network model, the error of this method in prediction of is greater than that of the neural network model. Nevertheless, in field applications, it can be useful to use Equation (28) in the initial estimate of the maximum horizontal stress using the breakdown pressure.

#### 11. Conclusion

In this article, by presenting an analytical–numerical model, first, the breakdown pressure of hydraulic fracturing and the parameters affecting it were examined. In this model, the stress intensity factor at the crack tip is calculated by the *J*-integral technique, which is more accurate than the other methods. The proposed model is solved by finite element numerical method and using the code written in MATLAB software. Although any form of pressure distribution can be considered on the crack surfaces, in this article, the pressure distribution on the crack surfaces was considered uniform and equivalent to the borehole pressure so that the result can be compared with the criterion provided by Haimson and Fairhurst [15].

According to the result obtained from the first part of the article; first, the numerical model showed a significant dependence of the breakdown pressure on the length of the preexisting initial crack, while the effect of initial crack length and in general, the size effect, is not present in classical models based on the tensile strength of rock. Second, by comparing the results of the validated numerical model with the model based on tensile strength presented by Haimson and Fairhurst [15], it is observed that the Haimson and Fairhurst [15] classical model for has a good prediction of the breakdown pressure. But as increases, the breakdown pressure obtained from the model based on the tensile strength differs from the breakdown pressure obtained from the numerical model. In other words, the range of validity of models based on tensile strength is for . This result is very important in field applications, since classical criteria are usually used in the field to estimate the in situ stresses using hydraulic fracture breakdown pressure.

In the second part of the paper, using 1,456 numerical analysis, complete datasets with various input parameters were prepared. Then an ANN model was trained, whose input parameters include minimum horizontal stress, breakdown pressure, rock tensile strength, and dimensionless crack length, and whose output parameter is maximum horizontal stress. The obtained neural network model has 25 neurons and 100 synaptic connections, which has good accuracy in predicting the maximum horizontal stress. The coefficient of determination of this model is 0.9997 and the RMSPE is 0.00002. In this way, by using the neural network model for each set of given input parameters, the maximum horizontal in situ stress can be obtained.

In summary, hydraulic fracturing breakdown pressure is dependent on several factors including initial crack length, and on the other hand, it is a key factor in determining in situ stresses. Hence, breakdown pressure criteria that are based only on the tensile strength of the rock cannot accurately determine the maximum horizontal in situ stress. Criteria based on fracture mechanics, although they are more accurate and consider all the factors involved in hydraulic fracturing, but due to their complexity, they cannot be used simply in practice. Therefore, techniques based on the machine learning and neural networks, if they are trained using accurate numerical or field data, can be a suitable alternative to simple the classical models or complex models of the fracture mechanics.

#### Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Authors’ Contributions

All authors contributed to the study conception and design. All authors read and approved the final manuscript.

#### Acknowledgments

The authors would like to acknowledge the University of Zanjan for providing the necessary resources and facilities for this research.