Application of improved data-driven diagnostics workflow based on coupled hydro-mechanical model for stimulation candidate selection in a complex carbonate reservoir

The focus of the present paper is on utilizing a comprehensive diagnostics workflow that combines coupled hydro-mechanical modeling with production data-driven diagnostics for optimization of stimulation candidate selection process. Reservoir fluids production and production-induced depletion affect reservoir mechanical environment due to reduced pore pressure and increased effective stress and the consequent porosity and permeability reduction, depending on rock sensitivity to stress changes. The coupled formation damage index (CFDI) is implemented in the traditional stimulation candidate selection workflow to capture the effect of production-induced stress changes on the near-wellbore permeability over time. The top potential stimulation candidates are recognized based on heterogeneity index, static formation damage index, and the CFDI parameter. CFDI is a dynamic parameter for stimulation candidate selection through estimating time-dependent permeability changes induced by stress state and fluid pressure. Probabilistic type curves and decline curve analysis of candidate wells versus reservoir unit are also applied to complement the stimulation and production enhancement candidate selection.


Introduction
Carbonate reservoirs contain more than 50% of world's conventional hydrocarbon reserves and on average have relatively low recovery factors of 20-30%, due to their inherent heterogeneity, complex pore system (matrix, fracture and vug), and wettability characteristics [1]. Due to a declining trend in the exploration of new super-giant fields and production decline of existing carbonate fields categorized as easy oil, better understanding of the productibility issues surrounding complex and tight carbonates such as permeability and fracture netwrok evolution has become increasingly important [2,3]. Also in recent years, exploration and development companies have gained more interest in development of deep carbonates with more complexities in reservoir rock and fluid behavior and application of improved oil recovery (IOR) techniques such as stimulation to tackle productivity issues [4]. Uncertainties in microfracture network distribution and fluid behavior heterogeneities are some of the complexity elements that impose challenges to the IOR candidate selection and make the use of simulation models as the main way of data integration limited [5]. Therefore, the importance of surveillance data-driven analytics integrated with the hydro-mechanical modeling for the IOR decision-making process such as formation damage removal is highlighted [4]. The surveillance data-driven analytics integrated with the hydro-mechanical modeling to identify the stimulation candidate selection in a complex carbonate oil field is discussed in this paper.
In its basic definition, near-wellbore formation damage is any unintended impedance to the flow of fluids into or out of a wellbore, caused by a reduction in permeability in the near-wellbore region, the changes in relative permeability to the hydrocarbon phase and the unintended flow restrictions in the completion itself [6]. The flow restrictions in tubing or imposed by well completion geometry are not included in this definition. The sources of near-wellbore formation damage can be internal or external sources. Internal sources include organic or inorganic deposits from the produced fluids such as scales, asphaltene, and fines. External sources include emulsions, precipitates or sludges caused by acid reactions, clay swelling or wettability changes caused by injected fluids and entrained particles such as solids in injected fluids [7]. All downhole operations including drilling, cementing, perforating, completion, work-over, gravel packing, production/injection, and even stimulation process can induce the nearwellbore formation damage [8]. In some cases of acidizing, stimulation can cause damage due to use of incompatible fluids, not properly mixed fluids, sludge formation, water blocking, wettability alterations or post-treatment fines migration. A more general definition, however, considers the formation damage in the rock matrix caused by the changes to the chemical-biological-mechanical-thermal environment during processes of enhanced oil/gas and geothermal recovery [9]. The reservoir-based enhanced oil recovery (EOR) methods such as CO 2 flooding, low salinity water (LSW), Alkali-surfactant-polymer (ASP) flooding, steam injection, and water injection can lead to reservoir formation damage.
One phenomenon that affects the mechanical environment of both conventional and unconventional reservoirs is the depletion due to production of oil and gas from the reservoir and consequently increased effective stress across the reservoir and compaction. The effective stress ( � ) , which is defined as the difference between total stress ( ), and pore pressure (p), is the average normal force per unit area transmitted from grain to grain of the rock. Figure 1 shows the relationship between total stress, effective stress, and pore pressure in a rock buried deep into the reservoir and acted upon by the overburden stresses.
Depending on how sensitive is the reservoir to the stress changes, compaction reduces the porosity and permeability of the reservoir formation, with permeability to be a more stress-sensitive property than porosity. The compaction-induced permeability reduction hinders the production of the well and can act as a near-wellbore formation damage [10]. Although it is difficult to experimentally predict the depletion-induced stress changes due to scale effects, stress arching, and faulting contributions in the field, it was observed that most of stress changes during production are irrecoverable on repressurization [11]. This observation highlights the importance of depletion rate management to reduce irrecoverable formation damage.
Various permeability relations from field and laboratory experimental results have been suggested for carbonate reservoirs, starting with the work of Dake [12]. Some models derive permeability relations via correlating a big data of the porosity and the permeability pairs for worldwide carbonate reservoirs classified in terms of the reservoir location and the rock types [13]. However, there exist stress-based models [14], in which a power law is proposed for the variation of the permeability with the effective stress. Various conditions were applied on permeability relations [15], variable stress conditions [14], or triaxial strain/ stress conditions [16]. Some models were based on complex distribution features of reservoir types [17], or were derived for carbonate reservoirs with extremely limited logging data [18].

