Regular Article
Stability analysis and uncertainty modeling of vertical and inclined wellbore drilling through heterogeneous field
^{1}
DEP/FEM, University of Campinas (UNICAMP), PO Box 6122, 13083860 Campinas, São Paulo, Brazil
^{2}
DES/FEC, University of Campinas (UNICAMP), PO Box 6143, 13083889 Campinas, São Paulo, Brazil
^{3}
CERMICS (ENPC), Université ParisEst, 77455, MarnelaVallée cedex 2, France
^{*} Corresponding author: nbatalha7@gmail.com
Received:
10
May
2019
Accepted:
20
January
2020
A stochastic twodimensional geomechanical model developed by the authors and presented herein is used to analyze wellbore stability in heterogeneous formations. It consists of a finite element model and assumes linear elastic and isotropic material behavior under the plane strain state. The model simulates the stress state around vertical and inclined wellbore when a formation is submitted to internal drilling fluid pressure. This new state of stress may lead to rock failure, which is analyzed through a failure criteria. Since the exact variation of formation’s mechanical properties is not known, a spatially correlated field is used to evaluate the variability of a rock mechanical material property. The correlation between each pair of finite elements is determined by a covariance function. A twodimensional spatially correlated field is used to verify the correlation between the elements of a vertical wellbore. A different approach, however, is necessary to model inclined wellbores. Once the direction and inclination of a wellbore are defined, a threedimensional spatially correlated field becomes necessary to best simulate the formation field. Simulations using the stochastic model proposed herein and considering constant elastic modulus have been compared. It is observed that when considering the same elastic modulus along the field, the area of the plastic zone is symmetric at the borehole wall; when using the stochastic model, however, the plastic zone area surrounding the well is not symmetric; thus, the most vulnerable plastic zone may lead to premature failure. Stochastic simulations have been carried out with different heterogeneous fields for vertical and inclined wells, and distributions of the total plastic zone areas are obtained for each case. Based on the stochastic field probabilistic analysis, a distribution function is presented, which aims to define a stability framework analysis to best assist a decision making process in determining if a mud pressure is operationally acceptable or not.
© N.A. Batalha et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In order to predict the formation behavior in drilling operations, numerical models have been developed to simulate the stress state around vertical and inclined wellbores [1–10]. These models intend to define the best drilling fluid density (or mud weight) to be used. The mud weight must be designed to keep the well stable and not cause rock failure. Insitu stresses magnitude and orientation, as well as rock strength and pore pressure, are some of the parameters required to a geomechanical approach in wellbore stability, as presented by [11]. The precision and variation of these parameters over a domain are generally not considered, but average values are usually obtained from multiple laboratories or onsite tests. Because of that, most models assume that the properties of the material are constant in the entire simulated domain. Thus, wellbore stability analysis may involve many uncertainties; to list a few: definition of formation’s mechanical properties, pore pressure, mud weight, system temperature, and insitu stresses. Therefore, analyzing the influence of these variables requires probabilistic methods and a consistent framework.
A probabilistic analysis considering the uncertainties of insitu stresses, pore pressure, and rock strength as input data in a geomechanical model was evaluated by Udegbunam et al. [12]. Monte Carlo simulations were used, and triangular and constant distribution types were applied for the random variable distributions in order to predict critical fracturing and collapse pressures. It was observed in this study that the minimum horizontal stress has a high impact on the cumulative uncertainty in the fracture pressure prediction. Cohesive rock strength was pointed as the variable that most influences the collapse pressure.
Tran and Abousleiman [13] studied how the mechanical properties, insitu stresses, and formation pore pressure conditions affect wellbore stability when using MohrCoulomb, DruckerPrager, Modified Lade, and tensile failure criterias. Monte Carlo method was used in this study, and different farfield regimes were also evaluated. Regarding the mechanical properties, cohesion and friction angle presented the most impact on collapse mud weight. It was noted in this work that the insitu stresses and formation pore pressure had more influence in tensile failure than rock tensile strength.
Fontoura et al. [14] analyzed the influence of the uncertainties in wellbore stability using Statistical Error Analysis Method (SEAM), First Order Second Moment (FOSM) and First Order Reliability Model (FORM) methods, which were compared with Monte Carlo method. Based on a sensitivity analysis, it was noted that formation cohesion and angle of internal friction are essential parameters for collapse failure, and, on the other hand, tensile strength is relevant for fracturing failure. Also, insitu stresses and formation pore pressure are considered important parameters in wellbore stability analysis, as well as Biot’s constant. In the case of vertical wellbore stability, the vertical stress was not relevant.
A numerical model for wellbore stability coupled to hydromechanical process was presented by Muller et al. [15], in which stochastic and deterministic analyses were compared; the spatial variability of initial stresses, pore pressure, and hydraulic and mechanical properties were taken into account. A reliability analysis was carried out limited to a pressure range, where lower limit, compressive failure, is based on a user defined plastic region, while the upper limit is based on the rock tensile failure stress. This study was later analyzed for threedimensional scenarios [16].
Quantitative risk assessment has also been used in wellbore stability uncertainties analysis [11, 17, 18]. The latter reference presented an analytical approach for wellbore stability uncertainties analysis, which takes into account the uncertainties of insitu stresses, operational activities, and rock mechanical properties. According to the author, their results have shown that cohesion, internal friction, and pump pressure were not significant parameters in the total risk of the drilling operation.
Probabilistic models have led to high impact in wellbore stability analysis, especially when sensitivity analysis is taken into account to define input parameters. The influence of elastic modulus heterogeneity on insitu stresses and damage in shale gas reservoirs was recently studied by Wang et al. [19], in which variability of von Mises stress, strain energy density, and damage factor was observed as heterogeneity increased.
Therefore, to obtain a more realistic wellbore stability analysis, a heterogeneous field shall be considered. Since the variation of rock mechanical properties along the field is seldom known due to difficulties to precisely define these properties, a random but correlated field can be assumed to carry out a probabilistic analysis. The present work introduces a geomechanical model coupled to a geostatistics algorithm to predict heterogeneous formation behavior in drilling operations.
The geomechanical model presented herein assumes linear elastic and isotropic mechanical properties and planestrain state. A yield criterion is selected to identify the formation’s failure area, plastic region, around the borehole wall. The stochastic field is based on a covariance function, which specifies the correlation between the mechanical characteristics of neighboring elements of a domain. This stochastic analysis approach enables a more realistic evaluation of the wellbore failure; the model is also extended to analyze vertical, horizontal, and inclined wellbore.
2 Rock and drilling fluid: a mechanical interaction
A stable and welldesigned wellbore guarantees operation safety and prevents several issues such as stuck pipe, cavity enlargement, fracturing formation, operation delays, and substantial expenditures. Note that the stability of a wellbore depends on equilibrium of stresses (insitu stresses), material properties, and drilling fluid’s density, also called mud weight.
Insitu stresses are naturally developed during a rock formation, and they are a consequence of gravity interactions, tectonic processes, and others. These stresses are composed by the overburden pressure σ_{ov}, defined by the weight of all rock layers over a rock element, and, in response to the overburden stress (or vertical stress σ_{ov} = σ_{V}), there are the horizontal stresses, Minor Horizontal Stress (σ_{h}) and Major Horizontal Stress (σ_{H}), which avoid the lateral deformation induced by vertical loading. Estimating the ratio between in situ rock stresses, it is assumed that vertical and horizontal stresses are principal values as reported and discussed in [20]. In this research, it is assumed the same considerations with a fixed stress ratio. Insitu stresses are subjected to the static equilibrium equation:
where, σ is the stress tensor and b is the vector of the body forces.
A new state of stress is created around the wellbore wall by the drilling process. This new stress state may generate significant change in formation integrity, leading to rock failure and operational problems [6, 21]. The formation around the borehole wall experiences deformation due to the relief of stresses, once the rock is removed, and then, replaced by the drilling fluid. The drilling fluid (or drilling mud) is a mixture of a fluid (water or oil based) and solid [22], used in the drilling operation to avoid the borehole closure and to fill the open space once occupied by the rock. The drilling fluid temporarily supports the borehole wall while drilling. The drilling fluid’s density must be high enough to prevent flow of formation fluids into the wellbore (a phenomenon named as “kick”) as well as wellbore closure, but low enough to prevent hydraulic fracture, formation failure and circulation loss [23]. Thus, the mud weight (or drilling fluid density) is responsible for keeping the wellbore walls stable. In order to minimize the cited operational risks, geomechanical models have been used in the design and drilling operations stages. Understanding the new stress state of vertical, inclined, and horizontal wellbore is essential to ensure a safe operation. The stress state simulation allows the definition of the optimal drilling fluid density to be used in a certain formation.
3 Methodology
A geomechanical model is introduced in order to simulate the formation behavior when mud pressure is applied to the borehole wall. The numerical model presented herein is developed using the Finite Element Method (FEM) and assumes linear elastic and isotropic material behavior under the planestrain state. The method is implemented in NeoPZ, which is an opensource library for the development of finite element simulations [24]. A yield criterion is selected to describe the formation’s plastic region at the borehole wall. The geostatistics analysis is developed based on a stochastic field model. The methodology used to implement this model is presented in this section.
3.1 The geomechanical model
Once the wellbore length is very long (kilometers long sometimes) compared to the size of the studied domain in x and y directions, a plain strain hypothesis is relevant considering null displacement in the zdirection. The simulation is, therefore, twodimensional, which significantly reduces computational time. As depicted in Figure 1, the finite element mesh is selected as a crosssection of the wellbore. The borehole cavity is the central circumference, and the remaining area refers to the original formation being drilled.
Fig. 1 Twodimensional mesh surrounding the wellbore. 
The FEM algorithm starts at finding the displacement u (x,y) from the equilibrium equation,
where,

