## Abstract

Shale gas has recently gained significant attention as one of the most important unconventional gas resources. Shales are fine-grained rocks formed from the compaction of silt- and clay-sized particles and are characterised by their fissured texture and very low permeability. Gas exists in an adsorbed state on the surface of the organic content of the rock and is freely available within the primary and secondary porosity. Geomechanical studies have indicated that, depending on the clay content of the rock, shales can exhibit a brittle failure mechanism. Brittle failure leads to the reduced strength of the plastic zone around a wellbore, which can potentially result in wellbore instability problems.

#### Author:

Mohsen S. Masoudian, Mir Amid Hashemi, Ali Tasalloti, Alec M. Marshall

Nottingham Centre for Geomechanics, University of Nottingham, University Park, Nottingham NG7 2RD, UK

Received: 19 July 2017 / Accepted: 22 December 2017 / Published online: 19 January 2018

© The Author(s) 2018.

Desorption of gas during production can cause shrinkage of the organic content of the rock. This becomes more important when considering the use of shales for CO_{2} sequestration purposes, where CO_{2} adsorption-induced swelling can play an important role. These phenomena lead to changes in the stress state within the rock mass, which then influence the permeability of the reservoir. Thus, rigorous simulation of material failure within coupled hydro-mechanical analyses is needed to achieve a more systematic and accurate representation of the wellbore. Despite numerous modelling efforts related to permeability, an adequate representation of the geomechanical behaviour of shale and its impact on permeability and gas production has not been achieved.

In order to achieve this aim, novel coupled poro-elastoplastic analytical solutions are developed in this paper which take into account the sorption-induced swelling and the brittle failure mechanism. These models employ linear elasticity and a Mohr–Coulomb failure criterion in a plane-strain condition with boundary conditions corresponding to both open-hole and cased-hole completions. The post-failure brittle behaviour of the rock is defined using residual strength parameters and a non-associated flow rule.

Swelling and shrinkage are considered to be elastic and are defined using a Langmuir-like curve, which is directly related to the reservoir pressure. The models are used to evaluate the stress distribution and the induced change in permeability within a reservoir. Results show that development of a plastic zone near the wellbore can significantly impact fracture permeability and gas production. The capabilities and limitations of the models are discussed and potential future developments related to modelling of permeability in brittle shales under elastoplastic deformations are identified.

## 1 Introduction

Unconventional gas reservoirs are increasingly being considered to provide a relatively clean energy source to meet burgeoning global demands. Unconventional gas resources, such as shale, represent around 40% of the remaining resource of technically recoverable natural gas; hence, they play an important role in the future global energy market (McGlade et al. 2013). In addition, unconventional gas reservoirs are being considered as a host rock for CO_{2} geo-sequestration, which may also provide an enhanced gas recovery (EGR) technique. Shales are fine-grained rocks formed from the compaction of silt- and clay-sized particles and are characterised by their fissured texture and very low permeability (10−18–10−23 m^{2}) (Alexander et al. 2011).

Gas shales can be considered as a dual porosity rock consisting of fractures and matrix blocks (Fig. 1), with a porosity of up to 15% (Wang and Reed 2009). Fracture sets can be subvertical or parallel to bedding, may be filled with mineral gauges, have openings ranging in scale from micrometres to centimetres, and spacing ranging in scale from a few centimetres to several metres (Gale et al. 2014). The matrix usually contains an organic content of 2–10%, while those with higher organic contents are usually considered too immature for production development (Alexander et al. 2011). The matrix contains porosity in the scale of µm to nm in size within different organic (fibrous) and nonorganic (clastic) structures (Mehmani et al. 2013).

Natural gas can be free within both the fracture (primary porosity) and clastic structure of the matrix (secondary porosity), or in an adsorbed state on the surface of the porous structure of organic content. The adsorption of gas in organic geomaterials can be described using an isotherm, among which Langmuir’s is most widely used (Gensterblum et al. 2015). For example, the gas content of five reservoirs in the USA, with in situ pressures between 2 and 28 MPa (600–8500 m deep), ranged from 1 to 10 m^{3}/ton, of which 20–85% was adsorbed in the organic content (Curtis 2002; Hill and Nelson 2000).

**Fig. 1 a** Fractures in shale core samples (Gale et al. 2014), b a subvertical calcite-cemented fracture (Soeder 1988), c simplified dual porosity concept for shale

