## Abstract

A competent low permeability and chemically inert geological barrier is an essential component of any strategy for the deep geological disposal of fluidized hazardous material and greenhouse gases. While the processes of injection are important to the assessment of the sequestration potential of the storage formation, the performance of the *caprock* is important to the containment potential, which can be compromised by the development of cracks and other defects that might be activated during and after injection. This paper presents a mathematical modelling approach that can be used to assess the state of stress in a surficial caprock during injection of a fluid to the interior of a poroelastic storage formation. Important information related to time-dependent evolution of the stress state and displacements of the surficial caprock with injection rates, and the stress state in the storage formation can be obtained from the theoretical developments. Most importantly, numerical results illustrate the influence of poromechanics on the development of adverse stress states in the geological barrier. The results obtained from the mathematical analysis illustrate that the surface heave increases as the hydraulic conductivity of the caprock decreases, whereas the surface heave decreases as the shear modulus of the caprock increases. The results also illustrate the influence of poromechanics on the development of adverse stress states in the caprock.

## 1. Introduction

The injection of fluids into porous formations is regarded as a procedure for managing contaminants and other hazardous waste [1–3]. Deep injection of fluidized forms of greenhouse gases into storage formations has been advocated as a means of mitigating the effects of climate change [4–7]. The process of injection into a porous storage formation pre-supposes that the fluid can be injected without detriment to the fabric of the porous medium and the geological setting can accommodate hydrodynamic trapping that is provided by a stable geological barrier (figure 1). How the caprock performs is therefore of critical importance to the longevity and effectiveness of the storage strategy. In the context of geological sequestration of greenhouse gases, the factors controlling the mechanics of the caprock during the injection process has to take into consideration a complex set of thermo-hydro-mechanical–chemical (THMC) processes. These can result from the geochemistry of the storage formation, the temperatures of the fluids being injected, which can perturb the *in situ* geothermal gradients, the buoyancy anomalies created by the injected fluid, the development of migration plumes that correspond to moving boundaries and chemomechanical processes that can lead to deterioration of the fabric of the storage formation. A comprehensive mathematical treatment of all these aspects presents a challenging problem in environmental geomechanics and such investigations are best addressed through investigations of site-specific situations. If attention is restricted to the investigation of caprock performance, then the hydromechanical (HM) processes that are created during injection come to the forefront. For example, the injection rates that will not cause distress to the caprock in terms of the potential for fracturing of a given caprock–storage formation can be examined by considering the appropriate HM problem. Even in a simplified HM investigation of the caprock–storage formation interaction, the actual geological setting is important and real-life problems are best examined through the adoption of a computational procedure. Some progress towards understanding caprock–storage formation interaction during internal injection can be made by examining the poroelasticity effects of the interaction. The idealization of the storage formation as a poroelastic medium without the influences of chemical interaction is considered to be appropriate for situations where the materials encountered are competent sandstones. This has been demonstrated by the selection of many storage settings in sandstone formations (e.g. Sleipner, Snøhvit [8], In Salah [9], Nagaoka [10]). When the storage medium has a dominant carbonate content, the reactive nature of carbon dioxide acidized water can lead to dissolution of the carbonates causing material erosion and the formation of features such as *wormholes* that can compromise the integrity of the storage effort (Selvadurai APS, Couture C-B. 2016 Wormhole generation in carbonate zones during reactive flows. McGill University, unpublished). Furthermore, the dissolved carbonates can accumulate in remote regions of the storage formation, which can result in pore clogging and ultimately lead to hydraulic fracture during continued injection into the storage formation, thus limiting its trapping potential. In addition, if extensive dissolution takes place in the storage horizon, then this can lead to collapse of the *wormholes* under geostatic stresses that can lead to caprock distress owing to loss of geological structural support provided by the storage formation.

