On the Derivation of Closed‐Form Expressions for Displacements, Strains, and Stresses Inside Poroelastic Reservoirs

We critically review the derivation of closed‐form analytical expressions for elastic displacements, strains, and stresses inside a subsurface reservoir undergoing pore pressure changes using inclusion theory. Although developed decades ago, inclusion theory has been used recently by various authors to obtain fast estimates of depletion‐induced and injection‐induced fault stresses in relation to induced seismicity. We therefore briefly address the current geomechanical relevance of this method, and provide a numerical example to demonstrate its use to compute induced fault stresses. However, the main goal of our paper is to correct some erroneous assumptions that were made in earlier publications. While the final expressions for the poroelastic stresses in these publications were correct, their derivation contained conceptual mistakes due to the mathematical subtleties that arise because of singularities in the Green's functions. The aim of our paper is therefore to present the correct derivation of expressions for the strains and stresses inside an inclusion and to clarify some of the results of the aforementioned studies. Furthermore, we present two conditions that the strain field must satisfy, which can be used to verify the analytical expressions.


Current Relevance
Pore pressure changes in subsurface reservoirs yield a poro-elastic response in the reservoir and its surroundings (Biot, 1941;Geertsma, 1973;Segall, 1985;Wang, 2000).Numerical computation of the induced stresses can be demanding, as a fine mesh resolution is required to capture the sharp stress variations at the reservoir boundaries, while singularities in the underlying mathematical structure may frustrate convergence.Therefore, analytical methods, developed several decades ago, have recently regained popularity to provide a fast estimate of injectioninduced or depletion-induced displacements, strains, or stresses (Jansen et al., 2019;Lehner, 2019;Wu et al., 2021).Moreover, these methods can nowadays be coupled with semi-analytical or numerical models to efficiently compute induced fault stresses as a basis for subsequent simulation of fault slip and seismic hazard assessment (Acosta et al., 2023;Barbosa et al., 2022;Cornelissen & Jansen, 2023b;Jansen & Meulenbroek, 2022;Smith et al., 2022;van Wees et al., 2018van Wees et al., , 2019)).

Inclusion Theory and the Nucleus of Strain Concept
Analytical expressions for the linear elastic displacements, strains, and stresses induced by pore pressure changes in a subsurface reservoir can be derived under some simplifying assumptions.Such expressions are mainly derived from either the nucleus of strain concept, based on potential theory, or inclusion theory.The nucleus of strain concept (Love, 1927) was introduced by Goodier (1937) for thermoelasticity and was later applied to poroelasticity by Geertsma (1966Geertsma ( , 1973)).In this concept, the displacement field is equated to the gradient of a scalar potential, denoted as the displacement potential.Substituting the displacement potential into the mechanical equilibrium equations yields the Poisson equation, for which a solution exists in integral form.Inclusion theory was introduced by Eshelby (1957) to compute the stresses around elliptical inclusions.Eshelby (1957) computes the displacement field due to a region, the inclusion, undergoing a change in shape and size through a series of imaginary steps that involve cutting, transforming, and re-inserting the inclusion.While the derivations of the governing equations in the two approaches are conceptually different, they yield identical expressions for the displacements, strains, and stresses.
Inclusion theory and the nucleus of strain concept were later adopted to estimate subsidence and stress fields outside elliptical subsurface reservoirs (Geertsma, 1973;Segall, 1985Segall, , 1989Segall, , 1992;;Segall et al., 1994).Segall and Fitzgerald (1998) considered the stresses inside ellipsoidal reservoirs, where the strains are uniform inside the inclusion as per Eshelby's result (Eshelby, 1957).More recently, Soltanzadeh and Hawkes (2008), Jansen et al. (2019), Lehner (2019), and Wu et al. (2021) also considered the non-uniform stresses inside non-elliptical reservoirs, such as rectangular, triangular, and trapezoidal reservoirs.Soltanzadeh and Hawkes (2008) did not present analytical expressions for the stresses, while the other three studies did.Although all these studies presented correct final solutions for the stress field, the derivations of Jansen et al. (2019) and Wu et al. (2021) contained errors.Furthermore, mathematical subtleties arise in all three studies.The aim of this paper is to clarify and correct the derivation of closed-form expressions for the stresses inside an inclusion and in particular address the derivations by Jansen et al. (2019) and Wu et al. (2021), but also to highlight the mathematical subtleties that arise when considering strains and stresses inside an inclusion.

