Large poroelasto-plastic deformations due to radially outward fluid injection

TL;DR

This paper develops a steady-state model for large poroelasto-plastic deformations caused by fluid injection into a borehole, highlighting the effects of permeability and plasticity on deformation and stress.

Contribution

It introduces a novel model accounting for pre-existing stresses, permeability transition, and compares poroelastic and poroelasto-plastic behaviors in geomechanical applications.

Findings

Plasticity enables larger deformations with smaller stresses.

Deformations and stresses increase as permeability transitions from impermeable to permeable.

Elastic region remains insensitive to model nonlinearities.

Abstract

Flow-induced failure of granular materials is relevant to a broad range of geomechanical applications. Plasticity, which is the inherent failure mechanism of most granular materials, enables large deformations that can invalidate linearised models. Motivated by fluid injection into a borehole, we develop a steady-state model for the large deformation of a thick-walled, partially-permeable, elastic--perfectly-plastic annulus with a pressurised inner cavity. We account for pre-existing compressive stresses, as would be present in the subsurface, by subtracting a compressed initial state from our solutions to provide the additional disturbance due to fluid injection. We also introduce a simple parameter that allows for a smooth transition from an impermeable material (i.e., subject to mechanical loading at the inner wall) to a fully permeable material (i.e., subject to an internal pore-pressure gradient), which would be relevant to coated boreholes and very-low-permeability materials. We focus on the difference between poroelastic and poroelasto-plastic deformations, the role of kinematic and constitutive nonlinearity, and the transition from impermeable to fully permeable. We find that plasticity can enable much larger deformations while predicting much smaller stresses. The former makes model choice increasingly important in the plastic region, while the elastic region remains insensitive to these changes. We also find that, for a fixed total stress at the inner wall, materials experience larger deformations and generally larger stresses as they transition from impermeable to fully permeable.

Equations

$$\tilde{\Gamma }:=\frac{\Lambda }{\mathcal{M}},{\tilde{a}}^{\mathrm{ref}}:=\frac{a^{\mathrm{ref}}}{b^{\mathrm{ref}}},{\tilde{\sigma }}_b^{\prime }:=\frac{{{\sigma }_b}^{\prime }}{\mathcal{M}},\tilde{c}:=\frac{c}{\mathcal{M}}\text{and either}\tilde{q}:=\frac{\mu Q}{2\pi k\mathcal{M}}\text{or}{\tilde{\sigma }}_a:=\frac{{\sigma }_a}{\mathcal{M}},$$

$$u_s(r)=r-R(r),$$

$$\frac{{\varphi }_f-{\varphi }_f^{\mathrm{ref}}}{1-{\varphi }_f^{\mathrm{ref}}}=\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left(r{u}_s-\frac{1}{2}{u}_s^2\right),$$

$$\frac{{\varphi }_f-{\varphi }_f^{\mathrm{ref}}}{1-{\varphi }_f^{\mathrm{ref}}}\approx \frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left(r{u}_s\right).$$

$$\frac{\mathrm{d}{\sigma }_r^{\prime }}{\mathrm{d}r}+\frac{{\sigma }_r^{\prime }-{\sigma }_{\theta }^{\prime }}{r}=\frac{\mathrm{d}p}{\mathrm{d}r},$$

$$\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left(r{\varphi }_f{v}_f\right)=0⟹{\varphi }_f{v}_f=\frac{q}{r},$$

$${\varphi }_f{v}_f=-\frac{\mathrm{d}p}{\mathrm{d}r}⟹\frac{\mathrm{d}p}{\mathrm{d}r}=-\frac{q}{r}.$$

$$\Delta p:=p(a)-p(b)=q\ln \left(\frac{b}{a}\right).$$

$${\sigma }_r^{\prime }={\epsilon }_r+\Gamma {\epsilon }_{\theta },{\sigma }_{\theta }^{\prime }=\Gamma {\epsilon }_r+{\epsilon }_{\theta },\mathrm{and}{\sigma }_z^{\prime }=\Gamma ({\epsilon }_r+{\epsilon }_{\theta })=\frac{\Gamma ({\sigma }_r^{\prime }+{\sigma }_{\theta }^{\prime })}{1+\Gamma },$$

$${\epsilon }_r=\frac{\mathrm{d}{u}_s}{\mathrm{d}r},{\epsilon }_{\theta }=\frac{u_s}{r},\mathrm{and}{\epsilon }_z\equiv 0,$$

$${\mathcal{F}}_{1,3}:=\alpha {\sigma }_1^{\prime }-{\sigma }_3^{\prime }-y,$$

$$\alpha :=\frac{1+\sin (\phi )}{1-\sin (\phi )}\text{and}y:=\frac{2c\cos (\phi )}{1-\sin (\phi )}.$$

$${\mathcal{F}}_{\theta ,r}:=\alpha {{\sigma }^{\prime }}_{\theta }-{{\sigma }^{\prime }}_r-y.$$

$$\epsilon ={\epsilon }^{\mathrm{e}}+{\epsilon }^{\mathrm{p}},$$

$$\beta {\dot{\epsilon }}_r^{\mathrm{p}}+{\dot{\epsilon }}_{\theta }^{\mathrm{p}}=0,$$

$$\beta :=\frac{1+\sin (\psi )}{1-\sin (\psi )}.$$

$${\sigma }_r^{\prime }(b)={\sigma }_b^{\prime }.$$

$$u_s(b)=b-1.$$

$$u_s(s^{+})=u_s(s^{-})\mathrm{and}{\sigma }_r^{\prime }(s^{+})={\sigma }_r^{\prime }(s^{-}),$$

$$\alpha {\sigma }_{\theta }^{\prime }(s^{+})-{\sigma }_r^{\prime }(s^{+})=y.$$

$${\sigma }_r^{\prime }(a)=(1-\zeta ){\sigma }_a$$

$$p(a)\equiv \Delta p=-\zeta {\sigma }_a,$$

$$q=-\frac{\zeta {\sigma }_a}{\ln \left(\frac{b}{a}\right)}.$$

$$u_s(a)=a-a^{\mathrm{ref}},$$

$$\frac{{\mathrm{d}}^2{u}_{s,0}}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}{u}_{s,0}}{\mathrm{d}r}-\frac{u_{s,0}}{r^2}=0$$

$${\sigma }_{r,0}^{\prime }(b)={\sigma }_b^{\prime }$$

$$u_{s,0}(r)=\left(\frac{{\sigma }_b^{\prime }}{1+\Gamma }\right)r,$$

$${\varphi }_{f,0}^Q={\varphi }_f^{\mathrm{ref}}+\frac{{\sigma }_b^{\prime }(1-{\varphi }_f^{\mathrm{ref}})[2(1+\Gamma )-{\sigma }_b^{\prime }]}{(1+\Gamma {)}^2}.$$

$$u_{s,0}^Q=\left(\frac{{\sigma }_b^{\prime }}{1+\Gamma -{\sigma }_b^{\prime }}\right)R.$$

$$a_0^Q=\left(\frac{1+\Gamma }{1+\Gamma -{\sigma }_b^{\prime }}\right)a^{\mathrm{ref}}.$$

$${\varphi }_{f,0}^L={\varphi }_f^{\mathrm{ref}}+\frac{2(1-{\varphi }_f^{\mathrm{ref}}){\sigma }_b^{\prime }}{1+\Gamma },u_{s,0}^L(R)=\left(\frac{{\sigma }_b^{\prime }}{1+\Gamma }\right)R,\text{and}{a}_{0,L}=\left(\frac{{\sigma }_b^{\prime }+1+\Gamma }{1+\Gamma }\right)a^{\mathrm{ref}}.$$

$$\frac{{\mathrm{d}}^2{u}_s}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}{u}_s}{\mathrm{d}r}-\frac{u_s}{r^2}=\frac{\zeta {\sigma }_a}{r\ln \left(\frac{b}{a}\right)}s<r<b.$$

Full text

I Introduction

A clear understanding of the failure of granular materials due to fluid injection has direct relevance to a variety of applications in geomechanics, including borehole stability, cavity expansion, and ‘fracking’ for the recovery of oil or natural gas from shales. These problems involve radially outward flow and deformation surrounding a long cylindrical hole through soil or rock. Most soils and shallow sedimentary rocks have a granular microstructure, and are intrinsically porous, water-saturated, and susceptible to plastic failure. These features suggest that these problems should be approached within the framework of large-deformation poroelasto-plasticity, but essentially all previous work has neglected at least one of these three ingredients (large deformations, permeability, or plasticity) and no previous study has assessed their relative importance.

A borehole is a long cylindrical cavity in the subsurface, from which the original material has been removed by drilling. The walls of the borehole can be supported by a metal casing (with or without perforations), or can be uncased and self-supporting. The stability of a borehole depends on many factors, including the state of stress in the material before drilling, the orientation of the borehole, and the loading to which the borehole is subjected. Problems in borehole integrity are concerned with predicting and preventing borehole collapse after drilling, or during subsequent operations (e.g., during fluid injection or extraction). Borehole integrity has been studied extensively (e.g, Wang, 1996; Wang and Dusseault, 1994; Detournay and Cheng, 1988; Risnes et al., 1982; Wang and Dusseault, 1991a, b; Wang et al., 1994; Detournay and Fairhurst, 1987), but exclusively under the assumption of infinitesimal deformations. Detournay and Fairhurst (1987) and Wang and Dusseault (1991a) neglect the permeability of the material while allowing for plastic failure, whereas Detournay and Cheng (1988) consider nonzero permeability while neglecting plasticity.

In cavity expansion, a borehole-like cavity is created via radially-outward mechanical displacement of material (typically soil) around an insertion point. Cavity expansion is a classical problem in geotechnical engineering, with relevance to pile-driving and penetrometer testing (Vesic, 1972; Carter et al., 1986; Yu and Houlsby, 1991; Yu, 2000; Davis and Selvadurai, 2002; Howell et al., 2009). Large plastic deformations are central to this process; hence, all of these studies consider rigorous large-deformation elasto-plasticity. However, for simplicity, they focus on drained or undrained limiting behaviours and thus neglect the role of fluid flow.

