4D Finite element modeling of stress distribution in depleted reservoir of south Iraq oilfield

The harvest of hydrocarbon from the depleted reservoir is crucial during field development. Therefore, drilling operations in the depleted reservoir faced several problems like partial and total lost circulation. Continuing production without an active water drive or water injection to support reservoir pressure will decrease the pore and fracture pressure. Moreover, this depletion will affect the distribution of stress and change the mud weight window. This study focused on vertical stress, maximum and minimum horizontal stress redistributions in the depleted reservoirs due to decreases in pore pressure and, consequently, the effect on the mud weight window. 1D and 4D robust geomechanical models are built based on all available data in a mature oil field. The 1D model was used to estimate all mechanical rock properties, stress, and pore pressure. The minimum and maximum horizontal stress were determined using the poroelastic horizontal strain model. Furthermore, the mechanical properties were calibrated using drained triaxial and uniaxial compression tests. The pore pressure was tested using modular dynamic tester log MDT. The Mohr–Coulomb model was applied in the 4D model to calculate the stress distribution in the depleted reservoir. According to study wells, the target area has been classified into four main groups in Mishrif reservoir based on depletion: highly, moderately, slightly, and no depleted region. Also, the results showed that the units had been classified into three main categories based on depletion state (from above to low depleted): L1.1, L1.2, and M1. The mean average reduction in minimum horizontal stress magnitude was 322 psi for L1.1, 183.86 psi for L1.2, and 115.56 psi for M1. Thus, the lower limit of fracture pressure dropped to a high value in L1.1, which is considered a weak point. As a result of changing horizontal stress, the mud weight window became narrow.


Introduction
The depleted reservoir is a popular feature of several oilfields in southern Iraq. Many production years without active water drive or water injection could cause depleted reservoir, which leads to change in in-situ stress, wellbore stability, formation pressure, and mud weight window. In hydrocarbon reservoirs, stress measurements appeared that the magnitude of horizontal stress related to decreasing in reservoir pressure (Addis 1997;Dahm et al. 2019;Neeamy & Selman 2020). As a result, the mud weight window becomes narrow, and mud weight should be an adjustment (Dutta & Kumar 2015;Li & Gray 2015). Therefore, lost circulation problems were faced during drilling operations when the mud weight was exceeded the fracture pressure, especially in depleted reservoir units (Abbas, Alhussainy, Abdul Hussien, & Flori, 2019;Allawi & Almahdawi 2019). Designing a safe mud window is necessary to prevent the hoop stress from causing a collapse, any change in vertical stress is neglecting (Akl 2015). Several problematic outcomes of a depletedcollapsed reservoir have been reported as the changes in stress associated with a depleted reservoir. Therefore, optimizes wellbore stability may be unachievable because the constraints on mud weight will appear (Addis, Cauley, & Kuyken 2001;Al Jawad & Tariq 2019;Alher et al. 2018;Fjar, Holt, Raaen, & Horsrud, 2008;Hettema et al. 2000;Johnson 1991;Neeamy & Selman 2020).