Organization of the Paper
In Section 2 we recapitulate the derivation of expressions for poroelastic displacements, strains, and stresses using inclusion theory, and in Section 3 we present the resulting expressions for rectangular and triangular elements.A critical review in Section 4 forms the heart of the paper in which we clarify the erroneous assumptions in earlier derivations.In Section 5 we provide a numerical example to illustrate the use of the theory to compute induced stresses in a faulted reservoir, whereafter Section 6 concludes the paper.Two appendices provide additional mathematical details.

General Procedure
Here, we will use inclusion theory to derive closed-form expressions for the displacements, strains, and stresses induced by pore pressure changes, but identical expressions can be obtained from potential theory and the related 'nucleus of strain concept' (Barbosa et al., 2022;Geertsma, 1973;Goodier, 1937;Lehner, 2019;Lehner et al., 2005;Rudnicki, 2002Rudnicki, , 2011;;Smith et al., 2022;Timoshenko & Goodier, 1970;Vasco & Johnson, 1987;Wang, 2000).For more background information and a thorough derivation of the governing equations in inclusion theory, readers are referred to Mura (1987) and Rudnicki (2011).The description in this section is partly based on the Supplementary Information belonging to Jansen et al. (2019).
We are interested in the elastic strain and stress fields caused by a change of pore pressure p in a region, denoted as the inclusion Ω, inside an infinite three-dimensional elastic domain.The boundary of the inclusion is defined by Γ.We consider linear elasticity theory, as is the norm in inclusion theory.Both the inclusion and the surroundings are assumed to have homogeneous and isotropic elastic properties, which are the same inside and outside the inclusion.An increase in pore pressure in the inclusion causes a reduction in effective stress in its matrix and consequently an elastic expansion of the inclusion.The induced stresses are assumed to vanish at infinity.As first described by Eshelby (1957), the displacements, strains, and stresses in and around the expanding inclusion can be computed by following a series of imaginary steps: 4. Re-attach the inclusion to its surroundings and remove the forces (i.e., apply the forces in opposite direction to the entire infinite solid).
These steps are visualized in Figure 1 for a rectangular inclusion, but the approach is applicable to inclusions of all shapes.
Only the incremental stresses (i.e., the stresses induced by changes in pore pressure) are computed with this method.These incremental stresses can be added to the background stresses (or initial stresses) to obtain the total stress field, as we later show in this paper.We assume that the incremental displacements, strains, and stresses are initially zero.
We assume that the change in pore pressure p is uniform in Ω.If the inclusion would be allowed to expand (for an increase in pore pressure) or shrink (for a drop in pore pressure) stress-free, as in step 2 above, it would experience a uniform dilational deformation.The corresponding strain tensor ϵ * ij , often referred to as the transformation strain, stress-free strain, or eigenstrain, has components (Wang, 2000, p. 38) where α = 1 c s /c m is Biot's coefficient, with c s the compressibility of the solid phase and c m the drained compressibility of the porous medium, K is the drained pore bulk modulus, δ ij is the regular Kronecker delta, and δ Ω is a modified Kronecker delta which equals 1 inside the inclusion and 0 outside.No strains would develop outside the inclusion as a result of step 2.
The inclusion can be brought back to its original shape and size by applying surface tractions with magnitude σ * ij n j to the boundary of the inclusion, where σ * ij is obtained from the transformation strain via Hooke's law where λ and μ are Lamé's first and second parameter, respectively, repeated indices imply summation, and n i denotes the i-component of the outward normal vector of the inclusion boundary.The stress tensor σ * ij is generally referred to as the eigenstress (not to be confused with the principal stresses, which are computed from the eigenvalues of the stress tensor).We employ the sign convention of continuum mechanics, such that positive normal stresses correspond to tension.At the end of step 3, the strain then equals zero everywhere, while the stress inside the inclusion equals σ * ij and vanishes outside of the inclusion.
Next, we reconnect the inclusion to the surrounding material and remove the distributed forces by applying their negatives.The resulting displacement field can be computed from where u i is the displacement in the i-direction and g ij are Green's functions, which denote the displacement in the idirection at point x due to an applied force in the j-direction at point x′ (Love, 1927).Expressions for the Green's functions will be given in the next section.Applying the divergence theorem allows us to transform the surface integral in Equation 3 into a volume integral The displacement field described by Equation 4 is due to a distribution of centers of dilatation (i.e., double forces without moment in all perpendicular directions), also referred to as nuclei of strain.
Once the displacement field is known from Equation 4, the total strains ϵ ij can be obtained from the compatibility equation The stresses can be obtained from the strains by Hooke's law.However, for points inside the inclusion Ω it is necessary to add the eigenstress defined in Equation 2. The stresses are then given by which is equivalent to the poroelastic form of Hooke's law (Wang, 2000).