σ is the stress tensor;

b is the body force vector;

u is the displacement vector;

u _{ D } is a known function;

σ _{ InSitu } is the insitu stress;

Ω is the domain, and applying:

δΩ_{D} and δΩ_{N} are the boundary conditions, Dirichlet and Neumann, respectively.
The space of testing functions is defined by:
where,
By multiplying the equilibrium equation by the testing function v , integrating over the domain Ω, and using the divergence theorem, we have that:
Assuming that between the elements,
and applying,
The forming function can be stated as:
Given that:
and by the stress–strain relation for linear elastic condition, we find the stress tensor:
with C being the elastic constant matrix, σ _{ 0 } the initial stress tensor and ε the strain tensor.
DiMaggioSandler yield criterion [25] is selected and implemented in the geomechanical model to describe the failure area/plastic region surrounding the wellbore. This failure criteria is chosen once it describes well the constitutive behavior of rocks [26–28]. The failure criteria yield function Φ is composed of two functions. In the present analysis, however, only the failure surface envelope function is considered. Therefore, the envelope which limits the elastic deformation under confinement is disregarded. Plastic deformation occurs when Φ > 0. In the present study, only the collapse failure is analyzed herein. For the sign convention, compressive stress is negative and tensile stress is positive.
It is important to mention that the present geomechanical model has been validated. An analytical approach is fully described in [29] to determine the stress state around wellbore for isotropic homogeneous rocks. This analytical approach, however, does not take into account rock elastic parameters such as Young’s Modulus and Poisson’s ratio, except by the axial stress that considers Poisson’s ratio in its equation. Thus, this work also aims to verify the influence of rock elastic modulus uncertainty as an input parameter. The geomechanical model solution presented herein was compared to this analytical approach; the comparison presented satisfactory agreement, as fully detailed in [30].
All analyses are based on a FE mesh of 750 elements in total, consisting of 25 radial elements and 30 circumferential elements; a quadratic approximation is used, and Dirichlet boundary condition is applied in the external part of the mesh to avoid rigid body displacement. The mud weight, or drilling fluid density pressure, is applied at the wellbore wall, and the formation’s initial stress state is defined by its insitu stresses.
3.2 Spatially correlated field for heterogeneous formation
Heterogeneous rock mechanics properties may significantly influence the formation behavior, and, as a consequence, also affect the mud weight to be used in the operation. Probabilistic evaluation of a given field allows a more realistic analysis of the rock stress state, and assist, with better accuracy, selecting the best density for the drilling fluid. The pore pressure and fracture gradients are very narrow, and a small change can affect wellbore stability [5]. Therefore, this stochastic field modeling represents an important issue to be addressed in wellbore stability analysis.
Since the model is a plane model, the simulation is set for a specific layer or rock formation. In other words, it assumes one average mechanical property for the whole field. Therefore, the stochastic model is built with one layer at a time.
The finite element mesh is built using a geometric progression, which leads to a greater mesh refinement close to the wellbore wall, where the greatest change in stress state is found. The stochastic grid follows the same geometry of the finite element mesh. Thus, each element receives a random but correlated value of Young’s Modulus. Following the geometric progression, the variation is greater closer to the wellbore wall than over the remaining domain. Nevertheless, it is also possible to not use the geometric progression, and then the elements would be equally distributed, and the Young’s Modulus would also be uniformly defined in the mesh.
A spatially correlated field is used to evaluate the variability of the formation’s material properties. The correlation between each pair of elements is defined by a covariance function. In this model, the random variable is the formation’s elastic modulus (Young’s Modulus). The Squared Exponential covariance function is chosen for this model, which is often used to describe soil properties autocovariance [31]. The covariance function, also known as “kernel,” is as follows, from [32]:
The selected kernel allows the user to control the distance for which the correlation between the elastic modulus of two elements is not relevant. The control is set by the constant γ, which is associated with the spatial frequency in the field [33]. Therefore, the scale γ defines the smoothness of the stochastic field.
The covariation function also depends on the distance between the centroids of the elements, defined by r. A matrix R_{m,m} holds all distances between each pair of elements, calculated by their elements centroids coordinates (X):
The correlation matrix K_{m,m}, holds the correlation between all pair of elements, and it is a function of the matrix R and the scale γ:
The matrix K is then decomposed by the singular value decomposition method, such that . The product of the left singular vector U and the square root of the diagonal matrix S, multiplies a Gaussian random and uncorrelated vector d _{ rand }, generating a correlated vector d _{ corr }.
The random but correlated vector d _{ corr } is then scaled multiplying it by a standard deviation (SD_{E}) and summed to an average Young’s Modulus (E_{avg}).
E_{stochastic field} refers to a stochastic but correlated field of Young’s Modulus. The variable SD_{E} is the standard deviation of the Young’s Modulus, and E_{avg} is the Young’s Modulus, both given as input. Figures 2 and 3 depict for the Young’s Modulus variable, the difference between a random correlated field, which will be considered in the analysis, and a simple random and not correlated field.
Fig. 2 Young’s Modulus (MPa) of a random correlated field. 
Fig. 3 Young’s Modulus (MPa) of a random but not correlated fied. 
3.3 A 3D spatially correlated field for inclined wellbore
A new approach must be considered for inclined wellbores. Inclined wellbores are obtained from drilling in a given direction and inclination angle. Inclination angle is the angle between the wellbore axis tangent and the local gravitational vector. Therefore, a vertical wellbore has an inclination angle of 0°, and a horizontal wellbore has an inclination angle of 90°. Direction (or azimuth) angle consists of the angle between the horizontal wellbore projection and the true geographic north. This angle is commonly referred to as “wellbore azimuth.” It varies from 0° to 360° and is measured in the clockwise direction starting at the geographic north [34].
Fjaer et al. [29] present a transformation of the insitu stresses to the local stress tensor of an inclined wellbore that can be obtained from a rotation between the z and y axes. This operation can be represented in matrix notation as follows:
where,