In this paper, we focus attention on the mathematical modelling of the hydromechanical effects that can materialize when a storage formation underlying a caprock layer of uniform thickness is pressurized by the injection of a fluid over a circular region located at a finite distance from the base of the caprock. The objective of the study was to determine the stress state in the caprock during the injection process. If injection is modelled as the rise and maintenance of the pressure in a specified region in the storage formation [11], then the state of stress in the caprock can be examined using an elastostatic model of the storage formation where the caprock is idealized as a structural element with a response that can be modelled as either a Poisson–Kirchhoff thin plate [5,6] or as a Mindlin–Reissner thick plate [12]. The elastostatic model can also provide quantitative estimates for the deflections of the caprock in the long term. The elastostatic model, however, *cannot* provide an assessment of the time-dependent effects of injection associated with geologic sequestration activities. In particular, the rates of injection need to be controlled in such a way that hydraulic fracturing will not be initiated within the storage horizon and the injection pressures will not compromise the sealing capacity of the caprock by initiation and extension of fractures. Because the primary focus of the paper is on the development of a mathematical solution to study the caprock–storage formation interaction during injection of fluids to a storage formation, we consider the situation where the caprock layer is located at the surface of a semi-infinite storage formation. In a practical situation, the caprock layer may be an embedded stratum (figure 1), but the consideration of a surface layer provides the opportunity to examine how the absence of confinement influences the performance of the caprock layer. The particular focus of this study is to provide an analytical result that can serve as a useful benchmark for examining the accuracy of computational approaches rather than to provide an analysis of a specific engineering situation associated with geological sequestration, which will aid the general applicability of the study. Even within the context of a simplification that incorporates a surficial caprock layer underlain by a semi-infinite storage horizon, the analytical solution of the poroelasticity problem can involve extensive mathematical manipulations. Neglecting the influence of a poroelastic surficial geological region is a reduction that makes the benchmark exercise manageable in a theoretical context. In this regard, it should be noted that the case of a more practical situation involving a caprock stratum embedded between a surficial geological medium and a storage horizon, the influence of the surficial geology can be approximately accounted for by specifying the geostatic loads associated with the surficial rocks. In such a treatment, the absence of poroelastic effects in the geological medium overlying the caprock will lead to an upper limit of caprock deflections during the injection process.

To preserve axial symmetry of the modelling, we consider the case where the fluid injection is carried out over a circular region located at a finite depth from the boundary of the caprock (figure 2).

The modelling adopts a poroelastic formulation for the mechanical behaviour of both the caprock and the storage formation. The poroelastic formulation of a defect-free geological medium is considered to be a satisfactory model for rocks that are encountered in competent storage horizons and with caprocks that are free of major defects. This allows the use of Biot’s classical theory of poroelasticity [13] for the study of the problem indicated in figure 2, by assigning different poroelasticity representations for both the storage horizon and the caprock. The governing equations of poroelasticity for mechanically and hydraulically isotropic fluid-saturated porous elastic media are given by several authors including Rice & Cleary [14], Atkinson & Craster [15], Cheng [16], Craster & Atkinson [17], Selvadurai [18–20], Selvadurai & Yue [21] and Selvadurai & Mahyari [22,23]. The mathematical formulation extends previous studies [24,25] that investigated, respectively, the axisymmetric problem of a poroelastic halfspace region with fluid extraction over a circular planar region located at the interior of a poroelastic halfspace and the three-dimensional problem of line-injection sources located in a poroelastic halfspace region. This study considers the more complex problem where the surface of the halfspace region contains a contiguous poroelastic layer with continuity of displacements, tractions and pore fluid pressures and fluxes. The mathematical analysis of the problem is developed through the application of Laplace and Hankel transform techniques, which reduces the initial boundary value problem to a set of ordinary differential equations that can be solved in exact closed form; numerical procedures are then used to invert both the Laplace and Hankel transforms. The paper presents numerical results for the surface deflections of the caprock and its stress state at critical locations as a function of the fluid injection rate and other mechanical and fluid flow characteristics of the geologic materials. The emphasis of the paper is on the study of the poroelastic response of geological media encountered in sequestration scenarios. The work can be extended to include more advanced nonlinear constitutive modelling of the rocks [26–30] to include thermoporoelastoplastic effects in sequestration activities. In the case of geological sequestration, however, the geological settings are usually chosen such that plasticity effects are negligible, but on occasions, such influences are unavoidable. The role of plasticity can have a notable influence on the mechanical behaviour of rocks in the vicinity of the injection location. The results of Selvadurai & Suvorov [30] for the thermoporoelastoplasticity problem for the internal pressurization of a cavity indicate that plasticity effects can contribute to an increase in displacements in the vicinity of the pressurized cavity.

## 2. Governing equations

