3D Geomechanical Model Construction for Wellbore Stability Analysis in Algerian Southeastern Petroleum Field

: The main objective of this research work was the wellbore stability evaluation of oil and gas wells based on a 3D geomechanical model, which as constructed using seismic inversion in a southeastern Algerian petroleum ﬁeld. The seismic inversion model was obtained by using an iterative method and Aki and Richards approximation. Since the correlation between the inversion model and the log data was high at the wells, the reservoir was efﬁciently characterized and its lithology carefully discriminated in order to build a reliable 3D geomechanical model. The latter was further used to suggest the drilling mud weight window for the ongoing wells (well 5) and to examine the stability of four previously drilled wells. The main contribution of this study is providing a 3D geomechanical model that allows the optimization of drilling mud weight parameters so that a wellbore’s stability is guaranteed, on the one hand, and, on the other hand, so that the reservoir damage brought about by excessive surfactant use can be prevented. Indeed, the mud parameters are not just important for the drilling process’s effectiveness but also for logging operations. Since the tools have limited investigation diameters, with excessive use of surfactant, the invaded zone can become larger than the tools’ investigation diameter, which makes their logs unreliable. Hence, the 3D geomechanical model presented here is highly recommendable for the proposition of new wells, entailing less exploration uncertainty and more controllable productivity.


Introduction
One of the most important steps in discovering new hydrocarbon fields is drilling wells. Once the field is discovered, the drilling system must be put in the appropriate location, and the well must then cross important geological layers to reach its targeted reservoir [1][2][3]. During drilling, the crossed formations will reveal different rock properties. Therefore, drillers must be cautious to ensure the well's stability and prevent its collapse as a result of differential pressure change. The drilling mud density is determined in a practical way that guarantees the underbalance status of the well [4], but high density values may cause damage to the drilled reservoir. In this context, this study surveyed the state of the art of pertinent methods in order to determine a compromise between ensuring the maximum wellbore stability and minimizing reservoir damages.
In this framework, the field of study selected has suffered from reservoir damage caused by the use of inappropriate drilling mud density values at the levels of wells 1, 2, 3, and 4. The study region in this research work was in the southeast of Algeria, and seismic inversion and a 3D geomechanical model were employed using seismic and well-logging

1D Geomechanical Model Construction
To construct a 1D geomechanical model, the following parameters were experimentally determined: the static Young's modulus, the static Poisson coefficient, the unconfined compressive strength (UCS), the resistance of the rock to the tension (TSTR), the friction angle, the pore pressure, and the vertical and horizontal constraints (minimum and maximum). Young's modulus is a mechanical property that quantifies the proportion between tensile stress σ and axial strain ε in the linear elastic region of a rock, extracted directly from the velocity of the P-wave (Vp) and the S-wave (Vs) in the seismic data [18]. The UCS is the maximum axial compressive stress that a cylindrical sample of rock can withstand under unconfined laboratory conditions. The UCS relationship and the static Poisson ratio are obtained from laboratory test results. We took into consideration the local mechanical properties of the rock formation in the reservoir under study (Turonian) [19]. Based on the laboratory test results, the following values were obtained: static Young's modulus = 0.8 × the dynamic Young's modulus, static Poisson coefficient = dynamic Poisson factor, and TSTR = 6~8% of UCS [20]. Further, the total horizontal stresses

1D Geomechanical Model Construction
To construct a 1D geomechanical model, the following parameters were experimentally determined: the static Young's modulus, the static Poisson coefficient, the unconfined compressive strength (UCS), the resistance of the rock to the tension (TSTR), the friction angle, the pore pressure, and the vertical and horizontal constraints (minimum and maximum). Young's modulus is a mechanical property that quantifies the proportion between tensile stress σ and axial strain ε in the linear elastic region of a rock, extracted directly from the velocity of the P-wave (Vp) and the S-wave (Vs) in the seismic data [18]. The UCS is the maximum axial compressive stress that a cylindrical sample of rock can withstand under unconfined laboratory conditions. The UCS relationship and the static Poisson ratio are obtained from laboratory test results. We took into consideration the local mechanical properties of the rock formation in the reservoir under study (Turonian) [19]. Based on the laboratory test results, the following values were obtained: static Young's modulus = 0.8 × the dynamic Young's modulus, static Poisson coefficient = dynamic Poisson factor, and TSTR = 6~8% of UCS [20]. Further, the total horizontal stresses were calculated using 1D poro-elastic theory. Since the rocks of this field have linear elastic and isotropic horizontal stresses, the calculation involved the Young's modulus E, the Poisson ratio ν, the Biot constant α, the vertical litho-static stress σ ν , and the minimum and maximum strains ε Energies 2022, 15, x FOR PEER REVIEW were calculated using 1D poro-elastic tic and isotropic horizontal stresses, th Poisson ratio ν, the Biot constant α, t and maximum strains εℎ and ε [3,15 mination of a certain window density mize the following damage: 1. Kicks: the kicks are represented Figure 2 [21]. When greater than will be a collapse of the wellbore and Energies 2022, 15,7455 ε H [3,15]. These different measurements allowed the determination of a certain window density for the mud used during drilling in order to minimize the following damage:

1.
Kicks: the kicks are represented by the pore pressure shown in grey in the track of Figure 2 [21]. When greater than the pressure generated by the mud density, there will be a collapse of the wellbore [4]; 2.
Shear failure: the shear failure is presented in red in the left track of Figure 2. It provides formation-breaking pressure due to the maximum in situ stress [21,22]; 3.
Breakdown: the breakdown is presented in dark blue in the left track of Figure 2. A break in traction occurs when the mud density is greater than the maximum horizontal stress [21]. were calculated using 1D poro-elastic theory. Since the rocks of this field have linear elastic and isotropic horizontal stresses, the calculation involved the Young's modulus E, the Poisson ratio ν, the Biot constant α, the vertical litho-static stress σ , and the minimum and maximum strains εℎ and ε [3,15]. These different measurements allowed the determination of a certain window density for the mud used during drilling in order to minimize the following damage: 1. Kicks: the kicks are represented by the pore pressure shown in grey in the track of Figure 2 [21]. When greater than the pressure generated by the mud density, there will be a collapse of the wellbore [4]; 2. Shear failure: the shear failure is presented in red in the left track of Figure 2. It provides formation-breaking pressure due to the maximum in situ stress [21,22]; 3. Breakdown: the breakdown is presented in dark blue in the left track of Figure 2. A break in traction occurs when the mud density is greater than the maximum horizontal stress [21].
The results of the 1D geomechanical modelling for well 1 and well 2 are presented in Figure 3a,b, respectively. These are examples of the results obtained from the four wells.  The results of the 1D geomechanical modelling for well 1 and well 2 are presented in Figure 3a   The different phases of well 1 and well 2 are oval along the Marly intervals, mainly in the Campanian and Santonian formations, where severe types of damage were observed. These intervals were characterized by a very low UCS and, thus, they required a higher mud weight to stabilize the wellbore during drilling. The total drilling mud losses reported in the Maastrichtian formation are related to the presence of natural fractures that are not open to the fracture gradient ( [24,25]). The Turonian and Coniacian formations, mainly in well 2, seemed to be more stable compared to the whole formation. We reduced the drilling mud weight in the Turonian formation to 1.05 sg in order to limit mud losses in the reservoir intervals. To confirm these conclusions and use them for the next planned wells, it was necessary to obtain reliable results. Therefore, we constructed a 3D geomechanical model based on seismic inversion. The model was then compared to 1D models for validation and generalization, as detailed in section four.  The different phases of well 1 and well 2 are oval along the Marly intervals, mainly in the Campanian and Santonian formations, where severe types of damage were observed. These intervals were characterized by a very low UCS and, thus, they required a higher mud weight to stabilize the wellbore during drilling. The total drilling mud losses reported in the Maastrichtian formation are related to the presence of natural fractures that are not open to the fracture gradient ( [24,25]). The Turonian and Coniacian formations, mainly in well 2, seemed to be more stable compared to the whole formation. We reduced the drilling mud weight in the Turonian formation to 1.05 sg in order to limit mud losses in the reservoir intervals. To confirm these conclusions and use them for the next planned wells, it was necessary to obtain reliable results. Therefore, we constructed a 3D geomechanical model based on seismic inversion. The model was then compared to 1D models for validation and generalization, as detailed in Section 4.

Generation of V p and V s Volumes
The relation between the V p and the acoustic impedance (AI) was extracted by using the cross plot v p = f (AI) from logging data from the four wells [26]. It can be seen from Figure 4 that the relationship is linear, and it is expressed by Equation (1):