Displacements
Up to this point, the equations that have been presented are valid for the general three-dimensional case.From this point onward, we will restrict ourselves to two-dimensional plane strain problems, similar to Jansen et al. (2019) and Wu et al. (2021).The same method is applicable to three-dimensional problems, but the expressions for the Green's functions and resulting integrals are different.For the plane strain case, the transformation strain is given by (Wang, 2000, p. 151) where ν is Poisson's ratio which is related to Lame's parameters as ν = λ 2(λ + μ) , such that the eigenstress still equals σ * ij = αpδ ij δ Ω (Wang, 2000).
Under plane strain conditions, the Green's functions are given by Mura (1987) as where √ is the distance between the points (x, y) and (x′, y′) (Figure 2).
with D = (1 2ν)αp 4π(1 ν)μ for ease of notation.These equations are also presented by Jansen et al. (2019), Lehner (2019), andWu et al. (2021).We define G x and G y as the double integrals of the Green's functions Outside of the inclusion, the Green's functions g x and g y can be integrated as usual.However, Wu et al. (2021) correctly pointed out that Fubini's theorem does not hold inside the inclusion.Fubini's theorem states that it is possible to compute a double integral by using an iterated integral, independent of the order of integration, if the double integral yields a finite answer when the integrand is replaced by its absolute value: For a proof, see, for example, Ghorpade and Limaye (2010), and for examples of the use of Fubini's theorem, see Stewart et al. (2021).Due to the singularity at (x, y) = (x′, y′), condition (Equation 15) is not met in Equations 13 and 14.Therefore, we must take special care of points inside the inclusion.
To solve the improper double integral over the inclusion domain Ω, we first remove a square region S δ of size 2δ × 2δ around the singularity such that the resulting domain Ω δ = Ω\S δ does not contain the singularity.The integral in Ω δ can then be computed directly with an iterated integral, as this domain does not contain the singularity.As an example, we consider the horizontal displacement for a point (x, y) inside the rectangular inclusion Ω The double integral over Ω δ is then given by (18) Since Fubini's theorem holds on the domain Ω δ , changing the order of integration does not affect the result.Furthermore, the result is independent of δ.Next, we consider the improper integral on the domain S δ .In Appendix A, we show that the improper double integral converges absolutely and equals zero This result is independent of δ.The solution to the integral in Equation 13is then given by The same approach can also be applied to the vertical displacement by replacing g x with g y in the equations above, as well as to shapes different from the rectangular inclusion considered in this example.

Strains and Stresses
Once the displacement field is known, the total strains ϵ ij are obtained by differentiation as Note that we take the derivative of the displacements outside of the double integral.Previous studies placed the derivative under the integral sign to obtain Green's functions for the strains and stresses (Barbosa et al., 2022;Jansen et al., 2019;Soltanzadeh & Hawkes, 2008;Wu et al., 2021).However, the resulting double integrals over the Green's functions for the stresses are not convergent (Ghorpade & Limaye, 2010).
Since the Green's functions for the displacements are discontinuous at the singularity, their derivatives are undefined inside the inclusion.Therefore, the order of differentiation and integration is not interchangeable for points inside the inclusion.In other words, Leibniz's integral rule cannot be applied for these points (Protter & Morrey, 1985).Indeed, it was already pointed out by Mura (1987, p. 12) that switching the order of differentiation and integration should generally be avoided.Placing the derivative under the integral sign is only valid for points outside the inclusion.The same conclusion follows from potential theory, in which the first derivatives of a potential (corresponding to displacements) obey Leibniz's integral rule for points inside the inclusion, whereas the second derivatives (corresponding to strains) do not (D'Urso, 2013).
Finally, we compute the stresses from the elastic strains using Hooke's law (Equation 6), which yields

Conditions to Verify Results
As obtaining closed-form expressions for the strains and stresses is prone to errors due to the mathematical subtleties surrounding the improper integrals, it is worthwhile to determine some requirements that the solution must adhere to.Here, we present two conditions that the strain field must satisfy.

