Development of a Fully Coupled Approach for Evaluation of Wellbore Stability in Hydrocarbon Reservoirs

Problem statement: When a well is drilled in a hydrocarbon reservoir, the original thermodynamic conditions are altered, the natural stresses are redistributed and a stress concentration occurs around the hole. The alteration of the original equilibrium can lead to wellbore damage, sometimes to its complete collapse. Loss of time associated with stability problems is estimated to account for 12-15% of drilling costs world-wide. Approach: The adoption of a reliable modeling approach to predict instability due to time-dependent alteration of natural equilibrium is fundamental for the optimization of drilling plans, completion design and production activities. Results: In this study the possibility of investigating instability phenomena in terms of both stress-strain and thermodynamic formation behavior through a fully coupled thermo-porous-elastic-plastic approach is demonstrated. According to the fully coupled approach, porous flow, temperature development and stress-strain calculations are performed together: the whole system is discretised on one grid domain and solved simultaneously for both the thermodynamic and the geomechanical variables. For the plastic analysis implementation, an iteratively coupled approach was adopted inside the fully coupled routine: the model basic equations (porous flow and rock deformation) and the plastic behavior equations were solved separately and sequentially at each non-linear iteration. The iterative coupling approach corresponds to an implicit treatment of the plastic variables, essential to preserve the stability of the elasto-plastic solution. The key points of the model analytical formulation, of the numerical formalization as well as of the implementation of the adopted solutions to make the thermo-porouselastic-plastic model applicable to assess wellbore stability are presented. Conclusion: The proposed model was first validated and then applied to several synthetic and real cases. In this study the effectiveness of the developed model to investigate the potential impact of instability phenomena on the well drilling design is demonstrated also by discussing the results from a case history.


INTRODUCTION
Wellbore stability related problems are reported to cost the petroleum industry more than 500 million US dollars each year [13] . In fact, alterations of natural equilibrium caused by drilling operations and extraction activities can lead to time-dependent wellbore instability phenomena, such as wall swelling and formation breakdown. Borehole instability also represents the main cause for loss of drilling fluids and consequent potential kick problems and can be so severe to compel the wellbore abandonment. In addition, sand production often associated with high production rate in poorly consolidated reservoirs can compromise oil production, increase completion costs and reduce the life of equipment.
Wellbore collapse and sand production are believed to be caused by large in-situ stresses and are the result of rock formation damage near a well. It is highly desirable to understand the main mechanisms for rock collapse and sand production and quantify these processes in order to optimize the drilling plan phase (i.e., mud characterization, drilling trajectory, rate of penetration) the design completion phase (casing placement and cementing operations) and the production activities.
In order to accurately analyze the complex behavior of many rock types encountered while drilling (i.e., shale minerals, salt formations), the effects of mechanical deformation, hydraulic diffusion and temperature gradient, as well as their combination, must be taken into account. However, as numerical modelling of coupled processes is extremely complex, it has been historically carried out in three separate areas: geomechanical modelling, which has the primary goal of computing stress/strain behaviour; formation simulation, which essentially models multiphase flow and heat transfer in porous media; and fracture mechanics, dealing with crack propagation and geometry. Usually, each of these disciplines makes simplifying assumptions about the part of the problem that is not of its own primary interest. However, such an approach is unacceptable in situations where there is a strong coupling, such as weak plastic reservoir rocks and unconsolidated porous media [2,11] . A coupled approach, combining constitutive equations with thermal and hydraulic laws, is then required to perform accurate time-dependent analyses of stress, pore pressure and temperature variation.
The technical literature shows several theoretic approaches to model the formation behavior with different degrees of coupling between rock deformation and fluid flow: from partially coupled to fully coupled methods. Although the couplings physically exist to some extent in all systems, the need for using more complex, fully coupled geomechanical modeling is generally acknowledged to be limited to the cases of compacting reservoirs or high-pressure injection operations or severe wellbore stability problem. Starting from Charlez's analytical thermo-porousmechanical model formulation [1] , the author of this study introduces the possibility of investigating the wellbore instability phenomena in terms of both the stress-strain and the thermo-dynamic formation behavior through a fully coupled approach. An iteratively coupled routine was adopted in the fully coupled approach in order to implement the plastic deformation analysis. The key points of both the numerical formulation and the implementation solutions developed by the author to make the thermoporo-elasto-plastic model applicable to assess wellbore stability are presented.