In fracking, fluid is injected into hydrocarbon-bearing rocks, usually shales, in order to open fractures around the injection point. These fractures provide hydraulic access deeper into the reservoir and allow gas to be collected from a larger region of the rock (Pye, 1973; Economides and Nolte, 2000; Haimson and Fairhurst, 1969). Fracking relies on the brittle failure of shale, but the mechanical properties of shale depend strongly on the composition (Economides and Nolte, 2000; Britt and Schoeffler, 2009; Rickman et al., 2008); many hydrocarbon-bearing shales have high clay content and may therefore have non-negligible ductility (Daigle et al., 2014; Swift et al., 2014; Vega et al., 2014; Vallejo, 1988). Both laboratory and field data suggest that the ductility of shales can have an important impact on the success of fracking (Haimson and Fairhurst, 1969; Britt and Schoeffler, 2009), but the ductility of shales is almost always neglected. Most studies also neglect fluid flow through the shale due to its extremely low permeability (e.g., Wang and Dusseault, 1991a; Wang et al., 1994; Detournay, 2004; Hubbert and Willis, 1957).

Here, motivated by these problems, we develop a kinematically rigorous, steady-state model for the large deformation of a porous, thick-walled, elastic–perfectly-plastic annulus with a pressurised inner cavity. To explore the role of fluid flow (i.e., the importance of permeability), we introduce a new parameter that allows us to transition smoothly from an impermeable material to a fully permeable material. This parameter is physically analogous to the presence of a thin, weak, low-permeability ‘skin’ near the cavity wall. In practise, such regions form naturally in boreholes when drilling and/or injection push fine grains into the pore space, or artificially due to the injection of ‘wall-building’ chemicals. The latter enables a fixed volume of fluid to support the walls of the borehole via hydrostatic pressure. To account for pre-existing stresses, as would be present around a borehole in the subsurface, we consider the impact of fluid injection relative to a pre-stressed initial state. We show that plastic failure enables large deformations that make it essential to account for rigorous, nonlinear kinematics in the plastic region. The elastic region is insensitive to these changes. We also show that, for a given applied total stress at the inner radius, a fully permeable material deforms much more than an impermeable one.

II Theoretical Model

We consider a simple model for the fluid-driven deformation of a thick-walled annulus with a pressurised inner cavity. We take the annulus to be made of a homogeneous cohesive granular material, having relaxed outer radius $b^{\mathrm{ref}}$ and relaxed inner radius $a^{\mathrm{ref}}$ (I in Figure 1). We focus on the final steady state during fluid injection; throughout, we assume axisymmetry.

II.1 Pre-stressed initial state

The subsurface exists in a state of compressive tectonic stress, the principal values of which typically vary from each other by less than one order of magnitude Hubbert and Willis (1957); Bazant et al. (2014). To mimic this compressed initial state, we consider a solid cylinder in plane strain in the plane orthogonal to its axis, with the two in-plane effective stresses equal and given by $\sigma_{b}^{\prime}\leq 0$. Imposing this stress on the relaxed cylinder decreases its outer radius, $b^{\mathrm{ref}}\mapsto{}b_{0}\leq{}b^{\mathrm{ref}}$ (I${}\mapsto{}$II in Figure 1); we assume that this compression is purely elastic.

We then assume that a cylinder of material of radius $a_{0}$ (relaxed radius $a^{\mathrm{ref}}$) is removed from this compressed state, and that the resulting cavity is supported with a permeable “casing” that preserves the radius $a_{0}$, such that the stresses in the remaining material are unchanged (II${}\mapsto{}$III in Figure 1). We assume that this casing has no impact on fluid exchange with the surrounding material, or on its outward deformation. Note that there will exist a minimum value of $q$ above which the casing is no longer supporting the cavity.

We then consider the elastic and subsequent elasto-plastic deformation of the resulting annulus due to radially outward fluid injection (III${}\mapsto{}$IV${}\mapsto{}$V in Figure 1). During injection, both radii expand such that $a_{0}\mapsto{}a>a_{0}$ and $b_{0}\mapsto{}b>b_{0}$. We assume that plastic yield initiates at the inner boundary and evolves outward to some intermediate radius $s$, such that the material has deformed elasto-plastically for $a\leq{}r<s$ and purely elastically for $s<r\leq{}b$, where $r$ is the radial coordinate. We assume that the effects of injection are localised around the cavity, such that the effective stresses and the pressure tend to their far-field values ($\sigma_{b}^{\prime}$ and [math], respectively) for finite $b$.

We seek the steady-state flow and deformation fields for this problem, for which the values of $a$, $s$, and $b$ are unknown a priori. We are interested in the deformation of the material relative to its compressed initial state (III in Figure 1), but we necessarily derive the steady state relative to the relaxed reference state (IV in Figure 1) due to the nonlinear nature of the model.

We next present our axisymmetric, plane-strain model for this problem in dimensionless form. We present our scalings and identify key dimensionless parameters in §II.2. We summarise the key aspects of the model, including kinematics, Darcy’s law, and mechanical equilibrium, in §II.3. We then discuss constitutive laws for the solid skeleton §II.4.

II.2 Scaling

To write the model in dimensionless form, we adopt characteristic scales for length, stress/pressure, and permeability. We take the dimensional relaxed (reference) outer radius ${b}^{\mathrm{ref}}$ as the characteristic length scale and the $p$-wave (oedometric) modulus ${\mathcal{M}}$ as the characteristic stress/pressure scale. We model fluid injection as a line source at the origin, characterised by either a fixed dimensionless flow rate $\tilde{q}$ or an applied dimensionless total stress at the inner cavity wall, $\tilde{\sigma}_{a}$. The dimensionless model is then characterised by the friction angle $\tilde{\varphi}$ and dilation angle $\tilde{\psi}$ (see §II.4.2), the reference (relaxed) porosity $\tilde{\phi}_{f}^{\mathrm{ref}}$, which we take to be uniform for simplicity, and five other dimensionless parameters:

[TABLE]

where ${\Lambda}$ is Lamé’s first parameter, ${k}$ is the permeability (assumed to be constant and uniform), ${\mu}$ is the dynamic viscosity of the fluid, ${{\sigma}}^{\prime}_{b}$ is the radial effective stress (see §II.3) at the outer boundary, ${c}$ is the cohesion between grains of the solid skeleton (see §II.4.2) and ${Q}$ is the volume injection rate per unit length into the page. Note that the dimensionless relaxed outer radius is $\tilde{b}^{\mathrm{ref}}\equiv 1$, and also that we take tension to be positive. Going forward, we work almost entirely in terms of dimensionless quantities and thus drop the tildes ($\tilde{\star}$) for convenience; we denote any further dimensional quantities with a breve ($\breve{\star}$).

II.3 Kinematics, Darcy’s law and mechanical equilibrium

Axisymmetry implies that the fluid velocity $\bm{v}_{f}$ and solid displacement $\bm{u}_{s}$ each have only one nontrivial component, such that $\bm{v}_{f}=v_{f}(r)\bm{\hat{e}}_{r}$ and $\bm{u}_{s}=u_{s}(r)\bm{\hat{e}}_{r}$, where $r$ is the Eulerian radial coordinate and $\bm{\hat{e}}_{r}$ is the radial unit vector. Steady state implies that $v_{f}$ and $u_{s}$ are independent of time. We work in an Eulerian (spatial) reference frame, such that the displacement is given by

[TABLE]

where the Lagrangian radial coordinate $R(r)$ denotes the original position of the material that is at position $r$ in the deformed state.

We assume that both the individual solid grains and the fluid are incompressible, such that deformation occurs through rearrangement of the grains and corresponding changes in the local porosity or void fraction $\phi_{f}$. We relate $u_{s}$ to $\phi_{f}$ via (Auton and MacMinn, 2017)

[TABLE]

which linearises under the assumption of infinitesimal strains to

[TABLE]

Additionally, mechanical equilibrium requires that $\nabla\cdot\bm{\sigma}=0$, where $\bm{\sigma}$ is the total Cauchy stress and where we have neglected inertia and the effect of gravity. The total stress can be decomposed as $\bm{\sigma}=\bm{\sigma}^{\prime}-p\bm{I}$, where Terzaghi’s effective stress $\bm{\sigma}^{\prime}$ is the portion of the total stress supported by deformation of the solid and $p$ is the fluid pressure. Combining these ideas for axisymmetric flow and deformation leads to

[TABLE]

where $\sigma^{\prime}_{r}$ and $\sigma^{\prime}_{\theta}$ are the radial and azimuthal (hoop) components of the effective stress, respectively. At steady state, conservation of mass for the fluid leads to

[TABLE]

by the definition of a line source of strength $q$.

Finally, we assume that the fluid flows through the solid according to Darcy’s law. For simplicity, we assume that the permeability remains constant. For steady axisymmetric flow, this leads to

[TABLE]

Combining Equations (6) and (7) and integrating leads to a direct relationship between $q$ and the pressure drop $\Delta p$:

[TABLE]

For a more detailed derivation and discussion of the above aspects of the model, see Auton and MacMinn (2017, 2018).

II.4 Elasticity and plasticity

II.4.1 Elasticity

Elastic deformations are quasi-static and reversible. Most ductile materials will yield before experiencing moderate or large elastic deformations, so we restrict our attention to infinitesimal deformations and therefore linear elasticity in the elastic region. In linear elasticity, the effective stress is related to the strain via

[TABLE]

where $\varepsilon_{r}$, $\varepsilon_{\theta}$, and $\varepsilon_{z}$ are the radial, azimuthal, and axial components of the strain, respectively, and the expressions for $\sigma_{z}^{\prime}$ and $\varepsilon_{z}$ are consequences of plane strain.

II.4.2 Plasticity

Plastic deformations are path-dependent and irreversible. Plastic failure implies that, once the state of stress exceeds a threshold, some fraction of any additional strain energy will be dissipated through irreversible rearrangements as opposed to being stored elastically. The simplest form of plasticity is ‘perfect plasticity’, in which the material properties are assumed to remain constant after yield and all additional strain energy is dissipated Hill (1950). The threshold that defines the transition from elastic to plastic behaviour is known as a yield condition. For granular materials, this is typically based on the idea of internal Coulomb-like friction. This implies that the effective shear stress $\tau^{\prime}$ anywhere within the material must be strictly less than a specified fraction of the corresponding effective normal stress $\sigma^{\prime}$ for the deformation to remain elastic. For perfect plasticity, $\tau^{\prime}$ is enforced to remain equal to this fraction of $\sigma^{\prime}$ after yield. For a cohesive granular material, it is standard to additionally include a yield strength $c$ due to the cohesion between the grains. The simplest cohesive-frictional yield condition is the cohesive Mohr-Coulomb condition, $|\tau^{\prime}|\leq-\sigma^{\prime}\tan\varphi+c$, where $\varphi$ is the friction angle Davis and Selvadurai (2002). We then express this condition in terms of the principal effective stresses111The principal effective stresses are the eigenvalues of the effective stress tensor, with $\sigma^{\prime}_{1}$ and $\sigma^{\prime}_{3}$ being largest (most tensile) and the smallest (least tensile), respectively. $\sigma_{1}^{\prime}\geq\sigma_{2}^{\prime}\geq\sigma_{3}^{\prime}$ by introducing a yield function $\mathcal{F}_{1,3}$,