The fully coupled equations governing the mechanics of a poroelastic medium are given in the references cited previously in this section, the relevant final forms of the equations governing the displacement and pore pressure fields are presented. Attention is restricted to a state of axial symmetry in the fluid injection process, which is applied with a prescribed time history over a circular disc-shaped region (figure 2). The dependent variables of the problem are the skeletal displacements *u*_{r}(*r*,*z*,*t*) and *u*_{z}(*r*,*z*,*t*) in the *r*- and *z*-directions, respectively, and the pore fluid pressure *p*(*r*,*z*,*t*), which are governed by
*Θ* is the volumetric strain. In addition, ∇^{2} is the axisymmetric form of Laplace’s operator given by
*G* and *ν* are, respectively, the shear modulus and Poisson’s ratio of the porous skeleton (i. e. the drained elastic parameters); *ν*_{u} is the undrained Poisson’s ratio of the fluid-saturated medium; *k* is the hydraulic conductivity; *γ*_{w} is the unit weight of the pore fluid; *B* is Skempton’s pore pressure parameter; *C*_{m} is the compressibility of the porous skeleton; *C*_{s} is the compressibility of the skeletal material; *C*_{f} is the compressibility of the pore fluid and *n* is the porosity. The constitutive parameters have to satisfy certain thermodynamic constraints to ensure positive definiteness of the strain energy potential; it has been shown [14] that these constraints can be expressed in the forms: *G*>0; 0≤*B*≤1; −1<*ν*<*ν*_{u}≤0.5. Alternative but equivalent descriptions are described by Cheng [16], Wang [31] and further references are given in Selvadurai [18,19].

### (a) Solution method of the governing equations

Several approaches for the solution of problems in poroelasticity have been proposed in the literature [19,21,32–41]. We use the approach proposed by McNamee & Gibson [34,35] and Gibson & McNamee [36] where the solution to the governing coupled partial differential equations (2.1)–(2.3) can be represented in terms of two scalar functions *S*(*r*,*z*,*t*) and *E*(*r*,*z*,*t*), which satisfy the partial differential equations
*S*(*r*,*z*,*t*) and *E*(*r*,*z*,*t*) as follows
*S*(*r*,*z*,*t*) and *E*(*r*,*z*,*t*) can be verified by back-substitution into the equations (2.1)–(2.3). Analytical solutions for the time-dependent axisymmetric poroelasticity problem can be obtained by using integral transform techniques. Laplace and zeroth-order Hankel transforms are used to remove, respectively, time-dependency and dependency on the radial coordinate *r*

## 3. Initial boundary value problem related to an internally located injection zone

We consider the problem of a poroelastic halfspace region that is overlain by a deformable poroelastic caprock (figure 2). The pore fluid pressures within the caprock layer and the halfspace are considered to be sessile and at hydrostatic values. The *in situ* stress state in the region is assumed to be geostatic with the effective stresses generated purely from the self-weight of the geological material and the pore fluid pressures. The role of tectonically induced stresses can also be considered in the analysis provided that the stress states can be represented by equivalent principal stress states that preserve axial symmetry. The developments in poroelasticity presented by (2.1)–(2.9) can therefore be regarded as the representations of displacement, stress and pore pressure fields associated purely with the effects of injection. The interface between the poroelastic caprock layer and the poroelastic halfspace region is assumed to be contiguous thereby allowing complete continuity conditions to be prescribed with respect to the displacements, effective stresses, pore fluid pressures and fluid fluxes. A planar circular area occupying the region 0≤*r*≤*a* at *z*=0 is subjected to fluid influx at a constant total flow rate *Q*_{0} (units L^{3}/T). (i.e. This represents the *total volume of fluid that is injected per unit time*, over the plane circular region and we note that as the radius of the injection region *a* increases the injection intensity measured as volume per unit area decreases.) The boundary conditions are prescribed in such a way that the injection rates are dimensionally consistent. The specification of the*total volume rate of injection* over the circular region of radius *a* also makes it possible to compare the results of this study with the *point injection* studies, such as the nuclei of fluid influx, that have been developed for both infinite and semi-infinite domains [5,6,12]. We note that the assumption of injection over a circular disc-shaped region enables the convenient formulation of the initial boundary value problem. In a practical context, if the injected fluids have different physical and mechanical characteristics to the resident fluids, then this will result in the development of a three-dimensional plume, the configuration of which needs to be determined from the solution of a moving boundary problem. Such problems are more suited to analysis via computational approaches. The assumption of a planar injection zone and the injection of fluids with identical properties make the resulting initial boundary value problem amenable to the mathematical formulation and solution. From a practical point of view, the specification of a uniform injection rate over a disc-shaped planar region is perhaps unrealistic, but the solution is intended to provide a fundamental result that can be used, with a superposition technique, to generate other non-uniform axisymmetric injection scenarios. The initial boundary value problem governing the fluid injection can be described by referring to the following poroelastic domains: (i) *a caprock poroelastic layer region*, identified by the superscript ()^{(C)} and occupying the region *z*∈(−*h*,−(*l*+*h*)), (ii) *a poroelastic layer region*, identified by superscript ()^{(L)} and occupying the region *z*∈(−*h*,0) and (iii) *a poroelastic halfspace region*, identified by superscript ()^{(H)}and occupying the region

