Numerical analysis of a caprock integrity during oil production by steam-assisted gravity drainage method

The work is devoted to the investigation of a caprock integrity during oil production by steam-assisted gravity drainage method. An originally proposed thermo-hydro-mechanical model was used for the evaluation of mechanical loading acting on the over-burden. The model includes mass conservation laws, the energy conservation law and the linear momentum balance. Filtration of each phase of the three-phase flow (steam, oil and water) is described by Darcy’s law. Effective stress concept is used to take into account the effect of pore pressure on the stress-strain state. Inelastic deformations are described by the phenomenological viscoplastic model based on Drucker-Prager yield criterion. Two qualitatively and quantitatively different scenarios of porosity evolution are considered. The first scenario corresponds to the case of the pore compression while the second one describes increase in porosity induced by the volumetric strains. The obtained stress-strain state in the reservoir was used for the assessment of the caprock integrity for both cases of porosity evolution. In addition, the effect of the porosity evolution on the oil production rate is studied. It has been shown that the oil production rate strongly depends on the prevailing physical mechanism of the porosity evolution.


INTRODUCTION
epletion of conventional oil resources requires the development of advanced techniques for crude oil recovery. Thermal methods are an effective way to reduce viscosity and enhance oil production rate. Steam-assisted gravity drainage (SAGD) method is characterized by a very high recovery factor (from 0.6 up to 0.8) [1]. This technique includes a number of horizontal parallel wells located one above the other. Steam injection into the upper wells D leads to the formation of a high-temperature steam chamber. Increase in the temperature reduces oil viscosity and enhances oil mobility. Oil together with the pore water and the condensed steam flow downwards into the production wells under the gravity. Optimization of SAGD requires formulation of mathematical models which takes into account multiphase flow, filtration, phase change processes and alteration of the stress-strain state due to the thermal and mechanical loading of the reservoir. Injection of the hot steam under the high pressure leads to the various geomechanical phenomena. For example, dilation and fracture significantly increase the porosity whereas compaction of soil grains and its thermal expansion reduces it [2]. Moreover, excessive stresses in the reservoir can break the integrity of the steam chamber. Loss of the caprock integrity is a serious problem for human safety and the environment which leads to the water and air pollution. Usually, numerical models of SAGD include momentum, mass and energy conservation laws which are supplemented by state and constitutive equations [3][4][5][6]. To describe the effect of fluid flow on a stress-strain state of porous media poroelastic Biot's theory is used. To provide a full coupling between fluid flow and the stress-strain state of the reservoir it is necessary to include a constitutive equation for porosity evolution into the model. The porous media is a complicated system which deformational behavior and properties are strongly depend on the level of the applied load, initial heterogeneity, anisotropy and many other factors. Therefore, there is a wide variety in the choice of relations for the description of the porosity evolution. In the simplest case, the porosity can be related to the pore pressure using the pore compressibility coefficient [7][8]. In [9][10] the incremental relation between the porosity and the mean stress with the coefficient which is linearly dependent on the current porosity value has been obtained. In case of the low compressibility, the porosity can be associated with the mean effective stresses using the composition of the linear fractional function and the exponent [11]. Another exponential dependence of the porosity on the mean effective stresses has been obtained in [12]. The incremental equation relating porosity to the linear combination of the volumetric strains and pore pressure has been proposed in [13][14]. In case of the low compressibility the dependence between the porosity and volumetric strains is also expressed as the exponential function [15][16]. At the same time, a linear fractional dependence of the porosity on volumetric strains has been proposed in [17][18][19]. As regards to the viscous oil production, the main mechanisms affecting the value of porosity is the shear dilation. In [1] it has been mentioned that propagation of the thermal front within the reservoir leads to the rise in the horizontal compressive stresses and substantial reduction of the vertical stresses. Such changes cause a significant shear dilation which enhances permeability and increases porosity of the soil. This effect can be described by the models which relate porosity to volumetric strains [4], [15][16][17][18][19]. However, when the effective confining pressures are high the shear dilation can be suppressed. In this case, the pore compression leads to the squeezing of the oil out of the pore space which can be described by relating of the porosity to the effective stresses [9][10][11][12]. Alteration of porosity during thermal oil production can significantly affect mechanical state of the reservoir near the caprock. Analysis of the existing works have shown that only few works are devoted to the study of the caprock integrity during SAGD. Rahmati et al. [20] investigated effect of the intrinsic anisotropy of caprock on maximum operating pressure of SAGD. For this purpose, they combined anisotropic Mohr-Coloumb failure criterion with transversely isotropic elastic model. They have shown that isotropic model gives a little bit higher values of maximum operating pressure than anisotropic. McLellan et al. [21] took into account statistical distribution of the natural fractures in caprock. Khan et al. [22] proposed an integrated geomechanics approach to evaluate caprock integrity which is based on the construction of the real three-dimensional formation taking into account petrophysical properties of the reservoir, strength parameters of the caprock and the results of mini-frac or leak-off tests. The obtained geological model is integrated with the coupled reservoir simulator which allows one to obtain stress-strain state of the considered formation. These results can be used for estimation of the caprock failure on the base of Mohr-Coulomb fracture criterion. Walters et al. [23] confirmed that the key factors determining maximum operating pressure in SAGD are the knowledge of the initial stress-strain state of the reservoir and its material properties. They have conducted analytical and numerical analysis of the stress change during SAGD. Analytical method was based on the theory of poro-thermo-elasticity and the Mohr-Circle analysis for evaluation of the shear failure. Numerical assessment included a coupled reservoir and geomechanical simulator. Tensile and shear stress ratios were used as failure indicators. Uwiera-Gartner et al. [24] have analyzed caprock integrity of SAGD project in northeast Alberta (Canada) using geological and geomechanical parameters. The drawback of this method is existence of the uncertainties in geomechanical parameters of the caprock. The key uncertainties affecting the analysis are heterogeneity, anisotropy, delamination and large fractures, which can reduce rock strength. This work is devoted to the investigation of SAGD operating parameters (the temperature and the pressure of the injected steam) on the caprock integrity near its base. The following algorithm has been used in order to estimate the integrity. Firstly, the coupled thermo-hydro-mechanical problem is solved to obtain mechanical state of the reservoir. An originally proposed model is applied to this purpose. According to the model, the process of steam injection is governed by mass conservation law, energy conservation law, Darcy's law for oil, steam and water filtration. Water-vapor phase change is determined by the additional heat source proposed in [25]. Inelastic deformations of the reservoir induced by propagation of the thermal front are described by the phenomenological viscoplastic model. On the second step of the algorithm, the values of mechanical stresses near the top of the reservoir are determined. The last step is the assessment of the caprock failure according to Drucker-Prager criterion. The specific feature of this work is that two qualitatively different scenarios of porosity evolution is considered. The first scenario corresponds to the case of the pore compression which can be described by the model [10]. The second scenario demonstrates increase in porosity induced by the volumetric strains. This effect can be described by the model [11] which is widely used in SAGD simulation.

