A Numerical Model for Caprock Analysis for Subsurface Gas Storage Applications

In considering a site for gas storage, it will be important to evaluate the effects of gas storage on the formation, so as to minimize the risk of a breach occurring in the system. Gas injection will result in an increase in formation fluid pressure, especially around the injection source, which in turn results in redistribution of the stress field. The induced deformations within the reservoir can potentially result in a damage zone within the caprock formation. This mechanical failure may involve shear along many of the existing fractures or creation of new fractures that reduce the sealing properties of the caprock system. The main objective of this paper is to develop a model to estimate the growth and extension of cracks in the caprock. In order to achieve this, the smeared crack approach is used to model the process of cracking in the caprock. Smeared cracking is a continuum approach for damage mechanics which is based on the idea that a crack is modeled by modifying the strength and stiffness of the material. The main model presented in this paper has three sub-models, which are the reservoir model, the caprock model and the smeared crack model. The reservoir model is a simplified coupled hydro-mechanical model that numerically simulates the radial fluid flow and analytically estimates the associated stress and strain within the reservoir. The results of the reservoir model are used as boundary conditions for the caprock model that estimates the stress and strain within the sealing caprock due to the deformation of the reservoir. Using the calculated stress and strain, the smeared crack model predicts the growth and extension of cracks within the caprock. The caprock is assumed to be initially crack free and impermeable. The developed model is then used to study the Yort-e-shah aquifer caprock in Iran to predict the growth and extension of cracks. Keywords-caprock integrity; smeared crack; reservoir geomechanic;Yort-e-shah aquifer


INTRODUCTION
Changing the pore fluid pressure and temperature in a reservoir will result in the generation of mechanical stresses in the vicinity of the reservoir [1,2].The potential exists for these stresses to induce a failure in the reservoir's bounding seal.These failure events may result in relatively permeable flow paths that enable leakage of the reservoir fluids into the surrounding geological formations, and potentially to ground surface or into shallow aquifers.As such, when evaluating the suitability of a reservoir for gas storage, it is important to recognize the types of failure mechanisms that may occur, and their likelihood.Authors in [3] have categorized typical geomechanical risks that may lead to leakage by influencing the hydraulic integrity of caprock during and after gas injection in three different sections: storage-induced, storage-activated and tectonic activity.In their research, storage-induced section includes capillary leakage, hydraulic fracturing and shear fracturing.Storage-activated section includes fault reactivation and reactivation of pre-existing fractures.Tectonically active regions are considered separately, and may be avoided by locating the project site in non-seismic regions.There are a great number of studies investigating the caprock integrity from different perspectives such as rock-fluid interactions [4,5], geomechanical failure [6,7], induced seismicity [8] and fault slip [9].A summary of the research activities in the area of caprock integrity has been provided in a number of review papers [10][11][12].Among all parameters influencing the hydromechanical integrity of caprocks, a number of studies have highlighted the importance of crack growth and extension on sealing performance of the caprock [10].Despite a long history of studies on crack initiation and propagation in rock mechanics, especially in the area of hydraulic fracturing, it remains a challenging topic where there is no consensus on an optimized unified approach and different mechanisms and methodologies have been employed [13][14][15][16][17].
To analyze the initiation and propagation of cracks in several structural problems, numerical studies such as finite element and boundary element methods have been used.There are generally three approaches in solid mechanics when modeling cracks: discrete crack, interface elements and smeared crack.Although, discrete crack approach reflects the fracture development phenomenon most closely, it does not fit the nature of numerical methods and it can be computationally expensive [18].With the development of cracks, each node is replaced by more nodes, which entails re-definition of finite element mesh and hence augments the need for computational resources.On the other hand, the smeared crack approach assumes that the cracked solid is a continuum and permits the description of the medium in terms of conventional stress-strain equations.Although this may conflict with the nature of the cracking phenomenon, it simplifies significantly the computational process.The smeared crack approach can model the crack development in any direction and it is independent of finite element discretization.Within the smeared crack framework, small cracks formed in the band are gradually connected to each other and also one or more cracks may be defined for each gauss point of each element.There are a number of comparative studies investigating the advantages and disadvantages of different approaches.In this paper, we employ the smeared crack approach to avoid the need for remeshing the finite elements, and to enhance the computational efficiency of the model.The idea of smeared cracks was first introduced by in [19] and was later named smeared crack approach in [20].
The results of this study are advantageous to many projects involving gas injection into subsurface formations, namely CO 2 sequestration, enhanced oil and gas recovery.In order to numerically study the crack development in an initially intact caprock following gas injection in an aquifer, an integrated reservoir-caprock model was developed.The reservoir modeling is a simplified single phase flow model coded in FORTRAN which also estimates the geomechanical response of the reservoir rock to gas injection.The spatial and temporal deformation distribution of the reservoir is then used as a boundary condition for the bottom of the caprock model constructed in ABAQUS.The model includes the elastoplastic behavior which is coupled with a smeared crack model.The smeared crack model provides an efficient tool to estimate the density and orientation of developed crack in all elements of the caprock.The developed model is then used to simulate the crack development within the caprock of Yort-e-Shah aquifer.The results of this study can have a significant impact on improving the current methodologies for caprock integrity assessments in geological storage activities.