Strain Discontinuity
As the Green's functions are solutions to the mechanical equilibrium equations, the normal traction must be continuous everywhere.While stresses develop both inside and outside the inclusion due to the change in pore pressure, the eigenstress is non-zero inside the inclusion and zero outside.Therefore, there must be a jump in the strains to maintain continuity of the normal traction across the interface of the inclusion.The magnitude of the required strain jump can be derived purely from the geometry of the inclusion boundary, independently from the analytical solution for the strains (Mura, 1987, p. 39).This approach can be used to verify our expressions for G xx , G yy , and G xy .The normal traction should be continuous across the boundary of the inclusion where T x and T y are the components of the normal traction vector, n x and n y are the components of the unit vector normal to the inclusion boundary (pointing outwards), and Δ indicates the jump in the variable from just inside the inclusion to just outside the inclusion (e.g., Δσ xx = σ out xx σ in xx ).We use Hooke's law to write the stresses in Equation 27 in terms of the strains The jump in strains can be expressed as where β i represents the jump in the derivative of the displacements (Mura, 1987, p. 39).Substituting this into Equation 28 and rearranging yields For a given boundary with known normal vector, this yields two equations for the two unknowns β x and β y .The general solution is which can be used to verify the expressions for the strains.

Volumetric Strain
The displacements, strains, and stresses must satisfy the mechanical equilibrium equation, which is given in terms of displacements by Journal of Geophysical Research: Solid Earth 10.1029/2023JB027733 where the modified Kronecker delta δ Ω appears because the change in pore pressure equals p inside the inclusion and 0 outside.Following Geertsma (1973), we can take the divergence of this equation and use the definition of the volumetric strain ɛ = ϵ kk = u k,k to obtain Hence, the volumetric strain must satisfy the condition Equation 37 can be used as another method to verify the expressions for the strains.This condition can be derived more directly from potential theory (Lehner, 2019;Wang, 2000).

Rectangular Inclusion
We now present expressions for the scaled displacements and strains.We consider a rectangle with corners (p, r), (q, r), (q, s), and (p, s) as shown in Figure 3.All the expressions presented in this section are valid for points inside and outside the inclusion.
Computing the double integrals of Equations 13 and 14 yields where atan is the arctangent function.
The scaled strains are obtained by taking the spatial derivatives of G x and G y as defined in Equations 21-23.This yields Figure 3. Geometry of a rectangular inclusion with corner points ( p, r), (q, r), (q, s), and ( p, s).

Triangular Inclusion
We consider a right triangle with vertices (o, r), (p, r), and (p, s) with angle θ = atan ( s r p o ), as shown in Figure 4. We assume that the hypotenuse lies along the line y′ = x′ tan θ.Using the fact that o = r cot θ, p = s cot θ, r = o tan θ, and s = p tan θ, the integrated Green's functions for the displacements due to a triangular inclusion are Taking the spatial derivatives of the integrated Green's functions yields the scaled strains Analytical solutions for the displacements, strains, and stresses for various other inclusion shapes, including elliptical inclusions, are given by Mura (1987, p. 80).