3
For reservoir deformation modeling, the most computations are based on two theories that depend on the deformation value, including linear elasticity and plasticity (Gai, Dean, Wheeler, & Liu 2003;Settari & Mourits 1998;Settari & Walters 2001). The finite element analysis provides an excellent possibility to deal with complex geometry calculations of a complete state for stress magnitudes and pore pressure. Moreover, the finite element analysis process accurately describes the reservoir geometry in terms of matrices and stress distribution. A finite element simulator uses to determine subsurface stress cases honoring the stress-strain relations set over each element from mapping seismically came mechanical properties onto the computational mesh (Jörg Herwanger & Koutsabeloulis 2011;NC Koutsabeloulis & Hope 1998;Rodriguez-Herrera, Zdraveva, & Koutsabeloulis 2014). (Nick Koutsabeloulis & Zhang 2009) applied static and dynamic reservoir modeling to determine the change in reservoir pore pressure, stresses and strain. Their study was mainly to estimate the initial stresses and predict the changes in the stress/strain over time. Their results presented necessary support to the drilling operation and reservoir management. (Jorg Herwanger 2014) studied the impact of the depleted reservoir and the change in pore pressure on the stresses. He finds out that the overburden formations transfer to the rock matrix, where the reduction in the stresses within the reservoir led to increasing the stresses to the side burden. (Hamid, Omair, & Guizada, 2017) summarized that the changes in the stress/ strain would impact the reservoir physical properties. This study applied a Finite Element approach for modeling the mechanical rock properties affected by formation pressure changes. His results showed that the permeability reduced by 18% in the depleted formation. (Nasreldin et al. 2019) constructed the 3D geomechanical model for sandstone depleted reservoir in the south Caspian basin in Turkmenistan using the simultaneous amplitude variation with offset (AVO) seismic inversion techniques. The study results showed the improvement in future well designing and field development.
The Y oilfield situates in Southern Iraq. Mishrif reservoir is one of the major reservoirs in Iraq because of prolific oil. The average pore pressure for Mishrif reservoir was 4000 psi in 2011 and dropped to 2700 psi close to the bubble point pressure of 2646 psi in 2020. Therefore, numerous instability problems have been faced through drilling. Tight holes, sticking pipes, lost circulation, and caving are most of these problems, as shown in Fig. 1. Nonproductive time NPT is produced from these problems because of increasing circulation and reaming time. The seriousness of these problems can cause sidetracks or losing the well.
The stratigraphy scheme used in the Y field was based on the terminology presented in exploration and appraisal well reports (G-1, G-2 & G-3) and various literature on the geology of Iraq, such as (Jassim & Goff, 2006;Nairn & Alsharhan 1997), and Abdelkarim et al. (2009). In this field, hydrocarbons were encountered in Mishref, Zubair, Ratawi and Yamama reservoirs. All formations were formed during the Cretaceous. Mishrif and Yamama are limestone formations, while Zubair and Ratawi are clastic formations.
Mishrif reservoir is subdivided into U1, U2, M1, M2, L1.1, L1.2 and L2 subunits. Hydrocarbon is present in Middle and Lower Mishrif at M1, L1.1, L1.2 and L2. Mishrif correlation was established based on G-1, G-2, G-3, G-4, G-A1 and G-A2 wells and extended to the wells drilled at Wellpads (WP) (A, B, C, E, F, J, P, Q), as shown in Fig. 2. The marl layer is thin and continuous that used as an indicator for middle Mishrif. Underlying the marl layer is the M1 reservoir which is present at all wells. Separating the M1 and the L1.1 is M2, which is a thin non-reservoir layer. Below the M2 is the L1.1 reservoir, characterized by the interbedding good and poor-quality reservoirs, as indicated by porosity logs in G-A2, G-1, G-A1 and G-4. Below the L1.2 layer is L2 which continues down to Rumaila.
This study used 1D and 4D numerical modeling to estimate the change in in-suite stress magnitude of M1, L1.1, and L1.2 units. 1D modeling was applied to model the geomechanical properties of the Mishrif reservoir to understand the root causes of the existing problems faced during drilling. The finite element of 4D modeling was used to estimate horizontal stress magnitude redistribution in two intervals in 2011 (before production) and 2020 (depleted reservoirs).

Methodology
This study shows the development of a four-dimensional stress finite elements model, incorporating pore pressure and stress elements, for a depleted reservoir section. Therefore, an integrated workflow was applied to build 1D and 4D models that determine the impact of pore pressure decrease in horizontal stress distribution. In addition, these models were calibrated to test their accuracy. The first step was data collected and audit; then, a 1D model was applied to estimate all mechanical properties. Later 4D using was built by dimensions target area setting and create layers with ten boreholes (ten wells). The mechanical properties of the 1D model were used as input data to the layers of the 4D model. The pore pressure was fed to the 4D model in two stages, the first of 2011(before production) and 2020 (depleted reservoir), then creating sideloads, mesh, calculating the horizontal stress, and mapping the stress. The workflow steps have been applied, which can be brief, as presented in Fig. 3.

1D geomechanical model
This model is based on the data from offset wells with pore pressure measurements and distributed on the entire reservoir. This process involves calculating the mechanical and elastic properties using correlations and physical equations . The necessary data are wells logs including density, gamma-ray, master/mud, neutron, resistivity, sonic (compression and shear wave velocities), formation micro-imager (FMI), caliper, in addition, and laboratory measurements.