MATHEMATICAL MODEL
he developed model of SAGD is based on the following assumptions. Fluid within the pore space consists of the three immiscible components (water, steam and oil). The flow of each component is slow and obeyed to Darcy's law. Capillary pressure is negligible due to relatively high porosity. The change in the oil viscosity is the result of the change in the temperature only. The deformations of porous media are small. Elastic behavior of solid skeleton is described by Hook's law. Effective stress concept is used to take into account the effect of pore pressure on the stressstrain state. Viscoplastic behavior of the reservoir is described by Perzyna's theory and Drucker-Prager yield criterion [26][27]. Therefore, the complete system of equations can be written as: where subscript i takes values s , w , o , r and denotes steam, water, oil and solid skeleton respectively; n is the porosity;  i is the density; t is the time; w S is the water saturation; s S is the steam saturation; o S is the oil saturation; i v is Darcy's velocity; K is the absolute permeability; ri k is the relative permeability;  i is the dynamic viscosity; p is the pressure; g is the gravity; T is the temperature; i c is the heat capacity; Q Lq is the heat source due to the phase change; L is the latent heat; σ is the stress tensor; n is the effective density; C is the stiffness tensor which has two components in case of linear isotropic elasticity (Young's modulus and Poisson's ratio); ε is the strain tensor;  T is the thermal expansion coefficient; J is the second invariant of the deviatoric part of the stress tensor; 1 I is the first invariant of the stress tensor; a , b are material parameters; u is the displacement vector; r is the mass transfer intensity factor; sat T is the phase change temperature. The following functions were used to define relative phase permeabilities [28]: where 1 a , 2 a , 3 a , 1 m , 2 m , 3 m are the empirical parameters; rw S , rs S , ro S are the residual values of water steam and oil saturations. In this work the effect of volumetric strains and effective stresses on the porosity evolution of the reservoir is investigated. For this purpose two qualitatively and quantitatively different models are considered. The first model [10] relates porosity to mean effective stresses by means of the coefficient which is linearly dependent on the current value of porosity: where p c is the pore compressibility; r c is the rock compressibility; ;  xx ,  yy ,  zz are the diagonal stress tensor components. This model can be applied in case of the high confining pressure when the effect of shear dilation is small. The second model [11] uses exponential dependence of volumetric strain on porosity: where 0 n is the initial porosity;  vol is the volumetric strain. This model can be applied to the cases when the prevailing mechanism of deformation is the increase in volumetric strains. The following boundary and initial conditions were used to close the developed system of equations for SAGD simulation: where 0 p is the initial pressure; 0 w S is the initial water saturation;  1 denotes the boundary of the injection well; b p is the injection pressure;  2 denotes the boundary of the production well; w p is the pressure in the production well; n is the unit normal vector; ' q is the heat flux vector;  3 ,  4 denote left and right boundaries of the rectangular domain;  5 ,  6 denote upper and lower boundary of the rectangular domain.