[TABLE]

where

[TABLE]

Plastic yield now corresponds to the condition $\mathcal{F}_{1,3}=0$. For $\mathcal{F}_{1,3}<0$, the material remains elastic. As soon as equality is first achieved, we enforce $\mathcal{F}_{1,3}\equiv 0$ locally thereafter. Taking $\varphi=0$ ($\alpha=1$) provides the Tresca yield condition, which is commonly used to model cohesive soils during undrained deformation (i.e., at fixed pore volume) because the associated plastic flow law is volume conservative (see below). The simplicity of the Tresca model enables analytical solutions in many cases (Auton, 2018). Additionally, these solutions can be used to derive the corresponding solutions for a von Mises material. The von Mises yield condition is commonly used to model multi-axial loading in metals.

Here, the three principal stresses are $\sigma^{\prime}_{r}$, $\sigma^{\prime}_{\theta}$, and $\sigma^{\prime}_{z}$, but not necessarily in this order. To apply the correct yield condition, it is necessary to determine the correct ordering of these stresses. For a cylinder in plane strain, it is commonly assumed that $\sigma^{\prime}_{1}=\sigma^{\prime}_{\theta}$ and $\sigma^{\prime}_{3}=\sigma^{\prime}_{r}$ (i.e., $\sigma_{\theta}^{\prime}>\sigma_{z}^{\prime}>\sigma_{r}^{\prime}$) Wang and Dusseault (1991a, b), such that the appropriate yield function is

[TABLE]

This ordering has been justified heuristically for an impermeable cylinder Yu and Houlsby (1991). We show in Appendix A that Equation (12) is indeed the appropriate yield function if the cylinder is sufficiently “strong”—that is, if $y(b^{2}-a^{2})>-2b^{2}\sigma_{b}^{\prime}$ and $\alpha>(1+\Gamma)/\Gamma$, assuming that $\mathcal{F}_{1,3}$ is maximised at $r=a$ and at steady state, and that, once the material yields according to a given yield condition, the material subsequently yields exclusively according to that condition. These assumptions are consistent with previous work Wang and Dusseault (1991b, a); Yu and Houlsby (1991).

We assume that the strain everywhere can be additively decomposed into elastic and plastic components Davis and Selvadurai (2002),

[TABLE]

where the superscripts $\mathrm{e}$ and $\mathrm{p}$ denote the elastic and plastic components of the strain, respectively. Stored elastic energy is associated with the elastic component of the strain, which is related directly to the effective stress via the elasticity law. Dissipation is associated with the plastic component of the strain, which has no direct connection to the effective stress; instead, the evolution of the plastic strain is described by a plastic flow law. Note that, whereas the strain decomposes into elastic and plastic components, there is one unique stress field.

For a Mohr-Coulomb yield condition, the so-called non-associated flow law is given by

[TABLE]

where the overdots denote the material time derivative following the solid and

[TABLE]

A nonzero dilation angle $\psi$ incorporates the fact that most granular materials undergo volumetric expansion (dilation) during plastic flow. For $\psi=\varphi$ ($\beta=\alpha$), this flow law is said to be ‘associated’; however, associated flow typically overestimates dilation for granular materials Yu (2000); Zhang and Salgado (2010). For $\psi=0$, this flow law reduces to the associated flow law for the Tresca yield condition and becomes volume conservative (no dilation). For the scenarios considered here, we expect that $0\leq\psi\leq\varphi$ ($1\leq\beta\leq\alpha$).

II.5 Boundary conditions

II.5.1 Outer boundary

[TABLE]

Note that we assume that $a\leq s<b$, so that the cylinder is not entirely plastic and hence the outer boundary conditions will always be applied to the elastic region. Note that the problem becomes over-constrained, and thus ill-defined, once $s$ reaches $b$.

II.5.2 The elastic-plastic interface

[TABLE]

At the elastic-plastic interface, we enforce continuity of displacement and of radial effective stress,

[TABLE]

where $s^{-}$ denotes $r=s$ as approached from within the plastic region $(a\leq r<s)$ and, likewise, $s^{+}$ denotes $r=s$ as approached from within the elastic region $(s<r\leq b)$. In addition, we require that the effective stresses in the elastic region must be at the point of yield at $r=s^{+}$,

[TABLE]

Note that Equations (17a) and (17b), and the fact that $\mathcal{F}_{\theta,r}\equiv 0$ throughout the plastic region, together imply continuity of $\sigma_{\theta}^{\prime}$ across $r=s$.

II.5.3 Inner boundary

We suppose that the cavity wall is coated by a thin, weak, low-permeability skin; for a borehole, this skin could be caused by local damage, clogging, or wall-building chemicals. We then suppose that the cavity is pressurised to a pressure $-\sigma_{a}$. The presence of the skin leads to a partitioning of this pressure between the fluid and the solid, such that the fluid pressure drops by some amount across the skin and the skin then exerts an effective stress on the solid at the cavity wall. For an impermeable skin, this would lead to an imposed radial effective stress $\sigma^{\prime}_{r}(a)=\sigma_{a}$; we refer to this as an ‘impermeable material’. For an ‘infinitely’ permeable skin, this would lead to an imposed fluid pressure $\Delta{p}\equiv{}p(a)=-\sigma_{a}$; we refer to this as a ‘fully permeable material’. To transition smoothly between these two limiting cases, we introduce a new ‘permeability-load parameter’ $\zeta$. We define $\zeta\in[0,1]$ such that the material can vary continuously from impermeable ($\zeta\equiv 0$) to fully permeable ($\zeta\equiv 1$). Pressurisation of the cavity then leads to two conditions at the inner boundary,

[TABLE]

We drive the system by imposing either $\sigma_{a}$ or the total flow rate $q$. The latter appears explicitly in Equation (7). We express $q$ in terms of $\sigma_{a}$ using Equations (8) and (18b),

[TABLE]

Note that, for injection, we expect $\sigma_{a}<0$. For fully permeable materials ($\zeta\equiv{}1$), in particular, it is convenient to impose $q$ rather than $\sigma_{a}$ (§V.2).

The inner boundary is also a material boundary. As such, the displacement must also satisfy a kinematic condition given by

[TABLE]

where $a^{\mathrm{ref}}$ is the relaxed reference position of the material that eventually comprises the cavity wall.

III Governing equations and model summary

III.1 Relaxed reference state to initial state

We now consider the transition from the relaxed reference state (I, Figure 1) to the initial state prior to pressurisation (III, Figure 1), which we take to be purely elastic compression under an imposed radial effective stress $\sigma^{\prime}_{b}$ in plane strain. As with $a$ and $b$, we denote quantities in the initial state by a subscript ‘0’.

Using Equations (9) to eliminate the stresses from Equation (5) in favour of the displacement, and further setting $q\equiv{}0$, purely mechanical linear-elastic deformation is governed by

[TABLE]

subject to

[TABLE]

and the requirement of boundedness at the origin. This problem has solution

[TABLE]

which gives $\sigma_{r,0}=\sigma_{\theta,0}=\sigma_{b}^{\prime}$, $\sigma_{z,0}=2\Gamma\sigma^{\prime}_{b}/(1+\Gamma)$, and $p_{0}=0$.

For a rigorous treatment of the kinematics, we impose the outer boundary conditions at $b_{0}$ rather than at $b^{\mathrm{ref}}$, and we calculate the initial porosity field according to Equation (3). We denote this initial state with a superscript ‘Q’ to indicate that it combines linear elasticity with rigorous kinematics (‘Quasi-linear’). The ‘Q’ initial porosity is given by

[TABLE]

We eliminate the Eulerian coordinate $r$ from Equation (21) in favour of the Lagrangian coordinate $R$ via the kinematic relationship $u_{s,0}(r)=r-R$,

[TABLE]

Finally, Equation (18d) leads to

[TABLE]

We also derive a fully linearised version of the initial state by imposing the outer boundary conditions at $b^{\mathrm{ref}}\equiv{}1$ and by calculating the initial porosity field according to Equation (4). We denote this initial state with a superscript ‘L’ to indicate that it is fully linearised under the assumption of infinitesimal strains (‘Linear’). The ‘L’ initial state is given by

[TABLE]

Note that Equation (18d) is not needed in this case.

III.2 Initial state to steady state

We now seek solutions to the poroelasto-plastic problem at steady state (V in Figure 1). We denote the displacement field in this state by $u_{s}$.

III.2.1 Elastic region

We again use Equations (9) to eliminate the stresses in Equation (5) in favour of the displacement, and we now also use Equation (18c) to eliminate $q$ in favour of $\sigma_{a}$. This leads to an ordinary differential equation (ODE) in $u_{s}$ in the elastic region,

[TABLE]

The domain of this second-order ODE has two free boundaries—a material boundary at $r=b$ and a constitutive boundary at $r=s$—and the solution will involve two constants of integration. We again define a ‘Q’ class of models subject to rigorous kinematics, where porosity is calculated according to Equation (3) and which requires four boundary conditions:

[TABLE]

We also again define an ‘L’ class of models that are fully linearised under the assumption of infinitesimal strains, in which we calculate the porosity according to Equation (4) and apply Equation (27c) at $b^{\mathrm{ref}}\equiv{}1$. Equation (27d) is not needed in this case.

III.2.2 Plastic region

In the plastic region, combining Equation (12) with Equation (5) and using Equation (18c) to eliminate $q$ in favour of $\sigma_{a}$ leads to

[TABLE]

which has solution

[TABLE]

where

[TABLE]

To determine the displacement, we integrate Equation (14) with respect to time to arrive at

[TABLE]

where $f$ is a quantity whose material time derivative following the solid must vanish. Equation (14) is a statement that the quantity $\beta\varepsilon_{r}^{\mathrm{p}}+\varepsilon_{\theta}^{\mathrm{p}}=f$ must be conserved (i.e., $f$ is a conservative tracer advected with the solid). Prior to fluid injection, there is no plastic strain in any direction ($\varepsilon_{r}^{\mathrm{p}}=\varepsilon_{\theta}^{\mathrm{p}}=0$) and therefore $f=0$; hence, it must be the case that $f\equiv 0$. We then decompose the plastic strains according to Equation (13),

[TABLE]

The elastic strains are always related to the effective stresses via the elasticity law. As such, Equations (9a) and (29) allow us to rewrite the right-hand side of Equation (32) as

[TABLE]

and $C_{1}$ and $C_{2}$ are as defined in Equations (30).