σ _{ InSitu } is the insitu stresses;

α is the wellbore direction/azimuth;

β is the wellbore inclination;

Q is the rotation matrix, composed by the rotation between the z and y axis.
In which Q and σ _{ InSitu } are composed, in matrix notation, by:
A twodimensional correlated field is considered for vertical wellbore since the mesh is perpendicular to the well. Inclined wellbore’ mesh is also perpendicular to the wellbore orientation, but embedded in the layer that the geomechanical model is being analyzed. Figures 4 and 5 illustrate the difference between the two cases.
Fig. 4 Vertical wellbore mesh. 
Fig. 5 Inclined wellbore mesh. 
The geomechanical model also assumes a twodimensional approach for inclined wellbores. The simplification is possible due to the local stress tensor transformation, previously described. For the stochastic analysis, however, the mesh is inclined compared to the horizontal formation layer, according to the wellbore inclination and direction. In order to evaluate the heterogeneous field of a drilled layer and determine the correlation between elements, a threedimensional correlated field is developed.
The threedimensional approach follows five steps: (1) definition of the height and radius of a cylindrical layer, (2) selection of a twodimensional mesh as if it had been in the vertical model, (3) selection of the number of cross sections inside the cylindrical layer using the same twodimensional mesh, (4) rotation of an extra twodimensional mesh inside the cylindrical layer following the wellbore direction and inclination, and (5) definition of the coordinates of all elements centroids. All the steps can be seen in Figure 6. All meshes are discretized in the same way, as well as the inclined mesh.
Fig. 6 Five steps of the 3D spatially correlated field methodology. 
The correlation matrix K_{m,m} is, therefore, calculated for all pairs of elements centroids. The stochastic field over the inclined mesh (in red) is correlated with inclined neighboring elements, which allows the use of the twodimensional geomechanical model. The methodology described herein to create threedimensional correlated layers, allows the user to generate drilling fields expeditiously.
4 Results and discussions
This section presents the results of a case study carried out using the numerical model and analysis method proposed herein. The wellbore in question is a carbonatic rock, and it is already known that the selected drilling fluid pressure (P_{w}) would induce a plastic zone at the wellbore wall. The wellbore parameters and rock features for the yield criterion are reported in [26]; the model proposed herein is validated against the deterministic result in [30].
The input values used in the simulations are summarized in Table 1, where r_{w} and r_{e} are the wellbore radius and mesh’s external radius, respectively; the same values are used in the vertical and inclined wellbore. The carbonatic material parameters for DiMaggioSandler criterion are A: 152.52; B: 0.0015489, and C: 146.29. The plastic zone area is calculated using the determinant of the Jacobian matrix from each integration point.
Geomechanical model input data.
4.1 Plastic zone for vertical and inclined wellbores
The plastic zone (obtained based on Φ – yield criterion envelope) at the borehole wall are compared between vertical and inclined wellbores considering homogeneous media. The analyses input follow Table 1, and the inclination and azimuth angles analyzed are reported in Table 2.
Wellbores orientation.
Figures 7–10 depict the yield function for V1 and I1 wellbores according to Table 2 nomenclature. Note that for the evaluated mud pressure of P_{w} = 19.5 MPa, the plastic area for the inclined wellbore is smaller than for the vertical wellbore.
Fig. 7 Yield function of wellbores V1. 
Fig. 8 Plastic zone of wellbores V1. 
Fig. 9 Yield function of wellbores I1. 
Fig. 10 Plastic zone of wellbores I1. 
Figure 11 presents the yield functions and Figure 12 the plastic zones for wellbores with same direction (α = 30°), but different inclination, as reported in Table 2. The inclined meshes are plotted with the vertical mesh, and their yield functions and plastic zone are compared. It can be seen that the plastic area is gradually decreasing as the wellbore increase its inclination. In all cases, the same wellbore mud pressure is assumed.
Fig. 11 Yield functions of I1, I2, I3, and V1 wellbores. 
Fig. 12 Plastic zones of I1, I2, I3, and V1 wellbores. 
The plastic zone of a heterogeneous media is compared to a homogeneous one for a vertical wellbore. Both models assume the same input parameters, except for Young’s Modulus. Note that a homogeneous field simulation – with a single value of Young’s Modulus for the entire domain – is a common practice by petroleum engineers. This average Young’s Modulus value is used as the mean value necessary to define the stochastic field. The standard deviation considered in this study is based on a Coefficient of Variation (CoV) of 10%. For a homogeneous and a heterogeneous field, Figure 13 depicts Young’s Modulus distribution in the field, the plastic zone at the borehole wall, and the plastic region areas. The total area of the plastic zone for the homogeneous field is 13.901 × 10^{−4} m^{2} and 12.631 × 10^{−4} m^{2} for this heterogeneous field example case. It can be noted that the plastic region around the wellbore for a homogeneous field is symmetric. For the heterogeneous case, however, the plastic zone surrounding the wellbore is not symmetric, presenting 56% of the total area on the top side and 44% of the total plastic zone area on the downside of the third column of Figure 13.
Fig. 13 Plastic zone area of the homogeneous and a random heterogeneous case: vertical wellbore. 
The scale factor used in this work for all heterogeneous case analysis is γ = 0.5, which is the ratio of the first element at the wellbore wall. The asymmetry commented above may lead to operational issues during a drilling operation. Once mud properties are designed to stabilize the most vulnerable side of the well, a drilling fluid density considered safe if drilling a homogeneous formation may present a higher and asymmetric plastic zone when a heterogeneous analysis is performed.
4.2 Stochastic analysis: inclined wellbores
The analyzed inclined wellbores present a decrease in the plastic zone area as the inclination increases when considering homogeneous media. Thus, an analysis is made in order to evaluate the heterogeneity effects in wellbores’ inclination and direction. The analysis assumes the same input parameters as Table 1, but for this inclined study, the threedimensional stochastic field generation method is used. Therefore, a cylinder is built for the stochastic analysis. For the present case, the cylinder has eight equal crosssections, and each section includes a twodimensional and horizontal mesh, while the inclined mesh is embedded within the cylinder.
In order to reduce computational processing time, the mesh is changed to 2 m of the external radius. The number of elements and the wellbore radius are kept the same as in all analyses. The space between each mesh in the cylinder is 0.5 m, and the cylinder’s total height is 4 m. Ten thousand cases for each wellbore inclination is presented in Figure 14. The vertical wellbore considered in the analyses have the same direction as considered in the inclined wellbores, but no inclination angle (wellbore named V2). The direction was assumed in the stochastic analysis, in order to evaluate the influence of the direction in the yield function. One can note that as well as in the homogeneous domain, the heterogeneous field cases also reduce its total plastic zone area (A_{T}) as inclination increases. Figure 16 compares the total plastic zone area histogram for the same vertical well with and without direction angle. The histograms depict that the azimuth has little effect in the vertical distribution dispersion. The Complementary Cumulative Distribution Function (CCDF) of each histogram are depict in Figures 15 and 17, which allows the user to verify the failure probability of A_{T} > x, 0 < x < 20 × 10^{−4} m^{2}.
Fig. 14 Histograms of A_{T} of wellbores I1, I2, I3, and V2. 
Fig. 15 CCDF of wellbores I1, I2, I3, and V2. 
Fig. 16 Histograms of A_{T} of wellbores V1 and V2. 
Fig. 17 CCDF of wellbores V1 and V2. 
Table 3 presents the total plastic area mean () of each stochastic case, as well as their CoV and the deterministic total area for the homogeneous field (). The highest CoV is found in the most inclined wellbores. One, however, may associate the dispersion of the variable under analysis to the threedimensional approach for stochastic inclined fields; thus, the same study presented in Figure 14 is evaluated disregarding the 3D field model. The result gives no significant change in and CoV values, with differences lower than 2%. Therefore, the small influence on the plastic zone area under the analyzed mud pressure is likely due to the very small plastic zone area under analysis. Also, one of the covariance function main feature is the distance between elements, and at the borehole wall, the length of the elements is significantly small due to the geometric progression. Thus, around the borehole wall, Young’s Modulus variability is not as high as it is in the farfield, even if the mesh has an inclination, once the plastic zone area is still very close to the borehole wall with very small elements.
Mean plastic area () and CoV of stochastic cases, and homogeneous plastic area ().
4.3 Probabilistic analysis: maximum plastic deformation area
As reported by Muller et al. [15], probabilistic analyses allow the user to verify the failure probability of wellbores based on a maximum value of the plastic zone area, according to the simulated drilling fluid’s density. If the maximum allowed area (A_{y}) of the plastic region at the borehole wall is a percentage of the wellbore cavity area (A_{w}), the failure probability of reaching the maximum failure area could be verified by:
where, A_{y} = p% of A_{w}, and A_{T} is the total area of plastic region.
In order to analyze the uncertainty of the drilling fluid pressure, a probabilistic analysis is carried out for a pressure range from 10 MPa to 35 MPa for the vertical wellbore studied in this work, and considering the same failure criteria. Ten thousand simulations are carried out for each mud weight. Additionally, random normal distributions are assigned for the insitu stresses. The distribution mean values are the respective stresses from Table 1, with a CoV of 10%. The random distributions follow the criteria:
Figure 18 depicts the plastic area for each simulation of ten thousand simulations for each mud pressure, while Figure 19 represents the failure probability for a given mud pressure, failure herein is given in a performance base manner based on account of plastic area being higher than 10%, 5%, and 3% of A_{w}. Considering 3% of A_{w} as an expected performance for the considering wellbore, no collapse would occur if using mud drilling’s pressure over 31 MPa since failure probability is close to zero. On the other hand, if 10% of A_{w} is the expected performance of the wellbore, the drilling mud weight shall be greater than 24 MPa to avoid shear failure. Note that the mud weight’s upper limit is defined by rock tensile strength, which is not evaluated in this work.
Fig. 18 Monte Carlo analysis. 
Fig. 19 Failure probability of maximum allowed area. 
5 Conclusion
A stochastic geomechanical model is presented to predict formation behavior in drilling operations considering heterogeneous fields. The simulation method proposed herein predicts the total area of the plastic region surrounding the wellbore once it is submitted to a mud/drilling fluid’s pressure. The results show that, when considering a heterogeneous field, a nonsymmetric plastic region may occur. This lack of symmetry may lead to operational difficulties.
A new procedure is proposed to evaluate stochastic fields of the inclined wellbore; this procedure allows the modeling of thick formation layers based on twodimensional geomechanical models with a low computational cost if compared to full 3D random field analyses.
A wellbore case is evaluated, and its results considering homogeneous and heterogeneous field are compared. The studied wellbore, when subject to a specific mud weight, exhibits different total plastic zone area according to its inclination. It is observed that the entire plastic area becomes smaller as the wellbore inclination increases. Therefore, considering the wellbore features analyzed herein, the mud pressure evaluated is safer for a wellbore with greater inclination, because the total plastic zone area is more significant for the vertical well compared to the inclined wellbore. This behavior can be seen in the homogeneous and heterogeneous field. Wellbore azimuth depicted a small influence in the distribution dispersion for vertical wells.
The total plastic zone area distributions determined in the stochastic analysis lead to an evaluation of the probabilistic failure. Moreover, it leads to information regarding the safety of a drilling operation, once an acceptable plastic region area at the borehole wall is defined by the group responsible for the operation. Note that a mud pressure considered safe in a homogeneous field analysis, may not be indicated based on the probabilistic study framework detailed herein.
The failure probability for a mud weight range can help in decisionmaking during drilling operations in heterogeneous formations. Besides, one is able to define a target wellbore performance based on plastic deformation area and verify the failure probability when using a specific drilling fluid’s density.
Acknowledgments
The authors are indebted to National Agency of Petroleum, Natural Gas and Biofuels (ANP) and to Brazilian Petroleum Corporation – Petrobras, Grant: 2014/000902, for funding this research project. Luiz C.M. Vieira Jr. is indebted to national council for scientific and technological development (CNPq) for the research productivity fellowship grant 304005/20177. Philippe R. B. Devloo is indebted to CNPq for the research productivity fellowship grant 305823/20175, and to São Paulo Research Foundation (FAPESP) for the research grant 2017/157363. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors only and do not necessarily reflect the views of the sponsors or affiliates. The authors would like to thank Dr. Thiago Dias dos Santos and Mr. Patrick Kherlakian for helpful discussions and contributions.
References
 Gholami R., Moradzadeh A., Rasouli V., Hanachi J. (2014) Practical application of failure criteria in determining safe mud weight windows in drilling operations, J. Rock Mech. Geotech. Eng. 6, 1, 13–25. [CrossRef] [Google Scholar]
 Ostadhassan M., Jabbari H., Zamiran S., Osouli A., Oster B., Lentz N. (2014) Wellbore instability of inclined wells in highly layered rocks – Bakken case study, Society of Petroleum Engineers. 2123. SPE171026. [Google Scholar]
 Zamiran S., Osouli A., Ostadhassan M. (2014) Geomechanical modeling of inclined wellbore in anisotropic shale layers of Bakken Formation, American Rock Mechanics Association. 21–23. ARMA20147422. [Google Scholar]
 Aadnov B.S., Rogaland U., Chenevert M.E. (1987) Stability of highly inclined boreholes, Society of Petroleum Engineers. SPE16052. [Google Scholar]
 Zervos A., Papanastasiou P., Cook J. (1998) Elastoplastic finite element analysis of inclined wellbores, Society of Petroleum Engineers. SPE47322. [Google Scholar]
 Garrouch A.A., Ebrahim A.S. (2001) Assessment of the stability of inclined wells, Society of Petroleum Engineers. SPE68861. [Google Scholar]
 Aslannezhad M., Manshad A.K., Jalalifar H. (2015) Determination of a safe mud window and analysis of wellbore stability to minimize drilling challenges and nonproductive time, J. Pet. Explor. Prod. Technol. 6, 3, 493–503. [Google Scholar]
 Manshad A.K., Jalalifar H., Aslannejad M. (2014) Analysis of vertical, horizontal and deviated wellbores stability by analytical and numerical methods, J. Pet. Explor. Prod. Technol. 4, 3, 359–369. [Google Scholar]
 Ma T., Chen P. (2015) Wellbore stability analysis of inclined wells in the AY field, Electron. J. Geotech. Eng. 20, 1313–1329. [Google Scholar]
 Cui L., Cheng A.H.D., Leshchinsky D., Abousleiman Y., Roegiers J.C. (1995) Stability analysis of an inclined borehole in an isotropic poroelastic medium, Rock Mechanics: Proceedings of the 35th US Symposium on Rock Mechanics, 5–7 June, Reno, Nevada. ARMA950307. [Google Scholar]
 Moos D., Peska P., Finkbeiner T., Zoback M. (2003) Comprehensive wellbore stability analysis utilizing Quantitative Risk Assessment, J. Pet. Sci. Eng. 38, 3, 97–109. [Google Scholar]
 Udegbunam J.E., Aadnoy B.S., Fjelde K.K. (2014) Uncertainty evaluation of wellbore stability model predictions, J. Pet. Sci. Eng. 124, 254–263. [Google Scholar]
 Tran M.H., Abousleiman Y.N. (2010) The impacts of failure criteria and geological stress states on the sensitivity of parameters in wellbore stability analysis, American Rock Mechanics Association. ARMA10328. [Google Scholar]
 Fontoura S.A.B., Holzberg B.B., Teixeira É.C., Frydman M. (2002) Probabilistic analysis of wellbore stability during drilling, Society of Petroleum Engineers. SPE78179. [Google Scholar]
 Muller A.L., Vargas E.A., Vaz L.E., Gonçalves C.J. (2009) Borehole stability analysis considering spatial variability and poroelastoplasticity, Int. J. Rock Mech. Min. Sci. 46, 1, 90–96. [CrossRef] [Google Scholar]
 Muller A.L., Vargas E.A., Vaz L.E., Gonçalves C.J. (2009) Threedimensional analysis of boreholes considering spatial variability of properties and poroelastoplasticity, J. Pet. Sci. Eng. 68, 3, 268–276. [Google Scholar]
 Ottesen S., Zheng R.H., McCann R.C. (1999) Borehole stability assessment using quantitative risk analysis, Society of Petroleum Engineers. SPE52864. [Google Scholar]
 Mostafavi V., Aadnoy B.S., Hareland G. (2011) Model – Based uncertainty assessment of wellbore stability analyses and downhole pressure estimations, American Rock Mechanics Association. ARMA11127. [Google Scholar]
 Wang D., Ge H., Yu B., Wen D., Zhou J., Han D., Liu L. (2018) Study on the influence of elastic modulus heterogeneity on insitu stress and its damage to gas shale reservoirs, J. Nat. Gas Geosci. 3, 3, 157–169. [CrossRef] [Google Scholar]
 Vallejo L.I.G., Hijazo T. (2008) A new method of estimating the ratio between in situ rock stresses and tectonics based on empirical and probabilistic analyses, Eng. Geol. 101, 3–4, 185–194. [Google Scholar]
 Roegiers J.C. (2002) Well modeling: An overview, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 57, 5, 569–577. [CrossRef] [Google Scholar]
 Charlez P.A., Roatesi S. (1999) A fully analytical solution of the wellbore stability problem under undrained conditions using a linearised camclay model, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 54, 5, 551–563. [CrossRef] [Google Scholar]
 Rocha L.A.S., Azevedo C.T. (2007) Projetos de Poços de Petróleo, 2nd edn., Interciência. [Google Scholar]
 Devloo P.R.B. (1997) PZ: An object environment for scientific programming, Comput. Methods Appl. Mech. Eng. 150, 133–153. [Google Scholar]
 DiMaggio F.L., Sandler I.S. (1971) Material model for granular soils, J. Eng. Mech. 97, 3, 935–950. [Google Scholar]
 Ceclio D.L. (2014) Elastoplastic modeling applied to the numerical simulation of wellbore stability analysis, PhD Thesis, University of Campinas, Brazil. [Google Scholar]
 Santos E.S.R. (2009) Numerical simulator of saturated elastoplastic porous media, PhD Thesis, University of Campinas, Brazil. [Google Scholar]
 Ceclio D.L., Devloo P.R.B., Gomes S.M., Santos E.R.S., Shauer N. (2015) An improved numerical integration algorithm for elastoplastic constitutive equations, Comput. Geotech. 64, 1–9. [Google Scholar]
 Fjaer E., Holt R.M., Raaen A.M., Risnes R., Horsrud P. (2008) Petroleum related rock mechanics, Elsevier. [Google Scholar]
 Batalha N.A. (2017) Linear elastic model for the simulation of the stress state in inclined wellbores, Master’s Thesis, University of Campinas, Brazil. [Google Scholar]
 DeGroot D.J., Baecher G.B. (1993) Estimating autocovariance of insitu soil properties, J. Geotech. Eng. 119, 1, 147–166. [CrossRef] [Google Scholar]
 Torquato S. (2013) Random heterogeneous materials: Microstructure and macroscopic properties, Springer, New York. [Google Scholar]
 Cranston P.G., Richie M.C., Vieira Jr L.C.M. (2017) Effects of soil support on the stability of corrugated metal pipe, Structural Stability Research Council. [Google Scholar]
 Rocha L.A.S., Aruaga D., Andrade R., Vieira J.L.B., Santos O.L.A. (2008) Perfuração Direcional, Interciência. [Google Scholar]