Critical Review
The expressions for the scaled displacements for a rectangular inclusion derived in this paper are equivalent to those presented by Jansen et al. (2019).Expressions for the scaled displacements for triangular inclusions were not given by Jansen et al. ( 2019).The expressions for the scaled strains are equivalent to the ones presented by Jansen et al. ( 2019) for both rectangular and triangular inclusions.Wu et al. (2021) did not consider the displacements.Their expression for the scaled shear strain is equal to the one presented in this paper, but their expressions for the scaled normal strains are different.They included an extra term πδ Ω for G xx and incorrectly used G yy = G xx .This latter equality is only valid for points outside the inclusion, but not for points inside the inclusion.It violates the condition for the volumetric strain (Equation 37) and is therefore incorrect.Note that the volumetric strain is constant inside the inclusion and equal to zero outside.This means that the inclusion experiences dilation, while the surroundings do not.The components of the strain tensor (i.e., the horizontal, vertical, and shear strain) are generally not constant in space.
The equations presented in this paper do satisfy the requirement for the volumetric strain.Writing Equation 37 in terms of G xx and G yy yields For a rectangular inclusion, it is relatively straightforward to prove that this holds for our expressions by taking the sum of Equations 40 and 41 and using the identity atan(x) + atan (1/x) = π 2 sign(x).This is a bit more cumbersome to prove for triangular inclusions, but the condition also holds in this case (Appendix B).The scaled volumetric strain is shown in Figure 5, which shows that the condition for the volumetric strain is met for both rectangular and triangular inclusions.
Additionally, we consider the jump in strains that arises from the condition of continuous normal traction across the inclusion boundary.Figure 6 shows the jump in the scaled strains G ij across a straight inclusion boundary for various boundary orientations, where we compare the jumps as computed from the closed-form expressions for the strains presented in Equations 40-42 and 45-47 and the closed-form expressions of Wu et al. (2021) with the jump required to satisfy continuity of normal traction as computed from Equations 29-31.The strains obtained from the closed-form expressions presented in this study satisfy the requirement of continuous normal traction given by Equations 29-31.The closed-form expressions for the strains given by Jansen et al. ( 2019) are equivalent to the expressions presented in Equations 40-42 and 45-47, and therefore also satisfy the condition of continuous normal traction.However, the normal strains as computed from the closed-form expressions of Wu et al. (2021) do not satisfy this condition.
Both Jansen et al. (2019) and Wu et al. (2021) considered Green's functions for the strains by placing the derivative under the integral sign in Equations 21-23.Moreover, Jansen et al. (2019) did not take into account the singularity at (x, y) = (x′, y′) and solved the double integral incorrectly, although this resulted in the correct expressions for the strains.However, changing the order of integration while disregarding the singularity would have given a different, incorrect result.Wu et al. (2021) did take into account the singularity and correctly solved the double integral.However, as we just showed, they did not obtain correct results for the strains.Because of the singularity at (x, y) = (x′, y′), we cannot interchange differentiation and integration for points inside the inclusion.Therefore, the Green's functions for the strains are only valid for points outside the inclusion.To circumvent this problem, we must not change the order of differentiation and integration.The correct way, as demonstrated in this paper, is to first solve the double integral for the displacements and then take the derivative to obtain the strains.
Nevertheless, both Jansen et al. (2019) and Wu et al. (2021) presented correct expressions for the stresses and their key results are therefore not affected by the mistakes in the derivation.We can show that our expressions for the normal stresses are equivalent to theirs.Using the expression for the volumetric strain (Equation 37) in the expression for the normal stresses (Equations 24 and 25) and using λ = 2μν 1 2ν , the normal stresses can be written as Defining the scaling factor C = 2 μD = (1 2ν)αp 2π(1 ν) allows Equations 49 and 50 to be written as which is identical to the expressions of the normal stresses given by Jansen et al. (2019) and Wu et al. (2021).
Their expressions for the stresses are therefore correct, although their explanations for the factor 2πδ Ω were incorrect.Jansen et al. (2019) had correct expressions for the strains, but incorrectly assumed that the volumetric strain equals zero everywhere, which is only valid for points outside the inclusion.They added the factor 2πδ Ω to comply with the condition of continuous normal traction.Wu et al. (2021) obtained incorrect expressions for the strains, but obtained correct expressions for the stresses by incorrectly using two different values for the eigenstress σ * ij throughout the solution process.Here, we have shown that the factor 2πδ Ω , instead, is the (scaled) sum of the contributions of the volumetric strain and the eigenstress to the normal stresses.
The jumps in the strains as required to satisfy continuous normal traction as computed by Equations 29-31 are represented by the solid lines, the jumps as computed from the closed-form expressions from Equations 40-42 and 45-47, which are equivalent to the expressions of Jansen et al. (2019), are represented by the circular markers, and the jumps as computed from the closed-form expressions of Wu et al. (2021) are represented by the asterisk markers.
The condition for the volumetric strain (Equation 37) can also be used to derive conditions for the stress arching ratios or the mean stress.The stress arching ratios γ ij are defined by Mulders (2003) and Soltanzadeh and Hawkes (2008) as the ratios of stress change over the pore pressure change multiplied by Biot's coefficient: Using Equations 37 and 53, the sum of the normal stress arching ratios can be expressed as which is identical to the expression as given by Soltanzadeh and Hawkes (2008), who derived this result for a cylindrical inclusion under plane strain conditions.Wu et al. (2021) correctly stated that this holds for any inclusion under plane strain conditions, regardless of its shape, as shown by Equation 54.
For a uniform eigenstrain as considered in this paper, the double integral for the displacements (Equations 13 and 14) can also be expressed as a line integral by applying Green's theorem.In that case, we can apply Leibniz's integral rule to obtain Green's functions for the strains, as long as we do not consider points (x, y) located on the boundaries of the inclusion.When only the strains and stresses are of interest, the line integral approach would save a computational step, as then it is not necessary to first derive expressions for the displacements.Additionally, analytical solutions for the integral over a straight line have already been given for both plane strain problems (Nozaki et al., 2001) and three-dimensional problems (D'Urso, 2013; Holstein & Ketteridge, 1996;Kuvshinov, 2008;Pohánka, 1988).These solutions can be used to determine the displacements, strains, and stresses for any polygonal or polyhedral inclusion.
While we focused on poroelastic stresses in this paper, the same approach is also applicable to determine thermoelastic stresses by simply replacing αp by 2μ(1 + ν) 1 2ν cT in all expressions presented here, where c is the linear thermal expansion coefficient and T is the change in temperature.In fact, the expressions for thermoelastic stresses were obtained much earlier than those for poroelastic ones (Goodier, 1937).The combined effects of poroelasticity and thermoelasticity can also be considered in a similar manner, as was done by Segall and Fitzgerald (1998) for the case of geothermal energy production.The mathematical subtleties discussed in this article are therefore not just limited to poroelasticity, but are also relevant for problems of (poro-)thermoelasticity.