The plastic strains can, in general, be large. Hence, it may be appropriate to adopt a nonlinear constitutive model for the total strains in the plastic region. Here, we consider the Hencky (logarithmic) model Yu (2000); Bažant (1998). For Hencky strains, the left-hand side of Equation (32) becomes

[TABLE]

which, on combining Equation (34) with Equations (32) and (33a), leads to

[TABLE]

where $D_{1}$ and $D_{2}$ are as defined in Equation (33b). We denote this nonlinear model by ‘N’ (‘Nonlinear’).

For comparison, we also consider a linear constitutive model for the total strains in the plastic region,

[TABLE]

which linearises to Equation (36) for infinitesimal strains. Combining Equation (36) with Equations (32) and (33a) yields

[TABLE]

We denote this model by ‘Q’. The ‘Q’ model can be further linearised by taking $a\approx{}a^{\mathrm{ref}}$; we denote this fully linearised model by ‘L’.

III.2.3 Model summary

For rigorous kinematics throughout and Hencky strains in the plastic region, the full boundary value problem (BVP) is

[TABLE]

where $\sigma_{r}^{\prime}(s^{+})$, $\sigma_{\theta}^{\prime}(s^{+})$, and $\sigma_{r}^{\prime}(b)$ can be expressed in terms of $u_{s}$ via Equations (9), and $\sigma_{r}^{\prime}(s^{-})$ is given in Equation (29).

We denote this model by ‘NQ’, where the first letter refers to rigorous kinematics and nonlinear strains in the plastic region (‘Nonlinear’) and the second refers to rigorous kinematics and linear elasticity in the elastic region (‘Quasi-linear’). Note that Equation (18a) has already been imposed in the above. For comparison, we also consider linear strains throughout the entire domain by coupling Equation (37) with Equations (38b)—(38h). We denote this model ‘QQ’ to indicate rigorous kinematics and linear strains throughout the domain.

Although deformations and strains in the plastic region may be large, we expect deformations and strains in the elastic region to remain small. Thus, we further linearise the kinematics in the elastic region under the assumption of small deformations to provide an intermediate model defined by Equation (37) and Equations (38b)—(38g), with Equation (38f) applied at $b\approx{}b^{\mathrm{ref}}\equiv{}1$. We denote this model by ‘QL’, where the ‘L’ refers to a fully linearised model in the elastic region. This problem involves only two free boundaries (at $a$ and $s$).

Finally, we consider a fully linearised model, which we denote ‘LL’. This model comprises Equation (37) with $a\approx{}a^{\mathrm{ref}}$, Equations (38b)—(38e), and Equation (38f) evaluated at $b\approx{}b^{\mathrm{ref}}\equiv{}1$. This problem involves only one free boundary (at $s$) and can be solved analytically, with $s$ determined via an implicit relation that is straightforward to solve numerically using conventional root-finding techniques. This solution provides a good first guess for the numerical solution of the other three models (see §IV.4).

IV Solutions for imposed $\sigma_{a}$ and general $\zeta$

For a given set of parameters, our poroelasto-plastic models will not be well posed for all values of $\sigma_{a}$. Specifically, the material will remain elastic if $\sigma_{a}$ is sufficiently small; this purely poroelastic problem was studied extensively by Auton and MacMinn (2017, 2018) for fully permeable materials ($\zeta\equiv 1$). Conversely, the material will yield completely if $\sigma_{a}$ is sufficiently large. We denote the value of $\sigma_{a}$ at which the material will first yield by $\sigma_{a}^{\mathrm{min}}$ and the value of $\sigma_{a}$ at which the material will yield completely by $\sigma_{a}^{\mathrm{max}}$. Note that $\sigma_{a}^{\mathrm{max}}<\sigma_{a}^{\mathrm{min}}<0$.

We next consider $\sigma_{a}^{\mathrm{min}}$ and $\sigma_{a}^{\mathrm{max}}$ for all models (§IV.1 and §IV.2, respectively). We then derive an implicit solution to the LL model (§IV.3). Finally, we discuss the numerical scheme used here for solving the QL, QQ, and NQ models (§IV.4).

IV.1 Value of $\sigma_{a}$ at which yield first occurs ($\sigma_{a}^{\mathrm{min}}$)

We determine $\sigma_{a}^{\mathrm{min}}$ by considering the purely poroelastic problem. Solving Equation (26), we obtain the general expression for the linear elastic displacement,

[TABLE]

For the poroelastic problem, we then apply the boundary conditions $\sigma_{r}^{\prime}(a)=(1-\zeta)\sigma_{a}$ and $\sigma_{r}^{\prime}(b)=\sigma_{b}^{\prime}$ to Equation (39b), to obtain

[TABLE]

As discussed above, we assume that yield first occurs at $r=s=a$ and that the yield function (Equation 12) is maximised at steady state. These assumptions require that

[TABLE]

where $a_{\mathrm{min}}$ is the value of $a$ associated with $\sigma_{a}^{\mathrm{min}}$. Combining Equations (39) and (41), we obtain

[TABLE]

where $b_{\mathrm{min}}$ is the value of $b$ associated with $\sigma_{a}^{\mathrm{min}}$ and $a_{\mathrm{min}}$. Evaluating Equation (40) at $a=a_{\mathrm{min}}$ and $b=b_{\mathrm{min}}$ and using the result in Equation (42), we arrive at

[TABLE]

For the LL model, we then take $a_{\mathrm{min}}=a^{\mathrm{ref}}$ and $b_{\mathrm{min}}=1$, at which point $\sigma_{a}^{\mathrm{min}}$ is fully determined. For the QL model, we take $b_{\mathrm{min}}=1$ and enforce the kinematic condition at the inner boundary to arrive at an implicit expression for $a_{\mathrm{min}}$,

[TABLE]

For the QQ and NQ models, we enforce kinematic conditions at both boundaries,

[TABLE]

and

[TABLE]

where $\mathcal{B}_{2}(a_{\mathrm{min}},b_{\mathrm{min}})$ and $\sigma_{a}^{\mathrm{min}}(a_{\mathrm{min}},b_{\mathrm{min}})$ are defined in Equations (40) and (43), respectively.

For $\zeta\equiv 1$, all of the above can be solved for $a_{\mathrm{min}}$ explicitly, whereas $b_{\mathrm{min}}$ must then be determined via numerical root-finding. This is not the case for $\zeta\not\equiv{}1$, in which case the QQ and NQ models require simultaneous root-finding for both $a_{\mathrm{min}}$ and $b_{\mathrm{min}}$.

IV.2 Value of $\sigma_{a}$ at which the material yields completely ($\sigma_{a}^{\mathrm{max}}$)

When the elastic-plastic interface reaches the outer radius, $r=s=b$, we enforce the constraint $\sigma_{r}^{\prime}(b)=\sigma_{b}^{\prime}$ on Equation (29) to arrive at

[TABLE]

where $a_{\mathrm{max}}$ and $b_{\mathrm{max}}$ are the values of $a$ and $b$, respectively, associated with $\sigma_{a}^{\mathrm{max}}$. The values of $a_{\mathrm{max}}$ and $b_{\mathrm{max}}$ are then determined by conditions at the inner and outer boundaries. For the LL model, we simply take $a_{\mathrm{max}}=a^{\mathrm{ref}}$ and $b_{\mathrm{max}}=1$ to arrive at an explicit expression for $\sigma_{a}^{\mathrm{max}}$. For the QL and QQ models222This expression also applies to the LL model, but is not needed to determine $\sigma_{a}^{\mathrm{max}}$., the solution of Equation (37) gives the displacement in the plastic region,

[TABLE]

where $u_{s}(s^{+})$ is the displacement from the elastic solution at $r=s$. For the QL model, we take $s=b_{\mathrm{max}}\equiv 1$. We can then calculate $u_{s}(s=1)$ from the elastic solution (Equation 39), subject to $\sigma_{r}^{\prime}(1)=\sigma_{b}^{\prime}$ and $\alpha\sigma_{\theta}^{\prime}(1)-\sigma_{b}^{\prime}=y$. This gives

[TABLE]

and hence

[TABLE]

so that $\sigma_{a}^{\mathrm{min}}(\zeta)$ is fully determined once $a_{\mathrm{max}}$ is found via the implicit relation

[TABLE]

For the QQ model, rigorous treatment of the kinematics gives $u_{s}(b)=b-1$ and hence

[TABLE]

To find $b_{\mathrm{max}}$, we once again appeal to the elastic problem at $r=s=b$. Solving Equation (39) subject to $\sigma_{r}^{\prime}(b_{\mathrm{max}})=\sigma_{b}^{\prime}$ and $\alpha\sigma_{\theta}^{\prime}(b_{\mathrm{max}})-\sigma_{b}^{\prime}=y$ and equating the result with $u_{s}(b_{\mathrm{max}})=b_{\mathrm{max}}-1$, we arrive at

[TABLE]

with $b_{\mathrm{max}}$ a known constant as defined in Equation (53a). For the NQ model, Equation (38) cannot be solved analytically, prohibiting the derivation of an algebraic expression for $a_{\mathrm{max}}$.

IV.3 Implicit analytical solution to the LL model

For the LL model, we derive analytical expressions for the displacement in both the elastic region ($s<r<b\equiv 1$) and the plastic region ($a\equiv a^{\mathrm{ref}}<r<s$), coupled with an implicit expression for the value of $s$. Solving the problem defined by Equations (37) and (38b) subject to boundary conditions (38c)–(38h) leads to

[TABLE]

The value of $s$ is then determined via the implicit expression

[TABLE]

where $C_{1}(s)$ and $C_{2}(s)$ are as defined in Equation (30) and $\mathcal{B}_{2}$ is as defined in Equation (40). Note that $s$ can take any value from $a^{\mathrm{ref}}$ to 1. We solve Equation (55) for $s$ via numerical root-finding using MATLAB’s fzero. The solution is then fully prescribed by Equation (54).

IV.4 Numerical method for the QL, QQ and NQ models

For all models, we have a closed, coupled, free-boundary BVP in terms of $u_{s}$, as presented in Equations (38) for the NQ model and Equations (37) and (38b)–(38h) for the other three models. For all models, the governing ODE in the elastic region (Equation 38b) can be solved analytically. For the LL, QL, and QQ models, the governing ODE in the plastic region (Equation 37) can also be solved analytically. However, doing so leads to a strongly nonlinear root-finding problem for $\{a,s,b\}$, in which we do not have a good initial guess for $s$. We avoid this by instead solving all of these models and the NQ model numerically by adapting the Chebyshev spectral collocation method of Auton and MacMinn (2017, 2018).

