Integrated mechanical earth model and quantitative risk assessment to successful drilling

Use of vital geomechanical parameters for determination of safe mud pressure window is generally associated with high level of uncertainty primarily because of absence of sufficient calibration data including laboratory and field test information. The traditional deterministic wellbore stability analysis methodologies usually overlooked the uncertainty of these key parameters. This paper exhibits implementing a quantitative risk assessment technique on the basis of Monte-Carlo modeling to consider uncertainty from input data so as to make it possible to survey not just the likelihood of accomplishing a desired level of wellbore stability at a particular mud weight, but also the impacts of the uncertainty in each single parameter on the wellbore stability. This methodology was implemented to a case study. The most important parameters have been recognized using a sensitivity analysis approach in which the outcome of this QRA procedure suggests the mud weight window with likelihood of well drilling success which can elude the wellbore collapse and lost circulation events. This sort of stochastic approach to deal with anticipated safe mud weight window can guarantee stable wellbore with considerable cost viability associated with drilling success. The technique built up in this paper can give the scientific foundation for assessment of wellbore stability under complicated geological circumstances. It was likewise noted that based on sensitivity analysis, uniaxial compressive strength and maximum horizontal stress are the most effective parameter in estimation of mud weight window. This accentuates the significance of trustworthy determinations of these two parameters for safe drilling of the future wells in the field.


Introduction
Proper mud weight plays a vital role in wellbore instability studies in oil and gas industry. As most of the current drilling projects including well construction and field development are challenging, safe mud weight windows have a deteriorate effect on the time schedule of the project and can drastically lessen non-productive time (NPT), issues related to instability of wellbores and unwanted costs on drilling operation (Plumb et al. 2000;Elyasi et al. 2013).
To minimize the risk of wellbore instability and following technical/economic consequences, a suitable pre-drilling geomechanical scheme needs to be developed. Building mechanical earth model (MEM) is a practical solution to guarantee well safety and reduce project costs especially during the well-planning phase. MEM is a numerical representation of the stress state and rock mechanical properties for a particular stratigraphic section in an oil field or a basin. Building a MEM requires integrating data from various sources to accurately describe the formations in terms of geomechanical attributes. Therefore, the more accurate data are available, and the more reliable model can be built (Noeth and Birchwood 2015; Afsari et al. 2009). However, in reality, not all geomechanical parameters are completely available, accurate or fully trustworthy due to lack of data sources and variability of information as a function of the depth. Furthermore, most of the inputs are derived from seismic data, wireline logs and core analysis results. The uncertainties arise in input dataset due to their own inherent limitations, limited amount of calibration with measured data and empirical equation for derivation of rock mechanical data. Accordingly, deterministic commonly utilized model generates results of calculating mud weights with considerable errors and wrong predictions; therefore, it is essential to quantify effects of input uncertainties on results of wellbore stability model via a stochastic approach. Ottesen et al. (1999), introduced a probabilistic approach based on quantitative risk assessment (QRA) for oil and gas drilling projects which further developed by Moos et al. (2003), to diagnose any uncertainty related to input data, identify wellbore collapse risk and finally minimize operational damages and consequent costs. In this method, a proper distribution function needs to be selected to evaluate and quantify errors. Then, an appropriate constitutive model must be identified to illustrate rock properties, deformation behavior and connect input parameters to purpose output. After identifying the constitutive model, the thresholds between success and failure should be determined through limit state function (LSF) (Gholami et al. 2015). This function makes a link between conventional stability models against operational failures via quantifying the risk of any operational failure and specifying proper mud densities range to lessen the likelihood of any instability related to wellbore projects. In the next step, several wellbore instability simulations need to be performed to build a response surface. The response surface is then used to determine the Likelihood of Success (LS) as a function of drilling fluid density to ensure reliable and safe drilling operation via quantifying uncertainties and errors associated with key parameters of wellbore, estimated for input and output data using probabilistic distribution functions. The last step is to perform an interactive numerical simulation method like finite element approach or Monte Carlo technique for each uncertainty and identify the thresholds to prevent wellbore collapse or lost circulation (Aadnoy and Looyeh 2010;Gholami et al. 2015). Monte-Carlo approach can be a confident alternative of traditional deterministic methods in oil and gas industry by sampling uncertainties and errors related to input datasets from the real distribution of measured parameters. Furthermore, it can be helpful to identify most critical parameters caused by utmost uncertainties and errors affecting the analysis and result to audit data optimally, maximize the likelihood of success and time-cost efficient reduction effort in any drilling or production operation.
In this study, QRA technique was undertaken utilizing log-derived rock elastic and strength parameters at a drilled well case to determine the disparities, errors and uncertainties associated with the obtained values, build limit state functions and response surfaces, recognize most influential parameters and their sensitivities, and finally evaluate how the likelihood of success can be maximized by carrying out the fundamental quantitative risk assessment technique. For this purpose, the reservoir in situ stress state and rock mechanical properties were derived by drilled well information and uncertainties associated with their calculations were incorporated. Once uncertainties for all the stresses and rock mechanical properties are defined, and random data samples were generated. Required mud weights from rock failure criteria were calculated for each of the input variable and stress regime and based on response surfaces, linear or polynomial equations were established to define required mud weight/pressure as a function of input variable. Then, Monte-Carlo simulations were run to generate cumulative probability curves for wellbore collapse and lost circulation. Based on these cumulative curves, probability of success of drilling is defined with optimum mud weight window.

