Large poroelasto-plastic deformations due to radially outward ﬂuid injection

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 ﬂuid 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 addi- tional disturbance due to ﬂuid 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 mate- rials. 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 ﬁnd 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 ﬁnd that, for a ﬁxed total stress at the inner wall, materials experience larger deformations and generally larger stresses as they transition from impermeable to fully permeable.


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, watersaturated, 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 selfsupporting. 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 , Risnes et al., 1982;Detournay and Fairhurst, 1987;Detournay and Cheng, 1988;Wang and Dusseault, 1991a;1991b;1994;Wang, 1996 ), 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 piledriving and penetrometer testing ( Vesic, 1972;Carter et al., 1986;Yu and Houlsby, 1991;Yu., 20 0 0;Davis and Selvadurai, 20 02;Howell et al., 20 09 ). 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 ( Haimson and Fairhurst, 1969;Pye, 1973;Economides et al., 20 0 0 ). Fracking relies on the brittle failure of shale, but the mechanical properties of shale depend strongly on the composition ( Economides et al., 20 0 0;Rickman et al., 20 08;Britt and Schoeffler, 2009 ); many hydrocarbon-bearing shales have high clay content and may therefore have non-negligible ductility ( Vallejo, 1988;Daigle et al., 2014;Swift et al., 2014;Vega et al., 2014 ). 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., Hubbert and Willis, 1957;Wang and Dusseault, 1991a;Detournay, 2004 ).
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 or injection pushes 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.

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 ref and relaxed inner radius a ref (I in Fig. 1 ). We focus on the final steady state during fluid injection; throughout, we assume axisymmetry .

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 σ b ≤ 0 . Imposing this stress on the relaxed cylinder decreases its outer radius, b ref → b 0 ≤ b ref (I → II in Fig. 1 ); we assume that this compression is purely elastic.
We then assume that a cylinder of material of radius a 0 (relaxed radius a 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 → III in Fig. 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 → IV → V in Fig. 1 ). During injection, both radii expand such that a 0 → a > a 0 and b 0 → 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 Fig. 1. Conceptual model for fluid-driven deformation of an annulus with a pressurised inner cavity: (I) We consider a relaxed, solid cylinder of materialthis is our relaxed reference state. (II) This cylinder is then subjected to a uniform radial compressive effective stress of magnitude σ b < 0 , as well as a compressive axial effective stress that enforces plane strain; we assume that the cylinder deforms purely elastically. (III) We then consider the removal of a concentric cylinder of material, followed by the insertion of a permeable "casing", such that the stresses and deformation are unchanged-this is our deformed initial state. (IV) Pressurisation of the cavity leads to elastic and then elasto-plastic expansion of the annulus; we assume that plastic failure initiates at the inner cavity wall and expands outwards, with the inner plastic and outer elastic regions separated by a sharp boundary (dark grey and light grey, respectively). (V) For steady injection, the flow and deformation will eventually reach a steady state. The inner and outer radii are material surfaces that move over the course of the deformation (dashed lines). The material comprising the inner radius has reference position a ref , initial position a 0 , and steady-state position a . The material comprising the outer radius has reference position b ref , initial position b 0 , and steady-state position b . The elastic-plastic interface also evolves, but is not a material surface. The material comprising the elastic-plastic interface at steady state is located at radius s ; the reference and initial positions of this material are s ref and s 0 , respectively (dotted lines). Note that we take tension to be positive. has deformed elasto-plastically for a ≤ r < s and purely elastically for s < r ≤ 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 ( σ b and 0, 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 Fig. 1 ), but we necessarily derive the steady state relative to the relaxed reference state (IV in Fig. 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 Section 2.2 . We summarise the key aspects of the model, including kinematics, Darcy's law, and mechanical equilibrium, in Section 2.3 . We then discuss constitutive laws for the solid skeleton Section 2.4 .

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 ref as the characteristic length scale and the p -wave (oedometric) modulus 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 ˜ q or an applied dimensionless total stress at the inner cavity wall, ˜ where is Lamé's first parameter, k is the permeability (assumed to be constant and uniform), μ is the dynamic viscosity of the fluid, σ b is the radial effective stress (see Section 2.3 ) at the outer boundary, c is the cohesion between grains of the solid skeleton (see Section 2.4.2 ) and Q is the volume injection rate per unit length into the page. Note that the dimensionless relaxed outer radius is ˜ b ref ≡ 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 ( ˜ ) for convenience; we denote any further dimensional quantities with a breve ( ˘ ).

Kinematics, Darcy's law and mechanical equilibrium
Axisymmetry implies that the fluid velocity v f and solid displacement u s each have only one nontrivial component, such where r is the Eulerian radial coordinate and ˆ 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 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 φ f . We relate u s to φ f via ( Auton and MacMinn, 2017 ) which linearises under the assumption of infinitesimal strains to Additionally, mechanical equilibrium requires that ∇ · σ = 0 , where σ is the total Cauchy stress and where we have neglected inertia and the effect of gravity. The total stress can be decomposed as σ = σ − p I , where Terzaghi's effective stress σ 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 where σ r and σ θ are the radial and azimuthal (hoop) components of the effective stress, respectively. At steady state, conservation of mass for the fluid leads to 1 r 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 Combining Eqs. (6) and (7) and integrating leads to a direct relationship between q and the pressure drop p : For a more detailed derivation and discussion of the above aspects of the model, see Auton andMacMinn (2017, 2018) .

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 small deformations and therefore linear elasticity in the elastic region. In linear elasticity, the effective stress is related to the strain via σ r = ε r + ε θ , σ θ = ε r + ε θ , and σ z = (ε r + ε θ ) = (σ r + σ θ ) and the strain is related to the displacement via where ε r , ε θ , and ε z are the radial, azimuthal, and axial components of the strain, respectively, and the expressions for σ z and ε z are consequences of plane strain.

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 τ anywhere within the material must be strictly less than a specified fraction of the corresponding effective normal stress σ for the deformation to remain elastic. For perfect plasticity, τ is enforced to remain equal to this fraction of σ 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, | τ | ≤ −σ tan ϕ + c, where ϕ is the friction angle ( Davis and Selvadurai, 2002 ). We then express this condition in terms of the principal effective stresses 1 σ 1 ≥ σ 2 ≥ σ 3 by introducing a yield function F 1 , 3 , Plastic yield now corresponds to the condition F 1 , 3 = 0 . For F 1 , 3 < 0 , the material remains elastic. As soon as equality is first achieved, we enforce F 1 , 3 ≡ 0 locally thereafter. Taking ϕ = 0 ( α = 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 σ r , σ θ , and σ 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 σ 1 = σ θ and σ 3 = σ r ( i.e. , σ θ > σ z > σ r ) Wang and Dusseault (1991a,b) , such that the appropriate yield function is This ordering has been justified heuristically for an impermeable cylinder ( Yu and Houlsby, 1991 ). We show in Appendix A that Eq. (12) is indeed the appropriate yield function if the cylinder is sufficiently "strong"-that is, if y (b 2 − a 2 ) > −2 b 2 σ b and α > (1 + ) / , assuming that 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 ( Yu and Houlsby, 1991;Wang and Dusseault, 1991a;1991b ). We assume that the strain everywhere can be additively decomposed into elastic and plastic components ( Davis and Selvadurai, 2002 ), where the superscripts e and 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 where the overdots denote the material time derivative following the solid and A nonzero dilation angle ψ incorporates the fact that most granular materials undergo volumetric expansion (dilation) during plastic flow. For ψ = ϕ ( β = α), this flow law is said to be 'associated'; however, associated flow typically overestimates dilation for granular materials ( Yu, 20 0 0;Zhang and Salgado, 2010 ). For ψ = 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 This means that the outer boundary is a free boundary, at which the relevant kinematic condition is Lastly, we take p(b) = 0 . Note that we assume that a ≤ 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 .

The elastic-plastic interface
At the elastic-plastic interface, we enforce continuity of displacement and of radial effective stress, where s − denotes r = s as approached from within the plastic region ( a ≤ r < s ) and, likewise, s + denotes r = s as approached from within the elastic region ( s < r ≤ b ). In addition, we require that the effective stresses in the elastic region must be at the point of yield at r = s + , Note that Eqs. (17a) and (17b) , and the fact that F θ ,r ≡ 0 throughout the plastic region, together imply continuity of σ θ across r = s .

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 −σ 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 σ r (a ) = σ a ; we refer to this as an 'impermeable material'. For an 'infinitely' permeable skin, this would lead to an imposed fluid pressure p ≡ p(a ) = −σ 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' ζ . We define ζ ∈ [0, 1] such that the material can vary continuously from impermeable ( ζ ≡ 0) to fully permeable ( ζ ≡ 1). Pressurisation of the cavity then leads to two conditions at the inner boundary, We drive the system by imposing either σ a or the total flow rate q . The latter appears explicitly in Equation (7) . We express q in terms of σ a using Eqs. (8) and (18b) , (18c) Note that, for injection, we expect σ a < 0. For fully permeable materials ( ζ ≡ 1), in particular, it is convenient to impose q rather than σ a ( Section 5.2 ).
The inner boundary is also a material boundary. As such, the displacement must also satisfy a kinematic condition given by where a ref is the relaxed reference position of the material that eventually comprises the cavity wall.

Relaxed reference state to initial state
We now consider the transition from the relaxed reference state (I, Fig. 1 ) to the initial state prior to pressurisation (III, Fig. 1 ), which we take to be purely elastic compression under an imposed radial effective stress σ b in plane strain. As with a and b , we denote quantities in the initial state by a subscript '0'.
Using Eq. (9) to eliminate the stresses from Eq. (5) in favour of the displacement, and further setting q ≡ 0, purely mechanical linear-elastic deformation is governed by and the requirement of boundedness at the origin. This problem has solution For a rigorous treatment of the kinematics, we impose the outer boundary conditions at b 0 rather than at b ref , and we calculate the initial porosity field according to Eq. (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 We eliminate the Eulerian coordinate r from Eq. (21) in favour of the Lagrangian coordinate R via the kinematic relationship Finally, Eq. (18d) leads to We also derive a fully linearised version of the initial state by imposing the outer boundary conditions at b ref ≡ 1 and by calculating the initial porosity field according to Eq. (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 Note that Equation (18d) is not needed in this case.

Initial state to steady state
We now seek solutions to the poroelasto-plastic problem at steady state (V in Fig. 1 ). We denote the displacement field in this state by u s .

Elastic region
We again use Eq. (9) to eliminate the stresses in Eq. (5) in favour of the displacement, and we now also use Eq. (18c) to eliminate q in favour of σ a . This leads to an ordinary differential equation (ODE) in u s in the elastic region, 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 Eq. (3) and which requires four boundary conditions: and 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 Eq. (4) and apply Eq. (27c) at b ref ≡ 1. Eq. (27d) is not needed in this case.

Plastic region
In the plastic region, combining Eq. (12) with Eq. (5) and using Eq. (18c) to eliminate q in favour of σ a leads to which has solution To determine the displacement, we integrate Eq. (14) with respect to time to arrive at where f is a quantity whose material time derivative following the solid must vanish. Eq. (14) is a statement that the quantity βε p r + ε 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 ( ε p r = ε p θ = 0 ) and therefore f = 0 ; hence, it must be the case that f ≡ 0. We then decompose the plastic strains according to Eq. (13) , The elastic strains are always related to the effective stresses via the elasticity law. As such, Eqs. (9a) and (29) allow us to rewrite the right-hand side of Eq. (32) as and C 1 and C 2 are as defined in Eq. (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 ( Bažant, 1998;Yu., 20 0 0 ). For Hencky strains, the left-hand side of Eq. (32) becomes where D 1 and D 2 are as defined in Eq. (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, We denote this model by 'Q'. The 'Q' model can be further linearised by taking a ≈ a ref ; we denote this fully linearised model by 'L'.

Model summary
For rigorous kinematics throughout and Hencky strains in the plastic region, the full boundary value problem (BVP) is where σ r (s + ) , σ θ (s + ) , and σ r (b) can be expressed in terms of u s via Eq. (9) , and σ r (s − ) is given in Eq. (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 Eq. (18a) has already been imposed in the above. For comparison, we also consider linear strains throughout the entire domain by coupling Eq. (37) with Eqs. (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 Eqs. (37) and (38b) . 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 Eq.
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 rootfinding techniques. This solution provides a good first guess for the numerical solution of the other three models (see Section 4.4 ).

Solutions for imposed σ a and general ζ
For a given set of parameters, our poroelasto-plastic models will not be well posed for all values of σ a . Specifically, the material will remain elastic if σ a is sufficiently small; this purely poroelastic problem was studied extensively by Auton andMacMinn (2017, 2018) for fully permeable materials ( ζ ≡ 1). Conversely, the material will yield completely if σ a is sufficiently large. We denote the value of σ a at which the material will first yield by σ min a and the value of σ a at which the material will yield completely by σ max a . Note that σ max a < σ min a < 0 . We next consider σ min a and σ max a for all models ( Section 4.1 and Section 4.2 , respectively). We then derive an implicit solution to the LL model ( Section 4.3 ). Finally, we discuss the numerical scheme used here for solving the QL, QQ, and NQ models ( Section 4.4 ).
where B 1 and B 2 are unknown functions of a and b . The general form of the stresses in the elastic region is then For the poroelastic problem, we then apply the boundary conditions σ r (a ) = (1 − ζ ) σ a and σ r (b) = σ b to Eq. (39b) , to As discussed above, we assume that yield first occurs at r = s = a and that the yield function ( Eq. 12 ) is maximised at steady state. These assumptions require that where a min is the value of a associated with σ min a . Combining Eqs. (39) and (41) , we obtain For the LL model, we then take a min = a ref and b min = 1 , at which point σ min a is fully determined. For the QL model, we take b min = 1 and enforce the kinematic condition at the inner boundary to arrive at an implicit expression for a min , For the QQ and NQ models, we enforce kinematic conditions at both boundaries, where B 2 (a min , b min ) and σ min a (a min , b min ) are defined in Eqs. (40) and (43) , respectively. For ζ ≡ 1, all of the above can be solved for a min explicitly, whereas b min must then be determined via numerical rootfinding. This is not the case for ζ ≡ 1 , in which case the QQ and NQ models require simultaneous root-finding for both a min and b min .
where D 1 and D 2 are as defined in Eq. (33b) and continuity of displacement at the elastic-plastic interface leads to an expression for D 1 , where u s (s + ) is the displacement from the elastic solution at r = s . For the QL model, we take s = b max ≡ 1 . We can then calculate u s (s = 1) from the elastic solution ( Eq. (39) ), subject to σ r (1) = σ b and ασ θ (1) − σ b = y . This gives and hence so that σ min a (ζ ) is fully determined once a max is found via the implicit relation For the QQ model, rigorous treatment of the kinematics gives u s (b) = b − 1 and hence To find b max , we once again appeal to the elastic problem at Note that this result also applies to the NQ model, as we have not used the plastic flow law in its calculation. Hence, σ min a (ζ ) is fully determined once a max is found via the implicit relation with b max a known constant as defined in Eq. (53a) . For the NQ model, Eq. (38) cannot be solved analytically, prohibiting the derivation of an algebraic expression for a max .

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 ≡ 1) and the plastic region ( a ≡ a ref < r < s ), coupled with an implicit expression for the value of s . Solving the problem defined by Eqs. (37) and (38b) subject to boundary conditions (38c) -(38h) leads to where D 1 (a(s),b( s) ) and D 2 (a(s),b( s) ) are as defined in Eq. (33b) and (33c) , and The value of s is then determined via the implicit expression where C 1 (a(s),b( s) ) and C 2 (a(s),b( s) ) are as defined in Eq. (30) and B 2 is as defined in Eq. (40) . Note that s can take any value from a ref to 1. We solve Eq. (55) for s via numerical root-finding using MATLAB 's fzero . The solution is then fully prescribed by Eq. (54) .

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 Eqs. (38) for the NQ model and Eqs. (37) and (38b) -(38h) for the other three models. For all models, the governing ODE in the elastic region ( Eq. (38b) ) can be solved analytically. For the LL, QL, and QQ models, the governing ODE in the plastic region ( Eq. (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 andMacMinn (2017, 2018) .
To do so, we map the plastic region ( r ∈ [ a, s ]) and the elastic region ( r ∈ [ 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 ∈ [ a, b ] is discretised by 2 N − 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 2 N unknowns, comprising the displacement at the 2 N − 1 nodes as well as the value of s ; we therefore require a system of 2 N constraints.
In the plastic region, the governing ODE ( Eq. (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 ∈ [ a, s )), having already used the boundary condition σ r (a ) = (1 − ζ ) σ a in the derivation this ODE. In the elastic region, the governing ODE ( Eq. (38b) ) is second order and must be enforced at exactly N − 2 nodes. We enforce this for the N − 2 interior nodes ( r ∈ ( s, b )) and then impose σ r (b) = σ b at the outer boundary ( Eq. (38f) ). Lastly, we enforce two conditions at r = s : Continuity of radial effective stress ( Eq. (38d) ) and the requirement that the elastic stresses must satisfy the yield condition ( Eq. (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 ( Eqs. (38g) and (38h) ) and the current value of s .
The main structural difference between the models presented here and the models presented in Auton andMacMinn (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 ref and b ref ≡ 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 ζ and a given, small value of σ a < σ min a . 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 σ a > σ max a using numerical continuation. Note that σ max a < σ min a < 0 . This approach allows for the fact that the model behaviours diverge from one another as the driving strength increases (see Fig. 3 ). The same approach can also be used for increasing ζ at fixed σ a , since deformation increases with ζ (see Figs. 4 and 5 ).

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 Section 3.2.3 and Section 4 , we have developed solutions for four distinct model combinations: LL, QL, QQ and NQ. We fix all parameters based on the values discussed in Section 5.1 below, except for σ a (or q ) and ζ , which we vary.
For comparison with our poroelasto-plastic models, we additionally extend the purely poroelastic models from Auton andMacMinn (2017, 2018) by modifying the boundary conditions to incorporate the permeability-load parameter ζ . 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 Figs. 2-5 below. This comparison illustrates the impact of plasticity across different values of ζ .
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 Fig. 1 ) from their values in the final steady state (V in Fig. 1 ). For the displacement, for example, we consider the disturbance δu s : = u s − u s, 0 . Recall that these compressed initial states are presented in Section 3.1 . We subtract the L initial state from the LL and QL models, and the Q initial state from the QQ and NQ models.

Parameters
We adopt a set of parameter values motivated by boreholes in the subsurface. As such, we adopt M ∼ 50 GPa and ˘ ∼ 27 GPa (Young modulus Ȇ ∼ 25 GPa and Poisson ratio ν ∼ 0.36) and a porosity of φ ref f ∼ 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;Rickman et al., 2008;Britt and Schoeffler, 2009;Bobko et al., 2011 ). We further assume moderate internal friction and cohesion, ϕ ∼ 35 • and c ∼ 120 MPa, but very little dilation, ψ ∼ 0.3 • ( Vermeer and De Borst, 1984;Mandl, 2005;Bobko and Ulm, 2008;Bobko et al., 2011 ). We take −σ b ∼ 50 MPa, as appropriate for a depth of ∼ 2.5 km ( Neuzil, 1994;Islam et al., 2010 ). These Two poroelastic models (left column) and four poroelasto-plastic models (right column), all fully permeable ( ζ ≡ 1), at steady state for fixed q ≈ 0.0012 ( s ≈ 0.03). We show the disturbances due to fluid injection relative to a compressed initial state (III → V in Fig. 1 ): δφ (first row), δu s (second row), δσ r (third row), δσ θ (fourth row), and δp (last row). For clarity, we plot these results against the Lagrangian coordinate R (r) = r − u s and on a logarithmic horizontal scale. Note that the top two rows are also on a logarithmic vertical scale. Plasticity enables large deformations in the plastic region that amplify the importance of model choice there; model choice is unimportant in the purely poroelastic case, and in the elastic region of the poroelasto-plastic case. For reference, we plot the true steady-state values of the same quantities ( i.e. , without subtracting the initial state) in Fig. 6 in Appendix B .

Fig. 3.
Four summary quantities against q for the fully permeable case ( ζ ≡ 1), with the same colours as in Fig. 2 . We plot the disturbance to the inner radius δa (top left), the maximum disturbance in porosity max( δφ f ) (top right), the injection pressure p (bottom left), and the maximum disturbance in effective stress max( δσ ) (bottom right). We focus on the poroelasto-plastic models in the main plots and on the poroelastic models in the insets. The shaded grey regions indicate the values of q for which the inner cavity contracts. These results again highlight that plasticity leads to a large deviation from poroelastic behaviour, and that model choice is unimportant for the poroelastic models but very important for the poroelasto-plastic models-particularly with regard to the kinematics in the plastic region. Note that yield first occurs for q ≈ 3 . 1 × 10 −4 (the second corner in max( δσ )).
αν ≡ α 1+ ∼ 1 . 4 > 1 and y (1 − a 2 ) + 2 σ b ∼ 0 . 008 > 0 , where we take a ≈ a ref ; as such, the relevant yield condition is F θ ,r = 0 (see Appendix A ). Below, we consider the response of this material to different driving strengths (varying σ a or q ), and we additionally consider the transition from impermeable to fully permeable by varying ζ from 0 to 1.

Fully permeable: fixed q
We begin comparing poroelasto-plastic behaviour with purely poroelastic behaviour in the context of a fully permeable material ( ζ ≡ 1; Figs. 2 and 3 ). Recall that, for this model, we drive the system with a fixed flow rate q instead of a fixed pressure difference p ≡ −σ a (see Section 2.5.3 ). For each model, a given value of q will correspond to a specific steadystate value of p ; however, this value will be different for each model. Note that, as with σ a , there exists a minimum value q min at which yield first occurs and a maximum value q max at which the material yields completely. These values can be calculated via q min = −σ min (1) and σ max a (1) are defined in Eqs.
(43) and (47) , respectively. In Fig. 2 , we consider q = 0 . 0012 . For these parameters, the LL model has q min ≈ 3 . 1 × 10 −4 and q max ≈ 0.0033, so that our chosen value of q is about 4 q min and about q 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) ≡ 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 Fig. 6 in Appendix B . Note that δp ≡ p because p 0 ≡ 0 for all models.
Although this value of q is substantially higher than q 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 ( Fig. 2 , left column).
The disturbances δφ f , δσ θ , and δp all have maxima at the cavity wall, with δσ θ falling off relatively steeply from this value as R increases. In contrast, δu s and δσ r have internal maxima near the outer boundary ( R ∼ 0.7) and near the inner boundary ( R ∼ 4 × 10 −4 ), respectively. The value of δσ r is pinned to −σ b at the inner boundary by construction since σ r (a ) = 0 and σ r, 0 ≡ σ b , and vanishes at the outer boundary since σ r (b) = σ r, 0 = σ b . Note that δφ f is strictly positive, meaning that the porosity everywhere has increased relative to the compressed initial state.
Upon introducing plastic yield ( Fig. 2 , right column), all quantities except δ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, max u s r ≈ 130 for the poroelasto-plastic LL model whereas max u s r ≈ 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, max ( u s r ) ≈ 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 Figs. 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 δφ f than any of the models. This discrepancy is due to the fact that d u s /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 δφ f > 1 near the inner boundary, highlighting the fact that linearising the kinematics in the plastic region can lead to nonphysical solutions 3 (see Fig. 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 δσ r 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 δσ θ 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 d d r (δσ θ ) is discontinuous across r = s . It is straightforward to show from the boundary conditions at r = s, the yield condition, and mechanical equilibrium that δu s , δσ r , d d r (δσ r ) , δσ θ , and δp , and d d r (δ p) must all be continuous across r = s . However, d d r (δu s ) and therefore also δφ f are weakly discontinuous across r = s (not readily visible in the figure). The discontinuity in d d r (δσ θ ) is particularly prominent because it occurs at the maximum value of σ θ , whereas those in δφ f and d d r (δu s ) occur in regions where δφ f and δu s are themselves small. For the poroelastic L and poroelasto-plastic LL models, δp decreases exactly logarithmically with r because d p d r = − 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 Fig. 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 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 Fig. 3 is q ≈ 0.0012, which is the value used in Fig. 2 .
For q 10 −4 ( Fig. 3 , grey band), the disturbance in the inner radius δa is negative ( Fig. 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 ≥ 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 poroelastoplastic models deform increasingly more than the poroelastic models, and the ordering of the models is the same as in Fig. 2 .
For small flow rates, the maximum porosity disturbance max( δφ f ) is small ( Fig. 3 , top right). Much like δa , max( δφ 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 max( δφ 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 p ≡ p ( a ) ≡ δp ( a ) increases linearly with q before yield, and continues to increase linearly with q for the LL model after yield ( Fig. 3 , bottom left). For the other poroelasto-plastic models, 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 ( δσ ) is constant and equal to −σ b = 10 −3 ( Fig. 3 , bottom  right). This is an artefact of the fact that σ a ≡ 0 , such that δσ a ≡ σ a − σ b = 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 10 −4 (outside the grey region). All models exhibit a corner in max ( δσ ) before yield due to a change in which stress component exhibits this maximum: max ( δσ ) is equal to δσ r (a ) for q 1 . 4 × 10 −4 and to δσ θ (a ) for q 1 . 4 × 10 −4 . The poroelasto-plastic models exhibit a second corner at yield ( q ≈ 3 . 1 × 10 −4 ) due to a change in the location of this maximum: max ( δσ ) is equal to δσ θ (a ) before yield and to δσ θ (s ) after yield.
In Fig. 8 in Appendix B , we plot the same four quantities as in Fig. 3 , but against −σ b 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 andMacMinn (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 −σ b varies. Additionally, all quantities converge smoothly towards poroelastic behaviour as −σ b increases. As a result, it seems reasonable to extrapolate the conclusions of Auton andMacMinn (2017, 2018) for a poroelastic cylinder with deformationdependent permeability to the poroelasto-plastic cylinders considered here. Fig. 2 a of Auton and MacMinn (2017) suggests that deformation-dependent permeability leads to a much smaller 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.

Impermeable to fully permeable: fixed σ a
In Figs. 4 and 5 , we vary the permeability-load parameter ζ , transitioning the model from impermeable ( ζ ≡ 0) to fully permeable ( ζ ≡ 1) for a fixed σ a = −0 . 0075 . This value of σ a is between σ min a and σ max a for all values of ζ , such that a < s < b . Recall that 0 ≤ ζ < 1 is analogous to a thin, weak, low-permeability membrane on the inner cavity wall, the permeability of which decreases as ζ decreases; varying ζ 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 Fig. 4 , we plot the same quantities as in Fig. 2 for four values of ζ ranging from nearly impermeable ( ζ ≈ 0) to nearly fully permeable ( ζ ≈ 1). For reference, we again plot the true steady-state values of the same quantities ( i.e. , without subtracting the initial state) in Fig. 9 in Appendix B .
Note that, because we now drive these models with fixed σ a , each model will produce a different value of q . We found earlier that increasing nonlinearity led to lower values of p for a given value of q ( Fig. 2 ). Here, we then expect that increasing nonlinearity will lead to a higher value of q for a given value of σ a , and therefore to larger deformations and stresses. This suggests that the relative ordering of the models should be reversed in Fig. 4 relative to Fig. 2 , which is indeed the case.
For the poroelastic models (left column), all quantities increase with ζ . This suggests that, for a given imposed cavity pressure ( −σ a ), a more permeable material will experience larger deformations and larger stresses than a less permeable material. The three quantities δφ f , δσ θ , and δp depend quantitatively but not qualitatively on ζ , varying monotonically from a maximum value at the inner radius to a minimum value at the outer radius; the former increases monotonically with ζ , whereas the latter is insensitive to ζ . The two quantities δu s and δσ r vary both quantitatively and qualitatively with ζ . The disturbance in displacement δu s is relatively insensitive to ζ at the cavity wall, implying that permeability makes very little difference to the size of the cavity. However, δu s is increasingly sensitive to ζ for larger values of R , such that the maximum displacement shifts from R = a 0 for small ζ to near R = 1 for moderate to large ζ . Although the δu s curves appear to get closer together as ζ increases, this is an artefact of the logarithmic vertical scale. The disturbance in radial effective stress δσ r vanishes at the outer boundary by construction, and is fixed to δσ r (a ) = (1 − ζ ) σ a − σ b at the inner boundary. The latter value increases with ζ , transitioning from compressive to tensile, and the rest of the curve essentially pivots about the constant outer boundary value. As a result, δσ r is monotonic in R for small ζ , with a minimum at the inner boundary and a maximum at the outer boundary, but transitions to being nonmonotonic in R for large ζ , 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 ζ increases (right column). This is because plasticity leads to larger deformations, as do larger values of ζ . The plastic region itself also grows larger as ζ increases (see δσ θ ). For this value of σ a and sufficiently small ζ , the deformations are small enough that the models are essentially indistinguishable. For this value of σ a and larger values of ζ , the deformations grow sufficiently large that significant differences emer ge between the models. For δφ f , δσ θ , and δp , the consequences of increasing ζ are again mostly quantitative: The deformations and stresses grow larger. For all models except the LL model, δσ r and δσ θ , exhibit weak non-monotonicity in ζ for some intermediate values of R as ζ approaches 1. Note that the QL and QQ models predict very large values of δφ f for ζ ≈ 1 relative to the other models and other values of ζ . For δu s and δσ r , the changes are again more complex. Plasticity leads to much larger displacements at the inner boundary as ζ increases, but to very little change in the displacement at the outer boundary. The maximum in δu s shifts from the inner boundary for all models for small ζ to the outer boundary for all models for intermediate ζ , and back to the inner boundary for the more nonlinear models for large ζ .
In Fig. 5 , we plot the same quantities as in Fig. 3 , but against ζ for fixed σ a = −0 . 0075 . As also illustrated in Fig. 4 , it is clear that model choice is unimportant for the poroelastic models for all ζ , and for the poroelasto-plastic models for sufficiently small ζ (for this set of parameters, ζ ࣠ 0.55). Model choice eventually becomes important for the poroelastoplastic models because plasticity enables large deformations and these deformations increase monotonically with ζ ( Fig. 5 ,  Fig. 4 . As in Fig. 2 , we plot two poroelastic models (left column) and four poroelasto-plastic models (right column) at steady state, here for fixed σa = −0 . 0 075 and ζ ∈ {0.0 02, 0.335, 0.665, 0.998}. Colours are the same as in Figs. 2 and 3 . For the poroelastic models, all quantities increase with ζ and there is again no discernible difference between the models. For the poroelasto-plastic models, all quantities increase with ζ to a point, after which the effective stresses from the more nonlinear models become weakly nonmonotonic in ζ . The difference between models grows with ζ as the deformations grow larger (except for the QL and QQ models, which are again indistinguishable).
top left). As also illustrated in Fig. 4 , the more nonlinear models again predict larger deformations. The two poroelastic models predict that max ( δφ f ) increases approximately linearly with ζ . These predictions are larger in magnitude than those of the poroelasto-plastic models until ζ ≈ 0.55 (QL and QQ) or ≈ 0.8 (LL and NQ), at which point the poroelasto-plastic models increase strongly following a corner where max ( δφ f ) shifts from an internal value to the value at the inner radius. The flow rate q increases monotonically with ζ for all models, with the more nonlinear models predicting larger flow rates  Fig. 3 , we plot four summary quantities, now against ζ for fixed σa = −0 . 0075 . We again focus on the poroelasto-plastic models in the main plots and on the poroelastic models in the insets. Colours are the same as in Figs. 2-4 . Note that plasticity enables much larger displacements, increasingly so as ζ increases. Whereas δa , max ( δφ f ), and q increase monotonically with ζ for all models, max ( δσ ) is nonmonotonic in ζ with an internal maximum near ζ = 1 for the QL, QQ, and NQ models. Note that, for q , the poroelastic models are behind the LL model (red) on the main plot and, for max ( δσ ), the poroelastic models predict much larger maximum stresses that exceed the vertical scale of the main plot. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) as ζ approaches 1. Note that q = − ζ σa ln b a for all models, and is thus pinned to q = 0 for an impermeable material ( ζ ≡ 0) and to q = p ln b a for a fully permeable material ( ζ ≡ 1). This dependence on a and b links q to the displacement field, and hence to model choice, for ζ > 0. Note that, for the L and LL models, q = ζ σa and is therefore independent of deformation. The maximum disturbance in effective stress max ( δσ ) increases monotonically with ζ for the L, Q, and LL models, and becomes weakly nonmonotonic in ζ for the QL, QQ, and NQ models as ζ approaches 1. This nonmonotonicity is more pronounced for the NQ model. The fact that max ( δσ ) is maximised for a particular value of ζ could be important in applications such as hydraulic fracturing in ductile shales or for borehole integrity. Note that the ordering of models for max ( δσ ) is reversed relative to the ordering of models for δ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 ζ ( Fig. 5 , top left and bottom right), whereas all models predict similar flow rates for most values of ζ .

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 thickwalled 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 ( Fig. 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' ζ that enables a smooth transition between these limiting states by (essentially) introducing a thin, weak, low-permeability skin at the cavity wall. In Figs. 2 and 3 , we considered a fully permeable material ( ζ ≡ 1) for a range of injection rates q . In Figs. 4 and 5 , we considered a fixed total stress at the inner radius σ a as ζ transitions from 0 to 1. We showed that the amount of deformation increases with ζ for a given value of σ a . The maximum tensile effective stress exhibits a maximum at an intermediate value of ζ 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 ζ , the choice of poroelasto-plastic model becomes increasingly significant as ζ 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 ζ . 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 ).
The research materials supporting this publication can be accessed by contacting the corresponding author.
σ z = 1 + σ r + σ θ ≡ ν σ r + σ θ , where ν ≡ 1+ is the Poisson ratio ( cf. Section 5.1 ). Recall that ν ∈ (0 , 1 2 ) for most physical materials, and that σ 1 ≥ σ 2 ≥ σ 3 by definition. We then have three possible constraints based on the signs of σ r and σ θ : If σ r , σ θ > 0 , then σ 1 = σ z ; if σ r , σ θ < 0 , then σ 3 = σ z ; and if σ r and σ θ have different signs, or if one of them is zero, then σ 2 ≡ σ z . The six possible yield functions are then F θ ,r : = ασ θ − σ r − y, where the material remains elastic for F 1 , 3 < 0 and yields when F 1 , 3 = 0 . We must then determine which of these yield functions first reaches zero. As stated in Section 2.4.2 , we assume that yield first occurs at r = a min ; that, prior to yield, the yield function 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 σ i (a min ) , where i ∈ { r, θ , z } ( i.e. , the effective stresses at the inner boundary at the point of first yield).

A1.1. F r,θ = 0
For F r,θ = 0 to be the appropriate yield condition, it must be the case that σ θ (a ) ≤ σ z (a ) ≤ σ r (a ) = 0 . For the L and Q models with ζ = 1 , it is straightforward to show that (see Section 4.1 ) where B 2 is given by .

(A.4)
We then rearrange the condition F r,θ < 0 , for which deformation remains elastic, to give . (A.5) 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 ≥ 0 (no suction/extraction). As a result, ensuring that q < 0 would then imply that F r,θ < 0 for q ≥ 0. Rearranging the requirement that q < 0 leads to the constraint that negative yield cannot occur for all relevant values of q if (A.6) meaning that the cylinder must be sufficiently 'strong' relative to the compressive far-field stress σ b . This constraint is satisfied for the parameter values used here (see Section 5.1 ). Hence, we can state conclusively that, for a sufficiently strong and fully permeable cylinder, F θ ,r is the correct yield condition during fluid injection ( q ≥ 0).
Note that, for σ r (a ) < 0 , the axial effective stress σ z (a ) can never be the minimum principal stress (see discussion above Eqs. A.2 ). Thus, we only have three alternative yield conditions to consider.

A2.1. F r,θ = 0
Ensuring that F r,θ is always strictly negative will prevent yield according to F r,θ . Using Eq. (39) leads to where B 2 is defined in Eq. (40) . Following the same procedure that leads to Eq. (43) , we find that yield according to F r,θ will not occur for σ a < σ min a ;r,θ = (A.8) meaning that σ min a ;r,θ is the least compressive value of σ a that would prevent negative yield according to F r,θ . Note that the σ min a presented in Eq. (43) is σ min a ;θ ,r , corresponding to positive yield according to F θ ,r .

ln
We focus on injection, so it must be the case that σ a < 0. As a result, ensuring that σ min ,r,θ a > 0 would then imply that F r,θ < 0 for all relevant σ a . This constraint simply requires that Eq. (A.6) must be satisfied since the denominator of Eq.

A3. Summary
In summary, for fixed q and ζ ≡ 1, there are only two possible yield conditions: F θ ,r = 0 and F r,θ = 0 . The latter corresponds to collapse or negative yield, and there is a value of q = q above which this cannot occur. We showed above that the cylinder will yield exclusively according to the former for any q ≥ 0, provided that the cylinder is sufficiently 'strong' relative to the compressive far-field stress-that is, if Eq. (A.6) is satisfied.
For a fixed total stress at the inner radius, ζ ∈ [0, 1), there are more possibilities; however, if αν > 1 and Eq. (A.6) is satisfied, then yield will occur exclusively according to F θ ,r for all σ a ≤ 0. Note that F r,θ and F z,θ both model negative yield (cavity collapse). Even if αν < 1, F θ ,r = 0 remains the appropriate yield condition for all ζ and a , provided that | σ a | > 2 | σ b | 1+ . Figures  Fig. 6. As Fig. 2 , but without subtracting the compressed initial state. Note that, unlike in Fig. 2 , the top two rows are no longer on a logarithmic vertical scale. Fig. 7. For all four poroelasto-plastic models-LL (top left), QL (top right), QQ (bottom left), and NQ (bottom right)-we plot the total thickness of the annulus b − a (magenta), the thickness of the elastic region b − s (green), and the thickness of the plastic region s − a (blue) against flow rate q for a fully permeable material ( ζ ≡ 1). Note that for the LL, QL, and QQ models, q is between q min ≈ 3 × 10 −4 , the smallest q that induces yield in all models, and q max ≈ 3 . 3 × 10 −3 , the first value of q for which the material is entirely yielded. Note that the range of q is reduced for the NQ model relative to the others. The LL model is the only one of these that linearises the kinematics in the plastic region, where deformations are largest. As a result, the LL model predicts non-physical behaviour for sufficiently large q ( i.e., s < a and b < a ); the other predict physical results for all q until the material has yielded entirely. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 8. As Figs. 3 and 5 , but plotted against varying −σ b for a fully permeable material ( ζ ≡ 1). This illustrates the transition from an unconstrained cylinder ( σ b ≡ 0 ) to a highly constrained cylinder, providing a link to our previous work on fully permeable unconstrained poroelastic cylinders ( Auton and MacMinn, 2017;2018 ).