Traditional stimulation candidate selection
The candidate selection is the most important step of stimulation process, as applying the best treatment design and field procedures to the wrong candidate will result in failure [19]. According to the traditional stimulation candidate selection workflows, the candidates are selected through a process of production profile analysis, nodal analysis, perforation and completion review to determine damage mechanisms and logs to evaluate mineralogy [8]. The damage mechanisms that contribute to total pressure where S total is total skin, S dam is skin due to formation damage, S c is completion skin due to sand control or partial penetration, S p is perforation skin due to crush zone and flow convergence, S dev is skin due to wellbore deviation that can be negative if perforated height is more than reservoir height, S stim is skin due to stimulation that can be zero by acid stimulation or negative for hydraulic fracturing, S non-darcy is rate-dependent skin that is prominent in gas well and S tp is skin due to two-phase flow, due to the competition of two flowing phases for the available flow channels [6]. The production-induced stress changes associated with depleting within reservoir and the consequent evolution of porosity, permeability and potential normal faulting are important to be studied in candidate selection process for further production enhancement and stimulation jobs. By using the coupled hydro-mechanical model that is described in the next section, one can estimate the total porosity loss in depleting reservoir and the associated permeability changes as a function of time.

Coupling hydro-mechanical model (HM) approach
In reservoir simulator approach, Darcy's law is solved for a case of single-phase, incompressible, isothermal flow, which in the absence of gravity can be written as Eq. (2).
with the porosity of the porous medium, the seepage velocity of the fluid in the porous medium, k the rock permeability, the fluid dynamic viscosity, and p the fluid pressure. In conventional simulation approach, rock permeability is assumed to be constant or depend only on fluid with uniaxial or isotropic strain for each mesh of the reservoir. In actual reservoir conditions, pressure, temperature and saturation variations generated by exploitation modify the stress distribution in the reservoir. Consequently, the porosity and compressibility of the rock vary during the reservoir production. According to Gutierrex and Lewis [20], the pressure dependence of porosity imposed via the rock compressibility introduced in reservoir simulators cannot alone describe the behavior of the reservoir rock during production and depends on the stress path followed by reservoir during its depletion. The stress path is the change in stress due to the change in pore pressure at a given point in the reservoir. It is defined as Δ h ∕Δp where Δ h is the change in horizontal stress and Δp is the change in pore pressure. If the depleting reservoir is infinite, i.e., the lateral length exceeds 10 times the reservoir thickness, using instantaneous application of force and pressure with no lateral strain, yields to a simple equation [Eq.
(3)] proposed to minimum horizontal stress in an isotropic basin [21]: where h is minimum horizontal stress, ν is Poisson's ratio, v is overburden stress, b is Biot's coefficient and P p is the pore pressure [22]. Taking the derivative of both sides of Eq. (3), and after simplification: In Eq. (4), it is assumed that overburden stress is not coupled with pore pressure, and therefore, the overburden stress vanishes during the derivation of Eq. (3) and Δ = 0 . Such an assumption was systematically studied in the way that minimum horizontal stress is fully coupled with pore pressure and is calibrated with hydraulic fracture stress data [23]. However, as researchers measured the effect of reservoir geometry and boundaries on the stress path [20,24], the necessity of using an empirically determined effective Poisson's ration which is calibrated against the least principal stress measurements obtained from leak-off tests is recognized. Therefore, the use of laboratory-defined Poisson's ratio for a reservoir in Eq. (4) is based on strong assumptions of uniaxial reservoir compaction, elastic behavior of the rock during depletion and homogenous drained Poisson's ratio at reservoir scale, identified to that measured in laboratory scale [10].
Poroelastic theory is applied in oil and gas industry mainly to understand subsidence, estimate reservoir productivity, and predict stresses around wellbores [25]. The carbonate reservoir is considered to be homogeneous. The elastic behavior of the reservoir is considered to be linear and isotropic elastic behavior. Only small strains are considered. The pore space is saturated with oil. Molecules in the pore space are assumed to be in a bulk state. Uniaxial strain behavior is considered for boundary conditions [11]. We consider a nonporous matrix in a representative elementary volume (REV) of the carbonate reservoir, i.e., the carbonate reservoir could be considered as a regular porous medium made of one pore network (i.e., the network of fractures). Therefore, the energy balance for the nonporous matrix in a REV of carbonate reservoir would be as Eq. (5) [26]: where f is defined as the Helmholtz free energy of the matrix per unit volume of reservoir. σ, s ij , and p are the volumetric stress, the deviatoric stresses, and the pore pressure, respectively. ϵ, e ij , and ϕ are the volumetric strain, the deviatoric strains, and the Lagrangian porosity, respectively. Energy can be added to such a system either by straining it with volumetric confining stress σ (the corresponding work being 'σdε'), by straining it with a deviatoric confining stress s ij (the corresponding work being 's ij de ij ') and by deforming the porosity ϕ with the pressure of the fluid in the matrix, p (the corresponding work being 'pdϕ'). Based on this energy balance, one can write the constitutive equations of the carbonate reservoir as [26]: where K, b, N, and G are poroleastic constants defined as the drained compression modulus, the Biot's coefficient, the Biot's modulus, and the shear modulus, respectively, all associated to the reservoir. Based on the derived constitutive equations, Eqs.
(2)-(4), we performed carbonate reservoir simulation. The performed simulations were asymmetric plane-strain two-dimensional. A Kozeny-Carman-type equation for permeability in the reservoir is considered, in which only the porosity intervened. Hence, the permeability k is given via the following relation [27]: where k 0 and 0 are the permeability and porosity in the state of reference and is the porosity at the actual state. Equation (9) shows that permeability (k) varies non-trivially with the porosity; with increase of permeability, porosity increases. The derivations performed in Eqs. (5-7 and 9) help to understand and predict how parameters such as porosity or permeability change for a representative elementary volume of the case carbonate reservoir saturated by oil. Moreover, a significant feature of our model is that twoway coupling between pore pressures and stresses or strains is captured: not only does the two-way behavior make it possible to predict how stresses or strains affect pore pressure, but also does it make it possible to predict how stresses or strains change in the presence of pore pressure. In other words, our model is able to capture the effect of stresses or strains on the permeability of the reservoir over time.