Geomechanics data
The first step in the geomechanical study is data collection. Building a mechanical earth model (MEM) requires integrating data from various sources to precisely describe the formations in terms of geomechanical characteristics. In reality not every single essential data can be obtained, so we need to take steps to obtain the most information from accessible data. This geomechanical study is aimed to accomplish the objectives by maximizing the utilization of available data by selecting appropriate techniques for computations.
Available data that were used to determine the mud weight window (MWW) are tabulated in Table 1.
Wireline logs are the key source of geomechanical characterization in this research study. GEOLOG software was used to determine the petrophysical properties, in situ stress state, rock elastic and rock mechanical strength properties from the well logs. The conventional and advanced well logs used in this study to determine the MWW are shown in Fig. 1. These include caliper, gamma ray, density, porosity together with compressional and shear sonic slowness logs. The applied procedure to determine the well mud weight window is presented in Sect. 4.

Assessment of key physical parameters and their probability distribution functions
To set up the QRA for wellbore stability investigation, key input parameters are required to be estimated and quantified by fitting a distribution function to their variations. The key input parameters, ordinarily considered by quantitative risk assessment method for wellbore stability study, are rock strength properties (Young's Modulus and Poisson's ratio), in situ stresses (overburden stress, minimum and maximum horizontal stress), pore pressure, unconfined compressive strength and internal friction angle. An evaluation result of rock properties and stress determination is shown in Fig. 2. The estimation procedures of these input parameters for QRA are discussed in 3.1-3.3.
The stress model indicated that the stress regime in interval between 3395 and 3470 m is primarily reverse faulting ( H > h > V ) and the interval between 3470 and 3700 m have mainly a strike-slip ( H > v > h ) stress regime as it can be seen in Fig. 2. It is well proven that the necessary mud weight for wellbore stability is significantly affected by formation stress regime. Hence, we perform QRA method for both sections independently.
Quantification of uncertainties is important since log data, which are utilized as the sources of data for estimation of input parameters, are exposed to different errors because of common mechanical or electrical functioning in a defective way. Moreover, the assessed properties from log data should be calibrated against field or core data which are subjected to their own uncertainties during measurements. The effect of inconsistency in geological setting of the filed may similarly escalate the level of uncertainty when necessary calibration points are not sufficient.
Taking a gander at the distribution functions of key parameters, it is seen that none of them can foresee the whole variety of the log information. Seven distribution functions, including Log Logistic, Inverse Gaussian, Weibull, Beta, Cauchy, Logistic and Chi-Squared were utilized for describing the probable changes in the logs data. For each individual parameter, the best function was chosen and utilized for quantitative risk assessment. For every one of the fitted distributions, goodness of fit tests was conducted using Kolmogorov-Smirnov, Anderson-Darling and Chi-Squared tests (D'Agostino and Stephens 1986). At that point, the fitted distributions were ranked and the best one was picked for the subsequent investigations. An exemplary of ranking for each distribution function is presented in Table 2 which demonstrates that Beta is the best likelihood distribution function for Young's modulus data in this case. In this regards, the Anderson-Darling does a good job at fitting the tails of the distribution and Kolmogorov-Smirnov is good at fitting the mid-range of the distribution. According to the tests, both the tails and mid-range of distribution are fitted good by Beta. In addition, Chi-square test is a nonparametric test. It is based upon the assumption that the test statistics will tend to follow a Chi-square distribution. Therefore, it is less reliable as parametric functions. The best distribution function for each individual parameter is presented in Table 3. The selected distribution functions are utilized for quantitative risk assessment.
The Beta-General distribution specifies a beta distribution with a characterized least and most value utilizing shape parameters α1, α2 which are constantly greater than zero. Weibull function is described by a shape (λ) and a mean (μ) parameters. On the other hand, Log-logistic distribution function is described by two parameters. The parameter α is a scale parameter and is additionally the middle of the distribution. The parameter b > 0 is a shape parameter. These two parameters must have a positive value. Additionally, Cushy distribution function which is the best fitted for density log of interval 3395-3470 m has only one parameter: k, a positive integer that determined the quantity of degrees of freedom.
From the above discussion, it tends to be inferred that the sources of uncertainty would be very high in estimation of geomechanical parameters since the logs which are displayed in Fig. 1 ought to be used. This uncertainty will even be enlarged when the stated geomechanical parameters are utilized for determination of safe mud weight window.
Figures 3, 4 and 5 demonstrate the uncertainties of the input parameters which they are given by probability density functions that are determined by means of the minimum, the maximum and the most likely values of each parameter. These quantify the uncertainties in the input parameters required to compute the mud weight limits necessary to avoid wellbore instabilities.

Estimation of rock properties
Wireline log data were utilized to derive dynamic properties and then applied established correlation (in study area) to derive static properties.  Fig. 3 Histogram and probability distribution functions fitted to density, porosity, P and S waves transit time for RF (left) and SS (right) stress regimes

Elastic properties
Supposing elastic isotropy, sonic travel time assessment (i.e., compressional slowness ∆t C and shear slowness ∆t S ) and bulk density ρ b were utilized together to compute the dynamic elastic properties. Young's Modulus of rock layers was predicted utilizing dynamic elastic equation (Fjaer et al. 2008) and then converted to their static values utilizing correlation proposed by Ameen et al. (2009).
There are numerous uncertainties engaged with computations of static Young's modulus. As logs information utilized for estimation of dynamic Young's modulus suffer from uncertainties because of malfunctioning of tools and complexity of information recorded. Aside from this, the correlations utilized for dynamic to static conversion are field specific and in many cases extracted within adequate representative data. The lab data that are utilized to calibrate the log based extracted properties are normally constrained to few samples which may not be well representative of the whole log area. Furthermore, it is noteworthy that the lab tests themselves are exposed to uncertainties associated with sampling and preparation approaches, test system and  Histogram and probability distribution functions fitted and used for quantifying the variation pore pressure and in situ stresses for RF (left) and SS (right) stress regimes 1 3 equipment errors. These will cause intrinsic uncertainties with the predicted magnitude for the static Young's modulus.
Although according to the published studies, Normal, Triangle, Uniform and Lognormal density functions have been suggested as the best distribution functions for rock mechanical analysis (Moos et al. 2003), Weibull and Loglogistic probability distribution functions were found to be progressively proper for describing much of data variation corresponding to Young's Modulus. Figure 4 (top) demonstrates the distribution function or Young's modulus.
However, dissimilar to Young's modulus it is not necessary to convert dynamic Poisson's ratio to its equivalent static parameter (Wang 2001). Despite the fact that the uncertainty due to converting dynamic to static parameters is excluded in determination of Poisson's ratio, the risk of using uncertain log information will make the estimation of this parameter uncertain. The best distribution functions found to cover the variation of Poisson's ratio was Log-Logistic function as shown in Fig. 4.

Uniaxial compressive strength (UCS)
Numerous relations have been proposed by many authors that relate UCS to wireline log data. Most of these equations utilize compressional slowness or velocity, Young's Modulus and porosity. The obtained UCS can only be utilized in specific fields with similarities in geological settings.
To estimate the UCS of the formations intersecting Well X, correlation proposed by Militzer and Stoll (1973) was used. The related distribution function is shown in Fig. 4.

Internal friction angle
There have been relatively few endeavors to find relationships between the angle of internal friction angle (Φ) and geophysical measurements in light of the way that even week rocks have generally high Φ, and there are comparatively complicated relations between Φ and mechanical characteristics of rock such as rock stiffness . Therefore, internal friction angle is numbered among the crucial parameters in wellbore stability analysis with the highest level of uncertainty in its prediction.
In this study, to estimate the friction angle, the correlation proposed by Ameen et al. (2009) was applied and Beta and Log-logistic functions were used to capture most variation of Φ of interval 3395-3470 m and 3470-3700 m, respectively, as shown in Fig. 4 (bottom).

Computation of pore pressure
In this study, Eaton Method (1975) was used to obtain the pore pressure in the intervals of overburden and reservoir section those contain clay rich sedimentary rocks. The calculated pore pressure profile was calibrated using drilling mud weight and wellbore instability events such as mud loss and well kick (influx of reservoir fluid into the well). For reservoir section pore pressure gradient obtained from XPT and MDT tools are used to calibrate pore pressure.
However, the MDT data are not available for the whole interval of the borehole and the measured values are restricted to the some parts of the reservoir interval only. Hence, uncertainties anticipated in approximation of pore pressure are possibly due to acquirement of P-wave transit time log data and measurement of the XPT/MDT information as calibration points.

Vertical stress (σ v )
Vertical stress (lithostatic pressure) is the pressure or stress produced by the weight of the overlaying materials. Rock density is generally obtained using wireline logs or from core density. In this study, only bulk density logs were used. Vertical stress (σ v ) was computed by integrating formation bulk density ρ b from surface to TD using Eq. (1) (Aadnoy and Looyeh 2010).
where g is the gravity (gravitational acceleration).
Because of the uncertainties corresponding to acquisition of density log, the estimated effective stress is subjected to uncertainties. The Beta distribution function was applied to present the variation of vertical stress as is shown in Fig. 5.

Horizontal stresses (σ hmin & σ Hmax )
Comparing methods for determination of two horizontal stresses, the value of σ h is more straightforward to determine if it happens to be the minimum in situ stress σ 3 and is less than the vertical stress (lithostatic pressure). With leak-off test (LOT), mini-frac or extended leak-off test (XLOT), indirect measurements of σ 3 (and therefore σ h ) can be obtained with reasonable accuracy. However, the magnitude of the maximum horizontal stress σ H is more difficult to determine, and it is best quantified by an integrated approach that applies the above mentioned constraints, employs some appropriate modeling and involves rigorous calibration to observed features such as wellbore breakout or drillinginduced tensile fractures.
In addition to the uncertainties corresponding to acquisition of each individual log, the operational and data interpretation errors associated with the downhole stress measurements, which are utilized for calibration of minimum horizontal stress, will add more uncertainties to the determination of horizontal stresses. The difficulty in direct calibration of maximum horizontal stress is a major origin of uncertainty in estimation of this principal stress.
Poroelastic method is one the well proven developed approaches which successfully utilized in this study to determine of horizontal stresses in this tectonically active field (Maleki et al. 2014).
Probability density function of the Weibull distribution was utilized to monitor the data variation associated with minimum and maximum horizontal stresses corresponding to Well X at the section with RF Stress regime. This distribution generates a distribution with the shape parameter alpha and a scale parameter beta. The results are shown in Fig. 5. Although there might be a little uncertainty in using this distribution function, the fitted function seems to be very sophisticated. The only problem is the fact that nobody can be sure about the accuracy of the values obtained for the maximum and minimum horizontal stresses in the first place. On the other hand, Log-logistic probability distribution function was used to capture most variation of horizontal stresses corresponding to Well X at the section with mainly SS Stress regime (3470-3700 m) as shown in Fig. 5.
Although the uncertainty in vertical stress of both intervals and pore pressure of 3470-3700 m interval contribute relatively small amounts which results in negligible change in the required mud weight, these parameters must be considered in the model because they are one of the effective parameters on the geomechanical effects of a well and safe drilling mud windows.

Quantitative risk assessment and mud weight determination
Wellbore instability due to rock failure is initiated by two major types of failure, tensile (cause by high mud weight) or shear (caused by low mud weight). Several techniques exist for forecasting rock failure (and wellbore instability). The most commonly used failure criteria which also are implemented in this study include Mohr-Coulomb criteria to determine shear failure, and Maximum Tensile Stress criteria to determine tensile failure.
The design of mud weight ought to pursue the accompanying rule: safe mud weight should be higher than the lowest safe mud weight and lower than the highest safe mud weight. By and large, the lowest value of safe mud weight required to be determined by utilizing the pore pressure and collapse (breakout) pressure, whereas the highest value of safe mud weight needs to be determined using the mud loss and tensile fracture pressure.
Quantitative risk assessment was employed to examine the impact of uncertainties in the parameters used to assess this design on the required mud weight. The minimum, most likely, and maximum amounts of each parameter of each interval are given by probability density functions. The mentioned values for interval 3470-3700 m are tabulated in Table 4 as follows.
Succeeding determination of each parameter distribution function, it is time to build a response surface for the critical drilling fluid density. Once the input uncertainties have been quantified, response surfaces for the well breakout (Fig. 6) and the mud loss (not shown here) can be defined. The stochastic result produced by these simulations has been utilized to build response surface for theoretical mud weight shown in Fig. 6. As such, once required mud weights for individual input variability have been calculated, it was plotted for these input parameters. In this regard, all the data in the range between the minimum value and the maximum value of each parameter were used in building the response surfaces. Accordingly, response surfaces reveal the sensitivity of the predicted mud weight to the uncertainty of each input parameter.
The data cluster has been fitted by linear or non-linear regression analysis that is used to fit the surfaces to theoretical values, where required mud weights are function of input parameters.
It should be noted that response surfaces like those shown in Fig. 6 are very case specific. As can be seen, rock UCS is an important factor in this case which has significant effect on the well mud weight, while the friction angle has less influence on it. From the other point of view, the uncertainty in vertical stress contributes a relatively small amount between 88.2 and 88.7 MPa which results in negligible change in mud weight required to stabilize the well whereas the magnitude of the principal horizontal stresses and pore pressure are quite important. Also, in contract to S Hmax and pore pressure, increase in Shmin decreases the amount of the required mud weight.
Accordingly, in this case, the uncertainty in UCS and S Hmax contributes the most to the uncertainty in the required mud weight (0.3 gr/cc and 0.58 gr/cc, respectively). In other words, the mud pressure in this interval is controlled mainly by fairly large uncertainty in the uniaxial compressive strength and the least principle stress. In general, sensitivity analyses such as those are shown in Fig. 6 are useful for guiding data collection (in future wells or in real-time) for reducing uncertainty in geomechanical models, evaluating the importance of laboratory rock strength determinations, etc. It is worthy to note that geomechanical core measurements are very important to calibrate the log-driven mechanical properties and deliver reliable results. In order to obtain static elastic properties from the dynamic properties, correlation from core data should to be used. Moreover, core data can be used to obtain rock strength properties. Using the trends (inclinations) of response surfaces, Monte-Carlo simulation (MCS) was implemented which creates cumulative probability density curves for wellbore instabilities related Fig. 6 Response surfaces for the key drilling parameters illustrating the sensitivity of the mud weight prediction and its degree of dependency on each parameter pressures such as collapse and lost circulation pressure. The probability of drilling success for both intervals is defined by mud weight window, considering all the uncertainties related to input parameters. Higher the probability of success indicates more narrow safe mud weight window which trends toward deterministic answer. Figure 7 shows the predicted safe mud window (left) and the cumulative likelihood of preventing wellbore collapse and lost circulation for both intervals as a function of the mud weight (right). The right-hand curves display the cumulative probability of avoiding borehole collapse and lost circulation for intervals with RF (upper right) and SS (lower right) stress regimes. The green horizontal line between the curves of both intervals indicate the range of mud weights that will simultaneously provide at least a 90% chance of avoiding both wellbore collapse and lost circulation (the operational mud window).
As it can be seen in Fig. 7, according to the implemented failure criterion, the mud weight shall be kept in the range between 2.33 and 2.46 gr/cc for the depth between 3395 and 3370 m and 1.82-1.97 gr/cc for the depth between 3470 and 3700 m to have a safe drilling operation with at least a 90% chance of success without exposing the well to breakouts and mud losses. It is obvious that due to the difference between the safe mud ranges of these two sections, before commencing the drilling of the second interval (3470-3700 m), the first one shall be cased.
It is worthy to note that the mud weight used to drill Well X in the first interval (3395-3470 m) and the second interval (3470-3700 m) were 2.25 gr/cc MPa and 1.95 g/cc, respectively. Since some collapses were observed while drilling Fig. 7 Safe predicted mud weight window (left) and cumulative likelihood of success (avoiding wellbore failure) for RF (upper right) and SS (lower right) stress regimes the first interval. On the other hand, according to the well daily drilling reports, no significant problem happened while drilling the second interval. Therefore, it seems that the mud weight needs to be increased in the first section which fairly matches the results obtained from QRA.

Conclusions
The aim of this study was to implement the quantitative risk assessment technique utilizing derived rock elastic and strength properties at a drilled well case to determine the disparities, errors and uncertainties associated with the deterministic outcomes, build limit state functions and response surfaces, recognize most influential parameters and their sensitivities, and evaluate how the likelihood of success can be maximized by carrying out the fundamental quantitative risk assessment technique. This uncertainty investigation can donate to decrease the non-productive times by diminishing issues related kick, tight hole, packing off, stuck pipe, well collapse and circulation losses.
The results pointed out that safe drilling requires a mud weight window of between 2.33 and 2.46 gr/cc in the interval with RF stress regime (MD between 3395 and 3370 m) and 1.82-1.97 gr/cc in the interval with SS stress regime (MD between 3470 and 3700 m). This achievement fairly calibrated via the real drilling mud weight and data from caliper log. Furthermore, the response surfaces of each parameter indicated that uniaxial compressive strength of the material and maximum horizontal stress are the most influencing parameters in wellbore stability analysis and the utmost effect on predicted breakout mud weight is a result of uncertainty in these two parameters. The consequence of the discrepancy in the maximum horizontal stress originates from the fact that the bound on this stress is enormous. Regarding UCS, it is necessary to perform laboratory testing to be able to calibrate the obtained amount from the empirical relations and reduce related uncertainty.
This methodology can definitely add considerable values in successful drilling and completion of future wells and corresponding cost efficiency.
Funding The authors received no financial support for the research, authorship, and/or publication of this article.