RESULTS AND DISCUSSIONS
he system of Eqns. (1)-(18) together with the initial and boundary conditions (19)-(32) was solved numerically using finite-element software Comsol Multiphysics®. For this purpose, the following algorithm has been proposed. Firstly, Eqns. (1)-(7) were reformulated using the total velocity of the three-phase flow. After this, the obtained equations were represented in the weak form. Each of the equations was multiplied by the test function. Then, the integration by parts was applied to reduce the order of differentiation. In the considered equations, the convective term prevails over the diffusive one. To smooth the oscillations caused by the convective terms the artificial diffusion was added to each equation. The resulting equations were implemented to Comsol Multiphysics® by Weak Form PDE interface. To solve the Eqns. (8)-(12) Heat Transfer and Structural Mechanics modules were used. The more detailed description of the algorithm can be found in [29]. The developed numerical algorithm has been applied to the three-dimensional numerical simulation of the rectangular reservoir containing two horizontal wells (injection and production). Schematic representation of the simulated domain is given in fig. 1 (a). The thickness of the considered area is 1 m. Dimensions in the X and Z directions are equal to 30 m and 24 m. The wells have diameter of 0.178 m. The production well is located at a distance of 10 m from the bottom of the reservoir. The distance between the wells is 5 m. The considered area has been discretized by triangular prismatic finite elements with the minimal element size of 0.03 m (near the wells) and the maximum element size of 1 m (near the boundary of the simulated domain). The reservoir and oil properties correspond to Yarega oil deposit (Russian crude oil deposit in Komi Republic). Physical and mechanical properties of water, steam, oil and reservoir are given in Tabs. 1-2. Fig. 1 (b) describes the temperature dependence of oil dynamic viscosity used in simulation. The main focus of this work is the assessment of the caprock integrity during oil production by SAGD with two sets of operational parameters: pb=2.5 MPa, Tb=213 K and pb=3.5 MPa, Tb=220 K. Values of other initial and boundary conditions are given in Tab. 3. The loads acting on a caprock are fully determined by the stress-strain state of the reservoir which depends on the structural changes induced by propagation of the steam chamber.

