Wellbore stability analysis during drilling is of special importance in the gas and oil industry. Loss of wellbore stability during drilling often results in increase of drilling cost and drop of drilling efficiency, even collapse of wellbore. Regardless of the drilling techniques (e.g., overbalanced or underbalanced) employed, if a wellbore is kept stable throughout the drilling process, integrity issues of the wellbore are unlikely to occur during oil and gas production^[1,2]. Planning and executing drilling operations with wellbore stability in mind is the focus of many studies and a priority for most drilling projects, therefore, understanding the factors affecting borehole stability during drilling is paramount^[3,4,5,6,7].

During drilling, as the rock on the wellbore track is drilled and brought out of the hole, the drilling fluid would exert a corresponding pressure on the wellbore wall, the stress of the rock surrounding the wellbore would be redistributed^[8], generating induced stress. To maintain wellbore stability, it is essential to drill with a drilling fluid of appropriate density (mud weight) to control induced wellbore stress^[9]. Shear and tensile failures are the major causes of mechanical instability in boreholes^[10]. If above the appropriate upper limit, drilling-fluid pressure would cause tensile failure in the wellbore wall; if lower than the appropriate lower limit, drilling fluid pressure would cause shear failure in the wellbore wall^[11,12,13]. The upper limit of safe mud weight is typically defined by mud loss to the formation and tensile failure of the wellbore wall, whereas the lower limit is typically defined by shear failure and sloughing^[14].

A number of approaches have been proposed to investigate wellbore stability. McLellan^[15] divided the prediction methods of wellbore instability into three categories: empirical, deterministic and probabilistic ones, of them, probabilistic methods were subdivided further into numerical, analytical and experimental analysis ones. Lee et al.^[16] assessed wellbore stability with the optimized mesh finite-element method. Wang and Sterling^[17] analyzed the stability of a horizontal well drilled into a loose sandstone formation, especially the effect of mud cake on wellbore stability, with finite-element method too. Salehi et al.^[18] evaluated wellbore stability by using finite-element and finite-difference methods and validated their results. Numerical models are widely used to research wellbore instability^{[19,20,21,22]}. Several numerical studies derived analytical solutions to express wellbore stability. For example, Manshad et al.^[23] adopted the three-dimensional analytical model proposed by Al-Ajmi and Zimmerman^[24] to work out appropriate mud weight, wellbore inclination and azimuth of a wellbore drilled in an oil field in southern Iran.

In this work, we employed FLAC3D (Fast Lagrangian Analysis of Continua in 3 Dimensions^[25]) Software to determine the safe mud weight window in a vertical wellbore penetrating the sandstone layers of the Asmari reservoir formation at one of the oilfields in south-western Iran. FLAC3D characterizes complex behaviors and relationships by continuously displaying limited volume, and is suitable to solve questions related to system failure and/or collapse of subsurface structures caused by continuous non-linear displacement and strain. In this paper, a vertical well encountering Asmari sandstone in an oilfield in southwest Iran was taken as an example, a finite volume model of the wellbore was built based on the geo-mechanical features of the formations drilled with the FLAC3D software, the upper and lower limits of safe mud weight for the Asmari sandstone were defined through simulating the formation of plastic deformation zone in the wellbore, and the effects of rock strength, main stresses around the wellbore and pore pressure on the safe mud window were evaluated.

Most issues during drilling are related with drilling fluid properties, especially with drilling fluid density directly or indirectly (Fig. 1). To analyze wellbore stability, well logs, laboratory test results and drilling reports of the wellbore were used to find out the geomechanical characteristics of the formations drilled. The elastic properties, rock strength characteristics, in situ stresses, and pore pressures are the most important geomechanical characteristics to be measured or calculated.

The strength and prevailing in situ stresses of a rock formation can be worked out from the rock elastic property^[26]. Poisson’s ratio and Young’s modulus are key parameters of rock elastic properties. Assuming elasticity is isotropic, the dynamic Young’s modulus (E_dyn)can be calculated by Eq. 1 with sonic well log data, and the dynamic Poisson’s ratio (υ_dyn) can be calculated by Eq. 2:

Testing core samples extracted from wellbores in the laboratory under reservoir pressure conditions is the typical method used to determine these two rock elastic parameters. However, only limited core samples are taken from commercially-drilled wells^[27]. Therefore, dynamic Poisson’s ratio and Young’s modulus are usually calculated by Eqs. 1 and 2 with sonic well-log data, and supplemented by static values obtained from few core samples taken from surrounding wellbores. Static Young’s modulus (E_sta) and static Poisson’s ratio (υ_sta) of the Asmari formation sandstone reservoir measured on cores have correlations with their dynamic counterparts calculated with well-log data, which can be expressed as Eqs. 3 and 4.