To do so, we map the plastic region ($r\in[a,s]$) and the elastic region ($r\in[s,b]$) to separate Chebyshev grids, each with $N$ Chebyshev nodes. Note that both domains include the elastic-plastic interface ($r=s$), meaning that the full domain $r\in[a,b]$ is discretised by $2N-1$ nodes (“collocation points”). We then discretise all derivatives using dense Chebyshev differentiation matrices, thus converting this coupled free-boundary BVP into a system of algebraic equations. We solve for $2N$ unknowns, comprising the displacement at the $2N-1$ nodes as well as the value of $s$; we therefore require a system of $2N$ constraints.

In the plastic region, the governing ODE (Equation 38a) is first order and thus must be enforced at exactly $N-1$ nodes. We enforce this for the first $N-1$ nodes of the discretised domain ($r\in[a,s)$), having already used the boundary condition $\sigma_{r}^{\prime}(a)=(1-\zeta)\sigma_{a}$ in the derivation this ODE. In the elastic region, the governing ODE (Equation 38b) is second order and must be enforced at exactly $N-2$ nodes. We enforce this for the $N-2$ interior nodes ($r\in(s,b)$) and then impose $\sigma_{r}^{\prime}(b)=\sigma_{b}$ at the outer boundary (Equation 38f). Lastly, we enforce two conditions at $r=s$: Continuity of radial effective stress (Equation 38d) and the requirement that the elastic stresses must satisfy the yield condition (Equation 38e). We solve this nonlinear system via Newton iteration. At each iteration, we update the domain via the kinematic conditions at the inner and outer boundaries (Equations 38g and 38h) and the current value of $s$.

The main structural difference between the models presented here and the models presented in Auton and MacMinn (2017, 2018) is the coupling of an elastic domain with a plastic domain at a free (but non-material) boundary. As a result, the derivation of an exact analytical Jacobian matrix is nontrivial and we instead approximate the Jacobian matrix numerically. In addition, whereas $a^{\mathrm{ref}}$ and $b^{\mathrm{ref}}\equiv{}1$ can be used as reasonable initial guesses for $a$ and $b$, no such guess exists for $s$. To accommodate the approximate nature of the Jacobian and the lack of a good initial guess for $s$, we begin by calculating the solution for the LL model for a given value of $\zeta$ and a given, small value of $\sigma_{a}<\sigma_{a}^{\mathrm{min}}$. We use this solution as an initial guess for the numerical solution of the QL model, which then provides an initial guess for the QQ model, which then provides an initial guess for the NQ model. We extend the solution for each model to larger values of $\sigma_{a}>\sigma_{a}^{\mathrm{max}}$ using numerical continuation. Note that $\sigma_{a}^{\mathrm{max}}<\sigma_{a}^{\mathrm{min}}<0$. This approach allows for the fact that the model behaviours diverge from one another as the driving strength increases (see Figure 3). The same approach can also be used for increasing $\zeta$ at fixed $\sigma_{a}$, since deformation increases with $\zeta$ (see Figures 4 and 5).

V Results

We now use our models and solutions to study the poroelasto-plastic deformation of a cohesive granular cylinder, subject to pressurisation of the inner cavity. As discussed in §III.2.3 and §IV, we have developed solutions for four distinct model combinations: LL, QL, QQ and NQ. We fix all parameters based on the values discussed in §V.1 below, except for $\sigma_{a}$ (or $q$) and $\zeta$, which we vary.

For comparison with our poroelasto-plastic models, we additionally extend the purely poroelastic models from Auton and MacMinn (2017, 2018) by modifying the boundary conditions to incorporate the permeability-load parameter $\zeta$. Specifically, we present solutions to the L-$k_{0}$ model (here denoted by ‘L’; linear elasticity, linearised kinematics, and constant permeability) and the Q-$k_{0}$ model (here denoted by ‘Q’; linear elasticity, rigorous kinematics, and constant permeability) alongside our poroelasto-plastic solutions in Figures 2–5 below. This comparison illustrates the impact of plasticity across different values of $\zeta$.

To investigate the impact of fluid injection into the pre-stressed subsurface, we wish to isolate the additional stress and deformation due to fluid injection from those due to background compression. We approximate the disturbance in all quantities due to fluid injection by subtracting their values in the compressed initial state (III in Figure 1) from their values in the final steady state (V in Figure 1). For the displacement, for example, we consider the disturbance $\delta{u_{s}}\vcentcolon=u_{s}-u_{s,0}$. Recall that these compressed initial states are presented in §III.1. We subtract the L initial state from the LL and QL models, and the Q initial state from the QQ and NQ models.

V.1 Parameters

We adopt a set of parameter values motivated by boreholes in the subsurface. As such, we adopt $\breve{\mathcal{M}}\sim 50$ GPa and $\breve{\Lambda}\sim 27$ GPa (Young modulus $\breve{E}\sim 25$ GPa and Poisson ratio $\nu\sim 0.36$) and a porosity of $\phi_{f}^{\mathrm{ref}}\sim 0.20$ as typical properties of sedimentary rocks such as sandstones and shales (e.g., Goodman, 1980; Hart and Wang, 1995; Bobko and Ulm, 2008; Bobko et al., 2011; Rickman et al., 2008; Britt and Schoeffler, 2009). We further assume moderate internal friction and cohesion, $\varphi\sim 35\text{\textdegree}$ and $\breve{c}\sim 120$ MPa, but very little dilation, $\psi\sim 0.3\text{\textdegree}$ Mandl (2005); Bobko and Ulm (2008); Bobko et al. (2011); Vermeer and Borst (1984). We take $-\breve{\sigma}_{b}^{\prime}\sim 50$ MPa, as appropriate for a depth of \sim$$2.5\,\mathrm{km} (Neuzil, 1994; Islam et al., 2010). These dimensional values correspond to $\Gamma=0.55$, $\alpha=4$, $\beta=1.01$, $y=10^{-2}$, $a^{\mathrm{ref}}=10^{-4}$, $\phi_{f}^{\mathrm{ref}}=0.2$ and $\sigma_{b}^{\prime}=-1\times 10^{-3}$. Note that $\alpha\nu\equiv\frac{\alpha\Gamma}{1+\Gamma}\sim 1.4>1$ and $y(1-a^{2})+2\sigma_{b}^{\prime}\sim 0.008>0$, where we take $a\approx{}a^{\mathrm{ref}}$; as such, the relevant yield condition is $\mathcal{F}_{\theta,r}=0$ (see Appendix A). Below, we consider the response of this material to different driving strengths (varying $\sigma_{a}$ or $q$), and we additionally consider the transition from impermeable to fully permeable by varying $\zeta$ from 0 to 1.

V.2 Fully permeable: fixed $q$

We begin comparing poroelasto-plastic behaviour with purely poroelastic behaviour in the context of a fully permeable material ($\zeta\equiv{}1$; Figures 2 and 3). Recall that, for this model, we drive the system with a fixed flow rate $q$ instead of a fixed pressure difference $\Delta{p}\equiv-\sigma_{a}$ (see §II.5.3). For each model, a given value of $q$ will correspond to a specific steady-state value of $\Delta{p}$; however, this value will be different for each model. Note that, as with $\sigma_{a}$, there exists a minimum value $q_{\mathrm{min}}$ at which yield first occurs and a maximum value $q_{\mathrm{max}}$ at which the material yields completely. These values can be calculated via $q_{\mathrm{min}}=-\sigma_{a}^{\mathrm{min}}(1)/\ln\left(\frac{b}{a}\right)$ and $q_{\mathrm{max}}=-\sigma_{a}^{\mathrm{max}}(1)/\ln\left(\frac{b}{a}\right)$, where $\sigma_{a}^{\mathrm{min}}(1)$ and $\sigma_{a}^{\mathrm{max}}(1)$ are defined in Equations (43) and (47), respectively.

In Figure 2, we consider $q=0.0012$. For these parameters, the LL model has $q_{\mathrm{min}}\approx 3.1\times 10^{-4}$ and $q_{\mathrm{max}}\approx 0.0033$, so that our chosen value of $q$ is about $4q_{\mathrm{min}}$ and about $q_{\mathrm{max}}/3$. For this value of $q$, we plot the disturbance due to fluid injection for various key quantities for the poroelastic models (left column) and the poroelasto-plastic models (right column) against the Lagrangian radial coordinate $R(r)\equiv{}r-u_{s}(r)$. For reference, we plot the true steady-state values of the same quantities (i.e., without subtracting the initial state) in Figure 6 in Appendix B. Note that $\delta{p}\equiv{}p$ because $p_{0}\equiv{}0$ for all models.

Although this value of $q$ is substantially higher than $q_{\mathrm{min}}$, it leads to relatively small elastic deformations in the absence of yield. As a result, the predictions of the two purely poroelastic models are essentially indistinguishable (Figure 2, left column). The disturbances $\delta{\phi_{f}}$, $\delta{\sigma_{\theta}^{\prime}}$, and $\delta{p}$ all have maxima at the cavity wall, with $\delta{\sigma_{\theta}^{\prime}}$ falling off relatively steeply from this value as $R$ increases. In contrast, $\delta{u_{s}}$ and $\delta{\sigma_{r}^{\prime}}$ have internal maxima near the outer boundary ($R\sim 0.7$) and near the inner boundary ($R\sim{}4\times{}10^{-4}$), respectively. The value of $\delta{\sigma_{r}^{\prime}}$ is pinned to $-\sigma_{b}^{\prime}$ at the inner boundary by construction since $\sigma_{r}^{\prime}(a)=0$ and $\sigma_{r,0}^{\prime}\equiv\sigma_{b}^{\prime}$, and vanishes at the outer boundary since $\sigma_{r}^{\prime}(b)=\sigma_{r,0}^{\prime}=\sigma_{b}^{\prime}$. Note that $\delta{\phi_{f}}$ is strictly positive, meaning that the porosity everywhere has increased relative to the compressed initial state.

Upon introducing plastic yield (Figure 2, right column), all quantities except $\delta{p}$ behave drastically differently than in the purely poroelastic case. Most importantly, plasticity enables much larger displacements in the plastic region (near the cavity wall), with much smaller associated stresses. The strains in the plastic region are no longer infinitesimal—for example, $\mathrm{max}\left(\frac{u_{s}}{r}\right)\approx 130$ for the poroelasto-plastic LL model whereas $\mathrm{max}\left(\frac{u_{s}}{r}\right)\approx 0.02$ for the poroelastic L and Q models, with this maximum occurring at the inner boundary in both cases. These large strains make model choice much more important in the plastic region, as evidenced by the substantial differences between the models there. In the elastic region, however, the maximum strain remains relatively small—for example, $\mathrm{max}\left(\frac{u_{s}^{\mathrm{e}}}{r}\right)\approx 0.003$ for the poroelasto-plastic LL model. This suggests that it is reasonable to linearise the elastic part of the solution, even in the presence of large displacements in the plastic region. This idea is further supported by the other results in Figures 2–5, in which the QL and QQ models are essentially indistinguishable. We conclude that model choice in the elastic region is relatively unimportant to the overall behaviour. Note also that although model choice in the plastic region is important for the behaviour in the plastic region, the resulting (substantial) differences in behaviour have relatively little impact on behaviour in the elastic region.