Applied stimulation candidate selection workflow
The proposed workflow in the current study is presented in Fig. 2, to screen the potential candidate wells for stimulation within a particular producing carbonate reservoir.
The presented diagnostics workflow is based on the monitored production data-driven analytics, well-test data acquired at the first day of production and coupled hydromechanical simulation accounting stress changes around the wellbore area. Figure 3 represents the schematic diagram of the established 2D coupled finite element model that used for investigating the effect of near-wellbore stress changes at perforation depth of 3500-4000 (m) for all the wells in the producing reservoir layer.
As described in Fig. 2, the heterogeneity index (HI) is a quick screening process to identify the preliminary candidate wells with anomalous behavior by comparing the relative production performance of each individual well to the performance of the field average [28]. The crosshair HI plot provides a view of individual well performances as better or worse than the field average in four quadrants. According to this plot, the candidates with oil and water production history lower than the average field performance proceed to the stimulation candidates' pool. The candidates with water production higher than the average field performance and oil production lower than the average field performance proceed to water shut-off (WSO) candidate selection workflow. Candidates with water cut higher than 50% may cause risk to the stimulation gain. Finally, the candidates with oil and water production history higher than the average field performance according to this plot are followed by continuous well performance monitoring. The formation damage indicator (FDI) plot is a scatter plot of reservoir capacity vs. the maximum oil rate for each well. The reservoir capacity is defined a permeability of the perforation zone obtained during well test multiplies the perforation interval (k 0 h). This scatter plot is to detect the underperforming well possibly due to skin, even though the formation has a high storage capacity to deliver. Thus, the reserves may remain trapped in a high percentage of the potentially productive zone. Considering the k 0 h values are not changing with time, it's a static plot.
Coupled formation damage index plot (CFDI) addresses the permeability (k) from Eq. (9) which is governed by porosity coupled with stresses or strains [Eqs. (5)(6)(7)(8)]. Thus, CFDI plot is a dynamic scatter plot that takes into account the stress changes during the reservoir production and the associated permeability changes of each well during its production time. According to the procedure developed by Zoback [10], by using Eq. (4) and setting b = 1 and ν = 0.25, one can find stress path as Δ h ∕Δp = 0.67. When the stress path is above 0.67, the production can induce normal faulting and the depletion-induced deformation should be checked. However, when the stress path is below 0.67, it would be a stable stress path.
Once the candidates have been screened based on HI plot, FDI and CFDI plots, the wells are analyzed individually to investigate the current status of the well and reservoir properties. Production performance of individual wells needs to be monitored to eliminate the inappropriate candidates such as the ones with high water cut percentage, i.e., above 70% [29]. If the produced sand concentration in a candidate well is higher than the field critical value, i.e., 15 lb per thousand barrels (pptb) [30], stimulation can hinder erosion risk and sand control options are the priority for the well, especially in gas wells. As the stimulation candidates are highly recommended if the reserves are sufficiently high enough to satisfy the fast payout time, the well decline curve is compared against the reservoir decline curve as the next step of screening process. If the reservoir pressure is available, average reservoir pressure decline is checked to be less than 5% to ensure the well decline is not due to the reservoir pressure decline. Before shortlisted as stimulation candidate, production performance for each candidate has to be investigated for the events such as shut-in well due to mechanical issue,