Bulk modulus (K) and shear modulus (G) can then be calculated from Poisson’s ratio and Young’s modulus by Eqs. 5 and 6.

A number of relationships involving rock strength have been proposed to quantify geomechanical characteristics of wellbore when there is no test data from core samples available. Generally, these relationships involve metrics that have a direct influence on rock strength, such as the elastic modulus and porosity^[28]. Of them, the relation between uniaxial compressive strength (UCS) and elastic modulus (E) can be expressed as Eq. 7.

In modeling with FLAC3D, friction angle (φ) (based on the Mohr-Coulomb criterion) and cohesion (C) are often used to characterize compressive strength of rock formation^[29,30]. Of them, the friction angle (φ) is worked out by Eq. 8^[31].

where

${{V}_{shale}}=\frac{GR-G{{R}_{min}}}{G{{R}_{max}}-G{{R}_{min}}}~~$

Cohesion (C) can be calculated from UCS and internal friction angle by Eq. 9^{[32, 33]}.

In order to evaluate tensile failure in a wellbore wall related to stress concentration, it is necessary to establish a criterion for the tensile strength. Generally, rock tensile strength is $\frac{1}{12}$ to $\frac{1}{8}$ of the UCS^[34]. Correlation analysis results show the best tensile strength of reservoir rock is about 1/10 of its UCS^[35,36]. Hence, the tensile strength was taken as $\frac{1}{10}$ of the UCS in this study.

The most definitive and accurate pore pressure is measured during drill stem tests (DST). Repeat formation tests (RFT) and modular formation dynamics tests (MDT) can also measure formation pressure quickly and accurately. These tools can measure pore pressure at specific points on the wellbore wall. RFT and MDT tools are operated on wireline and can provide accurate pore pressure data, particularly when many tests are performed repeatedly in close distance in specific formations. By plotting the relationship curve of pore pressure and depth of many zones in many wells, the formation pore pressure of a specific formation in a specific field can be obtained from extrapolation. In this study, the measured pore pressure data was obtained from numerous RFT tests.

One vertical and two horizontal stresses are typically used to represent the three principal stresses in rock formations, among them, two representing the maximum and minimum stresses. The vertical stress indicates the weight the overburden rock exerts on the formation. The overburden stress increases with depth as the weight of overburden rock increases. If the overburden rock layers are homogenous in the vertical depth direction (z), then the vertical stress (σ_v) can be simply calculated by Eq. 10.

But if the density of overburden layers is not uniform and varies with depth, then the vertical stress at a specific depth D is:

The estimation of horizontal stresses is key to accurate modeling of the in-situ stress regime^[37]. In isotropic conditions, the horizontal stresses in all directions are approximately equal. This case usually occur in rock formations not affected by large-scale regional tectonic movements giving rise to regional folds and faults within the Earth’s crust^[38]. However, in large geologic structures such as anticlines, salt domes or other tectonically-active zones, the two principal horizontal stresses are different and maximum horizontal and minimum horizontal stress directions can be identified^[39]. In modeling, the above two cases all should be considered.

A poro-elastic horizontal strain model (PHSM) can be employed to calculate the maximum and minimum horizontal stresses^[40,41,42].

The Biot coefficient is calculated from bulk modulus^[43,44]. For the relatively homogeneous Asmari sandstone formation, it is reasonable to assume a Biot coefficient equal to 1. However, in heterogeneous formations such an assumption isn’t reasonable^[45].

The geomechanical characteristics of a Asmari sandstone layer at 4 090 to 4 095 m were worked out by Eqs. 1 to 14 and the PHSM, then the wellbore stability model was built with FLAC3D software^[25]. The basic data used to calculate the geomechanical characteristics are: volume modulus of 20.7 GPa, shear modulus of 13.0 GPa, internal cohesion of 34.7 MPa, internal friction angle of 35.7°, maximum horizontal stress of 91.6 MPa, minimum horizontal stress of 57.6 MPa, ration of maximum horizontal stress to minimum horizontal stress of 1.59, overburden rock pressure of 102.5 MPa, pore pressure of 40.0 MPa, wellbore radius of 10.7 cm, and permeability of 110×10^-3 μm².

The artificial boundaries of a volume for modeling an underground structure must be selected carefully. The ratio between length and width of the modeled volume should be close to 1. Also, the aspect ratio (vertical to horizontal) of the elements of the finite volume model (FLAC3D) should be less than five in order to avoid degradation of the model’s numerical simulation performance^[25]. As the near wellbore zone is more important, the grids of the finite volume model around the wellbore are finer than those further away from the wellbore. In this study, the grids in the near-wellbore zone were set at 1 cm by 1 cm on the vertical plane. As the distance from the wellbore increases, the grid size became larger. The artificial boundaries for the modeled volume were set at 2 m by 2 m in the horizontal plane, and the wellbore was in the center of the model, with a vertical height of 5 m.