MATERIALS AND MATHODS
Different degrees of coupling: A multiphysical approach, able to simulate the correlations existing between rock geomechanics and reservoir fluid dynamics, requires the characterization of both the fluid and solid system components: A fluid-flow model describing the motion of the pore fluid and a stressdeformation model reproducing the deformation of the rock solid skeleton. Mass and force conservation laws and constitutive relations, which represent the coupling effects between components, are used to obtain the fluid-flow and stress-deformation model. The nonlinear behaviour typical of reservoir rocks can require the adoption of specific constitutive models, such as elastic-plastic law, visco-plastic law.
The methodologies of coupling flow through porous media and solid deformation that are found in the technical literature can be generally classified into fully coupled and partially coupled approaches.
In spite of recent advances in numerical calculation, the interaction between reservoir fluid flow and solid deformation still present a number of problems, such as accuracy, convergence issues and computing efficiency [12] .
The different coupling techniques are briefly presented in the following, highlighting the main advantages and shortcomings of each methodology [3] .
Fully coupled approach: In the fully coupled approach, porous flow and displacement calculations are performed together and the software solver must handle both the thermodynamic variables, i.e., temperature and pressure and the geomechanical response, such as displacement. The method is also called implicit coupling since the whole system is discretized in one grid domain and solved simultaneously. This approach has the advantage of internal consistency. The full coupling is also the most stable technique and preserves second order convergence of nonlinear iterations. However, the adoption of this methodology may complicate the coupling of existing flow simulators with geomechanical ones; it also requires more code development than other techniques and it can be slower than the explicit and iterative methods in some situations.

Iteratively coupled approach:
In the iteratively coupled approach, the basic equations for multiphase porous flow and rock deformation are solved separately and sequentially and the coupling terms are iterated at each time step. The exchange of information between the reservoir simulator and the geomechanic module is generally performed through a coupling code that also verifies the convergence of the coupling iterations. The adopted convergence criterion is typically based on pressure or stress changes between the last two solution iterations [12] . The adopted coupling variables are usually related to the key reservoir characteristics in order to highlight the most important coupling phenomena, such as volume changes, stress-dependent permeability, saturation-dependent rock strength.
The iterative coupling is straightforward when coupling existing reservoir simulators and stressanalysis codes. However, in this approach the calculation may display a first order convergence in the nonlinear iterations and, therefore, may require a large number of iterations, with a negative impact on computational time.

Explicitly coupled approach:
The explicit coupling is a special case of the iteratively coupled method: a simulator performs computations for multiphase porous flow at each time step and updates the geomechanical parameters during specified time steps, which are selected on the basis of pore volume change magnitude. If the pore volumes vary slowly during the time steps, then few geomechanical updates are required.
The main attraction of the explicitly coupled approach is the potential decrease of computational time, which is often mostly spent in calculating displacement for fluid flow/geomechanical problems. On the other hand, the explicit nature of the coupling can impose time step restriction on simulation runs in order to preserve stability and accuracy.

Model development for wellbore stability:
The presented wellbore stability model was developed on the basis of the analytical approach formulated by Charlez [1] and adapted to wellbore stability issues by Sartori1 [10] . It is based on a fully coupled thermoporous-elastic-plastic philosophy, which combines the constitutive equations with the thermal and hydraulic laws.
The proposed theoretical formulation is developed for a vertical borehole in a homogeneous, isotropic, elastoplastic, saturated porous medium subject to mono-dimensional, axisymmetric stress conditions in plane state of strain. According to Charlez [1] , the governing model equations can be divided into a basic group describing flow in the porous medium and rock deformation and a complementary group describing the plastic behavior of the formation.