In the plastic region, the LL model generally predicts the most extreme behaviour, with additional facets of nonlinearity increasingly moderating this behaviour. The QL and QQ models are effectively indistinguishable from each other, again demonstrating that model choice in the elastic region is unimportant. The NQ model predicts very similar behaviour to the QL and QQ models except with regard to porosity, where the NQ predicts much smaller values of $\delta{\phi_{f}}$ than any of the models. This discrepancy is due to the fact that $\mathrm{d}u_{s}/\mathrm{d}r$ is much more negative near the inner radius for the NQ model than for the QQ or QL models (Eq. 3). The LL model predicts $\delta{\phi_{f}}>1$ near the inner boundary, highlighting the fact that linearising the kinematics in the plastic region can lead to nonphysical solutions333In fact, we need $\delta{\phi_{f}}\lesssim 0.8$ since $\delta{\phi_{f}}\equiv\phi_{f}-\phi_{f,0}\approx{}0.8$ for these parameters (cf. Figure 6). (see Figure 7 in Appendix B). Note that, despite the linearisation of the kinematics in the elastic region, the QL model does not have this feature.

The internal maximum in $\delta{\sigma_{r}^{\prime}}$ has shifted further into the material and decreased in magnitude relative to the poroelastic case, now occurring at the elastic-plastic interface and with a magnitude of about $1/2$ of the poroelastic value. The maximum value of $\delta{\sigma_{\theta}^{\prime}}$ is now internal rather than at the inner radius, also occurring at the elastic-plastic interface and with a magnitude of about $1/4$ of the poroelastic value. This quantity has a sharp transition across the elastic-plastic interface, meaning that $\frac{\mathrm{d}}{\mathrm{d}r}(\delta{\sigma_{\theta}^{\prime}})$ is discontinuous across $r=s$. It is straightforward to show from the boundary conditions at $r=s$, the yield condition, and mechanical equilibrium that $\delta{u_{s}}$, $\delta{\sigma_{r}^{\prime}}$, $\frac{\mathrm{d}}{\mathrm{d}r}(\delta{\sigma_{r}^{\prime}})$, $\delta{\sigma_{\theta}^{\prime}}$, and $\delta{p}$, and $\frac{\mathrm{d}}{\mathrm{d}r}(\delta{p})$ must all be continuous across $r=s$. However, $\frac{\mathrm{d}}{\mathrm{d}r}(\delta{u_{s}})$ and therefore also $\delta{\phi_{f}}$ are weakly discontinuous across $r=s$ (not readily visible in the figure). The discontinuity in $\frac{\mathrm{d}}{\mathrm{d}r}(\delta{\sigma_{\theta}^{\prime}})$ is particularly prominent because it occurs at the maximum value of $\sigma_{\theta}^{\prime}$, whereas those in $\delta{\phi_{f}}$ and $\frac{\mathrm{d}}{\mathrm{d}r}(\delta{u_{s}})$ occur in regions where $\delta{\phi_{f}}$ and $\delta{u_{s}}$ are themselves small.

For the poroelastic L and poroelasto-plastic LL models, $\delta{p}$ decreases exactly logarithmically with $r$ because $\frac{\mathrm{d}p}{\mathrm{d}r}=-\frac{q}{r}$, and these linearised models do not account for the moving boundaries. The poroelastic Q model exhibits essentially the same behaviour since the poroelastic displacements are small. The poroelasto-plastic QL, QQ, and NQ models do account for the substantial increase in $a$ and the comparatively small increase in $b$ due to fluid injection, leading to lower injection pressures at steady state.

In Figure 3, we present results for a wide range of $q$ in terms of four summary quantities. Note that, for values of $q$ for which yield does not occur (i.e., $q<q_{\mathrm{min}}$), the LL and QL models behave according to the L model, while the QQ and NQ models behave according to the Q model. The maximum value of $q$ shown in Figure 3 is $q\approx 0.0012$, which is the value used in Figure 2.

For $q\lesssim{}10^{-4}$ (Figure 3, grey band), the disturbance in the inner radius $\delta{a}$ is negative (Figure 3, top left). This contraction occurs because $q$ is not large enough to support the compressive confining stress, such that the hypothetical casing would need to partially support the material to enforce $a\geq{}a_{0}$ (II to III in Fig. 1). Since we do not model the casing, our results in this range produce a steady state that is inconsistent with state III in Fig. 1. For larger values of $q$, the inner radius increases monotonically with $q$ as injection increasingly pushes the material radially outwards. After yield occurs, the poroelasto-plastic models deform increasingly more than the poroelastic models, and the ordering of the models is the same as in Figure 2.

For small flow rates, the maximum porosity disturbance $\mathrm{max}(\delta{\phi_{f}})$ is small (Figure 3, top right). Much like $\delta{a}$, $\mathrm{max}(\delta{\phi_{f}})$ increases linearly with $q$ until yield and then departs strongly from linear behaviour for the poroelasto-plastic models. For a small range of $q$ after yield, the poroelasto-plastic models predict values of $\mathrm{max}(\delta{\phi_{f}})$ that are less than the corresponding poroelastic predictions, which occurs because yield reduces the maximum tensile stress and relatively little plastic flow has occurred. For larger values of $q$, the predictions of the poroelasto-plastic models are much larger than those of the poroelastic models as deformations grow larger and plastic flow leads to significant dilation. The NQ model predicts a much weaker departure from poroelasticity than the other three poroelasto-plastic models.

The injection pressure $\Delta p\equiv{}p(a)\equiv\delta{p(a)}$ increases linearly with $q$ before yield, and continues to increase linearly with $q$ for the LL model after yield (Figure 3, bottom left). For the other poroelasto-plastic models, $\Delta{p}$ increases somewhat slower than linearly as the inner radius increases. However, this is a weak effect because the changes in inner and outer radii are ultimately small for all $q$. The injection pressure is important because it is one of the few observables during subsurface operations (i.e., it can be measured from the surface in realtime).

For small $q$, the maximum effective stress disturbance $\max(\delta{\sigma^{\prime}})$ is constant and equal to $-\sigma_{b}^{\prime}=10^{-3}$ (Figure 3, bottom right). This is an artefact of the fact that $\sigma_{a}^{\prime}\equiv 0$, such that $\delta{\sigma_{a}^{\prime}}\equiv\sigma_{a}^{\prime}-\sigma_{b}^{\prime}=10^{-3}$. This is the maximum stress (and the maximum stress disturbance) for small $q$ because all of the other stresses remain compressive until $q$ becomes large enough to generate tensile stresses. Note again that these values are only physically meaningful for $q\gtrsim{}10^{-4}$ (outside the grey region). All models exhibit a corner in $\max(\delta{\sigma^{\prime}})$ before yield due to a change in which stress component exhibits this maximum: $\max(\delta{\sigma^{\prime}})$ is equal to $\delta{\sigma_{r}^{\prime}}(a)$ for $q\lesssim{}1.4\times{}10^{-4}$ and to $\delta{\sigma_{\theta}^{\prime}}(a)$ for $q\gtrsim{}1.4\times{}10^{-4}$. The poroelasto-plastic models exhibit a second corner at yield ($q\approx{}3.1\times{}10^{-4}$) due to a change in the location of this maximum: $\max(\delta{\sigma^{\prime}})$ is equal to $\delta{\sigma_{\theta}^{\prime}}(a)$ before yield and to $\delta{\sigma_{\theta}^{\prime}}(s)$ after yield.

In Figure 8 in Appendix B, we plot the same four quantities as in Figure 3, but against $-\sigma_{b}^{\prime}$ for fixed $q$. This shows the transition from an unconstrained cylinder (no confining stress at the outer boundary) to a strongly compressed cylinder, providing a link with results of Auton and MacMinn (2017, 2018) for an unconstrained, fully permeable poroelastic cylinder. Within each class of models (poroelastic and poroelasto-plastic), the ordering of the models is preserved for each quantity as $-\sigma_{b}^{\prime}$ varies. Additionally, all quantities converge smoothly towards poroelastic behaviour as $-\sigma_{b}^{\prime}$ increases. As a result, it seems reasonable to extrapolate the conclusions of Auton and MacMinn (2017, 2018) for a poroelastic cylinder with deformation-dependent permeability to the poroelasto-plastic cylinders considered here. Figure 2a of Auton and MacMinn (2017) suggests that deformation-dependent permeability leads to a much smaller $\Delta{p}$ for a given value of $q$ at steady state, and thus to less deformation and lower stresses. We expect that the same would be true for the poroelasto-plastic scenarios considered here.

V.3 Impermeable to fully permeable: fixed $\sigma_{a}$

In Figures 4 and 5, we vary the permeability-load parameter $\zeta$, transitioning the model from impermeable ($\zeta\equiv 0$) to fully permeable ($\zeta\equiv 1$) for a fixed $\sigma_{a}=-0.0075$. This value of $\sigma_{a}$ is between $\sigma_{a}^{\mathrm{min}}$ and $\sigma_{a}^{\mathrm{max}}$ for all values of $\zeta$, such that $a<s<b$. Recall that $0\leq\zeta<1$ is analogous to a thin, weak, low-permeability membrane on the inner cavity wall, the permeability of which decreases as $\zeta$ decreases; varying $\zeta$ then varies the partitioning of the fixed total radial stress at the inner cavity wall between fluid loading (injection pressure) and mechanical loading (effective radial stress).

In Figure 4, we plot the same quantities as in Figure 2 for four values of $\zeta$ ranging from nearly impermeable ($\zeta\approx 0$) to nearly fully permeable ($\zeta\approx 1$). For reference, we again plot the true steady-state values of the same quantities (i.e., without subtracting the initial state) in Figure 9 in Appendix B.

Note that, because we now drive these models with fixed $\sigma_{a}$, each model will produce a different value of $q$. We found earlier that increasing nonlinearity led to lower values of $\Delta{p}$ for a given value of $q$ (Figure 2). Here, we then expect that increasing nonlinearity will lead to a higher value of $q$ for a given value of $\sigma_{a}$, and therefore to larger deformations and stresses. This suggests that the relative ordering of the models should be reversed in Figure 4 relative to Figure 2, which is indeed the case.