(i) The boundary conditions applicable to the surface *z*=−(*l*+*h*) are derived from the assumption that the surface of the caprock layer is unloaded and the pore fluid pressure is zero, i.e.
*H*(*t*) is the Heaviside step function. For completeness of the formulation of the initial boundary value problem related to the poroelasticity problem, we specify the initial conditions
**x** are the spatial coordinates, **u**(**x**,*t*) is the displacement vector and ** σ**(

**x**,

*t*) is the total stress tensor. In addition, the solution to the poroelasticity problem should satisfy the regularity conditions applicable to regions where the spatial dimension(s) can extend to infinity.

Because semi-infinite and layer regions that extend to infinity are encountered in the problem formulation, and because the stress states that induce the poroelastic effects are restricted to finite regions, the solution to the fluid injection problem should also satisfy the regularity conditions
*z*∈(−*h*,−(*l*+*h*)) and for ∀ *t*≥0 and
*t*≥0.

The uniqueness of solution to the initial boundary value problem posed by (2.1)–(2.3) and the consistent boundary and regularity conditions (3.1), (3.11) and (3.12) and the initial conditions (3.10) can be assured by the general proof of uniqueness of solution to the linear poroelasticity problem [42] and the general proof of uniqueness to the problem in classical elasticity [43,44], which is applicable to poroelastic media both at *t*=0 and

## 4. Solution of an initial boundary value problem

The solution procedure commences with the development of solutions of the coupled system of ordinary differential equations (2.12) and (2.13) applicable to the transformed dependent variables ^{(C)}, the poroelastic layer region within the halfspace region ()^{(L)} and the poroelastic halfspace region ()^{(H)} can be expressed in the following forms, ensuring that the regularity conditions (3.11) and (3.12) are satisfied *a priori*.

(i) For the poroelastic caprock layer occupying the region *z*∈(−*h*,−(*l*+*h*)), the solutions of (2.12) and (2.13) take the forms
*A*_{1}, *B*_{1}, *C*_{1}, *D*_{1}, *E*_{1} and *F*_{1} are arbitrary constants.

(ii) For the poroelastic layer above the injection zone and occupying the region *z*∈(0,−*h*), the solutions of (2.12) and (2.13) take the forms
*A*_{2}, *B*_{2}, *C*_{2}, *D*_{2}, *E*_{2} and *F*_{2} are arbitrary constants.

(iii) For the poroelastic halfspace region below the injection zone and occupying the region *A*_{3}, *C*_{3} and *E*_{3} are arbitrary constants. The 15 arbitrary constants can be uniquely determined by making use of the boundary conditions and continuity conditions given by the 15 equations in (3.1)–(3.9). It is sufficient to note that explicit expressions for the unknown constants *A*_{1},*B*_{1}… etc., can be obtained through a Mathematica^{®}-based symbolic evaluation. The complete expressions for the unknown constants are given in electronic supplementary material, appendix A. The time-dependent displacement fields, stress components and pore pressure fields in the poroelastic caprock layer and within the poroelastic halfspace region can be obtained through a Laplace transform inversion and a Hankel transform inversion of the resulting expressions. At this juncture, it is worth commenting on the poroelastic modelling of the problem where a layer of poroelastic caprock is embedded between a geological overburden, which can be modelled as a poroelastic layer of finite thickness and a storage formation, which is modelled as a poroelastic halfspace. In this case, additional boundary and interface conditions need to be prescribed between the poroelastic overburden and the poroelastic caprock, which will give rise to six additional interface conditions similar to those defined by (3.2)–(3.5), giving rise to 21 arbitrary constants governing the poroelasticity formulation of the idealized analogue of a problem depicted in figure 1. Quite apart from the complexities associated with the mathematical solution of the problem, the illustration of numerical results will have to contend with an inordinate number of material parameter groups that will make the exercise unmanageable and of limited utility as a benchmark problem. For this reason, the approach adopted in the reduced formulation of the benchmark problem is favoured.