Principal basic equations:
The basic equation group consists of: conservative laws (i.e.: Equilibrium, fluid mass and entropy), diffusion laws (i.e.: Darcy's law and Fourier's law), strain-displacement correlation law, stress-strain constitutive law and linearized behavior laws for pressure and temperature.
The constitutive law was developed in accordance to the hypothesis that small perturbations occurred, which implied that the total deformation is equal to the sum of an elastic increment and a plastic one: Where: σ = The stress tensor ε = The total strain tensor ε p = The plastic strain K B = The bulk modulus α = The Biot's coefficient P = The fluid pressure p 0 = The reference pressure α B = The drained thermal expansion coefficient of saturated rock T = The temperature T 0 = The reference temperature δ ij = The Kronecker tensor and C 0 reads: where, G B is the bulk shear modulus. Sartori proposed a peculiar formulation for the hydraulic diffusion equation and for the thermal diffusion equation which, compared to the analogous formulations of the thermoporous-elastic-plastic theory, contain an additional term typical of the plastic model. Combining Darcy's law with the mass conservation principle and the hydraulic diffusion equation [1] , the linearized pressure behaviour law was derived: Where: η = The Biot modulus L = The latent heat ρ 0 = The fluid density β = Plastic coefficient (analogous to Biot coefficient α) m = The fluid mass per unit volume I = The identity matrix Analogously, combining the Fourier's law with the heat conservation principle and with the thermal diffusion equation [1] , the linearized temperature behavior law was derived: Where: C ε = The heat capacity at constant volumetric deformation m 0 = The saturated rock density S = The entropy 0 m s = The specific fluid entropy at temperature T 0 Complementary equations: The complementary equations, which describe the plastic behavior of the formation, were taken as in the Desai's approach. Desai and coworkers [3][4][5][6][7] introduced a hierarchical approach for the development of constitutive models to account for various factors that influence the geological material behavior. This approach allows for progressive development of models of higher grades, corresponding to different levels of complexities. The basic model, termed the "δ 0 model", involves initially isotropic material, hardening isotropically with associative plasticity. In the proposed analytic formulation, the δ 0 model is used to describe the rock plastic deformation because, compared with models of higher grades, it provides a particularly accurate description of shale behavior and requires a smaller number of plastic input parameters. Belonging to classic incremental plasticity, Desai's model is based on the concepts of the yield function, f, on the plastic potential function, g and on the hardening function, α d , which are specific for the adopted plastic approach.
Desai proposed the following formulation for the yield function based on a polynomial expansion in terms of the direct invariants of the stress tensor (J 1 , J 2 , J 3 ): Where: σ ' e = The Biot's effective stress χ = The dimensionless trajectory of plastic strain p a = The atmospheric pressure α d = The dimensionless hardening function J e 1 = First invariant of the elastic stress n = The dimensionless shape change parameter γ d , β d , m = Dimensionless material response constants associated with the ultimate envelopes (i.e., the locus of points corresponding to asymptotic stress to stress-strain curves for different tests [7] ) S r = Expressed by the following equation: As the δ 0 model is based on associative plasticity, the plastic potential function g is the same as the yield function.
The adopted expression for the dimensionless hardening function reads: Where: a 1 , η 1 , a 2 , η 2 = Dimensionless hardening constants χ = The dimensionless trajectory of plastic strain defined as follows: and χ D is the deviatoric part of χ.
Solution implementation: Workflow: According to the solution workflow developed by the author, the physical phenomena induced by drilling operations and/or production activities into rock formations were simulated through two computational phases: the initialization and the working phase ( Fig. 1). During initialization, the plastic properties corresponding to the undisturbed formation were derived from complementary equations. The formation is assumed to be a homogeneous medium, subject to known and uniform temperature, pore pressure and stress distribution. As the complementary equation presents no spatial derivatives and so the system is spatially decoupled and thanks to the assumption of uniform thermodynamic and geomechanic original conditions, the initial phase is solved as a zerodimensional model. Absolute values of the initial stress tensor could be quite important, thus a gradual stress imposition approach could help minimizing convergence problems. Starting with a reference temperature and setting pore pressure and stress equal to zero, the known initial thermodynamic equilibrium and stress-strain field are reached incrementally and the corresponding plastic parameter values are calculated for each incremental step. This was achieved by solving iteratively the nonlinear complementary equation system.
During the working phase, spatial and temporal evolution of pressure, temperature, stress and strain was addressed based on the assigned time-variant boundary conditions. Axial symmetry and plane strain assumptions lead to a 1-D model. Spatial evolutions of the phenomena were simulated through macro-steps, each of which was divided into isochronal intervals (micro-steps). Boundary conditions were assumed to be constant or linear with time on the imposed macro-steps In the first case the increment was instantaneously assigned to the first step of the macro-step, while in the second case the increment was uniformly divided for each step composing the macro-step. This macro-step boundary value imposition allowed the simulation of the main phenomena connected with wellbore stability issues: drilling, completion and production conditions.   Starting from the imposed initial condition and original plastic parameter values derived by the initialization phase, the two nonlinear systems (i.e., basic and complementary ones) were solved separately and sequentially at each iteration during the working phase. This was achieved by the solution of two nested Newton-Raphson cycles at each time-step.