For the poroelastic models (left column), all quantities increase with $\zeta$. This suggests that, for a given imposed cavity pressure ($-\sigma_{a}$), a more permeable material will experience larger deformations and larger stresses than a less permeable material. The three quantities $\delta{\phi_{f}}$, $\delta{\sigma_{\theta}^{\prime}}$, and $\delta{p}$ depend quantitatively but not qualitatively on $\zeta$, varying monotonically from a maximum value at the inner radius to a minimum value at the outer radius; the former increases monotonically with $\zeta$, whereas the latter is insensitive to $\zeta$. The two quantities $\delta{u_{s}}$ and $\delta{\sigma_{r}^{\prime}}$ vary both quantitatively and qualitatively with $\zeta$. The disturbance in displacement $\delta{u_{s}}$ is relatively insensitive to $\zeta$ at the cavity wall, implying that permeability makes very little difference to the size of the cavity. However, $\delta{u_{s}}$ is increasingly sensitive to $\zeta$ for larger values of $R$, such that the maximum displacement shifts from $R=a_{0}$ for small $\zeta$ to near $R=1$ for moderate to large $\zeta$. Although the $\delta{u_{s}}$ curves appear to get closer together as $\zeta$ increases, this is an artefact of the logarithmic vertical scale. The disturbance in radial effective stress $\delta{\sigma_{r}^{\prime}}$ vanishes at the outer boundary by construction, and is fixed to $\delta{\sigma_{r}^{\prime}}(a)=(1-\zeta)\sigma_{a}-\sigma_{b}^{\prime}$ at the inner boundary. The latter value increases with $\zeta$, transitioning from compressive to tensile, and the rest of the curve essentially pivots about the constant outer boundary value. As a result, $\delta{\sigma_{r}^{\prime}}$ is monotonic in $R$ for small $\zeta$, with a minimum at the inner boundary and a maximum at the outer boundary, but transitions to being nonmonotonic in $R$ for large $\zeta$, having an interior maximum and eventually a minimum at the outer boundary.

For the poroelasto-plastic case, the difference between the models grows larger as $\zeta$ increases (right column). This is because plasticity leads to larger deformations, as do larger values of $\zeta$. The plastic region itself also grows larger as $\zeta$ increases (see $\delta{\sigma_{\theta}^{\prime}}$). For this value of $\sigma_{a}$ and sufficiently small $\zeta$, the deformations are small enough that the models are essentially indistinguishable. For this value of $\sigma_{a}$ and larger values of $\zeta$, the deformations grow sufficiently large that significant differences emerge between the models. For $\delta{\phi_{f}}$, $\delta{\sigma_{\theta}^{\prime}}$, and $\delta{p}$, the consequences of increasing $\zeta$ are again mostly quantitative: The deformations and stresses grow larger. For all models except the LL model, $\delta\sigma_{r}^{\prime}$ and $\delta{\sigma_{\theta}^{\prime}}$, exhibit weak non-monotonicity in $\zeta$ for some intermediate values of $R$ as $\zeta$ approaches 1. Note that the QL and QQ models predict very large values of $\delta{\phi_{f}}$ for $\zeta\approx{}1$ relative to the other models and other values of $\zeta$. For $\delta{u_{s}}$ and $\delta\sigma_{r}^{\prime}$, the changes are again more complex. Plasticity leads to much larger displacements at the inner boundary as $\zeta$ increases, but to very little change in the displacement at the outer boundary. The maximum in $\delta{u_{s}}$ shifts from the inner boundary for all models for small $\zeta$ to the outer boundary for all models for intermediate $\zeta$, and back to the inner boundary for the more nonlinear models for large $\zeta$.

In Figure 5, we plot the same quantities as in Figure 3, but against $\zeta$ for fixed $\sigma_{a}=-0.0075$. As also illustrated in Figure 4, it is clear that model choice is unimportant for the poroelastic models for all $\zeta$, and for the poroelasto-plastic models for sufficiently small $\zeta$ (for this set of parameters, $\zeta\lesssim 0.55$). Model choice eventually becomes important for the poroelasto-plastic models because plasticity enables large deformations and these deformations increase monotonically with $\zeta$ (Figure 5, top left). As also illustrated in Figure 4, the more nonlinear models again predict larger deformations. The two poroelastic models predict that $\max\left(\delta{\phi_{f}}\right)$ increases approximately linearly with $\zeta$. These predictions are larger in magnitude than those of the poroelasto-plastic models until $\zeta\approx{}0.55$ (QL and QQ) or $\approx{}0.8$ (LL and NQ), at which point the poroelasto-plastic models increase strongly following a corner where $\max\left(\delta{\phi_{f}}\right)$ shifts from an internal value to the value at the inner radius. The flow rate $q$ increases monotonically with $\zeta$ for all models, with the more nonlinear models predicting larger flow rates as $\zeta$ approaches 1. Note that $q=-\frac{\zeta\sigma_{a}}{\ln\left(\frac{b}{a}\right)}$ for all models, and is thus pinned to $q=0$ for an impermeable material ($\zeta\equiv 0$) and to $q=\frac{\Delta p}{\ln\left(\frac{b}{a}\right)}$ for a fully permeable material ($\zeta\equiv 1$). This dependence on $a$ and $b$ links $q$ to the displacement field, and hence to model choice, for $\zeta>0$. Note that, for the L and LL models, $q=\frac{\zeta\sigma_{a}}{\ln(a^{\mathrm{ref}})}$ and is therefore independent of deformation. The maximum disturbance in effective stress $\max\left(\delta{\sigma^{\prime}}\right)$ increases monotonically with $\zeta$ for the L, Q, and LL models, and becomes weakly nonmonotonic in $\zeta$ for the QL, QQ, and NQ models as $\zeta$ approaches 1. This nonmonotonicity is more pronounced for the NQ model. The fact that $\max\left(\delta{\sigma^{\prime}}\right)$ is maximised for a particular value of $\zeta$ could be important in applications such as hydraulic fracturing in ductile shales or for borehole integrity. Note that the ordering of models for $\max\left(\delta{\sigma^{\prime}}\right)$ is reversed relative to the ordering of models for $\delta{a}$ and $q$, reinforcing the fact that plasticity dissipates elastic energy and leads to lower (less tensile) stresses. Note, finally, that the poroelasto-plastic models predict significantly larger displacements and significantly lower stresses than the poroelastic models for all values of $\zeta$ (Figure 5, top left and bottom right), whereas all models predict similar flow rates for most values of $\zeta$.

VI Conclusion

Fluid-driven deformation is relevant to applications in borehole integrity and cavity expansion; motivated by these problems, this study provides the first kinematically rigorous poroelasto–perfectly-plastic model for fluid injection into a thick-walled annulus. To assess the importance of plasticity, and of large deformations, we performed a detailed examination of four such models: Classical linear poroelasto-plasticity (i.e., linear elasticity with linearised kinematics; LL); linear elasticity with rigorous kinematics in just the plastic region (QL) and in both regions (QQ); and linear elasticity, rigorous kinematics, and logarithmic (Hencky) strains in the plastic region (NQ). For a set of parameter values motivated by sedimentary rocks such as sandstone or shale, we then compared the predictions of these models with each other, and with those of two poroelastic models: Classic linear poroelasticity (i.e., linear elasticity with linearised kinematics; L) and linear elasticity with rigorous kinematics (Q).

We showed that there was negligible difference between the poroelastic L and Q models for these parameters because the deformations remain small. In contrast, plasticity enables large deformations in the plastic region, making model choice much more important there. Accounting for rigorous kinematics in the plastic region, in particular, can have a significant impact on the predicted behaviour (e.g., compare the LL and QL models). However, this effect is isolated to the plastic region, where deformations are large; in the elastic region, in contrast, deformations remain small and linearised kinematics remain appropriate, even in the presence of large deformations in the plastic region (e.g., compare the QL and QQ models). Linearised kinematics in the plastic region can lead to non-physical predictions (Figure 7 in Appendix B).

Previous models have treated low-permeability materials such as shale as either fully permeable or fully impermeable. Here, we proposed a new ‘permeability-load parameter’ $\zeta$ that enables a smooth transition between these limiting states by (essentially) introducing a thin, weak, low-permeability skin at the cavity wall. In Figures 2 and 3, we considered a fully permeable material ($\zeta\equiv 1$) for a range of injection rates $q$. In Figures 4 and 5, we considered a fixed total stress at the inner radius $\sigma_{a}$ as $\zeta$ transitions from 0 to 1. We showed that the amount of deformation increases with $\zeta$ for a given value of $\sigma_{a}$. The maximum tensile effective stress exhibits a maximum at an intermediate value of $\zeta$ near 1, such that an annulus with a slight reduction in permeability at the cavity wall experiences the greatest effective stress. Since the deformation increases with $\zeta$, the choice of poroelasto-plastic model becomes increasingly significant as $\zeta$ increases.

Our results highlight the significant qualitative and quantitative differences between fully impermeable and fully permeable materials, and provide a mechanism for smoothly transitioning between these two end-member behaviours. As such, many practicals scenarios that are currently modelled as either fully impermeable or fully permeable, such as boreholes with clogged or damaged walls, boreholes that have been treated with wall-building chemicals, or boreholes in low-permeability rocks such as shales, are probably best modelled with intermediate values of $\zeta$. We have also shown that, although plastic failure leads to drastically different material behaviour, including much larger deformations and much smaller stresses, it has a relatively minor impact on injection pressure. This indicates that injection pressure is a relatively weak indicator of plastic failure, which may be problematic in practice because pressure is one of the primary observables during injection. In other words, the onset of ductility may be quite difficult to detect from the surface, despite its strong impact on displacements and stresses.

In addition to the above qualitative points about the influence of constitutive behaviour, kinematics, and key parameters, our work here also provides a reference solution that can be used as a rigorous benchmark for finite-element algorithms. It may also be useful for interpreting laboratory experiments in similar geometries (e.g., MacMinn et al., 2015).

Acknowledgements.

The authors are grateful to EPSRC for support in the form of a Doctoral Training Award to L.C.A.

Appendix A Which yield condition?

We now consider all six possible cohesive Mohr-Coulomb yield criteria for a cylinder in plane strain with principal stresses $\sigma_{r}^{\prime}$, $\sigma_{\theta}^{\prime}$ and $\sigma_{z}^{\prime}$. For plane strain,

[TABLE]