Example: Depletion-Induced Fault Slip
As an example of the use of inclusion theory to compute induced fault stresses, consider a two-block hydrocarbon reservoir with two bounding faults and one central fault as depicted in Figure 7.An x y coordinate system is chosen with its origin in the central fault at depth D 0 = 3,000 m.The reservoir is represented with two rectangular and four triangular elements with dimensions listed in Table 1 where w is element width, r and s are bottom and top element y coordinates, and θ is fault angle, measured anti-clockwise with respect to horizontal.The fault throw and reservoir height follow as t = 125 m and h = 250 m.The central fault is permeable and the reservoir experiences a spatially uniform drop in pore pressure.Table 2 lists the physical properties that are required to compute the resulting depletion-induced stresses with the aid of inclusion theory.
The initial pore pressure p 0 and the initial stresses σ 0 ij in the reservoir and the faults can be computed as indicated at the bottom of Table 2.The depletioninduced incremental stresses σ ij are obtained by adding the contributions of the six elements.The individual element contributions for the incremental normal stresses σ xx and σ yy follow from Equations 51 and 52 and from a  similar expression for the incremental shear stresses σ xy , although without subtraction of the factor 2π inside the reservoir: where C = (1 2ν)αp 2π(1 ν) .The scaled strains G ij can be obtained from Equations 40-42 for rectangular elements.For triangular elements, they can be obtained from Equations 45-47, but only for right triangles oriented as in Figure 4, that is, with the right angle at the bottom-right corner and with the extension of the hypotenuse intersecting the y axis in the origin.To allow for the use of these equations for more general situations that may have different locations of the right angle and/or a shifted intersection with the y axis, four different cases can be identified, as depicted in Figure 8.The corresponding expressions G k ij ,k = 1,…,4, can be obtained as where G ij (x, y, o, p, r, s) indicates the original Equations 45-47, and where {x, ȳ, ō, p, r, s} = while y 0 is the intersection with the y axis according to The factor F can be expressed in terms of (modified) Kronecker deltas as (Valid for reservoir, overburden and underburden).
which implies that for configurations 2 and 4 the scaled shear strain G xy changes sign.
Once initial and incremental stresses have been computed, combined stresses Σ ij follow as their sum: The combined normal stress Σ ⊥ and the combined shear stress Σ ‖ in an inclined fault are then obtained through a standard tensorial coordinate transformation (see, e.g., Barber, 2010) where the x ŷ coordinate system is rotated anti-clockwise over an angle θ, and where a positive fault shear stress Σ ‖ is defined to correspond to an upward movement of the rock to the right side of the fault in relation to the rock at the left side.
To describe slip we use Coulomb friction theory (Ohnaka, 2013;Scholz, 2019).This implies that slip-provoking conditions occur when where Σ sl is the slip stress, defined as Here μ is the friction coefficient, Σ′ ⊥ is the effective normal stress, p 0 is the initial pore pressure (variable over the height of the reservoir), p is the incremental pore pressure (uniform in the reservoir; negative for depletion), and α f is an effective fault stress coefficient which is sometimes given the same value as the poroelastic Biot coefficient but also often set equal to one (Scholz, 2019).Note that in a typical deep subsurface situation the effective normal stress Σ′ ⊥ is negative, corresponding to compression.The friction coefficient may be a function of position, time, temperature, slip, slip rate, and/or additional state variables that represent internal mechanisms influencing the friction properties (Ohnaka, 2013;Scholz, 2019).Here, we only consider the simplest possible formulation: constant friction with a static friction coefficient μ.The Kronecker delta δ (Ω∪Φ) ensures that the pore pressure is only added for regions inside the reservoir domain Ω and/or the fault segment Φ (Jansen et al., 2019).Here we assume that the incremental pressure is felt in those fault segments that are in direct contact with reservoir rock from either one or two sides.Note that, in theory, all the peaks in the shear stresses (two in the bounding faults and four in the central fault) go to plus or minus infinity for even an infinitesimally small amount of depletion.In reality, pore pressure diffusion and a realistic geometry (e.g., rounded corners) will result in finite stress peaks, while non-elastic rock deformation will also put an upper limit on the stresses.However, also for more realistic physical conditions the emergence of sharp stress peaks is a typical characteristic of induced stresses in reservoirs with displaced faults.Moreover, an important consequence is that elastic numerical simulations aimed at finding a threshold depletion below which there will be no fault slip are not very meaningful because they are mesh dependent; that is, an increase in the mesh density near the peaks will lead to a decrease in the threshold change in pore pressure without converging to a final minimum.Also note that the present approach only provides information about potential slip patches, that is, zones where the absolute value of the fault shear stress exceeds the slip stress.Once slip occurs, a redistribution of shear stress will occur, leading to an increase in the slip patch size.Moreover, the nature of the friction (e.g., slip weakening, velocity weakening or rate-and-state-type) will determine whether the slip will be seismic or aseismic.For simple fault geometries, a (semi-)analytical approach of induced slip is still feasible (see, e.g., Jansen & Meulenbroek, 2022), but for more realistic situations, in particular those involving fault rupture, numerical simulations will be the method of choice (see, e.g., Buijze et al., 2017Buijze et al., , 2019)).
The present example considers spatially uniform reservoir depletion.However, inclusion theory can also be used to compute induced stresses in reservoirs with impermeable faults that cause the reservoir blocks to have different pore pressures (Wu et al., 2021), or in reservoirs with other nonuniform pore pressure distributions in which the integrals in the formulation may either be approximated with numerical quadrature or by simply using a large number of elements, each with a different pore pressure;  (2000,2002)), and Cornelissen and Jansen (2023b) (who used the 2D expressions of the present paper which are also volume based).