Generation of Vp and Vs Volumes
The relation between the Vp and the acoustic impedance (AI) was extracted by using the cross plot = ( ) from logging data from the four wells [26]. It can be seen from Figure 4 that the relationship is linear, and it is expressed by Equation (1): As the acoustic impedance volume was obtained from the seismic inversion results, a simple calculation was used to deduce the volume of the acoustic velocities, as shown in Figure 5a. Since the volumes of the ratio / and of Vp were calculated, the volume of Vs could, therefore, be directly deduced, as demonstrated in Figure 5b [16].

Density Volume Generation
Using the acoustic impedance and the acoustic velocity volumes, the density volume was extracted by using the relationship = × , as shown in Figure 5c. A comparison between the obtained inverted density model and the density logs from well 1 and well 2 is shown in Figure 6. The latter shows density values from inverted volumes (coloured) compared to logging data (black curves). Based on Figure 6, it can be concluded that the inverted model had high correlations at well 1 and well 2 in the studied field, as well as good correlations at well 3 and well 4. Therefore, we were able to proceed to 3D geomechanical model construction. The workflow of this process is summarized in Figure 7 ( [6,[27][28][29][30]).

Digital Grid Creation
When building a 3D geomechanical model, the first step is to create the grid in which the computation of the constraints will be carried out [31]. The grid is guided by the geological model of the region and must contain the overlying layers in addition to the reservoir's layers ( [6,32]).  As the acoustic impedance volume was obtained from the seismic inversion results, a simple calculation was used to deduce the volume of the acoustic velocities, as shown in Figure 5a. Since the volumes of the ratio V p /V s and of V p were calculated, the volume of V s could, therefore, be directly deduced, as demonstrated in Figure 5b [16]. (c)

Density Volume Generation
Using the acoustic impedance and the acoustic velocity volumes, the density volume was extracted by using the relationship AI = v p × ρ, as shown in Figure 5c. A comparison between the obtained inverted density model and the density logs from well 1 and well 2 is shown in Figure 6. The latter shows density values from inverted volumes (coloured) compared to logging data (black curves). Based on Figure 6, it can be concluded that the inverted model had high correlations at well 1 and well 2 in the studied field, as well as good correlations at well 3 and well 4. Therefore, we were able to proceed to 3D geomechanical model construction. The workflow of this process is summarized in Figure 7 ( [6,[27][28][29][30]).

Digital Grid Creation
When building a 3D geomechanical model, the first step is to create the grid in which the computation of the constraints will be carried out [31]. The grid is guided by the geological model of the region and must contain the overlying layers in addition to the reservoir's layers ( [6,32]).
The geological model of the studied field consisted of 8 horizons, from the Maastrichtian to the Aptian. Forty-nine (49) faults, provided from the structural interpretation of the seismic data, were also included in the constructed model. Thus, from the geological point of view, a simple rectangular grid covering all horizons was created; containing 37 layers. Particular attention was given to the reservoir layers and the cover layer, increasing the resolution gradually to provide more details [33,34]. At the end of this step, the obtained grid was composed of 1,583,550 cells, as shown in Figure 8. It was used for reservoir characterization and rock typing [35]. The geological model of the studied field consisted of 8 horizons, from the Maastrichtian to the Aptian. Forty-nine (49) faults, provided from the structural interpretation of the seismic data, were also included in the constructed model. Thus, from the geological point of view, a simple rectangular grid covering all horizons was created; containing 37 layers. Particular attention was given to the reservoir layers and the cover layer, increasing the resolution gradually to provide more details [33,34]. At the end of this step, the obtained grid was composed of 1,583,550 cells, as shown in Figure 8. It was used for reservoir characterization and rock typing [35].

Rock Typing
Once the grids were created, they were filled with the mechanical properties that characterized the rock; namely, the density, the elastic properties (Young's modulus, Poisson ratio, and Biot coefficient), and the resistance properties (UCS, the resistance of the

Rock Typing
Once the grids were created, they were filled with the mechanical properties that characterized the rock; namely, the density, the elastic properties (Young's modulus, Poisson ratio, and Biot coefficient), and the resistance properties (UCS, the resistance of the rock to tension, and the friction angle) [10]. These properties were calculated using the same method as in the 1D model constructed earlier in this study and with the same mathematical relationships [17]. The inputs of the model were the density volumes, V p and V s , and the dynamic elastic properties, which were determined directly from the inputs through the Gassmann equations (Equations (2)-(5)) [36,37].
where G dyn is the saturation gain function, k dyn is the bulk modulus of the rock, K dyn is the bulk modulus of the fluid, and PR dyn is the Poisson ratio. These properties are called "dynamic" because of their high-frequency contents. Nevertheless, their calibration can be guaranteed using laboratory tests on cores: the UCS cube was calculated from a UCS vs. Young's modulus correlation based on the regional core tests in the four wells, and the TSTR was estimated at 6-8% of the UCS [28]. The obtained 3D rock properties are illustrated in Figure 9. can be guaranteed using laboratory tests on cores: the UCS cube was calculated from a UCS vs. Young's modulus correlation based on the regional core tests in the four wells, and the TSTR was estimated at 6-8% of the UCS [28]. The obtained 3D rock properties are illustrated in Figure 9. It can be noted from Figure 10 that the properties calculated from the 3D cubes, which resulted from the inversion, remained faithful to the 1D models at wells 1, 2, and 3 [38]. Some inaccuracies were observed at the level of well 4 since the inputs (Vp, Vs, and Rhob cubes) were generated from logging data, which presented some unreliability problems for density logs at this well due to the presence of washout and doglegs (bad quality of the borehole).  It can be noted from Figure 10 that the properties calculated from the 3D cubes, which resulted from the inversion, remained faithful to the 1D models at wells 1, 2, and 3 [38]. Some inaccuracies were observed at the level of well 4 since the inputs (Vp, Vs, and Rhob cubes) were generated from logging data, which presented some unreliability problems for density logs at this well due to the presence of washout and doglegs (bad quality of the borehole).

Calculation of Limiting Conditions and Constraints
The constructed 3D geomechanical model allowed for the inclusion of all structural and mechanical information. In the constraint computation phase, a simulation was carried out using the finite element technique [39]. The outputs of the simulation were extracted and used for reservoir analyses and the design of new wells [40]. The simulation was undertaken in three steps. In the first step, the gravitational load was estimated over the entire volume to generate a cube of vertical stresses, which included faults and different mechanical properties. In step two, the boundary conditions were used to initialize the state of the constraints, and the pore pressure was also estimated. It should be noted that the state of stress at this step represented the initial conditions [18]. In step three, the constraints and pore pressure were iteratively calculated until they correlated with data at a predefined error threshold [17]. The different constraints that affected the model were calculated as follows: 1.
The vertical stress was determined from the density volume, and it was not quite accurate due to the quality of the seismic data. Nevertheless, the estimated vertical stresses remained close to the real values obtained from core data; 2.
The horizontal stresses, caused by the effects of tectonics, were estimated from the results of the 1D geomechanical model, which was established at the level of the four wells [30].
Effectively, the regional constraints were adjusted iteratively until they corresponded to the profiles of the different wells. In this study, seven simulations were required to find the best correlation between the estimated parameters and the data. The simulations were aimed at modelling the distribution of the stresses in the reservoir by taking into account the structures resulting from the seismic data interpretation ( [21,40]). The values obtained from the simulation were the minimum horizontal deformation E h = 0 ∼ 0.1 × 10 −3 , the maximum horizontal deformation E H = 0.1 ∼ 0.4 × 10 −3 , and the azimuth equal to 45 • , as shown in Figure 11. The results, shown in Figure 12, demonstrated a good correlation between the 3D and the 1D models. Since the 1D model could not include structural effects, the comparison was significant only for wells that were not located near to faults or structural features that could induce major stress rotations ( [18,40]). the structures resulting from the seismic data interpretation ( [21,40]). The values obtained from the simulation were the minimum horizontal deformation ℎ = 0~0.1 × 10 −3 , the maximum horizontal deformation = 0.1~0.4 × 10 −3 , and the azimuth equal to 45°, as shown in Figure 11. The results, shown in Figure 12, demonstrated a good correlation between the 3D and the 1D models. Since the 1D model could not include structural effects, the comparison was significant only for wells that were not located near to faults or structural features that could induce major stress rotations ( [18,40]).

Determination of the Mud Density Window
As soon as the position of the new well (well 5) was defined, the 3D geomechanical model was used to analyse the wellbore stability and determine its ideal direction, in addition to establishing the optimal mud density window to minimize reservoir damage ( [6,17,41]). The analysis included several trajectories, but only the alternative trajectories were selected for well 5. The proposed approach involved the creation of a mud density cube for each cell in the model and the definition of a limit of rupture (shear failure gradient), a limit of mud loss related to an induced fracture (fracture gradient), and a rupture limit for breakdown [42]. The size of the mud density window depended on the difference between the limits of these different gradients [32]. The stability of the wellbore depended on the relative orientation of the main stresses and the trajectory of the well. Thus, the creation of a single mud density cube required few assumptions about the orientation of the well. Typically, an initial cube is created assuming a vertical well (zero inclination), and then alternate cubes are created with tilts and azimuths, which represent the general intended directions of the wells [21]. Based on these assumptions, three mud density cubes were generated: a cube for the vertical well (0° inclination and 0° azimuth), a cube for the horizontal well simulated in the direction of the minimum horizontal stress (90° inclination and 45° azimuth), and a cube for the horizontal well simulated in the direction of the maximum horizontal stress (90° inclination and 135° azimuth).
The time slices shown in Figures 13-15 represent the targeted reservoir (the Turonian) obtained from the three generated cubes, respectively. The simulation results show that the horizontal trajectories of the wells were slightly more stable than the vertical trajectories in the Turonian layer. There was a small difference in the mud density window between drilling in the direction of the minimum and maximum horizontal stresses. However, drilling in the direction of the minimum horizontal stress seemed better. As it is subject to uncertainties in the properties of the faults, it is strongly recommended to keep the wellbore at least 200 m far from the faults if possible. The stresses in these areas are

Well Design Determination of the Mud Density Window
As soon as the position of the new well (well 5) was defined, the 3D geomechanical model was used to analyse the wellbore stability and determine its ideal direction, in addition to establishing the optimal mud density window to minimize reservoir damage ( [6,17,41]). The analysis included several trajectories, but only the alternative trajectories were selected for well 5. The proposed approach involved the creation of a mud density cube for each cell in the model and the definition of a limit of rupture (shear failure gradient), a limit of mud loss related to an induced fracture (fracture gradient), and a rupture limit for breakdown [42]. The size of the mud density window depended on the difference between the limits of these different gradients [32]. The stability of the wellbore depended on the relative orientation of the main stresses and the trajectory of the well. Thus, the creation of a single mud density cube required few assumptions about the orientation of the well. Typically, an initial cube is created assuming a vertical well (zero inclination), and then alternate cubes are created with tilts and azimuths, which represent the general intended directions of the wells [21]. Based on these assumptions, three mud density cubes were generated: a cube for the vertical well (0 • inclination and 0 • azimuth), a cube for the horizontal well simulated in the direction of the minimum horizontal stress (90 • inclination and 45 • azimuth), and a cube for the horizontal well simulated in the direction of the maximum horizontal stress (90 • inclination and 135 • azimuth). The time slices shown in Figures 13-15 represent the targeted reservoir (the Turonian) obtained from the three generated cubes, respectively. The simulation results show that the horizontal trajectories of the wells were slightly more stable than the vertical trajectories in the Turonian layer. There was a small difference in the mud density window between drilling in the direction of the minimum and maximum horizontal stresses. However, drilling in the direction of the minimum horizontal stress seemed better. As it is subject to uncertainties in the properties of the faults, it is strongly recommended to keep the wellbore at least 200 m far from the faults if possible. The stresses in these areas are unpredictable and the optimal mud density windows that guarantee wellbore stability will be very narrow. unpredictable and the optimal mud density windows that guarantee wellbore stability will be very narrow.   unpredictable and the optimal mud density windows that guarantee wellbore stability will be very narrow.   A comparison between the mud weight windows (MWWs) for well 5 and the other four wells in the studied field is shown in Figure 16. The MWW was more favourable in well 5 than in the other wells. In the Turonian formation, the main challenge was to keep the rupture limit lower than the interstitial pressure at certain intervals. However, this could reduce the drilling mud weight in this section to its limit, and it could even lead to mud loss in the naturally fractured subintervals ( [25,39]). The main limitation of this study is that the reliability and precision of the 3D geomechanical model obtained were strongly related to the vertical resolution of the seismic data. Furthermore, in order to ensure better results, it would be necessary to ensure good acquisition and design of the seismic field, which is not always available in southeastern Algeria due to its mountainous characteristics. A comparison between the mud weight windows (MWWs) for well 5 and the other four wells in the studied field is shown in Figure 16. The MWW was more favourable in well 5 than in the other wells. In the Turonian formation, the main challenge was to keep the rupture limit lower than the interstitial pressure at certain intervals. However, this could reduce the drilling mud weight in this section to its limit, and it could even lead to mud loss in the naturally fractured subintervals ( [25,39]). The main limitation of this study is that the reliability and precision of the 3D geomechanical model obtained were strongly related to the vertical resolution of the seismic data. Furthermore, in order to ensure better results, it would be necessary to ensure good acquisition and design of the seismic field, which is not always available in southeastern Algeria due to its mountainous characteristics.

Discussion
The 3D geomechanical model constructed provided us with the tools to estimate the optimal mud density window to be used to drill the planned well (well 5). This mud window could guarantee the stability of the wellbore of well 5, and there was no risk of reservoir invasion or fracturing in the Turonian formation. The mud window results are shown in Figure 16, and the curves can be interpreted as follows: 1. Gary block: the pore pressure values. If they become higher than the pressure created by the mud density, a kick could occur in the Turonian formation at well 5. The kicks correspond to lower shelling of the walls that that applied by the geological formations around the well 2. Yellow block: the fracture limit (shear failure or breakout). This occurs when the mud density is below the maximum in situ stress, but it should be controllable while drilling well 5 [42]; 3. Light blue block: the fracture gradient limit, which is caused by the minimum hori- Figure 16. Mud weight windows for the four wells compared to that of well 5.

Discussion
The 3D geomechanical model constructed provided us with the tools to estimate the optimal mud density window to be used to drill the planned well (well 5). This mud window could guarantee the stability of the wellbore of well 5, and there was no risk of 1.
Gary block: the pore pressure values. If they become higher than the pressure created by the mud density, a kick could occur in the Turonian formation at well 5. The kicks correspond to lower shelling of the walls that that applied by the geological formations around the well 2.
Yellow block: the fracture limit (shear failure or breakout). This occurs when the mud density is below the maximum in situ stress, but it should be controllable while drilling well 5 [42]; 3.
Light blue block: the fracture gradient limit, which is caused by the minimum horizontal in situ stress. When the mud density is higher than the latter, this will open the existing fractures. This should also be controlled during the drilling of well 5 in order to avoid fracturing the Turonian formation; 4.
Dark blue block: the breakdown (tensile break). This occurs when the mud density is greater than the maximum horizontal stress. The latter should also be considered when drilling well 5.
In Figure 16, the drilling mud density window is represented by the white block located in the middle of the tracks. Hence, to avoid the various fractures previously mentioned for well 5, mud densities were estimated for each interval based on the results of the 3D geomechanical model [43]. For the interval 0-600 m MD (Maastrichtian), the obtained mud weight was MW = 1.05 g/cm 3 (WBM); for the interval 600-1200 m MD (Maastrichtian-Campanian-Santonian), MW = 1.23 g/cm 3 (OBM); for the interval 1200-1681 m MD (Coniacian), MW = 1.05 g/m 3 (OBM); and for the interval 1681-2000 m MD (Turonian-Cenomanian), the obtained mud weight was MW = 1.05 g/cm 3 (OBM). These mud density values can ensure wellbore stability with minimum damage to the targeted reservoir, as confirmed during the process of drilling well 5 and its reservoir evaluation.

•
In this paper, a 3D geomechanical model was constructed in order to help drillers ensure the stability of the wellbores of hydrocarbon wells without damaging the reservoir in a petroleum field located in the southeast of Algeria; • The seismic inversion cubes were used as the input data for the construction of the 3D geomechanical model. The 1D geomechanical profiles, generated at the wells' scales, were used as basic tools for the construction of the 3D geomechanical model; • The 3D geomechanical model provided information on the spatial distribution of the existing stresses in the reservoir and all subjacent layers; • As a result, this study provides a powerful and reliable tool for drilling new wells and to optimize profitability by avoiding drilling incidents. It allows the estimation of the different mechanical properties of rocks at the seismic scale; • This is an effective way to determine the mud density to be used when drilling new wells in a studied area, thus allowing the determination of the preferred direction for their boreholes; • Furthermore, by obtaining such information, simulations of on-going wells can be performed at any position in the volume, and in any direction, so that an estimate of the drilling mud density window can be obtained simultaneously; • In addition, this model can help drillers choose the location and the optimal direction of new wells with the maximum possible wellbore stability; • The 3D geomechanical model can be also used to guide hydraulic fracturing supervisors in undertaking fracturing jobs in tight reservoirs based on the robust information on the spatial distribution of stresses, their directions, and their amplitudes.