The results that are of interest to geological sequestration are the *time-dependent vertical displacement* of the surface of the caprock layer and the *time-dependent peak radial stress* at the surface of the caprock layer. These results are particularly important for monitoring the influence of subsurface injection on caprock heave and estimating the potential for caprock distress, in terms of the development of tensile fractures in the caprock. The expressions for the time-dependent axial displacement and time-dependent radial stress are given by
*ζ* is a real constant associated with the Bromwich contour integration encountered in Laplace transform inversion. Other expressions for the displacements, pore fluid pressure and effective stresses are presented in electronic supplementary material, appendix B. The complexity of the expressions of the type (4.7) and (4.8) makes it difficult to conduct analytical evaluations of the inversion of the integral transforms. Previous experience [22,25,41,45,46] suggests that the operations involving Hankel transform inversions and Laplace transform inversions have to be carried out numerically to generate expressions for the displacements and stresses that can be expressed in terms of *r*, *z* and *t*. In this study, the integral transforms were evaluated using the algorithms given in Matlab^{®}. In particular, the inverse Laplace transform is evaluated using the procedures developed by Crump [47]. The accuracy of the inversion technique is confirmed through comparisons with results obtained by McNamee & Gibson [35] for the axisymmetric surface loading, *f*, of a halfspace region. Table 1 presents the numerical results evaluated in this study, and the analytical results given in McNamee & Gibson [35]. The error is large when *t*<0.1 but substantially decreases as *t* increases. The numerical results have an error of less than 0.1% when *t*=10, and, therefore, the procedures are considered satisfactory for evaluating a full range of analytical results for the poroelastic problem.

## 5. COMSOL^{™} modelling

An objective of the mathematical treatments presented in the paper is also to develop an analytical result that can be used to test the accuracy of computational schemes, which can be applied to examine more complex geologic sequestration settings. In this regard, the computational solutions were obtained using the finite-element-based multi-physics software COMSOL^{™}. This code has been extensively scrutinized, and details of the validation exercises are given in [24,25,29,30,45,48–51]. Figure 3 shows the representative finite-element discretization model for the injection problem that was used to develop the computational simulation. Because the problem is axisymmetric, only a section of the domain (about *z*-axis) needs to be modelled. A finite domain is used instead of the infinite region that was used in the analytical development. The effect of the finite boundary is minimized by taking the boundary of the domain sufficiently remote from the fluid injection region, thereby ensuring that the far-field boundary conditions are approximately satisfied. An alternative approach is to incorporate infinite elements that can model the far field behaviour [46,52,53]. This option is, however, unavailable in the multi-physics code COMSOL^{™}. At the surface, the pore pressure is set to zero, and a flux discontinuity is applied to the interface between the halfspace layer above the injection site and the halfspace region. A symmetry/zero flux boundary condition is applied at the base and the right-hand side boundary of the domain. Figure 4 shows the mesh configuration used in the COMSOL^{™} modelling; mesh refinement is used in the vicinity of the injection region and coarser meshes were used for the rest of the domain. The mesh consists of 12 570 elements for *a*/*h*=0.1 and increases to 33 840 for *a*/*h*=10. Each nodal point has three degrees of freedom.

It should be noted that as *k*_{s}/*k*_{c} increases, the COMSOL^{™} model could not produce results that were comparable to the analytical results. The maximum value for the permeability mismatch between the caprock layer and storage formation used in the COMSOL^{™} model was *k*_{s}/*k*_{c}=100, and therefore, the values of *k*_{c} and *k*_{s} were chosen accordingly.

## 6. Results

The results for the vertical surface displacements, pore fluid pressure and effective stress distributions along the axis of the stratified region are presented in this section. The parameters with subscript ‘c’ are the caprock variables, and those with subscript ‘s’ are the variables in the storage formation. For the case where *G*_{c}=*G*_{s} and *k*_{c}≠*k*_{s} we used the following variables to evaluate the analytical and computational results:
*G*_{c}≠*G*_{s}, *k*_{c}≠*k*_{s}, except that *G*_{c} was set to 50 GPa. The results are presented for a non-dimensional time, *t**=*ct*/*h*^{2}. The formulation of *c* is given in equation (2.4) and is evaluated by using *G*=Min(*G*_{c},*G*_{s}) and *k*=Min(*k*_{c},*k*_{s}). Throughout this study, the total volume flow rate *Q*_{0} and the thickness *l* of the caprock are fixed. For ease of presentation, we neglect the negative sign in the displacement, with the understanding that the surface displacement in the caprock occurs in the negative direction of the *z*-axis as indicated in figure 2.

### (a) Case when the permeability of the caprock differs from that of the storage formation

#### (i) Surface displacement