Numerical solution:
The author based the numerical formalization for the described multi-physical theory on the fully coupled approach, which allows porous flow and rock deformation calculations to be performed together. According to this implicit coupling technique the whole system is discretised on one grid domain and solved simultaneously for both the thermodynamic variables (temperature and pressure) and the geomechanical response (displacements). Moreover, in order to implement the plastic analysis, the author adopted an iteratively coupled approach inside the fully coupled routine. According to the iterative coupling technique, the model basic equations and the plastic behaviour equations were solved separately and sequentially at each nonlinear iteration. The iterative coupling approach corresponds to an implicit treatment of the plastic variables, essential to preserve the stability of the elasto-plastic solution.
Following this approach, the model could be arranged in two systems: F 1 containing the basic equations: Where: λ = The plastic multiplier χ = The hardening parameter tensor σ = The stress tensor ε = The strain tensor The coupling between basic and complementary equations arises from the plastic strain (ε p ), which appears in the elastic-plastic formulation of the hydraulic and thermal diffusivity equations (Eq. 3 and 4 respectively).
Basic equations: As the basic equations present both spatial and temporal derivatives of the problem unknowns, they were treated adopting the Finite Element Method approach and the Crank Nicolson time discretization, which guarantees unconditional stability. The model was solved adopting 1-D finite elements with a geometric progression grid in order to reach better accuracy in the neighborhood of the wellbore, where the variable's gradients are larger. The finite element selected for the discretization is the 1D element with 3 equidistant nodes. Quadratic Lagrangian polynomials N i , (with i = 1,2,3) were used as shape functions. Thus the corresponding finite element approximate solution is C 0 in pressure, displacement and temperature; continuity of strain between elements is not guaranteed but it is stepwise linear.
Concerning the finite element formulation of basic equations, from the strain-displacement equation and constitutive relation (1), the finite element formulation of the equilibrium equation becomes: Analogously, using the linearized pressure behavior law (3) and Darcy's law, the finite element formulation of mass balance becomes: Finally, adopting the linearized temperature behavior law (4) and the Fourier law, the finite element formulation of entropic balance is: where, B i is the vector: and N i the Lagrangian polynomial shape functions with nodal values as coefficients.
For each node a system of three equations (11,12,13) in three unknowns (∆u, ∆p, ∆T) could be written. Thus, the global finite element system in matrix formulation reads: Where: K = The basic coefficient matrix Δx = (Δu, Δp, ΔT) is the basic unknown vector Δb = Δb 0 +Δb 1 +Δb p is the known term which is composed by a boundary conditions contribution (Δb 0 ), an initial conditions contribution (Δb 1 ) and a contribution (Δb p ) due to the increment of the plastic strain tensor Δε p Note that Δb p expresses the coupling between basic and plastic equations; moreover it contains the nonlinearity of the system.
At each micro-time step (n) the basic nonlinear equation system (1) was solved with the Newton-Raphson method, adopting a convergence criterion based on the variability of the displacement increment between iterations. Moreover, at each iteration (i) of the basic equations solution, plastic deformation was determined in turn by solving iteratively the complementary Eq. 2, with temperature, pressure and stress values computed in the current basic system iteration (Fig. 1a).