A. Elasto-plasticity in Caprock
In this paper, cracking is assumed to be the most important aspect of caprock behavior, and it dominates the modeling.Cracking is assumed to occur when the stresses reach a failure surface, which can be named as crack detection surface.This failure surface is taken to be a simple Coulomb line written in terms of the first and second stress invariants, p and q.The smeared crack model does not track individual macro-cracks but rather constitutive calculations are performed independently at each integration point of the finite element model.The presence of cracks enters into these calculations by the way the cracks affect the stress and material stiffness associated with the integration point.When the principal stress components are dominantly compressive, the response of the caprock is modeled by an elastic-plastic theory, using a simple form of yield surface written in terms of the first two stress invariants.Associated flow and isotropic hardening are used.This model significantly simplifies the cracking and compression responses of caprock incorporated in the model.
When caprock is loaded in compression, it initially exhibits elastic response.As stress increases, some inelastic straining occurs, and the response of the material softens.An ultimate stress is reached, after which the material softens until it can no longer carry any stress.Note that the hysteresis of loading and unloading is not considered in this paper.When a uniaxial specimen is loaded into tension, it responds elastically until, at a stress that is typically 7-10% of the ultimate compressive stress, cracks are initiated very quickly.In order to model these phenomena, we employ a softening mechanism by which the developments of cracks results in a damage effect.In this framework, the open cracks are represented by the loss of elastic stiffness.In other words, we neglect permanent strains associated with cracking.In multi-axial stress states these phenomena can be generalized through the concept of surfaces of failure and of ultimate strength in stress space, which are defined by fitting to the experimental data.
The crack model used in this paper is the built-in smeared crack model of ABAQUS that is originally used for the cracking of concrete.A brief summary of the mathematical implementation of the model is given in this section, while further theoretical background of smeared crack model can be found in the literature [21][22][23][24][25][26].This model considers a "compressive" yield/flow surface to model the caprock response in compressive stress states, along with damage elasticity to represent cracks that have occurred at a material calculation point.The occurrence of cracks is defined by a "crack detection" failure surface.Following the classical theory of plasticity, the mechanical strain rate dε can be decomposed into elastic and plastic components as below: where el d  is the elastic strain rate of un-cracked rock which can be estimated from classical elasticity theory.
pl c d  is the plastic strain rate associated with the "compression" surface.
The "compression" surface, f c is: where p is the effective compressional stress, defined as: and q is the Mises equivalent deviatoric stress: is the deviatoric stress, The colon (operator ":") is applied for matrix multiplication, 0 a is a constant defined using the ratio of triaxial peak strength, ) as below:

www.etasr.com Rajabi et al.: A Numerical Model for Caprock Analysis for Subsurface Gas Storage Applications
Also ( )   is a hardening parameter, where c  is the size of the yield surface on the q -axis at 0  p .Note that surface is a straight line in p-q space and provides a good match to experimental data over a fairly wide range of pressure stress values [25].
where σ bc is the magnitude of each nonzero principal stress.The model uses associated flow rule in where the plastic strain rate is defined as below: Note that, c 0 is a constant defined as: When the stress state is tensile, cracking is the dominant form of failure, and hence the "crack detection" surface in stress space is used to estimate the spatial distribution of cracks.Damage elasticity is then used to describe the postfailure behavior of the rock.The behavior of the cracked rock is defined using a brittle fracture concept derived from [22].Thus, the elastic strain rate in (1), Considering an associated flow rule, the plastic components of strain in cracked rock can be written as: In which f t is a Coulomb-type tensile failure envelope as: where t  and u t  are the tensile stress and uniaxial tensile strength, respectively.Also, b 0 is a constant defined as: where . Also the associated plastic multiplier (hardening parameter) is defined as: Following crack detection, the damage elastic theory is used to relate the stress and strain in the failed rock: where D is the stiffness matrix of the rock, in which the effect of Poisson's ratio is neglected.It should be noted that cracking is assumed to be irrecoverable and that no more than three cracks can occur at any given point.