The variation in the time-dependent surface displacement is shown in figure 5 for different values of *a*/*h* when *h*/*l*=1. The solution at *t**=1000. As indicated in equation (4.7), for a constant *Q*_{0}, the injection intensity will decrease as *a* increases. As is evident from figure 5, the surface heave decreases as *a*/*h* increases. The corresponding computational results are denoted by circles. The discrepancy between the analytical and computational solutions is approximately 11% around the central location.

Selvadurai & Kim [24] investigated the ground subsidence due to fluid extraction from a circular disc-shaped region, which can also be used as a solution for the fluid injection problem owing to the linearity of the fundamental equations governing Biot’s theory of poroelasticity. Figure 6 compares the heave at the interface between the caprock and the storage halfspace, *z*=−*h*, with the surface heave when the caprock is absent [24]. The surface heave is normalized by the steady-state value at *r*=0. It is observed that without the caprock, the surface displacement occurs only around the injection region and decreases quickly as *r* is remote from the injection region. However, in the presence of the caprock, the surface displacement occurs over a broader region and decreases gradually as the point of interest is remote from the injection region. In the case without the caprock, surface flattening and widening of the centre were observed in the injection region at *t**=1 and 10 when *a*/*h*=5 and 10. However, when the caprock is present, no flattening of the surface is observed for the cases *a*/*h*=5 and 10; only widening of the centre is observed in the present results. When the present results are compared with those obtained for the case without the caprock, it was found that the ratio of the surface heave to the steady-state displacement is larger when the caprock is present. When the caprock is present, the surface displacement when *a*/*h*=10 attains in excess of 60% of the steady-state response at *t**=10 and 30% at *t**=1; however, the surface displacement for the case *a*/*h*=10 at *t**=10 was only 40% of the steady-state displacement without the caprock.

Figure 7 illustrates the surface heave for different depths of location of the injection zone, i.e. *a*/*h*=0.1 and *h*/*l*=1 to 10. It can be seen from figure 7 that varying the depth of the injection region influences the surface heave in a similar manner to varying the injection size. As *h*/*l* increases, the surface heave decreases. This observation is analogous to a *Saint Venant*-type effect in classical elasticity, where remoteness of the injection zone ensures that the caprock displacements are reduced. In the limit as *h*/*l* becomes large (e.g. *h*/*l*>10), the injection pressures have only a marginal influence on the caprock displacements. The corresponding computational results are shown as circles in figure 7. The agreement between analytical and computational results is generally good; however, there is a discrepancy between the two solutions.

In this study, the maximum discrepancy is 20% at the central location.

Figure 8 compares the surface heave for different *k*_{s}/*k*_{c} at *t**=0 and *k*_{s}/*k*_{c}=10^{5} are larger than those obtained for *k*_{s}/*k*_{c}=10^{2}. No distinct change in the shape of the surface heave is observed. The maximum surface heave for different values of *k*_{s}/*k*_{c} is shown in figure 9; as *k*_{s}/*k*_{c} increases, the surface displacement increases. A large *k*_{s}/*k*_{c} indicates a low hydraulic conductivity in the caprock layer. When the hydraulic conductivity is low, the pore fluid pressure tends to diffuse slowly from the vicinity of the injection zone, which contributes to a greater surface heave.

#### (ii) Pore fluid pressure

Figure 10 presents the changes in the pore fluid pressure owing to fluid injection into the poroelastic medium when *a*/*h*=0.1 and *h*/*l*=1. The pore fluid pressure is at its maximum at the injection site, which is located at the interface between the layer (*L*) and the halfspace, and decreases very rapidly as the position *z*/*h* is remote from the injection location. However, in the caprock where the hydraulic conductivity is low compared with the layer (*L*) and the halfspace, the change in the pore fluid pressure is very small as *z* moves away from the injection site. The corresponding computational results are presented as circles and they show good agreement with the analytical solutions except at locations close to the injection region.

#### (iii) Effective stresses and caprock distress

It was observed that the pore fluid pressure increased owing to fluid injection, and therefore, the effective stresses are expected to decrease in the injection region. The change in the effective stress (*σ*′_{ij})^{Q0} is evaluated from the result
*σ*_{ij})^{Q0} and (*p*)^{Q0}are the changes in the total stress and pore fluid pressure owing to fluid injection. Figures 11 and 12 show the effective stresses in the *r*- and *z*-directions, respectively. The negative values indicate a decrease in the effective stresses in both directions owing to fluid injection. In addition, in the caprock layer, it was observed that negative effective stress (tension) developed in the *r*-direction. At the surface, traction free boundary conditions are invoked, and therefore, only radial stresses are present. The possibility of tensile failure can be checked by comparing this effective stress with the tensile strength (*T*) of the caprock medium
*k*_{c}=10^{−12} m s^{−1}, *k*_{s}=10^{−10} m s^{−1} *a*=50 and *h*=500, the effective stress at the surface of the caprock layer in the *r*-direction is approximately 4.289×*Q*_{0}×10^{2} MPa (figure 11). For tensile failure to occur at the surface, the following condition needs to be satisfied
*T* is the tensile strength of the shale caprock, which is estimated at 2 MPa. Possible tensile failure can occur when *Q*_{0}>0.466×10^{−2} m^{3} s^{−1} (≈0.1456 megatonnes per year (Mt yr^{−1})).