Property
Value Unit   Fig. 2 shows the form of the steam chamber after 100 days of the heating with two different regimes. Steam injection with operational parameters p b =2.5 MPa, T b =213 K results in the vertical growth of the steam chamber while operational parameters p b =3.5 MPa, T b =220 K give a completely different shape with substantial sideway growth. It is noticeable that fig. 2(a) corresponds to the early stage of the steam chamber growth whereas fig. 2(b) illustrates the spreading phase. Therefore, the second regime of the heating is characterized by the very high values of the steam propagation within the reservoir.  Pressure in the production well 1.7·10 6 Pa Table 3: Initial and boundary conditions.
(a) (b) Figure 1: (a) Schematic representation of the simulated domain (  1 is the boundary of the injection well,  2 is the boundary of the production well,  3 is the left boundary,  4 is the right boundary,  5 is the upper boundary,  6 is the lower boundary). All sizes are in mm. (b) Dependence of the dynamic oil viscosity on the temperature [30]. The effect of volumetric strains and effective stresses on the porosity evolution is shown in fig. 3. It can be seen that the considered models (17) and (18) give a qualitatively and quantitatively different distribution of the porosity. The porosity Eqn. (17) describes compression of the soil in the process of SAGD ( fig. 3(a)). Pore compression leads to the decrease in pore pressure and increase in effective compressive stresses in such a way that all domain is in the compressive state. The highest contraction is observed within the area which corresponds to the spreading steam chamber. By contrast, the model (18) demonstrates pore expansion due to the influence of the pore pressure and temperature. According to this model, the porosity changes only in the heating zone and does not affect the rest of the domain. Thus, the first model is more sensitive to the applied mechanical and thermal loading than the second one. To choose the more appropriate model of the porosity for the specific reservoir it is necessary to determine the prevailing mechanism of its evolution because the obtained results can significantly affect the oil production rate. The effect of porosity evolution on the oil production rate is studied using the Kozeny-Carman model which is widely applied to the problems of petroleum engineering. The model relates absolute permeability to the porosity and can be written in the following form [31]: where 0 K is the initial value of an absolute permeability; 0 n is the initial value of porosity. Fig. 4 shows specific oil production rate obtained by the proposed SAGD model after the establishment of the interwell communication using models (17)- (18) and two abovementioned sets of operational parameters. The first case of the operational parameters gives an upward trend of the production rate while the second case shows a stable production rate corresponding to the spreading phase of the steam chamber. As the one would expect, the highest value of the production rate predicts the Eqn. (18) and the second set of operational parameters. The first set of the parameters for the same porosity model gives almost threefold smaller rate. The Eqn. (17) predicts moderate value of the oil production rate due to the pore compression. In this case, the second set of the parameters are also preferable to the first. The high values of injection pressure and temperature result in a faster development of the steam chamber, so the more oil becomes mobile. However, there is a possibility to change the structure of the reservoir in such a way that there will be little pore space available for the oil filtration and the reservoir can become almost impermeable. Thus, if the prevailing physical mechanism of the porosity evolution for the specific formation is the pore compaction, the increase in the values of the operational parameter can lead to the drop in the oil production rate instead of its rise.  The main criteria for the optimal choice of the operational parameters is to ensure the safe conditions which will not lead to the caprock failure. Loss of a caprock integrity is the serious problem for a human safety and the environment. The overburden provides isolation of the injected material from penetration of it into the water. Fig. 5 shows assessment of the caprock integrity according to the Drucker-Prager fracture criterion for two sets of operational parameters. The strength parameters a and b are equal to 0.19 and 243623 Pa. The following cases have been considered. The first and the second cases correspond to the solution to the thermo-poroelastic problem (neglecting of the visco-plastic strains defined by Eqn. (11)) obtained for both sets of operational parameters ( fig. 5, (a)-(b)). In all other cases, visco-plastic strains were taken into account. The second and the third case describe porosity evolution according to the Eqn. (17) (fig. 5, (c)-(d)). As regards to the fourth and the fifth, the Eqn. (18) was used ( fig. 5 (e)-(f)). The obtained results show that neglecting of the plastic strains in the reservoir leads to the high values of mechanical stresses which induce caprock failure for both sets of operational parameters ( fig. 5, (a)-(b)). On the contrary, the more sophisticated model which takes into account plastic strains in the reservoir predicts integrity of the caprock for all considered cases ( fig. 5, (c)-(f)). Therefore, neglecting of the inelastic strains leads to the overestimation of the caprock strength even for the small values of the operational parameters. The potentially dangerous zones are located over the steam chamber which produces a Gaussian-type loading acting on the base of the caprock. It is interesting to note that higher values of operational parameters predict safer condition than lower in case of the porosity model (17). This effect can be attributed to the Drucker-Prager model which is a pressure-dependent fracture criterion. According to this criterion, stress state with the lower value of the hydrostatic pressure is closer to the fracture surface than the state with the higher value under the same value of the equivalent (von Mises) stress.

CONCLUSION
n this study, an assessment of the caprock integrity during oil production by SAGD is carried out. The evaluation of the loadings acting on a base of the caprock was calculated with the use of the originally proposed thermo-hydromechanical model of steam injection into the porous media. The model describes convective heat transfer, phase change process, thermal expansion of solid grains and influence of the pore pressure on the stress-strain state of the reservoir. The process of steam injection is governed by mass conservation law, energy conservation law, Darcy's law for oil, steam and water filtration. Water-vapor phase change is determined by the additional heat source. Inelastic deformations of the reservoir induced by propagation of the thermal front are described by the phenomenological viscoplastic model. Two qualitatively and quantitatively different scenarios of porosity evolution is considered. The first I scenario corresponds to the case of the pore compression while the second one describes increase in porosity induced by the volumetric strains. The following conclusions can be made:  Operational parameters p b =3.6 MPa, T b =220 K substantially accelerate steam propagation within the reservoir compared to the values p b =2.5 MPa, T b =213 K.
 The oil production rate strongly depends on the prevailing physical mechanism of the porosity evolution. Pore compression impairs oil production. The operational parameters pb=2.5 MPa, Tb=213 K give almost threefold smaller oil production rate for both porosity models.  The considered values of operational parameters do not induce caprock failure if plastic strains in the reservoir is taking into account. Therefore, neglecting of the inelastic strains leads to the overestimation of the caprock strength even for the small values of the operational parameters.