Due to the multi-scale and multiple physical phenomena associated with shale gas production, the mechanism of gas flow is a combination of viscous (Darcy and non-Darcy) flow (Huang et al. 2016), and gas-slippage (Klinkenberg) effect (Mehmani et al. 2013), as well as continuum (Fickian), free-molecule (Knudsen), and surface (Yuan et al. 2014) diffusions. On the other hand, the gas sorption (adsorption and/or desorption) also plays an important role in the gas transport mechanism in shale reservoirs, given that the amount of adsorbed gas in the organic content is comparable with the free gas.

Despite the large number of studies on primary and CO_{2}-EGR, many aspects of shale gas development are still not well understood, mainly due to the complex interactions of fluid flow with different thermal, chemical, and mechanical phenomena. For instance, withdrawal of water and gas from a subsurface formation results in significant reservoir pressure alteration, which in turn leads to a redistribution of stress and strain within the shale formation.

Gas and water sorption-induced swelling/shrinkage (Lyu et al. 2015; Yuan et al. 2014), thermal phenomena (e.g. during thermal stimulation (Carroll et al. 2011)), and chemical reactions of water or CO_{2} (in case of CO2-EGR) with shale (Carroll et al. 2011) can additionally impact the stress–strain relations. Redistribution of stress and strain within any reservoir can subsequently lead to other changes, such as reservoir permeability alterations, subsidence/uplift of the ground surface, fault-reactivation mechanisms, and crack development in the caprock (Masoudian et al. 2016a).

Prediction of permeability has received considerable attention due to its significant effect on gas production from reservoirs. Numerous studies have related the permeability of gas shales to changes in the stress level (e.g. Amann-Hildenbrand et al. 2012; Chen et al. 2015a; Cho et al. 2013; Li et al. 2016; Spencer 1989; Zhou et al. 2016). However, these studies have only considered the elastic behaviour of the rock. Plastic deformations within the reservoir rock can have significant implications for production, injectivity, and stability of the wellbore.

Studies have demonstrated the significance of both elastic and plastic deformations on wellbore producibility and/or injectivity as well as permeability prediction in other types of reservoirs (e.g. Cui et al. 2007; Han and Dusseault 2003; Masoudian et al. 2016b). In addition, the stability of wellbores is an important issue in any oil and gas production project. Therefore, improved models are needed to properly estimate the deformation and stress distributions around the wellbore. Experimental studies have also revealed the brittle mechanism of failure in shale (e.g. Amann et al. 2011; Hull et al. 2015; Sone and Zoback 2013), which implies the need for improved models that consider this failure mechanism. The effect of gas and/or water sorption-induced swelling should also be taken into account when developing permeability models.

The most common geomechanical assumption when predicting permeability in gas reservoirs is of one-dimensional elastic deformation (vertical deformation). However, this assumption cannot provide an accurate estimation, especially near the wellbore where plastic deformations may occur. On the other hand, using fully numerical hydro-mechanical models leads to substantially more expensive computations. This is especially important in shale gas reservoirs where complex multiple porosity reservoir models are utilised.

This study aims to extend previous work (Masoudian and Hashemi 2016) by further developing brittle elastoplastic models for reservoirs coupled with sorption-induced swelling/shrinkage effects. In this paper, models are developed with boundary conditions that simulate different well completions and geological settings (e.g. open-hole and cased-hole completions). The models are then used to estimate the stress distribution within the reservoir, with which the changes in fracture permeability within both elastic and plastic zones are predicted. This is an important aspect of the work that has been largely neglected in other studies, mainly due to lack of computational resources for field scale fully coupled hydro-elastoplastic simulations. This work aims to highlight the role of plastic deformations in the prediction of permeability and its impact on gas production.

## 2 Theoretical Framework and Governing Equations

The reservoir is assumed to be a disc-shaped homogenous isotropic continuum. A vertical wellbore is assumed to be drilled at the centre of the disc, through which natural gas can be produced. It should be noted that these simplifying assumptions are employed in order to develop a fully analytical framework in this paper. Clearly, these assumptions cannot fully represent the heterogeneous nature of shale; however, they do not prevent evaluation of the effect of elastoplastic deformations on stress and permeability distributions within reservoirs, which is the focus of this paper. Here, a simplified multi-continuum approach is adopted, where the rock is considered as a continuum at the field scale, while the fractured nature of the rock is captured by the effect of stress on fracture porosity.

**Fig. 2** Schematics of a the geometry of the model, the multi-continuum approach, and the swelling/shrinkage effect on matrix size, and b the stress and strain relationships in the adopted elastoplastic framework

Figure 2a illustrates the geometry of the model and the multi-continuum approach adopted in this study. Figure 2b depicts the stress–strain relations within the brittle elastoplastic framework used in this paper. Note that the mechanical behaviour of brittle rock is idealised with a sudden drop of strength from peak to residual value. Similar to perfectly plastic behaviour, there is no further change in stress beyond the yield point as plastic strains continue to accumulate. The governing equations employed in this paper are explained below.