Mechanical rock properties
The determination of mechanical rock properties plays a significant role in any geomechanical analysis. The basic mechanical rock properties include elastic properties (Poisson's ratio and Young's modulus (E)), and rock strength properties (internal friction angle (φ), unconfined compressive strength (UCS) and cohesion). The shear and compression acoustic wave velocities and bulk density logs were used to calculate the dynamic values of elastic parameters as follows: The modified Morales's correlation was applied to compute the static Young's modulus, which (Morales & Marcinew 1993) proposed as modified by Tom Bratton. This correlation is used for grain-supported rocks (sandstone). It is based on the total porosity ( ∅ c ) and the dynamic Young's modulus, as follows: -Also, assuming the static Poisson's ratio equal to the dynamic Poisson's ratio by (Zhang & Bentley 2005). Furthermore, the static bulk and shear moduli were calculated based on the following expressions: -  (Chang et al. 2006) presented several correlations used to compute unconfined compressive strength (UCS) for sandstone and shale rocks. (McNally 1987) found a correlation that applied to calculate the UCS of the sandstone formation using Eq. 8, as shown in Figs. 4 and 5 (at track number six under the name (UCS YME)).
The friction angle was estimated based on correlation of (Plumb, 1994). Equation 9 was applied to estimate FANG, as presented in Figs. 4 and 5 (at track number eight under the name (FANG)).
where: GR ∶ Gamma ray log, GR max : Maximum values of gamma ray log, GR min : Minimum values of gamma ray log.
Laboratory tests (Triaxial and uniaxial) were implemented on retrieved rock specimens to measure the (8) UCS = 1200 exp (−0.03Δtc)  Figure 4 presents mechanical properties, pore pressure, all stress of well G-4 before production in 2011, whereas, Fig. 5 shows mechanical properties, pore pressure, all stress of well G-C43 at the depleted reservoir.
Laboratory tests are essential to calibrate models. Thus, there will be confidence that the outcomes of the models are accurate. Therefore, these tests are used to determine the accuracy of a geomechanical model-excellent matched between laboratory results and parameters measured from logs, as shown in Fig. 4.

4D geomechanical model
There is three orthogonal principal stress which effected any point beneath the surface. Stress can be classified based on orientation and magnitude of stress. Generally, it is compressive, nonhomogeneous, and anisotropic (Veatch Jr & Moschovidis, 1986). In-situ stress has an essential function in wellbore structure, drilling, completion, stimulation, production, planning of the wells. Therefore, before performing any rock stress analysis and failure assessment, one should know about the in-situ stress. Knowledge of in-situ stress is helpful to comprehend and determine the orientation and magnitude of the horizontal stress, the impact of stress on operations (drilling and production), and identify the directions of rock failure. (B. S. Aadnoy & Ong, 2003).
The maximum and minimum stress tend to be the same in value when no tectonic activities; only the vertical stress effect is present. Therefore, there are different values for the horizontal stress and should be estimated when the faulting and tectonic activities have existed. Many geomechanics problems can be solved by knowing the horizontal stress's orientation and magnitude (B. Aadnoy & Looyeh, 2019).
The minimum horizontal stress can be estimated and matching with a leak-off and extended leak-off test (LOT and XLOT), mini-frac test (i.e., direct methods) (Yamamoto, 2003;M. Zoback et al., 2003). Furthermore, the best way to determine the minimum horizontal stress is hydraulic fracturing because the mechanical properties of the rock is not needed.
The finite element analysis was used to calculate the stress distribution of the Mushrif reservoir before and after production. Thus, the finite element analysis calculations were applied to estimate the magnitude and orientation of horizontal and vertical stress. Mohr-Coulomb model was used because its outcome was the best.

Numerical project setting
The beginning of work was naming the project, choosing the system of units, and project dimensions. The dimensions of the target area were 46,000 m of the y-axis, 26,000 m of the x-axis, and 2500 m of the z-axis.

Layers and materials
Eight layers (units) have been created: Mishrif U1, U2, Marl, M1, M2, L1.1, L1.2, and L2. The ten wells defined the tops of all the layers. Many layers (marl, M1, L1.1, and L1.2) disappear in some places because they essentially did not appear in specific wells, especially the wells at the edge of the reservoir like well G-5.
This model is based mainly on the information that was calculated and verified in the 1D model. That information included Young's modulus, Poisson's ratio, cohesion, internal friction angle. In addition, the bulk density was from the wells log. Scale-up is essential for this modeling, which defines the property of each grid penetrated by the wells. The data from 1D and well logs were scale-up that used as input to the model. Thus, all layers were fed with the essential information that was mentioned.

Generating wells
The representation of the wells is significant because they are the information source of layers. Moreover, the tops of the layers were entered through the wells. Then the map's Fig. 6 The layers and wells contour was defined based on those tops of the layers. In addition, one of the most critical factors that were input through the wells was pore pressure. Thus, ten wells were created; that are represented in Fig. 6.

Mesh configuration
This part is essential because it determines the number and size of elements and the number of nodes used in the calculations. The results will be more accurate when the size elements are small. There are five options very coarse, coarse, medium, fine and very fine. Thus, the mesh construction was the medium. In this work, the number of elements in layers was 51,438, and the number of nodes was 74,161, as shown in Fig. 7.

The calculations
The calculations are the last part of the model, including vertical, maximum and minimum horizontal stress before and after depleted.