In a practical geological sequestration scenario, however, the caprock is most likely to be overlain by either surficial soils or rocks similar to the situation illustrated in figure 1. We consider the case where the caprock is overlain by a surficial soil, such as clay of saturated unit weight, *γ*_{s}=20 kN m^{−3}. The surficial soil layer is assumed to contribute little to the poromechanical interaction between the caprock and the storage formation but provides only a surcharge loading at the surface of the caprock layer by virtue of the self-weight. The self-weight of the surficial soils will introduce an additional effective radial compressive stress in the caprock, which is given by
*D* is the thickness of the surficial soil layer. Table 2 presents an estimate of how this surficial soil layer, by virtue of the additional self weight-induced loading, can contribute to the increase in the allowable injection rate that can be attained without the development of fracture at the upper boundary of the caprock layer. Table 3 provides a summary of the typical annual CO_{2} injection rates for the selected greenhouse gas storage sites.

The corresponding computational results for the *r*- and *z*-direction effective stresses are presented as circles in figures 11 and 12, respectively. The agreement is good for the effective stress in the *z*-direction; however, the COMSOL^{™} model underestimates the effective stress in the *r*-direction. In particular, the effective stress in the *r*-direction obtained by COMSOL^{™} at the surface is less than 40% of the analytical result. It must be emphasized that the analytical results presented in the paper are accurate and represent the solution to the initial boundary value problem defined by (3.1)–(3.12). The accuracy of the analytical technique has also been established through comparisons with known analytical results given by McNamee & Gibson [35], albeit for the surface loading of a poroelastic halfspace region. The issue of the discrepancy between the analytical results and the computational results obtained, particularly for the radial stress, can arise from a number of approximations including the representation of the halfspace region by a domain of finite extent and the application of specific boundary conditions at the far-field boundary of the computational domain as opposed to regularity conditions applicable to the semi-infinite domains. The mesh discretization in the vicinity of the injection zone can also contribute to the discrepancy. The trends in all computationally derived effective stress distributions are consistent with the results obtained from the analytical approach. We also note that the results presented for the poroelastic halfspace without a surficial layer is not meant to reflect any practical engineering situation but are included to illustrate the moderating influence of the surficial barrier in mitigating the boundary displacements of the storage formation.

#### (iv) Case when both elastic modulus and permeability of the caprock differ from those in the storage formation