Starting with the classical approach in plasticity, strains can be decomposed into elastic (ε_{e}) and plastic (ε_{p}) components, where the elastic strain can be further divided into elastic mechanical strain, (ε_{e},m), and elastic swelling (or shrinkage) strain, (ε_{e},s). Note that adsorption induces swelling while desorption leads to shrinkage of the matrix blocks. Assuming axial symmetry of the problem and considering a cylindrical system of coordinates (r, θ, z for radial, tangential and vertical directions, respectively), the strain compatibility equation iswhere r is the radial distance from the centre of a disc-shaped reservoir, and subscripts rr and θθ represent the radial and tangential coordinates, respectively. In the light of Biot’s definition of effective stress (σ′=σ−αP), the equation of equilibrium can be written as

where P is pore pressure, α is Biot’s coefficient, and σ′ is effective stress. Assuming that swelling strain is analogous to thermal strain, the elastic stress–strain constitutive relations for an isotropic rock can be written as

where *E*, *ν* and *ε*^{e,s}_{v} are the elastic modulus, Poisson’s ratio, and volumetric swelling strain, respectively. The linear Mohr-Coulomb yield criterion can be written as

where σ′_{θθ} and σ′_{rr} are the tangential and radial effective stresses, respectively. Poisson’s ratio, and volumetric swelling strain, respectively. The linear Mohr-Coulomb yield criterion can be written as

_{r}, and friction angle, ϕ

_{r}, into these equations, the residual parameters γ

_{r}and Y

_{r}can be defined for failed brittle rock. Considering a non-associated flow rule, the plastic components of radial and tangential strains are related as (Park and Kim 2006)

where u is the radial displacement, and β is a function of dilation angle, ψ, as

There are a large number of studies that have investigated modelling approaches for flow of gas (methane or carbon dioxide) in shale reservoirs, where the complex multi-porous nature of shale is carefully considered and different mass transport mechanisms (i.e. viscous and diffusive) are taken into account. However, for the comparative analyses in this paper, a simplified pressure solution is sufficient. To arrive at a simplified pressure solution, the reservoir is initially assumed to be under a uniform pore pressure regime. Gas production from the wellbore with a fixed pressure (p_{w}) leads to the development of a non-uniform pore pressure distribution, with pressure remaining equal to the initial pressure at the outer boundary only. Thus, a steady-state solution of radial flow can be found in the following logarithmic form (Cui et al. 2007)

_{0}and R

_{0}are the initial reservoir pressure and radius of the disc-shaped reservoir, respectively. The term Q is defined as

where *R*_{w} and *p*_{w} are the radius and pressure of the wellbore, respectively.

The volume of adsorbed gas, *V*s , can be estimated using Langmuir’s isotherm as a function of pressure (*P*)

where V^{s}_{L} is the maximum sorption capacity and b_{L} is the Langmuir isotherm constant. Initially, the matrix of the reservoir rock contains a specified amount of adsorbed gas at the initial pressure (p_{0}) which changes as the reservoir pressure decreases due to production. As a result, the volumetric swelling strain changes from the initial state of ε_{e},s_{ν},0 to ε^{e,s}_{ν} and is considered to have a linear relationship with adsorbed gas volume, which gives

where ε^{s}_{L} is the maximum swelling strain that can only take place when the gas within the matrix is equal to the sorption capacity. Note that ε^{e,s}_{ν,0} is used to account for the in situ equilibrium state of the reservoir and can be estimated directly using the equation above, replacing the pressure with p_{0}.

There are a number of fracture porosity/permeability models developed for reservoirs where adsorption-induced swelling or shrinkage is associated with gas injection or production (e.g. Gensterblum et al. 2015; Pan and Connell 2012). However, these models adopt two main approaches for porosity/permeability modelling: dependency on stress or strain. Both of these approaches are closely related and may lead to very similar results under certain conditions (Palmer 2009), as they both consider a cubic law to link permeability to fracture porosity (consistent with a Kozeny–Carman type relationship). In this paper, results are obtained with the most widely used approach where the permeability is related to the change in horizontal effective stress using an exponential equation as (Shi and Durucan 2005)

where k is permeability, the subscript 0 denotes an initial value, c_{f} is fracture compressibility, and Δσ′_{h} is the change in horizontal effective stress, which can be defined as the mean of the change in the radial and tangential components of effective stress. As illustrated in Fig. 2a, any element of rock (i.e. at a specific value of r) can be conceptualised as a hybrid of fractures and matrix blocks and any change in stress state (due to pore pressure reduction and swelling/shrinkage effects) results in a change in the aperture size of the fractures within that element. Equation 12 therefore provides an efficient method for relating the response of fractures and their permeability at fracture scale to the changes in effective stress at field scale. The use of mean horizontal effective stress implies a matchstick matrix model, which is consistent with the plane-strain assumption.