B. Reservoir Model
The reservoir model is a simplified coupled hydromechanical model that numerically simulates the radial fluid flow and analytically estimates the associated stress and strain within the reservoir [27].Gas is injected into an axisymmetric disc-shaped reservoir through a vertical wellbore at its center.The reservoir is treated as a homogeneous and deformable system whose porosity is filled with a single phase, single component and slightly compressible gas.It is a rather simplified model, but this paper is focused on the applicability of smeared cracking approach for caprock integrity and hence, using a simple reservoir model would be sufficient.Following the approach employed in [27,28], the coupled flow equation can be expressed as below: where t is time, r f is the radial distance from the wellbore, P is the gas pressure within the reservoir, K f is the intrinsic permeability, f  is the gas viscosity, f  is the porosity of the reservoir, c f is the compressibility of the fluid, 0 f  is the density of gas at reference pressure of 0 P ,  is the Biot's coefficient, and v  is the volumetric strain of the reservoir that is determined from the elastic constitutive relations.The term Q is the mass rate of gas diffusing into and/or out of matrix blocks, i.e. matrix-fracture transfer function, defined as below: where q ave is the mean value of the gas concentration inside the reservoir matrix and is determined by integration over the volume of the matrix block as below: where V m is the volume of matrix block and q is the gas concentration which is a function of time, t, global fracture coordinate, X f , and local matrix coordinate X m .

www.etasr.com Rajabi et al.: A Numerical Model for Caprock Analysis for Subsurface Gas Storage Applications
Assuming an isotropic poroelastic reservoir, the incremental change in horizontal effective stress becomes: Where υ is the Poisson's ratio, σ΄ is the effective stress and subscripts of r and θ represent the radial and tangential directions in global coordinates, respectively.The overbars in this equation and throughout this paper denote incremental values of parameters with respect to time.Also, note that the Biot's definition of effective stress is employed.The assumption of uniaxial vertical deformation has been widely used in reservoir geomechanics studies as a helpful simplification for reservoir and geomechanical coupling.Using the classical elasticity theory, the incremental strain can be determined as: where E is the Young's modulus of reservoir, m zz  and s zz  are the mechanical and chemical strains in vertical direction, respectively.Also, assuming that the displacement is zero at the bottom of the reservoir, the displacement at the top of the reservoir, z u can be calculated as below: where h r is the thickness of reservoir.A backward implicit finite difference scheme was applied to the flow equation ( 15), while the incremental mechanical stresses and deformations are calculated at each time step.Further discussions on features of the employed coupled hydro-mechanical model of the reservoir can be found in [27][28][29].

C. Computational Workflow
As explained earlier, the integrated reservoir-caprock model used in this paper is constructed by coupling a reservoir model with an ABAQUS model for elasto-plastic deformation of the caprock which includes the built-in smeared crack subroutine.The reservoir model calculates distributions of pore pressure and stress within the reservoir, using which the displacement of the top of the reservoir is estimated.The reservoir model uses a finite difference scheme coded in FORTRAN, the results of which are validated against an analytical solution as presented in [27].The estimated displacement is then used as a condition for the lower boundary of the caprock model, where a bonded (frictional) interface is assumed, i.e. the displacement of the top of the reservoir is the same as the displacement of the bottom of the caprock.The caprock model constructed in ABAQUS, will then estimate distributions of stress and strain within caprock using a finite element scheme.Using these estimated stress and strain, the crack model checks the crack criterion within all elements of the caprock.Upon crack initiation in each element, the constitutive equations for the cracked rock are employed to estimate the stress and strain within the rock.Results of the caprock model have been validated against the ones obtained from the analytical method of solution [30,31].The flowchart of the integrated model is given in Figure 1.

D. Model construction
In order to set up the model, a disk shape reservoir is considered underlying a caprock of the same shape, while their thickness can be independent of each other.Since the reservoir model only provides the displacement of the top of the reservoir as a boundary condition to the caprock model, the constructed model in ABAQUS includes a three dimensional representation of the overlying caprock.Considering isotropic reservoir and caprock, an axial symmetry around the wellbore can be assumed and the problem of disk can be reduced to a quarter.The spatial discretization of the caprock model in ABAQUS is illustrated in Figure 2a.The boundary conditions for the model are depicted in Figure 2b.As explained earlier, the displacement calculated for the top of the reservoir is used as the displacement of the bottom of the caprock which defines the condition for the bottom side of the model.The side faces are fixed in the normal direction, which means there would be only vertical and radial displacement on these surfaces.For the wellbore in the center of the disk, the pressure is assumed to be equal to the injection pressure to simulate an uncased wellbore.The total vertical stress applied at the top of the models is assumed to be constant during the whole injection time, and considered to be equivalent to the weight of the overburden.Also, the initial horizontal stress is assumed to be equal to the total vertical stress.
The model is a representation of the Yort-e-Shah aquifer in Iran which is located 70 km southeast of Tehran and 40 km south of Varamin.The geological studies of the formations in the region revealed that the aquifer is formed of limestone and gypsum, while the caprock is formed of anhydrite [32].The reservoir and caprock properties for Yort-e-Shah aquifer are mainly available from the site investigation reports as listed in Tables I-III.The radius of the caprock is considered to be 500 m, and the thickness is 60 m from depth of 240 m to 300 m.Data from [33] is used when relevant parameters are not available for Yort-e-Shah region.