The surface displacement was also investigated for the case where the caprock has a different shear modulus and hydraulic conductivity than those in the layer and halfspace. Figure 13 shows the non-dimensional surface displacement for different cases considered in this study when *a*/*h*=0.1 and *h*/*l*=1. It is observed that the magnitude of the surface heave is greater when *G*_{c}=*G*_{s} and *k*_{c}≠*k*_{s} compared with the cases when *G*_{c}=*G*_{s} and *k*_{c}=*k*_{s}. The surface heave obtained when *G*_{c}=*G*_{s} and *k*_{c}=*k*_{s} has a maximum at the centre and quickly converges to 0 as *r*/*h* moves away from the centre, whereas the results obtained from the other cases show a more gradual decrease. When the hydraulic conductivities differ, i.e. *k*_{c}≠*k*_{s}, the magnitude of the surface displacement obtained for the case where *G*_{c}=*G*_{s} is larger than that obtained when *G*_{c}≠*G*_{s}.

Figure 14 presents the maximum surface heave for different *G*_{s}/*G*_{c}. Two cases, (i) *k*_{c}=*k*_{s} and (ii) *k*_{c}≠*k*_{s}, are considered, and both cases show exactly the same trend except for the magnitude of the surface heave. In both cases, the surface heave increases as *G*_{s}/*G*_{c} increases.

## 7. Concluding remarks

The paper develops analytical solutions for the poroelastic modelling of the interaction between a poroelastic caprock and a poroelastic storage aquifer, for the case where steady injection takes place over a planar circular injection zone. The mathematical analysis can be performed to develop exact closed results in integral representations that can be evaluated using standard algorithms available for Laplace and Hankel transform inversions. Two cases were examined: in the first case, the elastic shear modulus for the caprock is the same as that of the storage horizon, but the hydraulic conductivities are different (*G*_{c}=*G*_{s}, *k*_{c}≠*k*_{s}). In the second case, both the shear modulus and the hydraulic conductivities are different (*G*_{c}≠*G*_{s}, *k*_{c}≠*k*_{s}). We investigated how the radius and the depth of the planar injection region influence the surface displacement for the first case. When the total flow rate to the injection zone is constant, as expected, an increase in the radius of the injection region leads to a reduction in the surface heave. Similarly, surface heave decreased as the injection depth increased. The zone of the storage formation above the injection area acts as a space for diffusion of the injected fluid, so that an increase in depth of the injection zone results in a decreased surface heave. In this study, it has been shown that both the shape and magnitude of the surface heave change significantly when there is a caprock layer present. As the hydraulic conductivity of the caprock decreases, the magnitude of the surface heave increased. We observe that in the hypothetical case where injection takes place into a poroelastic formation without a caprock, both the magnitude and profile of the surface heave can be distinctly different from the case where a caprock layer with different poroelasticity properties is present. The stiffening action of the caprock diffuses the spatial distribution of heave. These effects have also been observed in the purely elastic modelling of caprock–storage horizon interaction [5–7,12], the surface heave becomes flattened and localized to the region of the injection region for the case when the caprock is not present [24]. It was also found that the ratio of the surface heave to the steady-state displacement was larger when a caprock was present. Owing to the fluid injection, it was observed that the pore fluid pressure increased in the area surrounding the injection region and this caused a decrease in the effective stresses (both in the *r*- and *z*-directions). Because of the low hydraulic conductivities, the changes in the pore fluid pressure and the effective stresses were very small in the caprock compared with those within the halfspace storage region. One important observation is that adverse stresses, which lead to tensile failure, are developed at the caprock surface during fluid injection. For the case considered in this study, which is applicable to a shale similar to that encountered at the geological sequestration site at In Salah, it is observed that tensile failure can occur when *Q*_{0}>0.1456 Mt yr^{−1}. Assuming the presence of a surficial geological zone of thickness 1.5 km, which merely provides a surcharge load, the injection rate can be increased up to 0.517 Mt yr^{−1}, without distress to the caprock. The injection rates obtained through theoretical modelling of the caprock–storage formation interaction in this study are comparable with those associated with typical CO_{2} sequestration projects.

For the second case when *G*_{c}≠*G*_{s} and *k*_{c}≠*k*_{s}, the caprock experienced less surface heave compared with those obtained in the first case. As *G*_{c} increased, the surface heave decreased and the rate of decrease was larger when *k*_{c}=*k*_{s}. The magnitude of the surface heave decreased by 45% when *k*_{c}=*k*_{s}, whereas the magnitude of the heave decreased by 27% when *k*_{c}≠*k*_{s}. The analytical solutions were used to calibrate the results obtained using the finite-element-based code COMSOL^{™}. The computational results compare favourably with the exact analytical solutions.

## Ethics

The research conducted is unrelated ethical requirements indicated in the author guidelines.

## Data accessibility

The data files, computer codes and subroutines can be accessed at the following sites:On page: https://www.mcgill.ca/civil/people/selvadurai/research-repository File 1: https://www.mcgill.ca/civil/files/civil/instruction_for_source_file.pdf File 2: https://www.mcgill.ca/civil/files/civil/surficial_caprock_model.mph_.txt

## Authors' contributions

The conceptual development of the paper and the solution approach was developed by A.P.S.S. The mathematical formulation of the problem was prepared by A.P.S.S. and J.K. The computations were performed by J.K.; the accuracy of the computations was verified by A.P.S.S. The paper was written by A.P.S.S. and J.K.; the responses to the reviewer’s questions was written by A.P.S.S. and J.K.

## Competing interests

No competing interests.

## Funding

The work performed was supported by an NSERC Discovery Grant and the James McGill Professorship Grant awarded to A.P.S.S.

## Disclaimer

The use of the computational code COMSOL^{™} is for demonstration purposes only. The authors neither advocate nor recommend the use of the code without conducting suitable validation procedures to test the accuracy of the code in a rigorous fashion.

## Acknowledgements

The work described in this paper was supported by a *NSERC Discovery Grant*, the *James McGill Research Chairs* support and the *Carbon Management Canada Grant* awarded to the first author. J.K. acknowledges the partial research support provided by the Brace Center for Water Resources Management at McGill University and a McGill Engineering Doctoral Award. The authors are indebted to the reviewers for their constructive comments on the manuscript.

- Received June 20, 2015.
- Accepted February 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.