Closed-form solutions of permeability and porosity can then be obtained by substituting the stress and strain equations for different loading and boundary conditions into the above equations. It should be noted, however, that experimental results have shown that the relationship between post-failure permeability and the stress–strain state can be much more complex than what these models offer (e.g. Carey et al. 2015). However, due to the lack of dedicated post-failure models, this widely used simplified approach was adopted in this paper. Limitations of the model are discussed later in the paper.

## 3 Analytical Solutions for Elastic–Brittle–Plastic Rock

Following the approach employed by Masoudian and Hashemi (2016), analytical solutions for stress and strain around a wellbore can be found by assuming two distinct concentric zones; a plastic zone near the wellbore and an elastic zone towards the outer boundary. The derivation and development of these solutions are partially presented in “Appendix 1”. The interface of elastic and plastic zones is at a radial distance of R_{ep} from the wellbore centre. The general solution within the elastic zone (r>R_{ep}) can be written as

where C_{1} and C_{2} are the integration constants that can be determined from the boundary conditions, and

where E_{i}(x) is the exponential integral function.

The stress and strain relations within the plastic zone near the wellbore are governed by a Mohr–Coulomb failure envelope and a non-associated flow rule. Note that the brittle failure mechanism is implemented using residual strength parameters (i.e. cohesion and friction angle); perfectly plastic behaviour could be modelled by using peak strength values. The general solution in the plastic zone can then be written as

where Y_{r} and γ_{r} are obtained with residual values of cohesion and friction angle in Eq. 5, and uep is the displacement at the elastic–plastic interface, which can be found considering continuity of displacement at the interface, i.e. by replacing the value of Rep in the elastic solution of displacement (Eq. 13) once the value Rep is found. To obtain Rep, the integration constants C_{1}, C_{2}, and C_{3} are found by substituting different sets of boundary conditions into the general solution, as presented in Table 1. Each set of boundary conditions describes a particular geological and well completion condition.

**Table 1** Integration constants for different sets of boundary conditions

The constant stress at the wellbore represents an uncased hole; hence, the total stress applied to the inner boundary is equal to the bottomhole pressure (pw); i.e. zero effective stress on wellbore wall. On the other hand, zero displacement on the wellbore wall represents an ideally cased hole that does not allow the wellbore wall to converge. For the outer boundary, the constant stress represents an ideal scenario where the reservoir is completely separated from the adjacent formation. Zero displacement at the outer boundary represents a scenario where the reservoir is large enough to prevent deformation of the rock mass at the far-field. For example, if the reservoir was cut by low-angle faults at the outer boundary, the reservoir may be best simulated by applying a constant stress, whereas a zero horizontal displacement boundary would be more suitable if the reservoir was cut by steep faults and was in contact with mechanically strong geological strata.

However, due to the geological complexity of reservoirs, actual boundary condition will fall somewhere between these two conditions; hence, different boundary conditions should be tested. The boundary condition acronyms consist of four letters, with the first pair describing the condition of the wellbore wall, and the second pair describing the condition of the outer boundary. As such, case CSCS represents constant stresses (CS) on both the wellbore wall and outer boundary, case ZDCS represents zero displacement (ZD) of the wellbore wall and constant stress on the outer boundary, case CSZD represents constant stress on the wellbore wall and zero displacement of the outer boundary, and case ZDZD represents zero displacement on both the wellbore wall and outer boundary.

The interface between elastic and plastic zones is found considering the continuity of radial stress at r=R_{ep} and a trial-and-error scheme. To achieve this, an equality is constructed using the solution of radial stress within both elastic and plastic zones, and then its root for Rep can be found. If this equality does not have a root in the [R_{w},R_{0}] domain, it means that the rock mass does not include a plastic zone and hence a fully elastic solution is used by replacing Rep by R_{w }in Eq. 13 with the corresponding elastic integration constants presented in “Appendix 2”. These formulations and the introduced methodology were implemented in a MATLAB code and used to obtain the results discussed in the following section.

## 4 Results and Discussion

In order to examine the stress change and permeability evolution around a gas production well, a set of values were chosen for the model input parameters, as listed in Table 2. These values mainly represent gas shale reservoirs (suitable for hydraulic fracturing due to their brittleness) and were chosen to be within the reported ranges from the corresponding references.

**Table 2** Input parameters for the base simulation case