All Tables
Mean plastic area () and CoV of stochastic cases, and homogeneous plastic area ().
All Figures
Fig. 1 Twodimensional mesh surrounding the wellbore. 

In the text 
Fig. 2 Young’s Modulus (MPa) of a random correlated field. 

In the text 
Fig. 3 Young’s Modulus (MPa) of a random but not correlated fied. 

In the text 
Fig. 4 Vertical wellbore mesh. 

In the text 
Fig. 5 Inclined wellbore mesh. 

In the text 
Fig. 6 Five steps of the 3D spatially correlated field methodology. 

In the text 
Fig. 7 Yield function of wellbores V1. 

In the text 
Fig. 8 Plastic zone of wellbores V1. 

In the text 
Fig. 9 Yield function of wellbores I1. 

In the text 
Fig. 10 Plastic zone of wellbores I1. 

In the text 
Fig. 11 Yield functions of I1, I2, I3, and V1 wellbores. 

In the text 
Fig. 12 Plastic zones of I1, I2, I3, and V1 wellbores. 

In the text 
Fig. 13 Plastic zone area of the homogeneous and a random heterogeneous case: vertical wellbore. 

In the text 
Fig. 14 Histograms of A_{T} of wellbores I1, I2, I3, and V2. 

In the text 
Fig. 15 CCDF of wellbores I1, I2, I3, and V2. 

In the text 
Fig. 16 Histograms of A_{T} of wellbores V1 and V2. 

In the text 
Fig. 17 CCDF of wellbores V1 and V2. 

In the text 
Fig. 18 Monte Carlo analysis. 

In the text 
Fig. 19 Failure probability of maximum allowed area. 

In the text 