Complementary equations:
The complementary equations have a spatially decoupled nature, thus they are point functions. As the coefficients of the basic equation matrix K (derived from (11,12,13)), contain integrals, which were solved numerically via quadrature formula, the complementary equations were conveniently computed directly in each Gaussian integration point.
For each time step (n), Fig. 1b shows the complementary equation solution algorithm on a generic integration Gaussian point and for a generic Newton iteration (i) of the basic system solution. First, strain was computed from displacement of the current outer iteration, then an iterative predictor-corrector technique was used to determine the stress increment. In the prediction phase elastic assumptions were made (null plastic strain increment and null plastic multiplier increment) and the elastic stress was computed from stress-strain relation (1). Then, the elastic assumption was verified by determining the yield function for the predictor elastic stress value. In case of plasticity, plastic stress was computed again via the nonlinear plastic system (2) solution. This was achieved by the application of the Newton-Raphson iterative method in which the Jacobian matrix was computed analytically in order to preserve the convergence properties. A convergence criterion based on plastic strain variability between iterations was adopted.

RESULTA AND DISCUSSION
Time-dependent wellbore stability analysis: The reliability of the proposed approach, under both a theoretical and a numerical point of view, was first validated on the basis of the results obtained for several synthetic and real cases with the FLAC software; afterwards, the methodology was tested on several representative case studies [8,9] . FLAC (Fast Lagrangian Analysis of Continua) is one of the most widely used geomechanical numerical simulators for performing geotechnical and mining engineering analyses worldwide. FLAC is a two-dimensional explicit finite difference program which allows execution of purely geomechanical as well as coupled thermodynamic and geomechanical analyses. The validation was performed within elastic domain analyses because the Desai's approach is not one of the plastic constitutive laws implemented into FLAC.
The proposed wellbore stability model was used to optimize a vertical well drilling design in a 2000 m oilbearing shaly-sand formation. The results of the performed sensitivity analyses allowed the identification of the best drilling approach which could assure borehole stability based on the particular petrophysical characterization and thermodynamic conditions of the investigated formation.
Two drilling scenarios were examined: An "overbalanced case" and an "underbalanced case".
According to the conventional drilling technique while drilling a well, the hydrostatic pressure in the circulating downhole fluid system is maintained within a pressure window ranging between the formation hydrostatic pressure (lower limit) and the formation fracture pressure (upper limit), thus higher than pore pressure ("overbalanced conditions"). The confining overpressure exerted by the drilling fluids on the borehole walls generally assures wellbore safety in terms of formation stability and potential kick phenomena. On the other hand, in some cases (such as in depleted reservoirs, low pressure formations and high-permeability porous media) borehole drilling according to this conventional technique can be subject to loss of fluids into the formation and severe formation damage (mainly, permeability reduction in the near wellbore zone) with consequent shortcomings in terms of productivity performances. In fact, the productivity of a well depends on the rock permeability, which is a measurement of the ease with which fluid can flow through a porous medium.
According to another drilling consolidated technique, the so called underbalanced drilling technique, the fluid system is maintained at a pressure level lower than the prevailing formation pressure at the target depth. The use of light-weight drilling fluids can reduce drilling time significantly, improve bit life and consequently cut down drilling costs. With a proper understanding of the reservoir behavior, the use of underbalanced drilling can also minimize formation damage, with consequent maximization of the natural potential well productivity. Added cost and complexity are the downsides of this application. Moreover, reducing the mud confining pressure, the underbalanced drilling technique can involve potentially higher wellbore instability risks.
It is therefore evident that an accurate wellbore stability analysis approach, which accounts for both the petrophysical characteristics and the thermodynamic conditions of the investigated formation, plays a crucial role in determining the final performance and the safety issues of any drilled well, in either an overbalanced or underbalanced mode. Field case-results: The developed stability analyses considered an "overbalanced case" characterized by a mud pressure applied to wellbore wall, p w , very close to formation pressure p o and an "underbalanced case" in which p w was set equal to 0.85 p o . The input data used to describe the thermodynamic, petrophysical and elastoplastic behavior of the shaly-sand formation are summarized in Table 1.
As far as the geometrical discretization is concerned, the spatial investigation domain extended from the borehole radius, set equal to 0.1 m, to the For the initial conditions the values of the natural in situ stress, pore pressure and temperature typical of a 2000 m deep formation, characterized by normal geostatic and geothermal gradients, were imposed. The stress domain was set isotropic on the horizontal plane and characterised by a horizontal vs. vertical stress ratio equal to 0.76. As far as the boundary conditions are concerned, the boundary values were imposed equal to the natural initial conditions at the external investigation radius, whereas the boundary values set at the wellbore wall reproduced the drilling thermodynamic and stress conditions. The solutions of the stability analyses were expressed in terms of spatial and temporal variations of displacement, pressure, temperature, deformation (radial, tangential and volumetric) and of both effective and total stress (radial, tangential and volumetric).  show that adopting an overbalance technique, the magnitude of radial and tangential stresses (and, obviously, of the relative radial and tangential strains) induced by drilling operations into the formation, was not so large as to compromise wellbore safety. The application of the Mohr-Coulomb yield criteria confirmed the stability of the investigated reservoir strata under overbalance drilling conditions. On the other hand, if the hydrostatic pressure in the circulating downhole fluid system during drilling was maintained at a pressure lower than the pressure of the target formation ("underbalanced case"), the stress increment induced into the rock would compromise wellbore safety. For the drilling approach which was identified as the best solution, i.e., the overbalanced case, the temporal evolution of phenomena was analysed not only immediately after the wellbore drilling (24 h), but also twenty days after drilling operations, when a new equilibrium was achieved into the formation. The results in Fig. 4 and 5 shows that the geomechanical and thermodynamic temporal variations are fully compatible with wellbore safety.