Conclusions
We presented the derivation of closed-form expressions for displacements, strains, and stresses inside a poroelastic inclusion.This derivation differs from both Jansen et al. ( 2019) and Wu et al. (2021), who incorrectly applied Leibniz's integral rule and integrated Green's functions for the strains.Jansen et al. (2019) did not take note of the singularity and solved the double integral incorrectly, but managed to obtain correct expressions for the strains.Wu et al. (2021) correctly solved the double integral, but their solution for the strains does not satisfy the condition of continuous normal traction.As Leibniz's integral rule is not valid inside the inclusion, the Green's function for the strains used by Jansen et al. (2019) and Wu et al. (2021) is only valid outside the inclusion.To obtain the correct expressions for the strains, the double integral must first be solved for the displacements, after which we can take the derivative to obtain the strains.Jansen et al. (2019) and Wu et al. (2021) do present correct expressions for the stresses, which are equivalent to the expressions presented in this study.Hence, while their derivations contained errors, their key results remain valid.Furthermore, we have clarified that the term 2πCδ Ω , used in the aforementioned studies, comprises of the contributions of the eigenstress and the volumetric strain to the normal stresses.
Additionally, we derived two conditions for the strains which can be used to verify the analytical solutions.First, continuous normal traction across the inclusion boundary requires a discontinuity in the strains.Second, the volumetric strain is proportional to the incremental pore pressure.For a uniform pore pressure change as considered here, this means that the volumetric strain is constant inside the inclusion and zero outside.These conditions can be used to verify expressions for the strains, also for inclusions of different shapes than considered here.
It is clear that the integral over the region S′ δ \D′ δ is absolutely convergent, since f(u, v) is continuous and the domain S′ δ \D′ δ is compact.Therefore, we can compute the double integral by an iterated integral, which yields zero due to the symmetry in f (u, v).This can easily be shown by separating the integral into four parts Next, we consider the improper integral over the domain D′ δ .In order to establish convergence of the integral over D′ δ we will use the following proposition proven by Ghorpade and Limaye (2010).
Proposition (from Ghorpade and Limaye (2010)) Let D be a bounded subset of R 2 such that δD is of content zero and let f : D → R be a nonnegative unbounded continuous function.Suppose that D admits an expanding sequence D n of subsets of D such that f is bounded on D n for each n ∈ N. Then ∬ D f is unconditionally convergent if and only if there is α ∈ R such that ∬ D n f ≤ α for all n ∈ N, and in this case ∬ D f = lim n→∞ ∬ D n f .
To make use of this proposition, we split D′ δ in two parts using polar coordinates such that f(u, v) ≥ 0 in D′ δ,r and f(u, v) ≤ 0 in D′ δ,l .For D′ δ,r we choose the sequence D n as follows and we note that D n and f(u, v) satisfy the conditions for absolute convergence.Since f(u, v) is bounded on the domain D n , we can use an iterated integral to compute the integral over this subdomain, which yields The integral over domain D′ δ,r is then given by