III. NUMERICAL RESULTS
Using the parameters introduced above, a numerical analysis was conducted.The result of reservoir simulation is shown in Figure 3, where it can be seen that the pressure increase as injection continues.The outer boundary of the reservoir is assumed impermeable which means the pressure within the reservoir is built up with time until the whole reservoir reaches a uniform pressure.The same trend is observable for displacement at the top of the reservoir, as expected.These results are used as a boundary condition for the bottom of the caprock at every time step.The results of the simulations using the assumed conditions and reservoir outputs show that some limited cracks are developed within the caprock around the injection wellbore after 330 days.The locations of these cracked elements in the model are shown in Figure 4.The distribution of the maximum principal stresses within the caprock is depicted in Figure 5.As shown, the caprock is cracked where the elements are under tensile stress regime.The profile of the displacement at the bottom of the caprock is depicted in Figure 6.It can be seen that the caprock moves upward as injection continues, while the absolute magnitude of the uplift is larger near the injection well (note that all upward vertical displacements in the caprock are considered negative in this paper).By gas injection, the pressure in the reservoir rises gradually.By increasing of gas pressure in the reservoir, the effective stress reduces, therefore, the upward moving will be seen in caprock.

www.etasr.com Rajabi et al.: A Numerical Model for Caprock Analysis for Subsurface Gas Storage Applications
The orientation of the cracks developed within caprock is illustrated in Figure 7.Note that these are collected at the integration points of the finite elements of the model.It can be seen that the crack planes are directed perpendicular to the minimum principal stress.This is in agreement with what is widely known [38].Also note that this is the case around the whole circumference of the wellbore, which is an indication of the fact that the major and minor principal stresses are following wellbore's geometry in an axisymmetric manner.This is the case because an isotropic caprock has been assumed in this model, while developing anisotropic model can lead to a non-axisymmetric stress distribution, but the orientation of cracks are still expected to be perpendicular to the minor principal stress.The results of the caprock model have been validated against the ones obtained from the analytical method.A simplified uniaxial vertical displacement model is used to estimate the vertical displacement of the caprock bottom due to gas injection by (19) [30,31]: where ΔP is the change in pore pressure and ∆h is the vertical displacement.Ε, υ, α, and h r parameters were defined in Section II.Equation ( 21) can be used to calculate vertical displacement with the change in reservoir pore pressure.The vertical displacement of the caprock bottom due to gas injection, computed from the analytical solution and the numerical scheme are shown in Figure 8, after 99 and 330 days.There are differences between analytical and numerical results.When the maximum pore pressure change at the reservoir-caprock interface near injection well zone from the numerical result was used, the calculated vertical displacement from (21) was about 45.4 mm after 99 days of injection, which is 4 mm higher than the numerically calculated value (41.4 mm).The uniaxial analytical model overestimates the vertical displacement in comparison with a real case or the numerical results because the pore pressure is not distributed uniformly in the reservoir, unlike the assumption used in (21).Another reason is that the overburden stiffness, which is apt to restrict vertical deformation in a real 3D situation, is not considered in the uniaxial strain model of (21) [7,31].The integrity of a caprock is an important feature of storage sites.This paper presents preliminary results of a model for caprock integrity where smeared crack approach is employed.Changes in reservoir pressure change the stress conditions of the reservoir and the surrounding rocks and result in redistribution of stress within the caprock.The paper presented the hydro-mechanical coupling of a reservoir model and the caprock model.The latter included the use of smeared crack module of ABAQUS which was determined to be a computationally effective tool in modeling crack development within the caprock and the algorithm of the developed model was presented.The model was employed to analyze the crack development within the caprock of Yort-e-Shah aquifer in Iran.Results suggest that the caprock in this aquifer can help a safe storage since there is only a limited zone around wellbore prone to crack development.However, this conclusion is limited to the conditions assumed in this paper.Further investigation is required to assess the integrity of caprock under different injection pressure, loading and deformation conditions and assumptions.Further studies are being undertaken to include the impact of different boundary and initial condition, as well as different injection plans to assess the suitability of Yort-e-Shah aquifer subsurface gas storage.
el d  can be further decomposed into elastic, el d d  and plastic strains rates in cracked rock, pl t d  as below:

Fig. 3 .
Fig. 3. Results of reservoir model at different injection times; (a) pore pressure distribution and (b) vertical displacement of the reservoir ceiling

TABLE II .
PROPERTIES OF CAPROCK OF YORT-E-SHAH-AQUIFER f  c  y 

TABLE III .
RESERVOIR PROPERTIES IN YORT-E-SHAH REGION f  f 