where $\nu\equiv\frac{\Gamma}{1+\Gamma}$ is the Poisson ratio (cf. §V.1). Recall that $\nu\in(0,\frac{1}{2})$ for most physical materials, and that $\sigma^{\prime}_{1}\geq\sigma^{\prime}_{2}\geq\sigma^{\prime}_{3}$ by definition. We then have three possible constraints based on the signs of $\sigma_{r}^{\prime}$ and $\sigma_{\theta}^{\prime}$: If $\sigma_{r}^{\prime},\sigma_{\theta}^{\prime}>0$, then $\sigma^{\prime}_{1}\neq{}\sigma_{z}^{\prime}$; if $\sigma^{\prime}_{r},\sigma^{\prime}_{\theta}<0$, then $\sigma_{3}^{\prime}\neq\sigma^{\prime}_{z}$; and if $\sigma^{\prime}_{r}$ and $\sigma^{\prime}_{\theta}$ have different signs, or if one of them is zero, then $\sigma^{\prime}_{2}\equiv\sigma^{\prime}_{z}$. The six possible yield functions are then

[TABLE]

where the material remains elastic for $\mathcal{F}_{1,3}<0$ and yields when $\mathcal{F}_{1,3}=0$. We must then determine which of these yield functions first reaches zero. As stated in §II.4.2, we assume that yield first occurs at $r=a_{\mathrm{min}}$; that, prior to yield, the yield function $\mathcal{F}_{1,3}$ is maximised at the poroelastic steady-state and not during the preceding transient evolution; and that, once yield occurs according to a particular yield condition, the material will fail exclusively according to this condition. We are therefore only concerned with the values of $\sigma^{\prime}_{i}(a_{\mathrm{min}})$, where $i\in\{r,\theta,z\}$ (i.e., the effective stresses at the inner boundary at the point of first yield).

A.1 Fully permeable ($\zeta\equiv 1$)

For $\zeta\equiv 1$, $\sigma^{\prime}_{r}(a)=0$ implies that $\sigma_{z}^{\prime}(a)=\frac{\Gamma}{1+\Gamma}\sigma^{\prime}_{\theta}(a)$ and, as such, that $\sigma^{\prime}_{2}\equiv\sigma^{\prime}_{z}$. Hence, the only possible yield conditions are $\mathcal{F}_{\theta,r}=0$ and $\mathcal{F}_{r,\theta}=0$. The former is appropriate if $\sigma_{\theta}^{\prime}(a)\geq\sigma^{\prime}_{z}(a)\geq\sigma^{\prime}_{r}(a)=0$, which we assumed throughout our analysis in the main text. We now consider the conditions under which the latter would be appropriate.

A.1.1 $\mathcal{F}_{r,\theta}=0$

For $\mathcal{F}_{r,\theta}=0$ to be the appropriate yield condition, it must be the case that $\sigma^{\prime}_{\theta}(a)\leq\sigma_{z}^{\prime}(a)\leq\sigma^{\prime}_{r}(a)=0$. For the L and Q models with $\zeta=1$, it is straightforward to show that (see §IV.1)

[TABLE]

where $B_{2}$ is given by

[TABLE]

We then rearrange the condition $\mathcal{F}_{r,\theta}<0$, for which deformation remains elastic, to give

[TABLE]

This means that the material will only yield if the injection rate is small enough, or sufficiently negative (suction), to trigger cavity collapse. In the following, we refer to this as ‘negative yield’ and the converse as ‘positive yield’.

We are only concerned here with fluid injection problems, so we only consider $q\geq 0$ (no suction/extraction). As a result, ensuring that $q_{\star}<0$ would then imply that $\mathcal{F}_{r,\theta}<0$ for $q\geq 0$. Rearranging the requirement that $q_{\star}<0$ leads to the constraint that negative yield cannot occur for all relevant values of $q$ if

[TABLE]

meaning that the cylinder must be sufficiently ‘strong’ relative to the compressive far-field stress $\sigma_{b}^{\prime}$. This constraint is satisfied for the parameter values used here (see §V.1). Hence, we can state conclusively that, for a sufficiently strong and fully permeable cylinder, $\mathcal{F}_{\theta,r}$ is the correct yield condition during fluid injection ($q\geq 0$).

A.2 Partially permeable materials ($\zeta\not\equiv 1$)

For $\zeta\in[0,1)$, the radial effective stress at the inner radius is $\sigma_{r}^{\prime}(a)=(1-\zeta)\sigma_{a}\in[-|\sigma_{a}|,0)$. Hence, we are left with two cases to consider:

$\sigma^{\prime}_{r}(a)<0$ and $\sigma_{\theta}^{\prime}(a)>0$, in which case $\sigma^{\prime}_{2}\equiv\sigma^{\prime}_{z}$, and $\mathcal{F}_{\theta,r}=0$ and $\mathcal{F}_{r,\theta}=0$ are the possible yield criteria. 2. 2.

$\sigma^{\prime}_{r}(a),\sigma_{\theta}^{\prime}(a)<0$, in which case $\mathcal{F}_{\theta,r}=0$, $\mathcal{F}_{r,\theta}=0$, $\mathcal{F}_{z,\theta}=0$ and $\mathcal{F}_{z,r}=0$ are the possible yield criteria.

Note that, for $\sigma_{r}^{\prime}(a)<0$, the axial effective stress $\sigma_{z}^{\prime}(a)$ can never be the minimum principal stress (see discussion above Eqs. 57). Thus, we only have three alternative yield conditions to consider.

A.2.1 $\mathcal{F}_{r,\theta}=0$

Ensuring that $\mathcal{F}_{r,\theta}$ is always strictly negative will prevent yield according to $\mathcal{F}_{r,\theta}$. Using Equation (39) leads to

[TABLE]

where $\mathcal{B}_{2}$ is defined in Equation (40). Following the same procedure that leads to Equation (43), we find that yield according to $\mathcal{F}_{r,\theta}$ will not occur for

[TABLE]

meaning that $\sigma_{a;r,\theta}^{\mathrm{min}}$ is the least compressive value of $\sigma_{a}$ that would prevent negative yield according to $\mathcal{F}_{r,\theta}$. Note that the $\sigma_{a}^{\mathrm{min}}$ presented in Equation (43) is $\sigma_{a;\theta,r}^{\mathrm{min}}$, corresponding to positive yield according to $\mathcal{F}_{\theta,r}$.

We focus on injection, so it must be the case that $\sigma_{a}<0$. As a result, ensuring that $\sigma_{a}^{\mathrm{min},r,\theta}>0$ would then imply that $\mathcal{F}_{r,\theta}<0$ for all relevant $\sigma_{a}$. This constraint simply requires that Equation (61) must be satisfied since the denominator of Equation (63) is negative.

A.2.2 $\mathcal{F}_{z,\theta}=0$

Ensuring that $\mathcal{F}_{z,\theta}$ is always strictly negative will prevent yield according to $\mathcal{F}_{z,\theta}$. This is the appropriate yield function if $\sigma_{\theta}^{\prime}(a)<\sigma_{r}^{\prime}(a)<\sigma_{z}^{\prime}(a)<0$. If $\alpha\nu\geq 1$, as is the case for the parameters used above (cf. §V.1), then $\mathcal{F}_{z,\theta}<0$ rearranges to

[TABLE]

which is always satisfied. If $\alpha\nu<1$, the requirement is more complicated. As stated above, this yield condition is only appropriate if $\sigma_{\theta}^{\prime}(a)<\sigma_{r}^{\prime}(a)$, and therefore cannot be satisfied if $\sigma_{r}^{\prime}(a)<\sigma_{\theta}^{\prime}(a)$. Substituting Equation (39c) into the constraint that $\sigma_{r}^{\prime}(a)<\sigma_{\theta}^{\prime}(a)$, we obtain

[TABLE]

Substituting for $\mathcal{B}_{2}$ from Equation (40) and rearranging, we find that yield according to $\mathcal{F}_{z,\theta}$ cannot occur if

[TABLE]

This is true for sufficiently large $|\sigma_{a}|$,

[TABLE]

The maximum value of the right hand side of the above occurs in the limit $a\to 0$, and is given by $\displaystyle\frac{2|\sigma_{b}^{\prime}|}{2-\zeta(1-\Gamma)}$. Thus, the above is always satisfied if

[TABLE]

This inequality is safely satisfied for all $\zeta$ if

[TABLE]

since $\mathrm{max}(\zeta)=1$.

A.2.3 $\mathcal{F}_{z,r}=0$

Ensuring that $\mathcal{F}_{z,r}$ is always strictly negative will prevent yield according to $\mathcal{F}_{z,r}$. This is the appropriate yield function if $\sigma_{r}^{\prime}(a)<\sigma_{\theta}^{\prime}(a)<\sigma_{z}^{\prime}(a)<0$. The constrain that $\mathcal{F}_{z,r}<0$ rearranges to

[TABLE]

which is again always satisfied if $\alpha\nu>1$. For $\alpha\nu<1$, we appeal to the fact that $\mathcal{F}_{z,r}=0$ will only be the appropriate yield condition if $\sigma_{\theta}^{\prime}(a)<\sigma_{z}^{\prime}(a)$, and therefore if $\sigma_{\theta}^{\prime}(a)<\Gamma\sigma_{r}^{\prime}(a)$. This constraint implies that yield according to $\mathcal{F}_{z,r}$ cannot occur if

[TABLE]

which rearranges to

[TABLE]

Again, the right-hand side is maximised for $a\to 0$ and $\zeta\to{}1$, and is therefore satisfied for all $\zeta$ if

[TABLE]

which is the same constraint derived above for $\mathcal{F}_{z,\theta}$.

A.3 Summary

In summary, for fixed $q$ and $\zeta\equiv 1$, there are only two possible yield conditions: $\mathcal{F}_{\theta,r}=0$ and $\mathcal{F}_{r,\theta}=0$. The latter corresponds to collapse or negative yield, and there is a value of $q=q^{\star}$ above which this cannot occur. We showed above that the cylinder will yield exclusively according to the former for any $q\geq 0$, provided that the cylinder is sufficiently ‘strong’ relative to the compressive far-field stress—that is, if Equation (61) is satisfied.

For a fixed total stress at the inner radius, $\zeta\in[0,1)$, there are more possibilities; however, if $\alpha\nu>1$ and Equation (61) is satisfied, then yield will occur exclusively according to $\mathcal{F}_{\theta,r}$ for all $\sigma_{a}\leq 0$. Note that $\mathcal{F}_{r,\theta}$ and $\mathcal{F}_{z,\theta}$ both model negative yield (cavity collapse). Even if $\alpha\nu<1$, $\mathcal{F}_{\theta,r}=0$ remains the appropriate yield condition for all $\zeta$ and $a$, provided that $|\sigma_{a}|>\displaystyle\frac{2|\sigma_{b}^{\prime}|}{1+\Gamma}$.

Appendix B Additional Figures