Figure 1 .
Figure1.Two-dimensional view of the four imaginary steps involved in computing the displacements in and around an expanding rectangular inclusion in an infinite elastic solid using Eshelby's method.

Figure 2 .
Figure2.The Green's function g ij gives the displacement in the i-direction at point (x, y) due to a point force in the j-direction at point (x′, y′), which is a function of the distance R between both points.Each point (x′, y′) inside the inclusion contributes to the displacement at point (x, y).The total displacement at point (x, y) is then obtained by integrating the Green's function over the entire inclusion (in this example, a rectangular inclusion).

Figure 4 .
Figure 4. Geometry of a triangular inclusion with corner points (o, r), ( p, r),and ( p, s).The extension of the hypotenuse goes through the origin of the coordinate system (i.e., the hypotenuse lies along the line y = x tan θ).

Figure 5 .
Figure 5.The scaled volumetric strain (G xx + G yy )/π for a rectangular inclusion (left) and triangular inclusion (right).Inside the inclusion, G xx + G yy = 2π, while outside the inclusion G xx + G yy = 0.

Figure 6 .
Figure 6.Scaled jumps in the strains across the inclusion boundary as a function of the boundary orientation.Here, θ indicates the angle with respect to a horizontal line, such that θ = 0°corresponds to a horizontal boundary and θ = 90°corresponds to a vertical boundary.The jumps are defined as ΔG ij = G out

Figure 7 .
Figure 7. Reservoir with two fault blocks divided in six elements (not to scale).

Figure 8 .
Figure 8.The four possible configurations of a rectangular-triangular element in a horizontal reservoir.The corresponding expressions G k ij ,k = 1,…,4 are given by Equation 56.

Figure 9
Figure 9 depicts the combined shear stress Σ ‖ and the two slip thresholds Σ sl and Σ sl in the three faults.Values Σ sl < Σ ‖ correspond to a positive slip; this situation occurs in the central fault for values of y between slightly less than 50 m and slightly more than 50 m, resulting in a slip patch showing a local increase of the initial normalfault offset.Opposed, values Σ ‖ < Σ sl correspond to a negative slip tendency.This occurs in the left fault for y values between approximately 45 and 150 m.

Figure 9 .
Figure 9. Shear stress Σ ‖ and slip thresholds Σ sl and Σ sl in the three faults.The horizontal dotted black lines indicate the top and bottom boundaries of the reservoir blocks.

Figure A1 .
Figure A1.We consider a square region with side lengths 2δ around the singularity, termed S δ (left).This region is further split into a disk with radius δ around the singularity D δ and the domain S δ \D δ (right).

Journal of Geophysical Research: Solid Earth
The initial pore pressure, initial vertical stress, and initial horizontal stress can be computed as: p 0 Barbosa et al. (2022) Wees et al. (2019) (who used 3D point-based expressions fromOkada (1992)for the subdomains in the numerical integration),Smith et al. (2022)(who used 3D volume-based expressions fromKuvshinov (2008)which, unlike point-based expressions, don't suffer from near-fault anomalies),Barbosa et al. (2022)(who used 3D volume-based expressions from Nagy et al.