By combining the defined well-bore geometry (Fig. 2) and rock mechanics characteristics of formation above (Table 1), the wellbore stability in the Asmari sandstone can be modeled in detail. This is realized by monitoring the initiation of the plastic condition around the wellbore wall of the primary grids 0.1 times of the wellbore radius. The lower and upper limits of the safe mud weight window (SMWW) for the Asmari sandstone obtained by the FLAC3D modeling were 1 634 kg/m³ and 1 009 kg/m³, respectively. In other words, plastic zone would be created in the wellbore wall when the drilling mud weight is less than 1 009 kg/m³. When the mud weight is no less than 1 009 kg/m³, the plastic zone would disappear and does not reappear until drilling mud weight exceeds 1 634 kg/m³ (Fig. 3). Fig. 4 shows the stress on the wellbore wall according to the finite volume model for drilling mud weights outside the safe mud weight window range.

The criterion used to define the SMWW for the Asmari sandstone in the studied well is in good agreement with information recorded in the drilling reports and observed fracture development in the well. The determined SMWW (1 009-1 634 kg/m³) is therefore reasonable. Through sensitivity analysis, the effects of geomechanical properties on safe mud weight for this formation were evaluated further. Sensitivity of rock strength properties, in-situ stresses around wellbore, and pore pressure in the determined SMWW were analyzed with the FLAC3D model. By adjusting the basic data used to calculate the geomechanical parameters of Asmari sandstone, the safe mud weight under each condition was estimated by the finite volume model.

As mentioned above, cohesion and internal friction angle are two important factors influencing rock strength. Cohesion is the interaction force between rock particles (mineral grains and rock fragments). It is the shear strength when the normal stress is zero. This is a static force that only exists just before rock failure occurs. The internal friction angle is related to the roughness of mutual sliding surfaces within the rock fabric, and is a dynamic force and that exists during the failure process^[46].

To study the effect of cohesion (C) on SMWW, besides the base case value of 34.7 MPa, a range of cohesion values from 20 to 50 MPa were evaluated. The SMWWs for this range of cohesion values based on the criterion of plastic condition initiation in the finite volume model are shown graphically in Fig. 5. The sensitivity analysis results show that when the cohesion decreases, the SMWW turns narrower (i.e., the difference between the maximum and minimum limits of safe mud weight turns smaller). This implies that if the rock constituting the wellbore wall has smaller cohesion, the risk of well bore losing stability is higher, and the mud weight window preventing wellbore instability is narrower.

By using the same sensitivity analysis approach, the impact of the internal friction angle (φ) on the SMWW was studied. The friction angle was 35.7 degrees in the base case, and varied from 20 to 50 degrees in the sensitivity analysis. The maximum and minimum safe mud weights for all φ values are shown in Fig. 6. Similar with cohesion, when the friction angle turns smaller, the SMWW turns narrower, and the risk of wellbore instability becomes higher.

The prevailing stress field impacting a formation prior to penetration is referred to as the in-situ stress. In-situ stress is typically divided into three components: vertical stress (overburden), maximum horizontal stress and minimum horizontal stress.

Once a well is drilled into the formation, the wellbore wall has to withstand the stress field previously supported by the rock but now removed in the form of drill cuttings. The stresses exerted on the wellbore wall depend upon the wellbore’s orientation and the prevailing in-situ stress field impacting the formation at the drilling location^[47]. Before drilling, the formation is typically in a stable stress equilibrium. However, once penetrated by the wellbore, the stress distribution around the wellbore wall will be disturbed and this can potentially damage the wellbore stability. The disturbed stress field at or close to the wellbore wall is expressed in three different stresses (Fig. 7): 1) the radial stress exerting outwards from the center of the wellbore; 2) hoop stress exerting in tangential direction around the circumference of the wellbore wall; and, 3) axial stress exerting along the length of the wellbore’s trajectory (i.e., vertical stress in a vertical well). These three wellbore stresses are perpendicularly with each other, making it easy to characterize them in a wellbore coordinate system.

The long-established Kirsch equations, based on elasticity theory, can provide closed-form solutions to determine the stresses around a circular, excavated area, such as a wellbore. Eqs. 14 to 16 can be used to derive the three stresses on the wellbore wall of a vertical well^[48].

For a vertical wellbore, the maximum tangential stress (σ_t,max) occurs in the direction of the minimum far-field horizontal stress (i.e., θ=90°); whereas the minimum tangential stress (σ_t,min)is in the direction of the maximum far-field horizontal stress (i.e., θ=0°, see Fig. 8)^[49]. These two horizontal stress components can be calculated by Eqs. 17 and 18.