Vertical stress ( SV = 1 )
This type of stress depends mainly on the overall density of the formations. Vertical stress magnitude is not affected by depletion of the reservoir. Figures 8, 9 , 10, 11, 1213 were showed the vertical stress magnitude in the reservoir was low. At the same time, its value was high in the area around the reservoir.

Maximum horizontal stress ( SH = 2 )
Knowing the magnitude of maximum horizontal stress is very important to know the value of rock breakdown. Figs. 14,15,16,17,18,19 show the maximum horizontal stress distribution and correspond precisely to the maximum horizontal stress of the wells from the 1D Model.

Minimum horizontal stress ( Sh = 3 )
The minimum horizontal stress was also calculated and compatible with the magnitude of the minimum horizontal stress of the 1D Model. As is known, minimum horizontal stress is perpendicular to maximum horizontal stress. Therefore, its orientation is toward the southeast, as shown in Figs. 20,21,22,23,24,25.

Depletion effect on horizontal stress
The horizontal stress magnitude significantly impacts the reservoir's characteristics, drilling operations, and production processes. These stresses are affected by many factors,     including pore pressure, geomechanical properties of formations and tectonic movements. The pore pressure is one of the essential factors that affect horizontal stress magnitude, especially in the stages of reservoir development. The phases of reservoir development go through depletion and possibly the injection phase to support reservoir pressure. Therefore,  the stress magnitude must be calculated at each stage, and the model should update to avoid problems and reduce costs.
Thus, formation pressure data were used for exploration and appraisal wells (before production) in 2011 as a first stage to determine the stress's magnitude. So, this pressure is  the initial pressure. The second stage used the pore pressure data after the pressure drop (depletion) in 2020 to determine the change in the stress's magnitude. The 4D model was used to calculate the difference in the horizontal stress magnitude with the pressure drop in the reservoir production units, including M1, L1.1, and L1.2. Only units produced are the target because of the depletion (pressure drop), whereas in other units, there is no change in pressure. The results showed that the difference in the magnitude of the horizontal stress was in the productive layer. The highest decrease in the maximum horizontal stress magnitude was in unit L1.1 and less than in L1.2, and the minor effect was in unit M1, as presented in Table 1 (2011) and after production (2020) (depletion time) that plotted against stress points (S.P).
Moreover, the magnitude of the minimum horizontal stress (Sh) was affected but less than the maximum horizontal stress, as shown in Table 2. Minimum horizontal stress is significant because it is related to the strength of the rock. Therefore, any decrease in minimum horizontal stress magnitude means that the unit is likely lost circulation. Thus, the unit more susceptible to the lost circulation is L1.1 because the Sh magnitude is the lowest. The average difference of magnitude minimum horizontal stress  (2020)

Conclusion
This study was started with data collected and audit; then a 1D model was built to estimate all rocks' mechanical properties and calibrate them. The outcomes of the 1D model were used to feed the 4D model, which considers necessary for calculating the stresses before and after the depleted reservoir. Therefore, the impact of the depleted reservoir was examined using a four-dimensional finite element model to investigate stress redistribution. The main goal is to present stress mapping before and after the depleted reservoir, considering geomechanical stress mechanisms and potential associated problems. Moreover, it could provide the knowledge and comprehension of the utilized geomechanics in depleted reservoir studies. The results showed many essential points that should be considered in the upcoming development wells of the depleted reservoirs.
• According to the outcomes of the 4D model, the Mishrif reservoir has been classified into four major depleted regions. • The highly depleted area includes WPJ, WPC, and WPF (G-J40, G-C43, and G-F23). • The moderately depleted area includes WPQ, WPK, and WPH (appraisal well G-4, G-Q45, G K117, and G-H74). • The slightly depleted area includes WPD and WPE (G-D85 and G-E35). • No depleted includes the appraisal well Ga-5 at the edge of the reservoir. • The units of the depleted reservoir were also classified into three units based on depletion, highly depleted unit (L1.1), moderately depleted unit (L1.2), slightly depleted unit (M1). • The minimum horizontal stress map shows two areas (highly and moderately depleted area) with low fracture pressure. So, the mud weight should be less than the fracture gradient to avoid lost circulation. • The operating mud weight window becomes narrow, dependent on the depleted area and depleted units of the Mishrif reservoir. Therefore, the weakest point in the Mishrif formation is in the highly depleted area and unit L1.1. • All operations, including wellbore structure, drilling, completion, stimulation, production, planning of the wells, should be redesigned based on the change and redistribution in stress.
• The wellbore stability analysis model must be updated according to the redistribution of the stress map, and the upper limit mud weight should not exceed in particular.