Field application and analysis
Parameters used for REV of the case carbonate reservoir are listed in Tables 1 and 2. A finite element model is used to perform the two-dimensional coupled hydro-mechanical model for a radius of 50 m of matrix rock around wellbore area at the mid-perforation depth of each individual well with a number of 40'000 grid blocks. A sensitivity analysis was performed on number of grid blocks from 1500 to 42'000, and as changes in results approached zero by addition of grid blocks above 40,000, the mesh-size of 40,000 for our modeling was selected. The Biot coefficient (b) of the reservoir is required to be assumed in the simulations. Using stress path criteria explained in previous section and by setting Biot coefficient b = 1 and Poisson's ratio ν = 0.3 from Table 1, the stress path would be 0.57. Since the stress path in this study (i.e., 0.57) is shallower that 0.67, the stress and pore pressure will evolve with flatter gradient and get further from the failure condition (i.e., 0.67), thus the depletion is stabilizing the reservoir and the faulting would not happen. The designed workflow as per Fig. 2 has been utilized in a complex carbonate field with microfracture network system for stimulation candidate selection process. A continuous hydro-mechanical modeling approach was selected for this case due to uncertainties in microfracture network that implies single-medium porous system. The single completion production wells penetrate multiple formation layers vertically, where only one continuous shale barrier separates the two main producing layers. The parameters for HI plot of oil and water production performance in y and x axis are presented in Eqs. (10) and (11), respectively.
Based on Eqs. (10) and (11), the production wells that fall into negative HISUMO and HISUMW regions during their production history, which in this case are W003 and  Initial porosity (%) Reported in Table 2 for each well k 0 Initial permeability (mD) Sp. Gr Specific gravity oil (API 0 ) 32.00 W011, are selected as the first candidate group in stimulation candidate pool. As can be seen in Fig. 4, W003 has been underperforming throughout its production history, while W011 has just recently been added in the negative region in terms of oil and water production. Figure 5a, b are presenting the scatter plot of k 0 h versus maximum oil rate of wells and kh map of wells across the reservoir unit of the studied case, respectively. As seen in Fig. 5a, although W012 has a high reservoir capacity based on the k 0 h value extracted from the initial welltest report, the maximum oil rate the well is comparatively lower than the other wells with similar formation capacity. Also as presented in Fig. 5b, the k 0 h values of well across the reservoir unit are categorized in three groups in terms of values. Apart from W009 that has shown a very high formation capacity at the initial welltest due to the microfracture network effect around the wellbore, the rests of the wells have almost same values with W012 and W006 showing  relatively higher capacity. However, the maximum oil rate production monitored in W012 is so much lower than the one observed in W006, which makes W012 a potential stimulation candidate to carry on to the candidate pool. Figures 6a-d are presenting the results of coupled HM model for the changes in kh and oil rate versus time for the four categories of wells as defined in HI plot. As can be seen in Fig. 6a, the changes in kh values for W011 is more highlighted than W003 in this group of wells as low productive wells. For the group of wells in low oil and highwater production as shown in Fig. 6b, W012 has presented higher kh changes during its production history. For the high productive group with high oil and high water production as seen in Fig. 6c, the kh change in W009 is nearly 1000 (mD.ft), which is high in three years of production with 3000 psi reduction in pressure. The comparison of kh changes in the high productive group wells with high oil production and low water production, defines changes in W005 is more highlighted than in W006, which makes W005 also a candidate in the stimulation pool.
Considering that according to Fig. 6b, all the wells are in the same range of kh, except W006, W009 and W012, the kh changes over time are plotted in Fig. 7 against the production time.
The grid maps of initial porosities and the flowing wellhead pressure of wells at the target reservoir unit are presented in Fig. 8a, b. The selected candidates at this step are considered to proceed to DCA analysis and comparison of each candidate well's DCA against the reservoir unit DCA to ascertain the reservoir depletion is not the reason for the well's underperforming. The DCA plots of each candidate are shown in Figs. 9, 10 and 11. The DCA plot of the target reservoir unit is illustrated in Fig. 12.
The candidate wells from the same reservoir are compared with the average reservoir decline to provide another indication of production impairment. Table 3  reports the decline ratio of the well to the reservoir for all three candidates, which is above 1.2 for all three candidates. Therefore, the stimulation candidate wells are recommended as the reserves are sufficiently high enough to ensure the job is economic.

Conclusions
A diagnostics workflow for stimulation candidate selection based on production data-driven analytics integrated with the coupled hydro-mechanical model has been successfully developed and implemented in a complex carbonate field case. The proposed screening workflow takes into account the dynamic behavior of formation damage index over the production history of each well, by utilizing the coupled hydro-mechanical model to enhance the screening process and improve opportunity detection for stimulation jobs. The novelty of the applied workflow is that wells which are experiencing severe changes in near-wellbore permeability during their production history are detected to be considered for stimulation jobs, along with underperforming wells in comparison to the filed average. The integration of dynamic parameters such as flowing wellhead pressure and static parameters such as porosity, through the coupled hydro-dynamic model, provides a more comprehensive and enhanced stimulation candidate selection for the complex carbonate field of study. Another outcome of the approach is the automaton of data diagnostics process for saving significant manhours for the engineers and equipping managers with smart decision-making tool.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 12
Decline curve analysis (DCA) for target reservoir unit