Machine Learning: A Useful Tool in Geomechanical Studies, a Case Study from an O ﬀ shore Gas Field

: For a safe drilling operation with the of minimum borehole instability challenges, building a mechanical earth model (MEM) has proven to be extremely valuable. However, the natural complexity of reservoirs along with the lack of reliable information leads to a poor prediction of geomechanical parameters. Shear wave velocity has many applications, such as in petrophysical and geophysical as well as geomechanical studies. However, occasionally, wells lack shear wave velocity (especially in old wells), and estimating this parameter using other well logs is the optimum solution. Generally, available empirical relationships are being used, while they can only describe similar formations and their validation needs calibration. In this study, machine learning approaches for shear sonic log prediction were used. The results were then compared with each other and the empirical Greenberg–Castagna method. Results showed that the artiﬁcial neural network has the highest accuracy of the predictions over the single and multiple linear regression models. This improvement is more highlighted in hydrocarbon-bearing intervals, which is considered as a limitation of the empirical or any linear method. In the next step, rock elastic properties and in-situ stresses were calculated. Afterwards, in-situ stresses were predicted and coupled with a failure criterion to yield safe mud weight windows for wells in the ﬁeld. Predicted drilling events matched quite well with the observed drilling reports.


Introduction
Geomechanics has shown its potential to cover a broad range of work during the lifespan of a field from exploration to production and then abandonment [1,2]. By examining the available data and comprehending the key issues observed while drilling previous wells, we can predict the optimum way of drilling and finding proper casing shoe locations in the development phase. By coupling fluid-flow with the stress-strain regime in the field, compaction and subsidence analyses can be performed [3,4]. The basic approach to geomechanics analysis is to process available data for predicting rock elastic properties, in-situ stresses, and pore pressures. Shear velocity (Vs) is a kind of data that is not always available for wells, especially old wells, and needs to be predicted.
The main methods for shear wave velocity prediction categorized in the following ones: heuristic models [5], effective medium models [6], and empirical methods, most importantly the Greenberg-Castagna model [7]. Empirical methods are easy and fast to apply [8] compared to other methods [9,10]. However, in most cases, results from the empirical relationships are not accurate enough due to the [10]: • simplicity in which other parameters affecting the shear wave velocity are not included, • the fact that such relationships derived from a particular formation, and • the fact there are few studies on the carbonate rocks (i.e., most of the empirical relationships are derived from sandstone reservoirs).
During the last decade, improved computing hardware and software has led to the prosperous application of machine learning (ML) in different areas of oil industry such as seismic data, petrophysical analysis including synthetic log generation or prediction [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], which has shown to be a promising tool to help address their problems in a rigorous, repeatable way. Such methods, by considering various available parameters, can give a better prediction of the missing data than simple linear methods [10,26,27]. ML techniques can be categorized into two main types, known as supervised and unsupervised learning. Supervised learning applies to cases in which we have a set of inputs and known responses, while in unsupervised learning, the response is not available, and the method tries to identify natural patterns or cluster in the data [28].
Linear regression as a common ML algorithm intends to find a linear relationship between inputs and a response based on minimizing the sum of squared difference of distance between actual data points and the predicted ones [29,30]. Depending on the number of inputs, this technique can be called simple (having only one input) or multiple (having more than one input) linear regression, and results in slope(s) and intercept. The artificial neural network (ANN) as a nonlinear statistical data modeling tool tries to simulate the behavior of a system composed of neurons and can model complex relationships between inputs and outputs or find patterns between them [10]. Neural networks may have the input layers, hidden layers, and output layers representing the raw information, processing units by finding some weight functions, and results, respectively [31]. Different researchers have used such methods successfully in different disciplines of the oil industry, e.g., [10,[32][33][34].
In this study, we used data from three wells located in an offshore field in the Persian Gulf to demonstrate the potential of high accuracy shear velocity prediction in hydrocarbon-bearing and non-bearing intervals. The field is located in the Persian Gulf and the structural trap (carbonate reservoir) is formed by basal forces and salt diapirism [35]. Then, we calculated the rocks elastic properties followed by in-situ stresses to make a 1D mechanical earth model (MEM) for each well. It was shown that the ANN method can predict shear wave velocity employing other available logs with high accuracy.
Numerical representation of the stress states and rock mechanical properties are known as the mechanical earth model (MEM). As mentioned before, shear wave velocity is one of the main components in geomechanical studies. In a gas field of the southern part of Iran, due to harsh drilling conditions, a geomechanical study should be performed to find the safe mud weight of future wells. There are three available wells and two of them (Well-2 and Well-3) have dipole sonic logs (including shear and compressional velocities), while the old well (Well-1) excludes shear wave velocity. Due to the necessity of including all wells in the study, machine learning methods along with conventional methods for shear wave velocity prediction were used and a comparison was drawn. It should be noted that Well-3 as the longest well was considered for fitting the models and Well-2 as the validation.
In order to investigate a relationship between a dependent variable and an independent variable (or a set of them), regression analyses are usually employed [36]. One of the most common correlations for shear wave velocity prediction is Castagna equations. Castagna et al. [37] proposed empirical equations for shear wave velocity prediction in sandstone, limestone, shale, and dolomite rocks as below [38]: Dolomite Vs = −0.078 + 0.583Vp Shale Vs = −0.867 + 0.770Vp (4) where Vp and Vs are compressional and shear velocity in km/s. Following the Castagna equation, a simple linear regression (SLR) model between available Vp and Vs data was also fitted. Equation (5) shows the linear regression model from the Well-3, which gives a better insight into the way shear velocity changes with compressional velocity in our filed, compared with the empirical models by Castagna. As seen in Figure 1, empirical equations do not present a good response as they are derived from a particular area or reservoir rocks comparing with the linear model from the data. Table 1 also corroborates this fact by comparing the standard error of the estimate from different models. present a good response as they are derived from a particular area or reservoir rocks comparing with the linear model from the data. Table 1 also corroborates this fact by comparing the standard error of the estimate from different models.
In the next step, along with compressional wave as the main predictor, other logged properties such as gamma ray and bulk density were also considered (Equation (6)). This method known as multiple linear regression (MLR), will help to consider different rock properties for predicting the shear velocity. where DTS, DTC, GR, and RHOZ are shear and compressional sonic (us/ft), Gamma ray (API), and density (gr/cc) logs, respectively. Note that, in Figure 1, Equation (5) is transformed into the velocity form (green line) rather than slowness to be compared with Castagna equations.  Finally, the artificial neural network (ANN) was used, which has the ability to learn and model non-linear and complex relationships [10] and is known to produce good fitting functions. ANNs constitute a set of artificial neurons designed in a manner that exchange signals on the communicational links [39]. The nodes of the input layer do not change the data but send them to the hidden layer. In the hidden layers, some sort of transformations is applied to the input data through  In the next step, along with compressional wave as the main predictor, other logged properties such as gamma ray and bulk density were also considered (Equation (6)). This method known as multiple linear regression (MLR), will help to consider different rock properties for predicting the shear velocity.
where DTS, DTC, GR, and RHOZ are shear and compressional sonic (us/ft), Gamma ray (API), and density (gr/cc) logs, respectively. Note that, in Figure 1, Equation (5) is transformed into the velocity form (green line) rather than slowness to be compared with Castagna equations. Finally, the artificial neural network (ANN) was used, which has the ability to learn and model non-linear and complex relationships [10] and is known to produce good fitting functions. ANNs constitute a set of artificial neurons designed in a manner that exchange signals on the communicational links [39]. The nodes of the input layer do not change the data but send them to the hidden layer. In the hidden layers, some sort of transformations is applied to the input data through some weight functions and then are linked to the output layer that corresponds to the results we are looking for. In this study, initially, to remove the effect of bad data (e.g., null values, bad log measurements) available data were processed. Afterward, data from wells 2 and 3 were selected as the training (25596 samples), test (5484 samples), and validation set (5484 samples). Then, a three-layer feed forward neural network ( Figure 2) with sigmoid hidden neurons and linear output neurons was used to build the model. The network was trained with the Levenberg-Marquardt backpropagation algorithm.
Energies 2020, 13, 3528 4 of 16 some weight functions and then are linked to the output layer that corresponds to the results we are looking for. In this study, initially, to remove the effect of bad data (e.g., null values, bad log measurements) available data were processed. Afterward, data from wells 2 and 3 were selected as the training (25596 samples), test (5484 samples), and validation set (5484 samples). Then, a three-layer feed forward neural network ( Figure 2) with sigmoid hidden neurons and linear output neurons was used to build the model. The network was trained with the Levenberg-Marquardt backpropagation algorithm. After gathering all required data, by assuming elastic isotropy and using the following equations, elastic moduli such as dynamic Young's modulus and Poisson's ratio were calculated [40], where ν and ED represent Poisson's ratio and dynamic Young's modulus. However, Equation (8) will provide a dynamic elastic modulus which needs to be converted to static modulus using available correlations such as the equation presented by Eissa and Kazi [41]. The reason for this conversion is due to the difference in the measurement condition [42,43]. Uniaxial compressive strength (UCS) and friction angle (φ) were also estimated based on Plumb's correlation [44]. Calculated rock properties for Well-1 are presented in Figure 3. The next step in analyzing wellbore stability is predicting pore pressure. In this study, MDT (modular formation dynamics tester) data in the reservoir section, and Eaton method in the clay-rich intervals were used to predict pore pressure and interpolated for other sections as well [45]. In order to find the clay-rich intervals, the cross plot of the Gamma ray versus density log, as in Figure 4, was used. The pore pressure profile was then calibrated by used mud weight and observed kicks while drilling. The Eaton equation is presented below [46], , (9) where Ppg is pore pressure gradient, OBG is overburden gradient, Ppn is normal pore pressure gradient, and DTCn is the transient time of compressional wave in the normally pressured zone.
Vertical stress (σv) was then calculated based on the weight of overburden [47], RHOB , where g is the gravitational acceleration, RHOB is bulk density log, and z refers to depth. Principal horizontal stresses should be predicted using indirect methods such as poroelastic horizontal strain method [48], by considering the tectonic strains: After gathering all required data, by assuming elastic isotropy and using the following equations, elastic moduli such as dynamic Young's modulus and Poisson's ratio were calculated [40], where ν and E D represent Poisson's ratio and dynamic Young's modulus. However, Equation (8) will provide a dynamic elastic modulus which needs to be converted to static modulus using available correlations such as the equation presented by Eissa and Kazi [41]. The reason for this conversion is due to the difference in the measurement condition [42,43]. Uniaxial compressive strength (UCS) and friction angle (ϕ) were also estimated based on Plumb's correlation [44]. Calculated rock properties for Well-1 are presented in Figure 3. The next step in analyzing wellbore stability is predicting pore pressure. In this study, MDT (modular formation dynamics tester) data in the reservoir section, and Eaton method in the clay-rich intervals were used to predict pore pressure and interpolated for other sections as well [45]. In order to find the clay-rich intervals, the cross plot of the Gamma ray versus density log, as in Figure 4, was used. The pore pressure profile was then calibrated by used mud weight and observed kicks while drilling. The Eaton equation is presented below [46], where P pg is pore pressure gradient, OBG is overburden gradient, P pn is normal pore pressure gradient, and DTC n is the transient time of compressional wave in the normally pressured zone. Vertical stress (σ v ) was then calculated based on the weight of overburden [47], where g is the gravitational acceleration, RHOB is bulk density log, and z refers to depth. Principal horizontal stresses should be predicted using indirect methods such as poroelastic horizontal strain method [48], by considering the tectonic strains: where σ h , σ H , ν, E, σ v , and P p are minimum horizontal stress, maximum horizontal stress, Poisson's ratio, Young's modulus, vertical stress and pore pressure, respectively. α as Biot's coefficient was set as 1.0, and ε x and ε y are two horizontal tectonic strains [35,49]. Trial and error on a range of values of ε x and ε y should be performed in a way that makes a balance between the value of σ h from Leak off test/Extended Leak off test [50] and observed drilling events [51]. Generally, ε x and ε y are tuned in a way to predict losses and breakouts, respectively, i.e., tectonic strains are calibration factors [48]. Figure 3 shows the in-situ stress for Well-1. As can be seen, a normal tectonic system (σ v > σ H > σ h ) is the dominant stress regime, which has been detected in all three wells and also corroborated with seismic data.
Energies 2020, 13, 3528 6 of 16       Safe mud weight window (SMWW) comprises of four mud weight limits by which different conditions are defined: pore (formation) pressure, breakout, minimum horizontal stress and breakdown. Drilling with lower mud pressure compared to formation pore pressure, wellbore washout and/or kick might be expected. On the other hand, drilling with higher mud pressure compared with the breakout limit of the formation, rock shearing and wellbore enlargement may occur. To calculate breakout (P w BO ) and breakdown (P w BD ) boundaries, the following equations were used [52]:

Results and Discussion
Different methods are available to predict shear wave velocity using available petrophysical logs. Initially, a simple linear regression using DTC as the sole independent variable was created; which works similar to a locally calibrated Greenberg-Castagna (GC) method. Table 2 shows the statistical properties of such a single linear regression model. As mentioned before, Well-3 as the longest well was considered for fitting the models and Well-2 as the validation. It is worth mentioning that to avoid data bias, we selected 70% of the data [53] randomly in a way that covers the whole well path. Figure 5 shows the division of the training set and testing set. As seen, both train and test data are covering the entire profile.

Results and Discussion
Different methods are available to predict shear wave velocity using available petrophysical logs. Initially, a simple linear regression using DTC as the sole independent variable was created; which works similar to a locally calibrated Greenberg-Castagna (GC) method. Table 2 shows the statistical properties of such a single linear regression model. As mentioned before, Well-3 as the longest well was considered for fitting the models and Well-2 as the validation. It is worth mentioning that to avoid data bias, we selected 70% of the data [53] randomly in a way that covers the whole well path. Figure 5 shows the division of the training set and testing set. As seen, both train and test data are covering the entire profile. Shear velocity is a function of different rock properties, and merely using DTC is not accurate enough. Therefore, in the next step, multiple linear regression was performed to incorporate other properties as well. Table 3

Depth (m)
Vs (km/s) Figure 5. Training (red dots) and test (blue squares) data cover the whole well path to avoid data bias.  Shear velocity is a function of different rock properties, and merely using DTC is not accurate enough. Therefore, in the next step, multiple linear regression was performed to incorporate other properties as well. Table 3 shows the statistical properties of the multiple linear regression model. T-value is the value of each coefficient divided by its standard error. High (absolute) T-values in the regression shows the importance of each parameter in DTS prediction. As expected, DTC, GR, and Density had the highest T-values and consequently the highest impact on the results, respectively. In addition to the overall improvement in the accuracy of predictions compared to the simple linear regression model (Figure 6), a further improvement was noticed in the hydrocarbon-bearing zone (Figure 7) backed up by the MSE (mean squared error) and R 2 (R-squared).   At the final step for Vs prediction, the neural network was used. Mean square error (MSE) curve plotted vs. the number of the training courses which is showing after 34 rounds model arrived to the best learning and the lowest error ( Figure 8). According to the MSE and R 2 results (Table 4), the efficacy of ANN in estimating shear wave velocity is very high. So, the ANN model is an accurate method to predict the shear wave velocity for Well-1 with no such log. Bukar et al. [9] and Akhundi et al. [10] also had similar results and mentioned ANN and multiple linear regression yield better results than empirical relations and single linear regression. By comparing the predicted wellbore instabilities with the observed drilling events, the robustness of the methodology is tested. Figure 9 shows the generated safe mud weight window for all three wells in the field. Borehole breakouts and breakdown in reservoir and overburden sections matched fairly well with caliper data. From the caliper data shown in Figure 9 for Well-3, breakouts are observed in the depth intervals of 850-1250 m, 1530-2080 m, and 2440-3130 m. The deepest breakout caused pipe sticking. The figure also shows that there is not a SMWW at the interval of 2515-2600 m, corroborated by the observed breakouts and complete mud loss. As can be seen in Figure 9, drilling events were predicted robustly. This infers that the rock elastic parameters were precited accurately enough using the proposed method.  (Table 4), the efficacy of ANN in estimating shear wave velocity is very high. So, the ANN model is an accurate method to predict the shear wave velocity for Well-1 with no such log. Bukar et al. [9] and Akhundi et al. [10] also had similar results and mentioned ANN and multiple linear regression yield better results than empirical relations and single linear regression. By comparing the predicted wellbore instabilities with the observed drilling events, the robustness of the methodology is tested. Figure 9 shows the generated safe mud weight window for all three wells in the field. Borehole breakouts and breakdown in reservoir and overburden sections matched fairly well with caliper data. From the caliper data shown in Figure 9 for Well-3, breakouts are observed in the depth intervals of 850-1250 m, 1530-2080 m, and 2440-3130 m. The deepest breakout caused pipe sticking. The figure also shows that there is not a SMWW at the interval of 2515-2600 m, corroborated by the observed breakouts and complete mud loss. As can be seen in Figure 9, drilling events were predicted robustly. This infers that the rock elastic parameters were precited accurately enough using the proposed method.

Method
Step    Table 4. Neural network properties.

Method
Step

Conclusions
Geomechanical studies have become a critical discipline in petroleum engineering and offer solutions for safe drilling, compaction, and subsidence in the field, CO2 sequestration, etc. Available data play a key role in such studies where shear velocity data are one of the necessary types of data.
In the current study, simple and multiple linear regression and ANN were utilized to estimate the shear wave velocity, and results showed ANN was very powerful. The neural network method can grasp the convoluted relationship between shear velocity and other petrophysical parameters,

Conclusions
Geomechanical studies have become a critical discipline in petroleum engineering and offer solutions for safe drilling, compaction, and subsidence in the field, CO 2 sequestration, etc. Available data play a key role in such studies where shear velocity data are one of the necessary types of data.
In the current study, simple and multiple linear regression and ANN were utilized to estimate the shear wave velocity, and results showed ANN was very powerful. The neural network method can grasp the convoluted relationship between shear velocity and other petrophysical parameters, specifically where there is heterogeneity in the reservoir, which can provide more accurate results than linear regressions. Moreover, results also corroborate the fact that empirical relationships are not accurate enough as they are derived from particular fields in the world.
In the next step, 1D MEM was built for the wells in the field, where drilling events such as wellbore breakouts were well predicted with the developed MEM.