Fig. 9 shows that a mud weight too low will cause shear failure in the direction of the minimum horizontal stress. On the contrary, a mud weight too high will cause tensile failure in the direction of maximum horizontal stress^[50].

As the object of this study was a vertical well, the in-situ stresses around the wellbore (i.e., the maximum and minimum horizontal stresses) have a greater impact on wellbore stability than the vertical stress. For this reason, it is necessary to assess the effect of horizontal stress ratio on the SMWW. By changes the value of minimum horizontal stress and then the stress ratio and calculating the SMWW for each case, the sensitivity of the SMWW was evaluated.

It can be seen from Fig. 10 that as the ratio of maximum to minimum horizontal stresses (k) decreases (i.e., as the maximum horizontal stress value moves closer towards the minimum horizontal stress value), the wellbore becomes stable over a much wider range of mud weights (i.e. the SWMM widens considerably). Moreover, as k decreases, the maximum and minimum values of the safe mud weight window change in quite different gradients. This difference can be explained by that the impact of different k values on the maximum and minimum tangential stress calculated by Eqs. 17 and 18.

The pore pressure is the pressure of formation fluid within the pores of a reservoir formation. The pore pressure affects the stability of a wellbore by creating an effective stress that influences the wellbore stress field and therefore the SMWW. To assess the effect of pore pressure on SMWW for the Asmari sandstone in the wellbore studied, sensitivity analysis was conducted with pore pressure varying from 30 to 45 MPa, and the SWMM calculated for each case. Fig. 11 shows that as pore pressure goes down, the SMWW widens, i.e., when pore pressure decreases, the wellbore is stable under a wider range of mud weights.

It can be seen from Figs. 10 and 11 that as k and pore pressure increase, the SMWW narrows and the risk of wellbore instability goes up.

The results of this study demonstrate that a finite-volume simulation model of wellbore instability developed by FLAC3D software with calculated geomechanical characteristics of the Asmari sandstone in an Iranian oilfield, is effective in defining the limits of the safe mud weight window (SMWW). By using the initiation of plastic state as the criterion to determine wellbore stability, the maximum and minimum safe mud weights for Asmari formation in a vertical well were obtained from the simulation, which are consistent with those from drilling observations. The developed model can be used in sensitivity evaluation, to find out the effects of rock strength characteristics, in situ stresses, and pore pressure on the limits of the SMWW. The results of the sensitivity analysis show that the reduction of cohesion and internal friction angle would lead to a significant narrowing of the SMWW and the rise of wellbore instability risk. Whereas, the reduction of pore pressure and the ratio of maximum horizontal stress to minimum horizontal stress (K) would make the SMWW widen significantly and the risk of wellbore instability reduce. The model, which can quantify the variation of safe mud weight window, can be taken as a drilling design and monitoring tool. The follow-up study will verify the results of this study by experiments, tri-axial stress experiments will be done on reservoir samples to find more information on formation collapse and fracturing pressure.

C—cohesion, MPa;

D—depth of a specific point, m;

E—elastic modulus, GPa;

E_dyn—dynamic elastic modulus, GPa;

E_sta—static elastic modulus, GPa;

G—shear modulus, GPa;

GR—Gamma log, API;

GR_max—maximum value of gamma log, API;

GR_min—minimum value of gamma log, API;

g—gravity acceleration, m/s²;

K—volume modulus, GPa;

k—ratio of maximum horizontal stress to minimum horizontal stress;

p_m—hydrostatic pressure of drilling fluid, MPa;

p_p—pore pressure, MPa;

UCS—uni-axial compressive strength, MPa;

V_shale—shale content, %;

z—formation depth, m;

α—Biot coefficient;

ε_x—strain in the direction of maximum horizontal stress, dimensionless;

ε_y—strain in the direction of minimum horizontal stress, dimensionless;

r, θ—polar coordinates, (°);

υ—Poisson’s ratio, dimensionless;

υ_dyn—dynamic Poisson’s ratio, dimensionless;

υ_sta—static Poisson’s ratio, dimensionless;

ρ—rock density, g/cm³;

ρ_b—bulk density, g/cm³;

σ_a—axial stress, MPa;

σ_H—maximum horizontal stress, MPa;

σ_h—minimum horizontal stress, MPa;

σ_r—radial stress, MPa;

σ_t—tangential stress, MPa;

σ_{t, max}—maximum tangential stress, MPa;

σ_{t, min}—minimum tangential stress, MPa;

σ_v—vertical stress, MPa;

φ—internal friction angle, (°);

ϕ_CNL—neutron porosity (rock porosity from neutron log), %;

Δt_c—interval transit time, μs/m;

Δt_s—shear wave time difference, μs/m.