CONCLUSION
The adoption of a borehole instability modelling approach is fundamental when addressing technical and safety issues related to any drilling, completion and extraction design phase because it permits a systematic consideration of the stress variations around the borehole and the associated rock deformations and pore pressure changes.
The importance of adopting of a coupled thermodynamical/geomechanical theory for an accurate investigation of the complex behavior of many rock formations encountered while drilling is widely accepted and the definition of a numerical solution suitable for such complex coupled multi-physical theory implementation is anything but trivial.
In this study the numerical solution developed for implementing Charlez's thermo-porous-elastic-plastic theory, where plasticity is described by the basic Desai's model, is presented.
The adoption of a fully coupled approach allows porous flow, temperature development and rock deformation calculations to be performed together. The method is also called implicit coupling since the whole system is discretised on one grid domain and solved simultaneously for both the thermodynamic variables and the geomechanical response. This approach has the advantage of internal consistency and stability; moreover, it preserves second order convergence of nonlinear iterations. Furthermore, in order to implement the plastic analysis, an iteratively coupled approach was adopted inside the fully coupled routine. According to the iterative coupling technique, the model basic equations and the plastic behavior equations are solved separately and sequentially at each non-linear iteration. This was achieved by solving two nested Newton-Raphson cycles at each time-step. The iterative coupling approach corresponds to an implicit treatment of the plastic variables, essential to preserve the stability of the elastic-plastic solution.
Application of the proposed approach to synthetic and real cases showed the value and potential impact on economical and safety issues of wellbore drilling.