# Analytical Solution to Assess the Induced Seismicity Potential of Faults in Pressurized and Depleted Reservoirs

Haiqing Wu, Victor Vilarrasa, Silvia de Simone, Maarten Saaltink, Francesco Parisio

<figure><img src="image_1.png" /></figure>

## To cite this version:

Haiqing Wu, Victor Vilarrasa, Silvia de Simone, Maarten Saaltink, Francesco Parisio. Analytical Solution to Assess the Induced Seismicity Potential of Faults in Pressurized and Depleted Reservoirs. Journal of Geophysical Research: Solid Earth, 2021, 126 (1), pp.e2020JB020436. 10.1029/2020JB020436. insu-03116199

HAL Id: insu-03116199

https://insu.hal.science/insu-03116199

Submitted on 20 Jan 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

10.1029/2020JB020436

**Key Points:**

* We propose an analytical solution for the stress response to pore pressure changes in the reservoir crossed by a fault with an arbitrary dip

* Induced seismicity potential of impermeable faults is always larger than that of permeable faults under any initial and injection conditions

* Induced seismicity potential is overestimated and the injectivity underestimated if the poromechanical coupling is neglected

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

H. Wu,
haiqing.wu@upc.edu

**Citation:**

Wu, H., Vilarrasa, V., De Simone, S., Saaltink, M., & Parisio, F. (2021). Analytical solution to assess the induced seismicity potential of faults in pressurized and depleted reservoirs. *Journal of Geophysical Research: Solid Earth*, 126, e2020JB020436. https://doi.org/10.1029/2020JB020436

Received 19 JUN 2020

Accepted 4 DEC 2020

<figure><img src="image_4.png" /></figure>

# Analytical Solution to Assess the Induced Seismicity Potential of Faults in Pressurized and Depleted Reservoirs

Haiqing Wu¹,² ID, Victor Vilarrasa²,³,⁴ ID, Silvia De Simone⁵ ID, Maarten Saaltink¹,² ID, and Francesco Parisio⁶ ID

¹Department of Civil and Environmental Engineering (DECA), Universitat Politècnica de Catalunya (UPC), Barcelona, Spain, ²Associated Unit: Hydrogeology Group (UPC-CSIC), Barcelona, Spain, ³Institute of Environmental Assessment and Water Research (IDAEA), Spanish National Research Council (CSIC), Barcelona, Spain, ⁴Mediterranean Institute for Advanced Studies (IMEDEA), Spanish National Research Council (CSIC), Esporles, Spain, ⁵Univ Rennes, CNRS, Géosciences Rennes-UMR 6118, Rennes, France, ⁶Chair of Soil Mechanics and Foundation Engineering, Institute of Geotechnics, Technische Universität Bergakademie Freiberg, Freiberg, Germany

**Abstract** Displaced faults crossing the reservoir could significantly increase the induced earthquake frequency in geo-energy projects. Understanding and predicting the stress variation in such cases is essential to minimize the risk of induced seismicity. Here, we adopt the inclusion theory to develop an analytical solution for the stress response to pore pressure variations within the reservoir for both permeable and impermeable faults with offset ranging from zero to the reservoir thickness. By analyzing fault stability changes due to reservoir pressurization/depletion under different scenarios, we find that (1) the induced seismicity potential of impermeable faults is always larger than that of permeable faults under any initial and injection conditions—the maximum size of the fault undergoing failure is 3–5 times larger for impermeable than for permeable faults; (2) stress concentration at the corners results in the occurrence of reversed slip in normal faults with a normal faulting stress regime; (3) while fault offset has no impact on the slip potential for impermeable faults, the slip potential increases with the offset for permeable faults, which indicates that non-displaced permeable faults constitute a safer choice for site selection; (4) an impermeable fault would rupture at a lower deviatoric stress, and at a smaller pressure buildup than a permeable one; and (5) the induced seismicity potential is overestimated and the injectivity underestimated if the stress arching (i.e., the poromechanical coupling) is neglected. This analytical solution is a useful tool for site selection and for supporting decision making during the lifetime of geo-energy projects.

## 1. Introduction

Induced seismicity has become a widespread issue as a result of the proliferation of geo-energy projects (Foulger et al., 2018). On one hand, geothermal energy production and geologic carbon storage are essential technologies to reach zero or negative net carbon emissions. On the other hand, the increased energy demand is boosting other operations, such as seasonal natural gas storage, subsurface energy storage, and disposal of wastewater from conventional and non-conventional oil and gas production. Injecting or pumping fluids at depth—a widespread practice in geo-energy operations—alters the in-situ stress field and may lead to fault rupture and induced seismicity (Buijze et al., 2017; Ellsworth, 2013; Grigoli et al., 2018). In several cases, the authorities have decided to cancel projects believed to be associated with large induced earthquakes and a non-exhaustive list includes the Deep Heat Mining Project in Basel, Switzerland (Deichmann & Giardini, 2009; Haring et al., 2008; Terakawa et al., 2012), the Castor Natural Storage Project, Spain (Del Potro & Diez 2015; Juanes et al., 2017; Villaseñor et al. 2020), and the Enhanced Geothermal System (EGS) Project at Pohang, South Korea (Ellsworth et al., 2019; Grigoli et al., 2018; Lee et al., 2019). To reduce the risks of induced seismicity and safely promote sustainable energy development, tools to predict and subsequently mitigate induced seismicity should be developed.

Despite considerable advancements in understanding the triggering mechanisms of induced seismicity in recent years, forecast and mitigation of induced seismicity remains challenging and some fundamental questions remain open (Ellsworth et al., 2016; Lee et al., 2019). Potential triggering mechanisms include pore
pressure diffusion caused by single phase (Shapiro & Dinske, 2009; Simpson et al., 1988) and multiphase flows (Zbinden et al., 2017), poroelastic and thermally induced stress changes (Chang & Segall, 2016; De Simone et al., 2017; Langenbruch & Zoback, 2016), and strength weakening due to geochemical reactions (Rohmer et al., 2016; Vilarrasa et al., 2019). These mechanisms, separately or acting jointly, can lead to fault slip (Lehner, 2019; Lele et al., 2016; Orlic et al., 2013; Orlic & Wassing, 2013; Rutqvist et al., 2016; Van den Bogert, 2015; Van Wees et al., 2017) and nucleation of dynamic rupture (Buijze et al., 2017, 2019; Galis et al., 2017, 2019; Garagash & Germanovich, 2012) on different geological settings (Bourne & Oates, 2017; Haug et al., 2018), even at very large distances (Goebel et al., 2017).

Faults intersecting the injection/pumping formation undergo pore pressure and stress changes, affecting their stability. Pore pressure changes are controlled by the hydraulic properties of faults, which are highly variable, ranging from conductive faults to flow barriers (Caine et al., 1996). For example, low-permeable faults are present at the Snohvit CO₂ storage site, Norway (Chiaramonte et al., 2013; Hansen et al., 2013), at Pohang EGS project, South Korea (Ellsworth et al., 2019; Kim et al., 2018), and at many compartmentalized reservoirs (e.g., Castelletto et al., 2013), while permeable faults are found at the Groningen gas field, the Netherlands (Jansen et al., 2019; Van Wees et al., 2014) and the Corinth rift, Greece (Duverger et al., 2015; Geraud et al., 2006). Stress changes arise when the reservoir deformation is restricted, as in the case of closed or compartmentalized reservoirs, and they are governed by the poromechanical properties of the rock—the stiffer the rock, the larger the induced stress—and by the fault offset, which generates an additional stress concentration (Buijze et al., 2017; Galis et al., 2017, 2019). Such generated stress could lead to an increase in induced earthquake frequency, as was observed in the Groningen gas field (NAM, 2016; Van Wees et al., 2014, 2017).

Numerical simulations can account for great physical and geometrical complexity, but the computational cost often prevents systematic explorations of the parametric space. Analytical methods offer an alternative to obtain fast estimations, but require more stringent hypotheses and simplifications on the geometry and physics of the problem when compared with numerical methods. Interestingly, their drawback turns into an advantage when the perspective is changed and the goal becomes a quick and efficient parametric space analysis, ultimately highlighting the factors controlling the problem. For the problem of reservoir pressurization/depletion, Eshelby's inclusion theory (Eshelby, 1957) is at the heart of several analytical solutions describing displacement, strain, and stress fields in an infinite half-space with an elliptic inclusion (Segall, 1985, 1992; Segall & Fitzgerald, 1998). The theory was applied to study subsidence and induced seismicity (Segall, 1985, 1989; Segall et al., 1994), recognizing the influence of reservoir geometry and orientation (Soltanzadeh & Hawkes, 2008, 2009) and the importance of including the contribution of crack-tip resistance to fault strength (Wang et al., 2016). Existing analytical solutions either assume non-displaced faults (Segall, 1985, 1992; Segall & Fitzgerald, 1998; Soltanzadeh & Hawkes, 2008, 2009; Wang et al., 2016) or displaced but permeable faults (Jansen et al., 2019). No solution currently exists for low-permeable faults that cross the reservoir with an offset: the aim of our contribution is to fill this knowledge gap and analyze the difference in terms of fault stability between permeable and low-permeable faults crossing a pressurized/depleted reservoir for both non-displaced and displaced faults.

In this paper, we propose an analytical solution for stress variations in response to injection/pumping into a reservoir crossed by a fault that could be either permeable or impermeable with offset ranging from zero to the reservoir thickness. Note that by stress variations we refer to the total stress changes, whereas for the effective stress, we explicitly mention effective in our terminology. The structure of the paper is as follows. In Section 2, we introduce the conceptual problem, develop the analytical solution, and show its validation. In Section 3, we present the methods to assess fault stability and fault slip potential based on our solution for both permeable and impermeable faults. In Section 4, we illustrate the effect of fault permeability on fault stability and perform systematic parametric space analyses of fault offset, fault dip, initial stress state, and pressure change. Finally, we provide extended discussion of the results and its principal implications in Section 5.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Geological model and its schematic geometry. Pore pressure changes within the reservoir on one or both sides of the fault depending on whether the fault is impermeable or permeable, respectively. The caprock and bedrock are assumed impermeable and thus, no pressure changes occur outside the reservoir.</figcaption></figure>

## 2. Analytical Solution for Stress Changes Around a Fault Crossing a Pressurized/Depleted Reservoir

### 2.1. Problem Formulation and Assumptions

We evaluate the induced stress arising in a deep reservoir crossed by a displaced permeable or impermeable fault as a consequence of fluid injection or production. The reservoir is treated as an inclusion that is hydraulically disconnected from the overlying caprock and underlying bedrock (Figure 1). A displaced fault with an arbitrary dip angle $\theta$ crosses the whole reservoir and extends to the surrounding rock, dividing the entire domain into two parts: the left part is the hanging wall and the right part is the footwall for a normal fault (Jansen et al., 2019). A non-displaced fault is a particular case in which the fault offset is zero. The fault offset, $ht$, is defined as $b-a$ (see Figure 1), and fault geometry is parameterized by four corner points (P1, P2, P3, and P4). The height (thickness) and width of the faulted reservoir are $a+b$ and $c+d$, respectively, where the width can be assumed as infinite by imposing $c = d = \infty$. The reservoir length is assumed as infinite in the out-of-plane direction.

We adopt the solid mechanics sign convention of stress and strain, that is, negative normal components denote compression, and a positive shear stress is assumed to rotate the material element in the counterclockwise direction, which indicates that the left part of the fault moves downward relative to the right part. For pore pressure, a negative pore pressure change refers to production and a positive one to injection.

We apply the following hypotheses: (1) the reservoir, assumed horizontal, elastic, homogeneous, and isotropic, and its surroundings have identical stiffness but different permeability and porosity—the latter implies that flow and pressure changes take place exclusively within the reservoir; (2) two-dimensional (2D) plane strain conditions apply based on the assumption that the reservoir extends infinitely in the out-of-plane direction; (3) quasi-steady-state uniform pore pressure changes occur in the reservoir as a consequence of injection/production, that is, the transient effect of flow is neglected; and (4) reservoir depth is great enough so that the effect of the free surface can be neglected (Lehner, 2019).

### 2.2. Analytical Solution

According to the inclusion theory (Eshelby, 1957; Mura, 1987; Rudnicki, 2011), pore pressure changes induce stress variation $\sigma_{ij}$ in the reservoir as (see Appendix A for the full mathematical development)

$$ \sigma_{ij}(x, y) = C \left[ \iint_{\Omega} g_{ij}(x, y, \zeta, \xi) d\Omega - \pi \delta_{ij} \delta_{\Omega} \right] = C \left[ G_{ij}(x, y) - \pi \delta_{ij} \delta_{\Omega} \right], \quad (1) $$

$$ C = \frac{(1 - 2\nu)\alpha\Delta p}{2\pi(1 - \nu)}, \quad (2) $$

where $\Omega$ is the inclusion domain, $g_{ij}$ and $G_{ij}$ represent the Green's function for stress and its surface integral, respectively, $x$ and $y$ are the Cartesian coordinates, $\zeta$ and $\xi$ denote the coordinate values within the domain $\Omega$, $\alpha$ is the Biot's coefficient, $\nu$ is the Poisson's ratio, and $\Delta p$ is the pore pressure change. $\delta_{ij}$ is the Kronecker delta, which equals 1 if $i = j$ or 0 if $i \neq j$ and $\delta_{\Omega}$ is the modified Kronecker delta, which equals 1 if $(x, y) \in \Omega$ or 0 if $(x, y) \notin \Omega$.

Green's function $g_{ij}$ gives the magnitude of the stress in the $i$-th direction at point $(x, y)$ in response to a body force in the $j$-th direction applied at point $(\zeta, \xi$; see Equations A17-A19). To perform their integration over the inclusion domain as in Equation 1, we set the origin to coincide with the midpoint of the fault (Figure 1). The integration domain is different, whether we are in the case of permeable or impermeable
faults. In the former case, pore pressure changes within the reservoir on both sides of the fault. For the latter case, pore pressure changes only on the side of the fault where injection or depletion takes place. In the permeable case, the entire inclusion consists of two trapezoids (Figure 1), each of which can be divided into two subdomains to simplify the integration. Thus, we apply the superposition principle of integral to combine the solutions for rectangular and triangular domains, which returns (see Appendix B for the full derivation) as

$$ 
\begin{aligned}
G_{xx}(x,y) = -G_{yy}(x,y) = & a \tan \frac{y+b}{x+c} - a \tan \frac{y-a}{x+c} + a \tan \frac{y-b}{x-d} - a \tan \frac{y+a}{x-d} \\
& - \left[ f_1(x,y,-b) - f_1(x,y,a) + f_1(x,y,b) - f_1(x,y,-a) \right] \sin^2 \theta \\
& - \frac{\sin \theta \cos \theta}{2} \ln \left( \frac{f_2(x,y,-b)f_2(x,y,b)}{f_2(x,y,a)f_2(x,y,-a)} \right) - \pi \delta_{\Omega}
\end{aligned}
\quad (3) $$

$$ G_{xy}(x,y) = \left[ f_1(x,y,-b) - f_1(x,y,a) + f_1(x,y,b) - f_1(x,y,-a) \right] \sin \theta \cos \theta - \frac{\sin^2 \theta}{2} \ln \frac{f_2(x,y,-b)f_2(x,y,b)}{f_2(x,y,a)f_2(x,y,-a)} + \frac{1}{2} \ln \frac{f_3(x+c,y+b)f_3(x-d,y-b)}{f_3(x+c,y-a)f_3(x-d,y+a)}, \quad (4) $$

where the functions $f_1$, $f_2$, and $f_3$ are

$$ f_1(x, y, \hat{y}) = a \tan \frac{(x - \hat{y} \cot \theta) \cot \theta + (y - \hat{y})}{x - y \cot \theta}, \quad (5) $$

$$ f_2(x, y, \hat{y}) = (x - \hat{y} \cot \theta)^2 + (y - \hat{y})^2, \quad (6) $$

$$ f_3(x - \hat{x}, y - \hat{y}) = (x - \hat{x})^2 + (y - \hat{y})^2, \quad (7) $$

and where a, b, c, and d are the geometrical parameters shown in Figure 1. The last term ($-\pi\delta_{\Omega}$) in Equation 3 results from the solution of improper integral (Courant & John, 1989) because the Green's function for stress becomes unbounded for points (x, y) located in the inclusion domain. The corners of two trapezoidal domains are singularities for the solutions (Equations 3 and 4; see Appendix B). The vertical fault is a special case of inclined fault, which is obtained by setting $\theta = 90^\circ$ (Equations B21-B22). And faults with no offset are also a special case in which $a = b$ (Equations B23-B24).

When substituting Equations 3 and 4 into Equation 1, we obtain the x-y planar solution for describing the distribution of induced stress in the pressurized or depleted reservoir and its surrounding rock. Our solution is consistent to the one developed by Jansen et al. (2019), for the case of a horizontal infinite reservoir crossed by a permeable fault. However, our current solution is also valid for any arbitrary reservoir width, with the solution for the infinite reservoir being a special case, that is, $c = d = \infty$.

Fault stability and its likelihood of rupture depend on the distribution of the normal and tangential stress components along the fault plane. Thus, the above x-y planar solution along the fault plane needs to be transformed into the coordinate system placed on the fault and oriented along it. We apply the stress transformation with axis rotation (Equations C4 and C5) to derive the closed expressions for such induced stress along an arbitrary fault plane with dip angle $\theta$, which yields (see Appendix C for the full derivation)
$$
\begin{align}
\bar{\sigma}_n(y \cot \theta, y) = & -\frac{\pi}{2} \left[ \operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) + \operatorname{sgn}(y-b) - \operatorname{sgn}(y+a) \right] \sin^2 \theta \nonumber \\
& -\cos 2\theta \left( a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} + a \tan \frac{y-b}{y \cot \theta - d} \right) \nonumber \\
& -a \tan \frac{y+a}{y \cot \theta - d} - \pi \delta_{\Omega} + \frac{\sin 2\theta}{4} \ln \frac{f_4(y,b)}{f_4(y,a)} \nonumber \\
& -\frac{\sin 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b) f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a) f_3(y \cot \theta - d, y+a)} - \pi \delta_{\Omega}, \tag{8}
\end{align}
$$

$$
\begin{align}
\bar{\tau}(y \cot \theta, y) = & -\frac{\pi}{4} \left[ \operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) + \operatorname{sgn}(y-b) - \operatorname{sgn}(y+a) \right] \sin 2\theta \nonumber \\
& + \sin 2\theta \left( a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} + a \tan \frac{y-b}{y \cot \theta - d} \right) \nonumber \\
& - a \tan \frac{y+a}{y \cot \theta - d} - \pi \delta_{\Omega} - \frac{\sin^2 \theta}{2} \ln \frac{f_4(y,b)}{f_4(y,a)} \tag{9} \\
& - \frac{\cos 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b) f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a) f_3(y \cot \theta - d, y+a)}, \nonumber
\end{align}
$$

where $\sigma_n(y \cot \theta, y)$ and $\tau(y \cot \theta, y)$ are the induced normal and tangential stress components along the fault plane, respectively; they will be shortened to $\sigma_n$ and $\tau$ for convenience hereafter. The stress components with an overbar denote the dimensionless stress components, which are normalized by the scaling parameter $C$ (Equation 2). The sgn $(\cdot)$ is the sign function defined as $1$ if $(\cdot) > 0$, $0$ if $(\cdot) = 0$ or $-1$ if $(\cdot) < 0$, and function $f_4$ is defined as

$$f_4(y, \hat{y}) = (y + \hat{y})^2 (y - \hat{y})^2. \qquad (10)$$

The four corners, P1, P2, P3, and P4, on the fault plane are singularities of Equations 8 and 9 (Figure 1). With such a general solution, one can easily find the solutions for the special cases of vertical faults (Equations C8 and C9) and zero offset faults (Equations C10 and C11). In particular, $\sigma_n$ and $\tau$ just correspond to $\sigma_{xx}$ and $\sigma_{xy}$ for the case of vertical faults, respectively. In the above equations, the segment P1-P2 of the fault belongs to the inclusion, while the segment P3-P4 belongs to the surroundings (see Appendix C).

For an impermeable fault, we assume that the pore pressure change is restricted to the side of the fault where injection/production takes place, while pore pressure on the other side remains unaltered. Thus, the integration of the Green's function for stress only entails one part of the inclusion domain, that is, one trapezoidal domain. Considering that fluid is injected into the left-hand side of the domain, such integrations are

$$
\begin{align}
G_{xx}(x,y) = -G_{yy}(x,y) = & a \tan \frac{y+b}{x+c} \nonumber \\
& -a \tan \frac{y-a}{x+c} - \frac{\sin \theta \cos \theta}{2} \ln \left( \frac{f_2(x,y,-b)}{f_2(x,y,a)} \right) \nonumber \\
& -[f_1(x,y,-b) - f_1(x,y,a)] \sin^2 \theta - \pi \delta_{\Omega}, \tag{11}
\end{align}
$$

$$
\begin{align}
G_{xy}(x,y) = & [f_1(x,y,-b) - f_1(x,y,a)] \sin \theta \cos \theta \nonumber \\
& - \frac{\sin^2 \theta}{2} \ln \frac{f_2(x,y,-b)}{f_2(x,y,a)} + \frac{1}{2} \ln \frac{f_3(x+c,y+b)}{f_3(x+c,y-a)}, \tag{12}
\end{align}
$$

and the dimensionless induced normal and tangential stress components along the fault plane are
<table><caption>**Table 1**

Geometrical Parameters of the Reservoir and Rock Properties Adopted for the Validation Example</caption><thead><tr><th>Parameter</th><th>Physical meaning</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>a</td><td>Geometrical parameters (refer to Figure 1)</td><td>100</td><td>m</td></tr><tr><td>b</td><td></td><td>200</td><td>m</td></tr><tr><td>c</td><td></td><td>2,000</td><td>m</td></tr><tr><td>d</td><td></td><td>2,000</td><td>m</td></tr><tr><td>θ</td><td>Fault dip</td><td>60</td><td>°</td></tr><tr><td>μ</td><td>Shear modulus</td><td>6,500</td><td>MPa</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>0.15</td><td>-</td></tr><tr><td>α</td><td>Biot's coefficient</td><td>0.9</td><td>-</td></tr><tr><td>Δp</td><td>Pressure buildup</td><td>20</td><td>MPa</td></tr></tbody></table>

$$ \bar{\sigma}_n(y \cot \theta, y) = -\cos 2\theta \left( a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} - \pi \delta_{\Omega} \right) \\ - \frac{\pi}{2} \left[ \operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) \right] \sin^2 \theta \\ + \frac{\sin 2\theta}{4} \ln \frac{(y+b)^2}{(y-a)^2} \\ - \frac{\sin 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b)}{f_3(y \cot \theta + c, y-a)} - \pi \delta_{\Omega}, \qquad (13) $$

$$ \bar{\tau}(y \cot \theta, y) = \sin 2\theta \left( a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} - \pi \delta_{\Omega} \right) \\ - \frac{\sin^2 \theta}{2} \ln \frac{(y+b)^2}{(y-a)^2} - \frac{\pi}{4} \left[ \operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) \right] \sin 2\theta \qquad (14) \\ - \frac{\cos 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b)}{f_3(y \cot \theta + c, y-a)}. $$

The corners of the left-hand trapezoidal domains of the fault are singularities of Equations 11 and 12, and the corners P2 and P4 on the fault plane are singularities of Equations 13 and 14 (Figure 1). The whole impermeable fault belongs to the surroundings for fluid injection into the left-hand side of the fault because we apply the right limit, that is, the limit that the argument approaches the fault from its right-hand side, as the value of the fault plane. Note that in Equation 11 the extra term originating from improper integral is always $-\pi\delta_{\Omega}$ because it depends only on the integrand (the limit of integration at improper points only depends on the integrand) and not on the geometry. The solution is similar in the case of injection into the right-hand side of the fault. To avoid confusion or repetitions, in the following we will always consider the case of injection into the left-hand side as an example to represent the impermeable case.

## 2.3. Validation Against Numerical Solution

To verify the accuracy and correctness of our analytical solution, we compare fluid injection-induced stress distribution along a permeable and an impermeable fault against numerical solutions. The numerical simulations are performed with the fully coupled finite element code CODE_BRIGHT (Olivella et al., 1994, 1996). The geometry is shown in Figure 1. We adopt dimensions and rock properties as in Jansen et al. (2019) in order to also compare our results with theirs (Table 1). In the numerical simulations, we mimic the impermeable rock, that is, caprock and bedrock, and the impermeable fault by assigning low values of intrinsic permeability, that is, $10^{-18}$ m$^2$. We impose mechanical boundary conditions of zero normal displacement to the lateral and lower boundaries and an overburden of $-70$ MPa on the upper boundary, corresponding to a depth of 3.5 km. We assume that the initial stress state is isotropic. The magnitude of the initial stress and pressure is irrelevant because we are interested in the stress changes induced by pore pressure changes and both the hydraulic and mechanical processes are linear. We impose a pressure buildup of 20 MPa in the reservoir, which is the entire reservoir in the case of permeable fault while half of it is in the case of impermeable fault. We make sure that the size of the reservoir is large enough to minimize the boundary effects.

We compare the numerical and analytical results for a permeable and an impermeable fault (Figure 2). Further results are presented in Section S1. The numerically computed induced stress on the fault plane is almost identical to the analytical one. Small discrepancies near the corners are a consequence of the existence of singularities for the analytical solution, which leads to an infinite stress, and of the discrete nature of the numerical solution. We also consider the case of vertical permeable fault crossing a horizontal infinite reservoir in order to compare our results with those of Jansen et al. (2019). Comparisons are shown in Section S2 and they exhibit good agreement.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Comparison between analytically and numerically evaluated induced (a) horizontal, (b) vertical, and (c) x-y planar shear stress components along the fault plane. The results along the fault are projected on the vertical axes y. The legend is shown in (a), indicating that solid lines represent the analytical results (AN) and dotted lines represent the numerical results with CODE_BRIGHT (CB), for both a permeable (per) and an impermeable (imp) fault. A schematic of the reservoir geometry, with the four corners, is indicated by the gray background.</figcaption></figure>

### 3. Induced Seismicity Potential Assessment

#### 3.1. Coulomb Failure Stress and Coulomb Failure Stress Change

We assess the fault stability by means of the Coulomb Failure Stress (CFS; King, 1994)

$$CFS = |\tau + \tau^0| + \eta_{st}(\sigma_n' + \sigma_n'^0), \quad (15)$$

where $\eta_{st}$ is the static friction coefficient, superscript 0 represents the initial state, and superscript ' denotes effective stress. Here, the initial normal and tangential stress components on the fault plane are also calculated according to the stress transformation (Equations C4 and C5). Shear stress always drives the fault to slip regardless of whether it is positive or negative. The second term on the right-hand side of Equation 15 denotes the fault slip resistance and increases as the effective normal stress becomes more compressive. A positive value of CFS indicates that slip is activated in the direction of the shear stress along the fault, that is, a positive shear stress represents the normal slip and a negative shear stress represents the reversed slip for a normal fault (Jha & Juanes, 2014).

For the assessment of induced seismicity and to identify whether a portion of the fault becomes more or less stable, we use the CFS change ($\Delta CFS$) as

$$\Delta CFS = |\tau + \tau^0| - |\tau^0| + \eta_{st}\sigma_n', \quad (16)$$

A positive $\Delta CFS$ implies that the induced stress is driving the fault toward failure and eventually co-seismic slip. CFS and $\Delta CFS$ can be normalized by the scaling parameter C (Equation 2) as
<table><caption>**Table 2**

Properties of Berea Sandstone and the Initial Stress State of the Reservoir Adopted for the Failure Potential Analysis</caption><thead><tr><th>Parameter</th><th>Physical meaning</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>θ</td><td>Fault dip</td><td>60</td><td>°</td></tr><tr><td>ht<sub>D</sub></td><td>Dimensionless fault offset</td><td>1/3</td><td>-</td></tr><tr><td>μ</td><td>Shear modulus</td><td>4,600</td><td>MPa</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>0.29</td><td>-</td></tr><tr><td>α</td><td>Biot's coefficient</td><td>0.7</td><td>-</td></tr><tr><td>Δp<sub>D</sub></td><td>Dimensionless pressure buildup</td><td>4/7</td><td>-</td></tr><tr><td>p<sup>0</sup></td><td>Initial pore pressure</td><td>35</td><td>MPa</td></tr><tr><td>σ<sup>0</sup><sub>xy</sub></td><td>Initial vertical stress</td><td>-70</td><td>MPa</td></tr><tr><td>k<sub>0</sub></td><td>Stress ratio of horizontal to vertical stress</td><td>0.6</td><td>-</td></tr><tr><td>σ<sup>0</sup><sub>xy</sub></td><td>Initial shear stress in the x-y plane</td><td>0</td><td>MPa</td></tr><tr><td>η<sub>st</sub></td><td>Static friction coefficient</td><td>0.6</td><td>-</td></tr><tr><td>C</td><td>Scaling parameter for stress (Equation 2)</td><td>1.318</td><td>MPa</td></tr><tr><td>CFS<sup>0</sup><sub>D</sub></td><td>Initial dimensionless CFS (Equation 17)</td><td>-1.954</td><td>-</td></tr></tbody></table>

$$CFS_D = \left| \bar{\tau} + \frac{\tau^0}{C} \right| + \eta_{st} \left( \bar{\sigma}_n + \frac{\alpha \Delta p}{C} + \frac{\sigma_n^0 + \alpha p^0}{C} \right), \quad (17)$$

$$\Delta CFS_D = \left| \frac{\bar{\tau} + \tau^0}{C} \right| - \left| \frac{\tau^0}{C} \right| + \eta_{st} \left( \bar{\sigma}_n + \frac{\alpha \Delta p}{C} \right), \quad (18)$$

where we make use of Equations 8 and 9 for impermeable faults or Equations 13 and 14 for impermeable faults, and the subscript D denotes dimensionless variables. Here the dimensionless effective pore pressure change, that is, the term $\alpha \Delta p / C$, is an initial physical property within the reservoir, which is independent of the process of integration. Thus, the segment P1-P4 undergoes the same pore pressure change as the reservoir for permeable faults, and the segment P2-P4 undergoes the same pore pressure change as the left-reservoir compartment for fluid injection into the left-hand side of the impermeable fault.

## 3.2. Fault Slip Size

To quantitatively evaluate the fault slip potential, and thus, the induced seismicity potential, we define the fault slip size as

$$S_i = \frac{\ell_i}{\sin \theta}, \quad (19)$$

where $\ell_i$ is a continuous interval in coordinate y with CFS > 0. The slipping area can be discontinuous, so more than one $S_i$ may exist. We assume that the greatest magnitude of induced earthquakes is proportional to the maximum fault slip size, defined as

$$S_{\max} = \max(S_i), \quad (20)$$

which can be expressed in dimensionless form as

$$S_{D \max} = \frac{S_{\max}}{L_s} = \max \left( \frac{\ell_i}{a+b} \right), \quad (21)$$

where $L_s$ is a characteristic length of fault, here assumed as the length of the fault intercepting the reservoir

$$L_s = \frac{a+b}{\sin \theta}, \quad (22)$$

Assuming that each grain is restricted by its surrounding grains, that is, the existence of cohesion between grains, the fault will not slide until the maximum unstable patch reaches a threshold. We set $S_{Dmax} = 0.01$ as the threshold of fault slip in this paper, that is, the fault is always regarded as stable for $S_{Dmax} < 0.01$.

## 3.3. Properties of the Base Case Scenario

We evaluate the stress variation and fault stability as well as the fault slip potential for a pressurized reservoir whose properties are derived from laboratory measurements on Berea sandstone (Makhnenko & Labuzet 2015; Vilarrasa et al., 2016; Table 2). We assume the same geometrical model as in Section 2.3 (Figure 1 and Table 1), with the center of the reservoir at 3.5 km depth, and the initial stress state (normal faulting stress regime) as shown in Table 2.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Dimensionless (a) shear and (b) normal stress components on the fault plane for the case of permeable (blue lines) and impermeable (red lines) fault. Dashed lines are the term of pressure buildup normalized by the scaling parameter C. Blue auxiliary line and label denote the zero point of $\bar{\tau}$ for a permeable fault, and the red ones for an impermeable fault. The results along the fault are projected on the vertical dimensionless axes $y_D$. A schematic of the reservoir geometry, with the four corners, is indicated by the gray background.</figcaption></figure>

To generalize the problem, we normalize the coordinate y and the fault offset by the reservoir thickness, and we scale the pressure buildup by the initial pore pressure, which yields the dimensionless variables

$$y_D = \frac{y}{(a + b)}, \qquad (23)$$

$$ht_D = \frac{(b - a)}{(a + b)}, \qquad (24)$$

$$\Delta p_D = \frac{\Delta p}{p^0}. \qquad (25)$$

We compare the results for the two scenarios of permeable and impermeable faults to understand the influence of the hydraulic properties of faults on the fault stability and fault slip potential. In addition to these base case scenarios, we perform a parametric space analyses to explore the effects of fault geometry, initial stress state, and operational aspects. We compare in all cases the difference between permeable and impermeable faults.

# 4. Results

## 4.1. Effect of Fault Permeability in the Base Case Scenario

We evaluate the dimensionless induced shear and normal stress components on the fault plane according to Equations 8 and 9 for a permeable fault as well as Equations 13 and 14 for an impermeable fault (Figure 3). For the permeable fault, the dimensionless induced stress is symmetrical with respect to $y_D = 0$, as it is reflected by Equations 8 and 9 when the geometric parameters are $c = d$. The corner points are singular, such that the induced shear stress tends to $+\infty$ at P1 and P4, and $-\infty$ at P2 and P3 (for representation purposes, the infinite shear stress is cut off to a finite value). The induced normal stress has a reverse behavior with respect to the induced shear stress (compare Figures 3a and 3b), that is, it tends to $-\infty$ at P1 and P4, and $+\infty$ at P2 and P3. The entire fault plane except for a tiny vicinity at corners P2 and P3 shows a negative induced normal stress implying an increase in slip resistance, which contributes to the fault stability.

Unlike the permeable case, the induced shear and normal stress components are not symmetrical with respect to $y_D = 0$ for the impermeable fault, but the feature of reverse behavior for shear and normal stress components still holds (compare Figures 3a and 3b). The induced shear and normal stress components tend to infinity at corners P2 and P4 as a consequence of injecting from the left-hand side. These stress singularities in both the permeable and impermeable faults correspond to the points of stress concentration. The infinite value is a theoretical consequence of the integration of the Green's function and it is unrealistic for faults in nature where the material will undergo nonlinear deformation bounding stress values. We also plot the initial shear and effective normal stress components as well as the pressure buildup in Figure 3 to identify the contribution of each term to fault slip.

The dimensionless $\Delta CFS$ ($\Delta CFS_D$, Equation (18) along the fault plane, reflecting variations in the fault stability, remains symmetrical with respect to $y_D = 0$ for the permeable fault (Figure 4a) because the arithmetic operations of the symmetrical stress does not alter its symmetry. The stability of the permeable fault decreases everywhere, except for a small region close to the internal corners P2 and P3. Around the external corners, $\Delta CFS_D$ reaches its maximum value because of the stress concentration, which will likely induce fault slip locally. Conversely, the impermeable fault (Figure 4a) becomes more stable above the internal corner P2 and less stable below it. To determine the actual fault stability and assess whether failure conditions occur, $CFS_D$ is computed as $CFS_D^0$ (Table 2) plus $\Delta CFS_D$. It results that $CFS_D$ has the same trend as
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Dimensionless (a) Coulomb Failure Stress Change ($\Delta CFS_D$) and (b) Coulomb Failure Stress ($CFS_D$) on the fault plane for the case of permeable (blue lines) and impermeable (red lines) faults. $I_i$ and $y_{D_i}$ are the subinterval and zero point of $CFS_D$, respectively, in which the blue auxiliary lines and labels denote the ones for a permeable fault, and the red ones for an impermeable fault. The results along the fault are projected on the vertical dimensionless axes $y_D$. A schematic of the reservoir geometry, with the four corners, is indicated by the gray background. We also show the case of no stress arching (denoted by NO SA), that is, stress changes induced by poromechanical effects are neglected. The legend is shown in (b).</figcaption></figure>

$\Delta CFS_D$ (compare Figures 4a and 4b), but shifted by the magnitude of $CFS_D^0$. Therefore, the size of the fault that potentially undergoes failure ($CFS_D > 0$) is smaller than that where stability decreases ($\Delta CFS_D > 0$).

To analyze the slip mechanism of permeable faults, given the symmetry of $CFS_D$, we divide the upper half part of the $CFS_D$ curve into six subintervals (I1, I2, I3, I4, I5, and I6) by the four zero-points of $CFS_D$ ($y_{D1}$, $y_{D3}$, $y_{D4}$, and $y_{D5}$) and the zero point of $\bar{\tau}$ ($y_{D2}$, see Figure 3a; blue symbols in Figures 3a and 4b). For impermeable faults, the whole $CFS_D$ curve, however, is divided into six subintervals (I1, I2, I3, I4, I5, and I6) by the four zero points of $CFS_D$ ($y_{D1}$, $y_{D2}$, $y_{D3}$, and $y_{D5}$) and the zero point of $\bar{\tau}$ ($y_{D4}$, see Figure 3a; red symbols in Figures 3a and 4b). Detailed fault state and mechanics for each subinterval are shown in Table 3, in which both the shear and normal stress components changes belong to the fluid injection-induced poroelastic response. Overall, pore pressure buildup, which mainly results in the decrease of slip resistance, induces fault slip in the reservoir or makes it less stable. The poroelastic response, however, represents a stabilizing effect on the fault within the reservoir except for a small vicinity around the corners because of the local stress concentration. While the poroelastic response has a small negative effect on fault stability both in the caprock and bedrock for permeable faults, it performs a positive effect on fault stability in the caprock but a negative effect in the bedrock for impermeable faults (Figure 4a).

Following such combined characteristics of pore pressure buildup and poroelastic response, a permeable fault has four disconnected unstable patches (two normal slip patches and two reversed slip patches, Figure 4b), and an impermeable one has two unstable patches (one normal slip patch and one reversed slip patch, Figure 4b). The unstable patches of permeable faults, located between the external and internal corners (i.e., between P1 and P2 and between P3 and P4), are symmetric with respect to $y_D = 0$ and are separated by the stable central portion of the reservoir (between the internal corners P2 and P3). In contrast, practically the whole section of impermeable faults located in the pressurized reservoir become unstable.
<table><caption>**Table 3**
Slip Mechanism of Permeable and Impermeable Faults in the Base Case Scenario</caption><thead><tr><th>Fault</th><th>Subinterval</th><th>State (see Figure 4b)</th><th>Mechanism (see Figure 3)</th></tr></thead><tbody><tr><td rowspan="6">Permeable</td><td>I1</td><td>Less stable</td><td>Increase in shear stress</td></tr><tr><td>I2</td><td>Normal slip</td><td>Increase both in shear stress and pore pressure</td></tr><tr><td>I3</td><td>Normal slip</td><td>Increase in pore pressure</td></tr><tr><td>I4</td><td>Stable</td><td>Decrease in shear stress</td></tr><tr><td>I5</td><td>Reversed slip</td><td>Reversed increase in shear stress</td></tr><tr><td>I6</td><td>Stable</td><td>Decrease in shear stress and increase in normal stress</td></tr><tr><td rowspan="6">Impermeable</td><td>I1</td><td>More stable</td><td>Decrease in shear stress</td></tr><tr><td>I2</td><td>Reversed slip</td><td>Reversed increase in shear stress</td></tr><tr><td>I3</td><td>Stable</td><td>Decrease in shear stress</td></tr><tr><td>I4</td><td>Normal slip</td><td>Increase in pore pressure</td></tr><tr><td>I5</td><td>Normal slip</td><td>Increase both in shear stress and pore pressure</td></tr><tr><td>I6</td><td>Less stable</td><td>Increase in shear stress</td></tr></tbody></table>

In this case, while $S_{Dmax}$ = 0.32 for the permeable fault, it reaches 0.99 for the impermeable case, so both of them slide but the slip size of the impermeable fault is more than three times that of the permeable fault. Multiplying $S_{Dmax}$ by the fault characteristic length (Equation 22) yields the dimensional maximum fault slip size, which is $S_{max}$ = 110.85 m for the permeable case and $S_{max}$ = 342.95 m for the impermeable case.

For illustrative purposes, we include in Figure 4 the case in which the stress arching (Rudnicki 2002; Segall, 1985; Soltanzadeh & Hawkes, 2008) is neglected, that is, stress changes both inside and outside the reservoir induced by poromechanical effects are neglected and the effective normal stress variation equals the pressure changes. Thus, only the pressure buildup in the reservoir induces the increase in CFS, that is, $\Delta CFS = \eta_{st}\alpha\Delta p$ = 8.4 MPa in the base case scenario, which is significantly larger than the $\Delta CFS$ for the case of including the stress arching, except for the infinite values at the corner points. Neglecting stress changes significantly overestimates the decrease in fault stability, because the compression induced in the rock in response to reservoir expansion caused by pressurization is not taken into account. We will discuss this further in Section 5.

## 4.2. Effect of Fault Offset and Fault Dip

Fault offset affects differently permeable and impermeable faults (Figure 5). While fault stability significantly varies with offset for permeable faults (Figure 5a), impermeable faults undergo the same stability changes, but shifted, coinciding with the center of the pressurized/depleted reservoir (Figure 5b). For a permeable fault, $\Delta CFS_D$ slightly increases (stability decreases) within the reservoir and it is barely altered in the surrounding rock when the offset is equal to zero (Figure 5a). The stability-decreasing section increases proportionally to the fault offset and it concentrates at the external corners, where it tends to infinity. The section of the fault where the reservoir is juxtaposed on both sides of the fault presents a slight increase in stability. The size of this stabilized section decreases with fault offset, becoming negligible when the fault offset equals the reservoir thickness. In contrast, the size of the symmetric stability-decreasing sections between the internal and external corners of the reservoir increases with fault offset. Furthermore, both the stability of the caprock and the bedrock also slightly decreases. For an impermeable fault, the size of the stability-decreasing section, which is mainly constrained by the reservoir thickness, is independent of fault offset (Figure 5b), because the effect of stress concentration, which is controlled by the horizontal boundaries of the reservoir and the fault plane, is always the same for the impermeable fault regardless of its offset.

We analyze the impact of the fault dip, $\theta$, (between 0°, horizontal fault, and 90°, vertical fault) on fault stability for both permeable and impermeable faults with no offset (Figure 6) and 1/3 of dimensionless offset (Figure 7). Figure 6a displays the schematic geometric model of no offset fault and the initial $CFS_D$ as a
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Dimensionless Coulomb Failure Stress Change ($\Delta CFS_D$) on the fault plane for several dimensionless fault offsets for (a) permeable and (b) impermeable faults. The numbers on the curves denote the dimensionless fault offset. The results along the fault are projected on the vertical dimensionless axes $y_D$. The gray background indicates the position of the hanging and foot walls, which move simultaneously as the offset increases.</figcaption></figure>

function of the dip angle. The initial $CFS_D$ shows that the fault is stable, with the most critical dip around 61°, as expected for a normal faulting stress regime with a fault friction coefficient of 0.6. As a result of reservoir pressurization, fault stability changes differ depending on the hydraulic nature of the fault. For a permeable fault, $\Delta CFS_D$ is constant in the reservoir for each value of the dip angle (Figure 6b). It exhibits the maximum value for $\theta = 0°$ (horizontal fault) and the minimum value for $\theta \approx 61°$, while it is close to zero in the surrounding rock for any value of the dip angle. For an impermeable fault, $\Delta CFS_D$ is not constant along the fault plane (see also Figure 5b), although its variation is only strongly relevant for $\theta > 45°$ (Figure 6c). Its maximum value is always located at the external corners (the horizontal boundary between the reservoir and its surrounding rock) because of the stress concentration (recall Figure 4a), especially for a high dip angle because the effect of stress concentration becomes maximum for the angle of corners at 90° (Ahmadi et al., 2012). $\Delta CFS_D$ in the caprock decreases for increasing dip angle except for $\theta > 85°$, being negative (more stable) in the range of $31° < \theta < 85°$. $\Delta CFS_D$ increases with the dip angle in the bedrock and is positive (less stable) for $\theta > 31°$. Note that to assess fault stability, $\Delta CFS_D$ for either a permeable or an impermeable fault has to be added to the initial $CFS_D$, which also changes with the dip angle (Figure 6a). The difference in $\Delta CFS_D$ between the permeable and impermeable faults is negative throughout the reservoir and baserock for almost all the dip angles, that is, the impermeable fault is less stable, and positive over a small area located in the overlying caprock as a consequence of the left-hand side fluid injection (Figure 6d). A right-hand side injection would yield symmetrical results, with a positive difference located in the bedrock. Overall, an impermeable fault is more likely to be reactive than a permeable one when there is no offset.

As for the case of 1/3 of dimensionless offset (Figure 7), $CFS_D^0$ is the same as for the fault with no offset (compare Figures 6a and 7a). $\Delta CFS_D$ for a permeable fault is not constant along its plane and is symmetrical with respect to $y_D = 0$ (Figure 7b). The fault has greater stability in the section between the internal corners P2 and P3, where the reservoir is juxtaposed on both sides of the fault, for $50° < \theta < 80°$ and lower stability for other dips. For an impermeable fault (Figure 7c), $\Delta CFS_D$ has the same distribution as for the case of
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Fault stability changes along the fault plane as a function of the dip angle, $\theta$, for the case of zero-offset fault. (a) Schematic geometry and $CFS_D^0$ (shown in a polar coordinate system, the polar angle and diameter denote the fault dip angle and $CFS_D^0$, respectively), (b) $\Delta CFS_D$ for the permeable fault, (c) $\Delta CFS_D$ for the impermeable fault (for comparison purposes, (b) and (c) have the same color scale but the range of the legend is adapted to the values shown in each case), and (d) the difference between $\Delta CFS_D$ for the permeable and impermeable faults (i.e., (b-c)), where negative values indicate that the impermeable fault is less stable.</figcaption></figure>

no offset (Figure 6c), with the discontinuity shifted downward as a consequence of the downward shift of the boundary between the reservoir and the surrounding rock (recall Figure 5b). The difference in $\Delta CFS_D$ between the permeable and impermeable faults shows that for the section above the internal corner P2, the permeable fault is more unstable than the impermeable one, but the impermeable fault is less stable in the rest—a similar result to the zero-offset case (Figure 7d).

We present the fault slip potential as a function of the fault dip for several fault offsets while keeping the other parameters as the base case scenario, expressed in terms of the dimensionless maximum fault slip size ($S_{Dmax}$), in Figure 8. For permeable faults, there is a clear onset value for fault dip ($\theta_o = 42^\circ$) corresponding to the threshold of fault slip, that is, the fault undergoes slip for dip angles above the onset dip. Once $\theta > \theta_o$, $S_{Dmax}$ increases rapidly and then gradually reaches its peak around $\theta_c = 56^\circ$, that is, close to the critical dip angle for a normal faulting stress regime. The general trend of $S_{Dmax}$ as well as the onset and critical fault dips for slip are almost independent of fault offset, but its maximum value ($\theta = \theta_c$) increases with fault offset. According to our initial and injection conditions (Table 2), rupture does not occur when the offset is zero. However, for more critical initial stress state or larger pressure buildup, failure would occur also in the case of zero offset as a consequence of pore pressure buildup (recall Figure 3). When the fault is close to be vertically oriented, there is an inflection point in $S_{Dmax}$ that is barely visible for a small-offset (1/3) fault, but becomes more evident with larger offset. The occurrence of this inflection reflects that $S_{Dmax}$ for
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Fault stability changes along the fault plane as a function of the dip angle, $\theta$, for the case with fault dimensionless offset equal to 1/3. (a) Schematic geometry and $CFS_D^0$ (shown in a polar coordinate system, the polar angle and diameter denote the fault dip angle and $CFS_D^0$, respectively), (b) $\Delta CFS_D$ for the permeable fault, (c) $\Delta CFS_D$ for the impermeable fault (for comparison purposes, (b) and (c) have the same color scale but the range of the legend is adapted to the values shown in each case), and (d) the difference between $\Delta CFS_D$ for the permeable and impermeable faults (i.e., (b)–(c)), where negative values indicate that the impermeable fault is less stable.</figcaption></figure>

reversed slip exceeds the one for normal slip, that is, the reversed slip becomes the primary slip form. For impermeable faults, the effect of fault dip on $S_{Dmax}$ is similar to the case of permeable fault, with an onset dip $\theta_0 = 41.3^\circ$ and a critical dip $\theta_c = 59.4^\circ$, but no inflection is observed and the maximum value of $S_{Dmax}$ is larger, approaching the reservoir thickness at the critical dip. The fault slip potential in impermeable faults is independent of the fault offset also because the effect of stress concentration is always the same, regardless of its offset—a result similar to the fault stability (recall Figure 5b).

### 4.3. Effect of Initial Stress and Pore Pressure Changes

The initial stress state determines the initial $CFS$, which significantly affects the fault stability and fault slip potential. We explore the influence of the initial stress state by applying different values of the horizontal to vertical stress ratio ($k_0$) while keeping the vertical stress constant (Figure 9). Given that we adopt a value of the static friction coefficient of 0.6, we set a minimum stress ratio of 0.563 to ensure that the initial conditions correspond to $CFS < 0$ (stable fault). We define a critical stress ratio $k_0^c$ (marked by dots in Figure 9) which corresponds to the threshold of fault slip ($S_{Dmax} = 0.01$), that is, $k_0 > k_0^c$ implies a stable fault (negligible rupture size), and $k_0 < k_0^c$ a ruptured one for the applied pore pressure change. The dimensionless maximum fault slip size increases with decreasing stress ratio for both the permeable and impermeable
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Dimensionless maximum fault slip size ($S_{Dmax}$) as a function of fault dip for both permeable (blue color-scale lines) and impermeable faults (red line) and for several values of the dimensionless offset. The number on the blue lines denotes the dimensionless offset of permeable fault, and the results of different offsets for impermeable fault coincide in one line.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Dimensionless maximum fault slip size ($S_{Dmax}$) as a function of the initial stress ratio for several pressure buildups. Dots represent the critical stress ratios. Blue and red color scales for curves and dots correspond to permeable and impermeable faults, respectively, and the number on the curves denotes the dimensionless pressure buildup.</figcaption></figure>

faults (Figure 9). In other words, the smaller the stress ratio, the larger the deviatoric stress and thus, the larger the fault slip size. The rate of increase in $S_{Dmax}$ with decreasing $k_0$ is not steady though and is controlled by the cusp-like shape of $CFS_D$ (recall Figure 4b). The sharp increase in $S_{Dmax}$ once $k_0$ becomes lower than $k_0^c$ is due to the progressive failure of the pressurized reservoir. Once the whole reservoir is in failure, the portion of the caprock or bedrock that undergoes failure increases slowly as $k_0$ decreases (see Figure 4b). As the initial CFS approaches 0, the rupture size sharply increases because the asymptotic increase in $\Delta CFS_D$ within the caprock or bedrock is reached. The maximum rupture size coincides with the minimum in $k_0$, and for an impermeable fault in the base case scenario (only changing the stress ratio to its minimum) the rupture size is three times greater than the reservoir thickness ($S_{Dmax} \approx 3$). The rupture size is a lower bound because our model does not incorporate frictional strength weakening (Buijze et al., 2017, 2019; Garagash & Germanovich, 2012), and stress redistribution (De Simone et al., 2017; Sacks et al., 1978) associated with shear slip activation.

The effect of operational aspects, expressed as pore pressure changes, is mainly controlled by the injected volume, injection rate, and reservoir boundaries, that is, compartmentalization (Mathias et al., 2009; Nordbotten et al., 2005; Wu et al., 2016, 2018), and it affects the magnitude of induced earthquakes. Therefore, we further explore the ($S_{Dmax}$, $k_0$) space for different values of pressure buildup (Figures 9 and 10). We find that both $S_{Dmax}$ and $k_0^c$ (corresponding to the contour of $S_{Dmax} = 0.01$ in Figure 10) linearly increase with pressure buildup, with the highest increments of $S_{Dmax}$ corresponding to the impermeable fault. The contour plots in $S_{Dmax}$ with $k_0$ and $\Delta p_D$ show that the rupture size and the critical stress ratio for an impermeable fault are larger than for a permeable one under any initial and injection conditions (Figure 10). Thus, impermeable faults would rupture at lower initial deviatoric stress and with larger earthquake magnitude. Generally, $S_{Dmax}$ for an impermeable fault is 3–5 times greater than for a permeable one under a given $k_0$ and $\Delta p_D$.

The previous analyses on fault dip show that the onset and critical dip angles ($\theta_o$ and $\theta_c$) are barely related to the fault offset (Figure 8). Nonetheless, Equation 17 suggests that $\theta_o$ and $\theta_c$ depend upon the initial stress ratio $k_0$ and pressure buildup $\Delta p_D$. Such dependencies are the object of our next analyses. For a given pressure buildup (Table 2), while $\theta_o$ monotonically increases with the stress ratio, $\theta_c$ exhibits a more complex behavior (Figure 11a). $\theta_c$ evolution can be divided into three phases, with an increasing trend as $k_0$ increases, except for a decreasing branch in the mid-valued range of $k_0$. The differential $\theta_c - \theta_o$ decreases with increasing $k_0$ and equals 0, that is, the onset dip angle coincides with the critical dip angle, at $k_0 = 0.674$ (i.e., the critical stress ratio in this case) for a permeable fault and at $k_0 = 0.694$ for an impermeable fault. A direct consequence is that the range of dip angles favorable to slip is reduced for increasing $k_0$. For $k_0 > 0.674$ (0.694) the permeable (impermeable) fault is always stable regardless of its inclination (recall Figures 8 and 9). This means that geological sites with a higher in-situ stress ratio, that is, lower initial deviatoric stress, are intrinsically less prone to fluid injection-induced seismicity. The onset dip angle is smaller for impermeable faults than for permeable faults. Thus, the range of dip angles favorable to slip is larger for impermeable than for permeable faults, and its difference increases with $k_0$. The critical stress ratio corresponding to such critical fault dip
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Dimensionless maximum fault slip size ($S_{Dmax}$) as a function of the initial stress ratio and dimensionless pressure buildup for (a) permeable and (b) impermeable faults (for comparison purposes, both figures have the same color scale). The numbers on the contours denote the values of $S_{Dmax}$.</figcaption></figure>

is also greater for impermeable faults than for permeable ones, similar to the case for an arbitrary fault dip as shown in Figures 9 and 10. The difference in critical stress ratio between Figures 9 and 11 indicates that it increases with dip angle (only for $\theta_o < \theta < \theta_c$) and peaks at $\theta_c$. Thus, the site characteristics significantly affect the results of the induced seismicity assessment.

In contrast, for a given $k_0$, the onset dip angle $\theta_o$ monotonically decreases with increasing dimensionless pressure buildup $\Delta p_D$, whereas $\theta_c$ decreases first and then increases with increasing $\Delta p_D$ (Figure 11b). Thus, the differential $\theta_c - \theta_o$ decreases with decreasing $\Delta p_D$ and equals 0 at $\Delta p_D = 0.213$ (defined as the critical pressure buildup for this case) for a permeable fault and at $\Delta p_D = 0.186$ for an impermeable fault. Similarly, for $\Delta p_D < 0.213$ (0.186) the permeable (impermeable) fault is always stable regardless of its inclination. This means that a larger pressure buildup is necessary to induce seismicity, which can be translated into a value of maximum sustainable injection pressure (Bai et al., 2017; Rutqvist et al., 2007; Zhou et al., 2008) to

<figure><img src="image_15.png" /><figcaption>**Figure 11.** The onset and critical dip angles ($\theta_o$ and $\theta_c$) for fault slip as a function of (a) the initial stress ratio with a dimensionless pressure buildup of 4/7 and (b) the dimensionless pressure buildup with a stress ratio of 0.6. The blue and red lines correspond to the permeable and impermeable faults, respectively. The solid and dashed lines denote $\theta_o$ and $\theta_c$, respectively, while the dotted lines mean the difference between $\theta_c$ and $\theta_o$.</figcaption></figure>
minimize the risk of inducing seismicity. The main difference between permeable and impermeable faults is that the range of dip angles favorable to slip is larger for impermeable than for permeable faults, and the critical $Δp_D$ related to $θ_c – θ_o = 0$ is smaller for impermeable than for permeable faults. It implies that an impermeable fault is more likely to induce seismicity in terms of pressure buildup and its maximum sustainable pressure is smaller compared with a permeable fault.

## 5. Discussions

We present an analytical solution to assess fault stability changes ($ΔCFS_D$) as a result of reservoir pressurization/depletion. The results are shown in dimensionless form, which generalizes the problem with respect to the pore pressure change. All the stress components, $ΔCFS$ and CFS, are normalized by the scaling parameter C, which strongly depends on the Poisson's ratio ($ν$). C monotonically decreases with the increase of $ν$ (Equation 2). Fault stability is obtained by adding $ΔCFS_D$ to the initial $CFS_D$, which depends on the initial stress state and pore pressure (Equation 17). The variation in $ν$ will change the final $CFS_D$ and its zero points. An extra calculation shows that the magnitude of $CFS_D$ and $S_{Dmax}$ increases with increasing $ν$ for both permeable and impermeable faults. Furthermore, the zero-offset permeable fault is stabilized when $ν$ is lower than 0.24 for the case considered in Table 2, that is, fault dip of 60° and pore pressure buildup of 20 MPa (it is destabilized for $ν = 0.29$ in Figure 4a).

The analytical solution is a useful tool to quickly evaluate the induced seismicity potential of geo-energy projects. Injection control strategies of the maximum sustainable pressure, which have been validated in CO₂ sequestration projects (Bai et al., 2017; Rutqvist et al., 2007; Zhou et al., 2008), could benefit from adding this solution in the decision-making process. A detailed site characterization is needed for its effective application, because the maximum rupture size, and thus, the magnitude of the induced earthquakes, not only depends on operational aspects, like pressure change, but also on the stress state and hydro-geomechanical characteristics of the fault, such as permeability, strength, offset, and dip. For given initial and injection conditions, the range of dip angles that may undergo failure can be defined, but knowing the hydraulic properties of faults is critical because all the parametric space analyses confirm the unfavorable effect in terms of induced seismicity potential of an impermeable fault, and an identical result can be predicted for fluid injection into the right-hand side of the fault (i.e., the footwall). In particular, low-permeable faults lead to larger rupture area at lower initial deviatoric stress and at smaller pressure changes (Figures 8–11). One factor for causing such unfavorable effect is that the differential deformation between the two sides of the impermeable fault distributes in the whole reservoir compartment, that is, the segment P2–P4 for injection into the hanging wall or the segment P1–P3 for injection into the footwall, while the centered segment P2–P3 of a permeable fault always keeps the same deformation.

Pore pressure buildup arising from fluid injection reduces the slip resistance, and thus drives the fault toward failure (Shapiro & Dinske, 2009). The subsequent poroelastic response to adapt to such change tends to balance the system (Figure 3). However, the presence of the corners P1, P2, P3, and P4 (see Figure 1) results in a strong stress concentration (Buijze et al., 2017; Galis et al., 2019), leading to excessive adjustment nearby the corners and fault sliding. Therefore, pore pressure buildup and stress concentration at the corners during reservoir pressurization are the main reasons for fault reactivation and induced seismicity when ignoring thermal effects and geochemical reactions. Stress concentration is mainly controlled by the geometry of corners (Buijze et al., 2017), which changes with fault offset for permeable faults, but is independent of offset for impermeable faults, resulting in different fault stability patterns between permeable and impermeable faults (Figures 5–7).

Stress concentration not only highlights the effect of the hydraulic nature of faults, but also results in the occurrence of reversed slip at the internal corners P2 and P3 for normal faults with a normal faulting stress regime. In particular, reversed slip becomes the primary slip form for the permeable faults with a large dip angle and offset (Figure 8). In contrast, normal slip will occur and become the primary slip form for a small dip angle and permeable fault in thrust faults with a thrust faulting stress regime, that is, $k_0 > 1$.

Regarding the criterion for fault slip, one option is to focus on points, that is, when the stress state of one point reaches its failure conditions, for example, $CFS > 0$, the fault slips. The disadvantage of such a
criterion is that the infinite-induced stress at the corners leads to a misjudgment that even small amounts of injection will induce fault slip (Jansen et al., 2019). Another option consists in considering a minimum unstable patch, like the fault slip size defined in this paper. We have considered a threshold (0.01) for the dimensionless maximum fault slip size, that is, fault slip will not occur unless $S_{Dmax} > 0.01$. Although further research is needed to determine how much the threshold should be, the parametric space analyses and the related conclusions that can be drawn from Section 4 are independent of the threshold value adopted.

Concerning fluid production, the results of $\Delta CFS_D$ are basically symmetrical to those of injection presented in this paper and are easily obtained by a sign change, because the pressure change is included in the scaling parameter $C$. Thus, the change in fault stability is the opposite to the case of injection and, while fault slip first occurs at the external corners (P1 and P4) during reservoir pressurization, it first happens at the internal corners (P2 and P3) during reservoir depletion for both the permeable and impermeable faults (Jansen et al., 2019). The reverse effect may represent a potential method to control or mitigate induced seismicity: short-term production followed by an injection phase could improve the stability of the reservoir. Analogously, such an operation could be performed before decommissioning. Moreover, we also observe that the fault slip tendency in the case of production is much smaller than that in injection because $\Delta CFS_D$ is mainly negative (see Figure 4). Thus, the risk of induced seismicity in the case of injection is higher than that of production in a normal faulting stress regime like the one considered in this study.

Our analytical solution provides an accurate (Figure 2) and fast estimation of the stress variation in the reservoir and its surroundings, which takes into account the increase in stress due to the geological constraints to deformation. The stress variation, which is often called stress arching (Rudnicki 2002; Segall, 1985; Soltanzadeh & Hawkes, 2008), is positive within the reservoir in response to fluid injection. This means that the effective stress reduction is smaller than the pressure buildup (Figure 3). However, many engineering applications in subsurface energy-related projects neglect the stress arching in assessing fault stability and the effective stress reduction is assumed as equivalent to the pore pressure buildup during injection (e.g., Karvounis et al., 2014). For example, simplified models for CO₂ sequestration calculate the maximum sustainable pressure as the fracture pressure (Bandilla & Celia, 2017; Mathias et al., 2009; Zhou et al., 2008), neglecting stress arching. We compare the slip tendency analysis estimated by means of our solution with that estimated by neglecting the stress arching in Figure 4. The results show that ignoring the stress arching means overestimating the slip potential for both permeable and impermeable fault in a normal faulting stress regime. This implies an underestimation of the maximum sustainable pressure (i.e., injectivity).

## 6. Conclusion

In recent decades, increasing interest in the subsurface as a source of carbon-free energy resources has led to an increasing number of induced earthquakes, with some of these earthquakes resulting in the cancellation of projects. To improve the prediction capability of induced seismicity, we have developed an analytical solution to compute the induced stress along both permeable and impermeable faults as a result of reservoir pressurization/depletion. The solution is based on the inclusion theory and has been validated by comparing it with a numerical solution. We have performed a comprehensive analysis on induced seismicity potential due to reservoir pressurization or depletion and obtained the following conclusions:

* The induced seismicity potential of impermeable faults is always larger than that of permeable faults under any initial and injection conditions. Generally, the maximum size of fault undergoing slip for the impermeable faults is 3–5 times greater than that for permeable ones under a given initial stress ratio and pressure buildup. Moreover, an impermeable fault would rupture at a higher stress ratio, that is, less deviatoric stress, and at a smaller pressure buildup than a permeable one.

* Pore pressure buildup and stress concentration at the corners during reservoir pressurization/depletion are the main reasons for fault reactivation and induced seismicity. Stress concentration not only amplifies the effect of the hydraulic properties of faults, but also results in the occurrence of reversed slip at the corners for normal faults with a normal faulting stress regime, and of normal slip for thrust faults with a thrust faulting stress regime.

* The slip potential of permeable faults resulting from reservoir pressurization/depletion increases with the fault offset because of the change in stress concentration, which implies that non-displaced
permeable faults constitute a safer choice for site selection. In contrast, the offset has no impact on the slip potential of impermeable faults, because the effect of stress concentration is always the same.

* For a given pressure buildup, the difference between the critical and onset dip angles, that is, the range of dip angles favorable to slip, reduces for increasing the initial stress ratio and equals to zero at its critical stress ratio. This means that geological sites with a higher in situ stress ratio (lower initial deviatoric stress) are intrinsically less prone to fluid injection-induced seismicity. This finding is useful for site selection in geo-energy projects.

* For a given stress ratio, the range of dip angles favorable to slip reduces for decreasing the pressure buildup and equals to zero at its critical pressure buildup as a larger pressure buildup is more likely to induce seismicity. Thus, the methodology of the maximum sustainable injection pressure to minimize the risk of inducing seismicity is feasible and should be available for designing and managing the injection parameters.

* The fault slip potential increases if we ignore the stress arching, that is, assuming the effective stress reduction is equivalent to the pore pressure buildup during injection, for both permeable and impermeable faults in a normal faulting stress regime, which implies that the induced seismicity potential is overestimated and the maximum sustainable pressure is underestimated.

* Our analytical solution includes two limitations resulting from our simplifying assumptions to solve this complex problem: (1) we assume a linear elastic material, which is physically unrealistic for rock materials; in reality, a nonlinear elastic or inelastic deformation cannot be avoided during the reservoir pressurization/depletion; (2) we assume a quasi-steady-state pore pressure change in the reservoir and neglect the transient effect of flow, so the calculated induced stress are still overestimated, particularly in the low-permeable fault zone because pore pressure will eventually diffuse into the portion of the caprock and baserock in contact with the reservoir. Such limitations are worthy to be investigated further and complemented by more detailed numerical solutions

## Appendix A: Inclusion Theory and Induced Stress

The fundamental concept of the inclusion theory lies in a series of imaginary steps involving cutting, transforming, and restoring the inclusion itself (Eshelby, 1957; Mura, 1987; Rudnicki, 2011). In the last step, restoring the inclusion to its original shape and size, that is, with zero strain, corresponds to the application of a stress field $\sigma^*$ (eigenstress) inside the inclusion, to neutralize the volumetric eigenstrain $\epsilon^*$ that it would undergo if unbounded, and of a simultaneous body force $f$ (restoring force) over the entire matrix, to keep the stress equilibrium (Jansen et al., 2019; Rudnicki, 2011). For the case of uniform pore pressure change inside the inclusion, $\Delta p$, then is

$$f(x, y, z) = -\alpha \Delta p \mathbf{n}_{\Gamma}(x, y, z), \qquad (A1)$$

where the vector $f$ has three components in the coordinate directions for a three-dimensional (3D) problem, $\mathbf{n}_{\Gamma}$ is the unit normal vector pointing outward from the boundary ($\Gamma$) of the inclusion, $x, y$, and $z$ are the Cartesian coordinates, and $\alpha$ is Biot's coefficient. The term $-\alpha \Delta p$ indicates the normal eigenstress for the 3D scenario.

Affected by the restoring force field, a displacement field $u_i(x, y, z)$ is provoked in the inclusion and its surrounding rock as

$$u_i(x, y, z) = -\alpha \Delta p \iint_{\Gamma} \tilde{\mathbf{g}}_i(x, y, z, \xi, \psi) \cdot \mathbf{n}_{\Gamma}(\xi, \xi, \psi) d\Gamma, \qquad (A2)$$

where

$$\tilde{\mathbf{g}}_i(x, y, z, \xi, \psi) = [\tilde{g}_{ix}(x, y, z, \xi, \psi), \tilde{g}_{iy}(x, y, z, \xi, \psi), \tilde{g}_{iz}(x, y, z, \xi, \psi)]^T, \qquad (A3)$$
and $\tilde{g}_{ij}(x, y, z, \zeta, \xi, \psi)$ is the Green's function describing the displacement at any point $(x, y, z)$ under a unit body force at point $(\zeta, \xi, \psi)$, $\zeta, \xi, \psi$ are the coordinate values on $\Gamma$, and subscripts $i$ and $j$ are free indexes with $i, j \in (x, y, z)$ in the 3D Cartesian space.

For an infinite elastic unbounded domain, $\tilde{g}_{ij}(x, y, z, \zeta, \xi, \psi)$ can be expressed as (Love, 1944; Mura, 1987)

$$ \tilde{g}_{ij}(x, y, z, \zeta, \xi, \psi) = \frac{\delta_{ij}}{4\pi\mu\tilde{R}} - \frac{1}{16\pi\mu(1-\nu)} \frac{\partial^2}{\partial x_i \partial x_j} \tilde{R}, \quad (A4) $$

where $\mu$ and $\nu$ are the shear modulus and Poisson's ratio, respectively, $\delta_{ij}$ is the Kronecker delta, which equals 1 if $i=j$ or 0 if $i \neq j$, and

$$ \tilde{R}^2 = (x - \zeta)^2 + (y - \xi)^2 + (z - \psi)^2. \quad (A5) $$

Equation A2 is a standard surface integral and can be transformed into a volume integral by applying the Gauss's divergence theorem

$$ u_i(x, y, z) = -\alpha \Delta p \iint_{\Omega} \nabla \tilde{g}_i(x, y, z, \zeta, \xi, \psi) d\Omega, \quad (A6) $$

where $\Omega$ means the inclusion volume. Note that now $\zeta, \xi, \psi$ in Equation A6 denote the coordinate values in the domain $\Omega$.

Under the assumption of plane strain and integrating Equation A4 along the out-of-plane dimension (z) yields (Jansen et al., 2019; Mura, 1987)

$$ \tilde{g}_{ij}(x, y, \zeta, \xi) = \frac{1}{8\pi\mu(1-\nu)} \left[ \frac{\bar{x}_i \cdot \bar{x}_j}{R^2} - (3-4\nu)\delta_{ij} \ln R \right] \quad i, j \in (x, y), \quad (A7) $$

where

$$ \bar{x}_i = |x - \zeta| \text{ or } |y - \xi|, \quad (A8) $$

$$ R^2 = (x - \zeta)^2 + (y - \xi)^2. \quad (A9) $$

Equation A6 can be simplified into

$$ u_i(x, y) = -\alpha \Delta p \iint_{\Omega} \frac{\partial \tilde{g}_{ix}}{\partial x} + \frac{\partial \tilde{g}_{iy}}{\partial y} d\Omega = D \iint_{\Omega} g_i(x, y, \zeta, \xi) d\Omega, \quad (A10) $$

where $g_i(x, y, \zeta, \xi)$ is

$$ g_x(x, y, \zeta, \xi) = \frac{x - \zeta}{2R^2}, \quad (A11) $$

$$ g_y(x, y, \zeta, \xi) = \frac{y - \xi}{2R^2}, \quad (A12) $$

and D is a dimensionless scaling parameter

$$ D = \frac{(1 - 2\nu)\alpha\Delta p}{2\pi(1 - \nu)\mu} = \frac{(1 + \nu)\alpha\Delta p}{3\pi(1 - \nu)K}, \quad (A13) $$

where K is the bulk modulus.

For an isotropic linear elastic material, the Hooke's equation relating the stress $\sigma_{ij}$ to the strain $\epsilon_{ij}$ tensor is
$$ \sigma_{ij} = 2\mu\epsilon_{ij} + \frac{2\mu\nu}{1-2\nu}\delta_{ij}\epsilon_{kk}, \qquad (A14) $$

where $\epsilon_{kk}$ is the volumetric strain and $\epsilon_{ij}$ equals the symmetric part of the displacement gradient

$$ \epsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right). \qquad (A15) $$

Substituting Equations A10 into A15, and taking the results into Equation A14, yields the expression of the stress field in the whole matrix induced by the restoring force **f**

$$ \sigma_{ij}(x, y) = C \iint_{\Omega} g_{ij}(x, y, \zeta, \xi) d\Omega, \qquad (A16) $$

where $C = \mu D$ is a scaling parameter and $g_{ij}(x, y, \zeta, \xi)$ is the Green's function for stress at $(x, y)$ given a unit point force at $(\zeta, \xi)$

$$ g_{xx}(x, y, \zeta, \xi) = 2 \frac{\partial g_x}{\partial x} = \frac{(y - \zeta)^2 - (x - \zeta)^2}{R^4}, \qquad (A17) $$

$$ g_{yy}(x, y, \zeta, \xi) = 2 \frac{\partial g_y}{\partial y} = \frac{(x - \zeta)^2 - (y - \zeta)^2}{R^4}, \qquad (A18) $$

$$ g_{xy}(x, y, \zeta, \xi) = 2 \frac{\partial g_x}{\partial y} = -2 \frac{(x - \zeta)(y - \zeta)}{R^4}. \qquad (A19) $$

For the inclusion, the stress field, however, is also affected by the eigenstress resulting from the application of the surface traction to bring the inclusion back to its initial configuration (Eshelby, 1957). Since the surface traction depends on the geometry of the inclusion, such eigenstress is different from that of the 3D scenario for the case of the plane strain problem. Regarding that the eigenstrain is a pure dilatational strain, only the normal components of eigenstress have a finite value (Eshelby, 1957; Soltanzadeh & Hawkes, 2008), thus

$$ \sigma_{ij}(x, y) = C \iint_{\Omega} g_{ij}(x, y, \zeta, \xi) d\Omega + \sigma^* \delta_{ij} \delta_{\Omega} = C G_{ij} x, y + \sigma^* \delta_{ij} \delta_{\Omega}, \qquad (A20) $$

where $G_{ij}(x, y)$ denotes the surface integral of Green's function $g_{ij}(x, y, \zeta, \xi)$, and $\delta_{\Omega}$ is the modified Kronecker delta

$$ \delta_{\Omega} = \begin{cases} 1 & \text{if } (x, y) \in \Omega, \\ 0 & \text{if } (x, y) \notin \Omega. \end{cases} \qquad (A21) $$

From the perspective of stress arching effect, the induced stress field caused by pore pressure change in the inclusion can be described as (Soltanzadeh & Hawkes, 2008)

$$ \sigma_{ij} = -\gamma_{ij} \alpha \Delta p, \qquad (A22) $$

where $\gamma_{ij}$ is the normalized stress arching ratio, which depends on the geometry of the inclusion. The minus means that injection corresponds to compression.

For the ellipsoidal inclusions, the stress and strain fields are uniform for all points inside the inclusion (Eshelby, 1957; Rudnicki, 1999). And for the case in which the inclusion extends infinite in one direction (implying that we can apply the plane strain assumption), such as the elliptic cylindrical inclusion, the sum of the stress arching ratio in the other two directions is a constant (Soltanzadeh, 2009; Soltanzadeh & Hawkes, 2008)
$$ \gamma_{xx} + \gamma_{yy} = \frac{(1 - 2\nu)}{1 - \nu}. \qquad (A23) $$

According to Equations A20–A23, we find that

$$ C[G_{xx}(x, y) + G_{yy}(x, y)] + 2\sigma^* = -\frac{1 - 2\nu}{1 - \nu}\alpha\Delta p. \qquad (A24) $$

Given that $g_{xx}(x, y, \zeta, \xi) = -g_{yy}(x, y, \zeta, \xi)$ (Equations A17 and A18), thus

$$ \sigma^* = -\frac{1 - 2\nu}{2(1 - \nu)}\alpha\Delta p. \qquad (A25) $$

Although Equation A25 is derived for an elliptic cylindrical inclusion, it is valid for all the geometric inclu- sions that can be regarded as plane strain problem, because the inclusion only undergoes a pure dilatational deformation. Introducing Equation A25 into Equation A20, we obtain the final expression of the induced stress tensor both in the inclusion and its surrounding rock, as given in Equation 1.

## Appendix B: Surface Integral of Green's Function for Stress

In this appendix, we solve the surface integrals of Equation 1. The entire inclusion is divided into two trap- ezoids by the fault. The surface integral of one function over a trapezoid can be regarded as the sum of the integrations over a triangle and a rectangle (Figure B1). The integrand function is the Green's function for stress given in Equations A17–A19. We need to consider the existence of singularities for values of $\zeta$ and $\xi$ equal x and y, respectively, where the Green's function becomes infinite. This only occurs for points (x, y)

<figure><img src="image_16.png" /><figcaption>**Figure B1.** Schematic geometry of triangular and rectangular domains in an arbitrary coordinate system.</figcaption></figure>

located inside the inclusion, for which the integral becomes improper. Thus, we perform the regular bound- ed integral for (x, y) located outside the inclusion, while for (x, y) located inside the inclusion we solve the improper integral by excluding a neighborhood of the singularity. To generalize the integration of Green's function, we apply an arbitrary coordinate system, as shown in Figure B1. After solving the surface integrals over the triangular and rectangular domains, the solutions are transformed into the coordinate system of Figure 1.
## Rectangular Inclusion

For $(x, y)$ located outside of the rectangle, we calculate the surface integral of the Green's function for horizontal stress $g_{xx}(x, y, \zeta, \xi)$ using standard techniques and we obtain

$$ 
\begin{aligned}
G_{xx}^{\text{out}}(x, y) &= \iint_{rp}^{s} g_{xx}(x, y, \zeta, \xi) d\zeta d\xi \\
&= \int_{r}^{s} \left[ -\frac{x - \zeta}{(x - \zeta)^2 + (y - \zeta)^2} \right]_{p}^{q} d\xi \\
&= \left[ a \tan \frac{y - \zeta}{x - q} - a \tan \frac{y - \zeta}{x - p} \right]_{r}^{s} \\
&= a \tan \frac{y - s}{x - q} - a \tan \frac{y - s}{x - p} - a \tan \frac{y - r}{x - q} + a \tan \frac{y - r}{x - p}.
\end{aligned} 
\quad (B1) $$

If $(x, y)$ is located in the rectangle, we need to consider the improper integration of Green's function. Note that if the standard integration (like in equation (B1)) is applied for such points, the result depends on the order of integration, because Fubini's theorem (DiBenedetto, 2016) does not hold in this case. We solve the improper integral excluding a neighborhood of the singular point $(\zeta, \xi) = (x, y)$, which makes the integrand function being bounded. We observe that $g_{xx}(x, y, \zeta, \xi)$ is antisymmetric with respect to the line $y-\xi = x-\zeta$, that is, $g_{xx}(x, y, \zeta, \xi) = -g_{xx}(y, x, \zeta, \xi)$. This means that the integral is zero for a square domain centered in the singular point, because the contribution from the triangle above the symmetry line cancels with the contribution from the triangle below. Therefore, we can exclude a square neighborhood of any size contained in the domain $\Omega$ and centered in the singular point. In the rest of the domain, Fubini's theorem holds and we can apply the standard sequential integration (Equation B1). For the square domain, the actual contribution is zero while standard integration technique gives $\pi$, which can be checked by substituting $y-r = x-p = s-y = q-x$ into Equation (B1). Therefore, we must remove the offending contribution, so that the integral for points inside the inclusion is

$$ G_{xx}^{\text{in}}(x, y) = G_{xx}^{\text{out}}(x, y) - \pi. \quad (B2) $$

Note that this result is independent on the order of integration.

Equations B1 and B2 only differ for the last term $-\pi$, which is a result of the improper integral for points located inside the reservoir. We therefore express the surface integral in the general form

$$ G_{xx}(x, y) = a \tan \frac{y - r}{x - p} - a \tan \frac{y - r}{x - q} - a \tan \frac{y - s}{x - p} + a \tan \frac{y - s}{x - q} - \pi \delta_{\Omega}, \quad (B3) $$

with $\delta_{\Omega}$ defined in Equation A21.

There are four singularities at the corners (A, B, C, and D) of the rectangular inclusion domain (Figure B1), where the arguments in Equation B3 become indefinite (i.e., 0/0). The integral is in fact not defined at these points, and it is discontinuous there (different values are obtained when approaching it from one side or the other).

The same procedure can be applied to integrate the Green's function for vertical stress $g_{yy}(x, y, \zeta, \xi)$. Since $g_{yy}(x, y, \zeta, \xi) = -g_{xx}(x, y, \zeta, \xi)$ (see Equations A17 and A18), it follows that

$$ G_{yy}(x, y) = \int_{rp}^{s} \int_{q}^{q} g_{yy}(x, y, \zeta, \xi) d\zeta d\xi = -G_{xx}(x, y). \quad (B4) $$

To integrate the Green's function for shear stress $g_{xy}(x, y, \zeta, \xi)$, we follow the same procedure used to integrate $g_{xx}(x, y, \zeta, \xi)$. We observe that $g_{xy}(x, y, \zeta, \xi)$ is antisymmetric with respect to $y-\xi = 0$ or $x-\zeta = 0$,
therefore the integral over a square domain centered in the singular point is 0. We also observe that the
standard sequential integration technique in this case gives

$$
\begin{align*}
G_{xy}(x,y) &= \int_{pr}^{qs} \int_{xy} g_{xy}(x,y,\xi,\xi') d\xi d\xi' \\
&= \int_{p}^{q} \left[ \frac{-(x-\xi)}{(x-\xi)^2 + (y-\xi)^2} \right]_{r}^{s} d\xi = \frac{1}{2} \left[ \ln \frac{(x-\xi)^2 + (y-s)^2}{(x-\xi)^2 + (y-r)^2} \right]_{p}^{q} \tag{B5} \\
&= \frac{1}{2} \ln \frac{\left[ (x-q)^2 + (y-s)^2 \right] \left[ (x-p)^2 + (y-r)^2 \right]}{\left[ (x-q)^2 + (y-r)^2 \right] \left[ (x-p)^2 + (y-s)^2 \right]}.
\end{align*}
$$

Application of Equation B5 for a square neighborhood of the singular point gives 0, which demonstrates
that Equation B5 is valid for both the cases of (x, y) located outside or inside the inclusion.

In Equation B5, there are four singularities at the corners A, B, C, and D.

Triangular Inclusion

For the integration over a triangular inclusion, we apply the same technique adopted for the case of the
rectangular inclusion. As an example, the integral of the Green's function for horizontal stress is

$$
\begin{align*}
G_{xx}(x,y) &= \int_{r}^{s} \int_{\xi \cot \theta}^{p} \frac{(y-\xi)^2 - (x-\xi)^2}{\left[(x-\xi)^2 + (y-\xi)^2\right]^2} d\xi d\xi - \pi \delta_{\Omega} \\
&= \int_{r}^{s} -\left[\frac{x-\xi}{(x-\xi)^2 + (y-\xi)^2}\right]_{\xi \cot \theta}^{p} d\xi - \pi \delta_{\Omega} \tag{B6} \\
&= \left[\frac{\tan \frac{y-\xi}{x-p}}{x-p}\right]_{r}^{s} + \int_{r}^{s} \frac{x-\xi \cot \theta}{(x-\xi \cot \theta)^2 + (y-\xi)^2} d\xi - \pi \delta_{\Omega},
\end{align*}
$$

where the last term $-\pi\delta_{\Omega}$ takes into account the effect of the improper integral for points located inside the inclusion. The integration of the second term on the right-hand side of Equation B6 is nontrivial and a potential solution can be obtained by the following transformation:

$$
x - \xi \cot \theta = \left[ \left( (x - \xi \cot \theta) \cos^2 \theta + (y - \xi) \sin \theta \cos \theta \right) + \left[ (x - \xi \cot \theta) \sin^2 \theta - (y - \xi) \sin \theta \cos \theta \right] \right]. \quad (B7)
$$

Substituting Equation B7 into Equation B6 and after some derivations, we obtain

$$
\int_{r}^{s} \frac{x - \xi \cot \theta}{(x - \xi \cot \theta)^2 + (y - \xi)^2} d\xi = -\frac{\sin \theta \cos \theta}{2} \left[ \ln \left[ (x - \xi \cot \theta)^2 + (y - \xi)^2 \right] \right]_{r}^{s} \\
-\sin^2 \theta \left[ \tan \frac{(x - \xi \cot \theta) \cot \theta + (y - \xi)}{x - y \cot \theta} \right]_{r}^{s}.
\quad (B8)
$$

The final expression of the integration of the Green's function for horizontal stress over a triangular inclusion domain is
$$G_{xx}(x,y) = a \tan \frac{y-s}{x-p} - a \tan \frac{y-r}{x-p} - \frac{\sin \theta \cos \theta}{2} \ln \frac{f_2(x,y,s)}{f_2(x,y,r)} \\ -[f_1(x,y,s) - f_1(x,y,r)]\sin^2 \theta - \pi \delta_{\Omega}, \quad (B9)$$

where functions $f_1$ and $f_2$ are

$$f_1(x, y, \hat{y}) = a \tan \frac{(x - \hat{y} \cot \theta) \cot \theta + (y - \hat{y})}{x - y \cot \theta}, \quad (B10)$$

$$f_2(x, y, \hat{y}) = (x - \hat{y} \cot \theta)^2 + (y - \hat{y})^2. \quad (B11)$$

Similarly, integrations of the Green's function for the vertical and shear stress components over a triangular inclusion domain are

$$G_{yy}(x,y) = \int_{r}^{s} \int_{\xi \cot \theta}^{p} g_{yy}(x,y,\xi) d\xi d\xi = -G_{xx}(x,y), \quad (B12)$$

$$\begin{aligned} G_{xy}(x,y) &= \int_{r}^{s} \int_{\xi \cot \theta}^{p} g_{xy}(x,y,\xi) d\xi d\xi \\ &= [f_1(x,y,s) - f_1(x,y,r)]\sin \theta \cos \theta \\ &\quad - \frac{\sin^2 \theta}{2} \ln \frac{f_2(x,y,s)}{f_2(x,y,r)} + \frac{1}{2} \ln \frac{f_3(x-p,y-s)}{f_3(x-p,y-r)}, \end{aligned} \quad (B13)$$

where function $f_3$ is

$$f_3(x - \hat{x}, y - \hat{y}) = (x - \hat{x})^2 + (y - \hat{y})^2, \quad (B14)$$

The singularities are located at the three corners A, D, and E in this case.

## Trapezoidal Inclusion

We apply the superposition principle and obtain the integrations of the Green's function for the stress components over a trapezoid as a combination of the above integrals for rectangular and triangular inclusions as

$$\begin{aligned} G_{xx}(x,y) &= \int_{r}^{s} \int_{\xi \cot \theta}^{p} g_{xx}(x,y,\xi) d\xi d\xi + \iint_{r,p}^{s,q} g_{xx}(x,y,\xi) d\xi d\xi \\ &= a \tan \frac{y-s}{x-q} - a \tan \frac{y-r}{x-q} - \frac{\sin \theta \cos \theta}{2} \ln \frac{f_2(x,y,s)}{f_2(x,y,r)} \\ &\quad - [f_1(x,y,s) - f_1(x,y,r)] \sin^2 \theta - \pi \delta_{\Omega}, \end{aligned} \quad (B15)$$

$$G_{yy}(x,y) = \int_{r}^{s} \int_{\xi \cot \theta}^{p} g_{yy}(x,y,\xi) d\xi d\xi + \iint_{r,p}^{s,q} g_{yy}(x,y,\xi) d\xi d\xi = -G_{xx}(x,y), \quad (B16)$$

$$\begin{aligned} G_{xy}(x,y) &= \int_{r}^{s} \int_{\xi \cot \theta}^{p} g_{xy}(x,y,\xi) d\xi d\xi + \iint_{r,p}^{s,q} g_{xx}(x,y,\xi) d\xi d\xi \\ &= [f_1(x,y,s) - f_1(x,y,r)]\sin \theta \cos \theta \\ &\quad - \frac{\sin^2 \theta}{2} \ln \frac{f_2(x,y,s)}{f_2(x,y,r)} + \frac{1}{2} \ln \frac{f_3(x-q,y-s)}{f_3(x-q,y-r)}. \end{aligned} \quad (B17)$$
Note that the last term in Equation B15 is a consequence of having improper integrals either for the rectangle or the triangle, i.e., (x, y) falls within the rectangle or the triangle. In Equations B15–B17, the singularities are located at the four corners A, B, C, and E of the trapezoid.

## Application into a Specific Coordinate System

For a permeable fault, pore pressure changes at both sides of the fault during fluid injection or production and thus, the inclusion is composed of two trapezoids. The analytical expressions of integration of the Green's function for the stress components in the general coordinate system is transformed into the coordinate system of Figure 1 as

$$ 
\begin{aligned}
G_{xx}(x, y) &= \int_{a}^{-b} \int_{\xi \cot \theta}^{-b-b \cot \theta} g_{xx}(x, y, \xi) d\xi d\xi + \int_{a}^{-b} \int_{-b \cot \theta}^{-c} g_{xx}(x, y, \xi) d\xi d\xi \\
&\quad + \int_{-a \xi \cot \theta}^{b} \int_{a}^{b} g_{xx}(x, y, \xi) d\xi d\xi + \int_{-ab \cot \theta}^{b} \int_{a}^{d} g_{xx}(x, y, \xi) d\xi d\xi \\
&= \int_{a}^{-b} \int_{\xi \cot \theta}^{-c} g_{xx}(x, y, \xi) d\xi d\xi + \int_{-a \xi \cot \theta}^{b} \int_{a}^{d} g_{xx}(x, y, \xi) d\xi d\xi \quad (B18) \\
&= a \tan \frac{y+b}{x+c} - a \tan \frac{y-a}{x+c} + a \tan \frac{y-b}{x-d} - a \tan \frac{y+a}{x-d} \\
&\quad - [f_1(x, y, -b) - f_1(x, y, a) + f_1(x, y, b) - f_1(x, y, -a)] \sin^2 \theta \\
&\quad - \frac{\sin \theta \cos \theta}{2} \ln \frac{f_2(x, y, -b) f_2(x, y, b)}{f_2(x, y, a) f_2(x, y, -a)} - \pi \delta_{\Omega},
\end{aligned}
$$

$$ G_{yy}(x, y) = \int_{a}^{-b} \int_{\xi \cot \theta}^{-c} g_{yy}(x, y, \xi) d\xi d\xi + \int_{-a \xi \cot \theta}^{b} \int_{a}^{d} g_{yy}(x, y, \xi) d\xi d\xi = -G_{xx}(x, y), \quad (B19) $$

$$ 
\begin{aligned}
G_{xy}(x, y) &= \int_{a}^{-b} \int_{\xi \cot \theta}^{-c} g_{xy}(x, y, \xi) d\xi d\xi + \int_{-a \xi \cot \theta}^{b} \int_{a}^{d} g_{xy}(x, y, \xi) d\xi d\xi \quad (B20) \\
&= [f_1(x, y, -b) - f_1(x, y, a) + f_1(x, y, b) - f_1(x, y, -a)] \sin \theta \cos \theta
\end{aligned}
$$

If c and d have finite value, the inclusion represents a finite reservoir and there are eight singularities for the solution: the four corners P1, P2, P3, and P4 located on the fault (Figure 1) and the four corners on the outer boundary of the faulted reservoir. If c = d = ∞, the reservoir is infinite and there are only the four singularities located on the fault.

For a vertical fault (θ = 90), such integral solutions read as

$$ 
\begin{aligned}
G_{xx}(x, y) &= -G_{yy}(x, y) \\
&= \int_{-b-c}^{a} \int_{-b-c}^{0} g_{xx}(x, y, \xi) d\xi d\xi + \int_{-a0}^{b} \int_{-a0}^{d} g_{xx}(x, y, \xi) d\xi d\xi \\
&= a \tan \frac{y-a}{x} - a \tan \frac{y-a}{x+c} - a \tan \frac{y+b}{x} + a \tan \frac{y+b}{x+c} \\
&\quad + a \tan \frac{y-b}{x-d} - a \tan \frac{y-b}{x} - a \tan \frac{y+a}{x-d} + a \tan \frac{y+a}{x} - \pi \delta_{\Omega},
\end{aligned}
$$

$$ 
\begin{aligned}
G_{xy}(x, y) &= \int_{-c-b}^{0} \int_{-a}^{a} g_{xy}(x, y, \xi) d\xi d\xi + \int_{0-a}^{b} \int_{-a}^{d} g_{xy}(x, y, \xi) d\xi d\xi \\
&= \frac{1}{2} \ln \frac{f_3(x, y-a) f_3(x+c, y+b) f_3(x-d, y-b) f_3(x, y+a)}{f_3(x, y+b) f_3(x+c, y-a) f_3(x-d, y+a) f_3(x, y-b)},
\end{aligned}
$$
and for a zero offset fault in which $a = b$, we have

$$ 
\begin{aligned}
G_{xx}(x,y) &= -G_{yy}(x,y) \\
&= a \tan \frac{y+b}{x+c} - a \tan \frac{y-a}{x+c} + a \tan \frac{y-b}{x-d} - a \tan \frac{y+a}{x-d} - \pi \delta_{\Omega},
\end{aligned}
\quad (B23) $$

$$ G_{xy}(x,y) = \frac{1}{2} \ln \frac{f_3(x+c, y+b) f_3(x-d, y-b)}{f_3(x+c, y-a) f_3(x-d, y+a)}. \quad (B24) $$

## Appendix C: Induced Stress on the Fault Plane

The expression for normal and tangential stress on the fault plane can be found by transforming the main coordinate system into a coordinate system placed on the fault and oriented along it. We first operate the translation, thus the horizontal and vertical integral solutions of the Green's function for the stress components on the fault plane are evaluated by setting $x = y \cot \theta$ in Equations B18-B20, such that

$$ 
\begin{aligned}
G_{xx}(y \cot \theta, y) &= -G_{yy}(y \cot \theta, y) \\
&= -\frac{\pi}{2} \left[ \operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) + \operatorname{sgn}(y-b) - \operatorname{sgn}(y+a) \right] \sin^2 \theta \\
&\quad + a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} + a \tan \frac{y-b}{y \cot \theta - d} \\
&\quad - a \tan \frac{y+a}{y \cot \theta - d} - \frac{\sin \theta \cos \theta}{2} \ln \frac{f_4(y,b)}{f_4(y,a)} - \pi \delta_{\circ},
\end{aligned}
\quad (C1) $$

$$ 
\begin{aligned}
G_{xy}(y \cot \theta, y) = & \frac{\pi}{2} \left[ \operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) + \operatorname{sgn}(y-b) - \operatorname{sgn}(y+a) \right] \sin \theta \cos \theta \\
& + \frac{1}{2} \ln \frac{f_3(y \cot \theta + c, y+b) f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a) f_3(y \cot \theta - d, y+a)} - \frac{\sin^2 \theta}{2} \ln \frac{f_4(y,b)}{f_4(y,a)}.
\end{aligned}
\quad (C2) $$

where sgn (•) is the sign function defined as 1 if (•) > 0, 0 if (•) = 0 or -1 if (•) <0, and function $f_4$ is

$$ f_4(y, \hat{y}) = (y + \hat{y})^2 (y - \hat{y})^2, \quad (C3) $$

Substituting Equations C1 and C2 into Equation 1 yields the x-y planar induced stress along the fault plane. Subsequently, applying the stress transformation with axis rotation to transform the coordinate system, one can obtain the expression of induced normal $\sigma_n(y \cot \theta, y)$ and tangential $\tau(y \cot \theta, y)$ stress components along a fault plane with an arbitrary dip angle. Considering the sign convention and geometry adopted here, such stress transformation equations are

$$ \sigma_n(y \cot \theta, y) = \frac{\sigma_{xx} + \sigma_{yy}}{2} + \frac{\sigma_{xx} - \sigma_{yy}}{2} \cos(2\theta + \pi) + \sigma_{xy} \sin(2\theta + \pi), \quad (C4) $$

$$ \tau(y \cot \theta, y) = -\frac{\sigma_{xx} - \sigma_{yy}}{2} \sin(2\theta + \pi) + \sigma_{xy} \cos(2\theta + \pi). \quad (C5) $$

Substituting Equations C1, C2, and 1 into Equations C4 and C5 and normalizing by the scaling parameter C (Equation 2), the dimensionless induced normal $\bar{\sigma}_n(y \cot \theta, y)$ and tangential $\bar{\tau}(y \cot \theta, y)$ stress components along the fault plane are
$$
\begin{align}
\bar{\sigma}_n(y \cot \theta, y) &= (G_{xx} - \pi \delta_{\Omega}) \sin^2 \theta + (G_{yy} - \pi \delta_{\Omega}) \cos^2 \theta - 2G_{xy} \sin \theta \cos \theta \nonumber \\
&= -\frac{\pi}{2} [\operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) + \operatorname{sgn}(y-b) - \operatorname{sgn}(y+a)] \sin^2 \theta \nonumber \\
&\quad - \cos 2\theta (a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} + a \tan \frac{y-b}{y \cot \theta - d}) \tag{C6} \\
&\quad -a \tan \frac{y+a}{y \cot \theta - d} - \pi \delta_{\Omega}) + \frac{\sin 2\theta}{4} \ln \frac{f_4(y,b)}{f_4(y,a)} \nonumber \\
&\quad -\frac{\sin 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b) f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a) f_3(y \cot \theta - d, y+a)} - \pi \delta_{\Omega}, \nonumber
\end{align}
$$

$$
\begin{align}
\bar{\tau}(y \cot \theta, y) &= \sin \theta \cos \theta [(G_{xx} - \pi \delta_{\Omega}) - (G_{yy} - \pi \delta_{\Omega})] - G_{xy} (\cos^2 \theta - \sin^2 \theta) \nonumber \\
&= -\frac{\pi}{4} [\operatorname{sgn}(y+b) - \operatorname{sgn}(y-a) + \operatorname{sgn}(y-b) - \operatorname{sgn}(y+a)] \sin 2\theta \nonumber \\
&\quad + \sin 2\theta (a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} + a \tan \frac{y-b}{y \cot \theta - d}) \tag{C7} \\
&\quad -a \tan \frac{y+a}{y \cot \theta - d} - \pi \delta_{\Omega}) - \frac{\sin^2 \theta}{2} \ln \frac{f_4(y,b)}{f_4(y,a)} \nonumber \\
&\quad -\frac{\cos 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b) f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a) f_3(y \cot \theta - d, y+a)}, \nonumber
\end{align}
$$

where $G_{ij}(y \cot \theta, y)$ is shortened into $G_{ij}$ for convenience, and the four corners P1, P2, P3, and P4 on the fault plane are singularities (Figure 1). To simplify Equations B18–B20 into Equations C1 and C2, we apply the general rule of taking the right limit, i.e., the limit that the argument approaches the fault from its right-hand side, as the value of the fault plane. Thus, the segment P1–P2 belongs to the inclusion, while the segment P3–P4 belongs to the surroundings for a permeable fault in the above equations.

In particular, for a vertical fault ($\theta = 90^\circ$), we obtain

$$
\begin{align}
\bar{\sigma}_n(0, y) &= G_{xx}(0, y) - \pi \delta_{\Omega} \nonumber \\
&= \frac{\pi}{2} [\operatorname{sgn}(y-a) - \operatorname{sgn}(y+b) - \operatorname{sgn}(y-b) + \operatorname{sgn}(y+a)] \tag{C8} \\
&\quad + a \tan \frac{y+b}{c} - a \tan \frac{y-a}{c} - a \tan \frac{y-b}{d} + a \tan \frac{y+a}{d} - 2\pi \delta_{\Omega}, \nonumber
\end{align}
$$

$$
\begin{align}
\bar{\tau}(0, y) &= G_{xy}(0, y) \nonumber \\
&= \frac{1}{2} \ln \frac{f_3(c, y+b) f_3(d, y-b)}{f_3(c, y-a) f_3(d, y+a)} - \frac{1}{2} \ln \frac{f_4(y, b)}{f_4(y, a)}, \tag{C9}
\end{align}
$$

and for a zero offset fault, we have

$$
\begin{align}
\bar{\sigma}_n(y \cot \theta, y) &= -\cos 2\theta (a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} \nonumber \\
&\quad + a \tan \frac{y-b}{y \cot \theta - d} - a \tan \frac{y+a}{y \cot \theta - d} - \pi \delta_{\Omega}) \tag{C10} \\
&\quad -\frac{\sin 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b) f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a) f_3(y \cot \theta - d, y+a)} - \pi \delta_{\Omega}, \nonumber
\end{align}
$$
$$ \begin{aligned} \bar{\tau}(y \cot \theta, y) = & \sin 2\theta (a \tan \frac{y+b}{y \cot \theta + c} - a \tan \frac{y-a}{y \cot \theta + c} \\ & + a \tan \frac{y-b}{y \cot \theta - d} - a \tan \frac{y+a}{y \cot \theta - d} - \pi \delta_{\Omega}) \\ & - \frac{\cos 2\theta}{2} \ln \frac{f_3(y \cot \theta + c, y+b)f_3(y \cot \theta - d, y-b)}{f_3(y \cot \theta + c, y-a)f_3(y \cot \theta - d, y+a)}. \end{aligned} \quad (C11) $$

## Data Availability Statement

The Python files used to produce most of the figures in this paper and the input files used for the numerical simulations in CODE_BRIGHT are available from the CSIC data repository (https://digital.csic.es/handle/10261/221111).

## Acknowledgments

H. Wu acknowledges the financial support received from the AGAUR (Generalitat de Catalunya) through the “grant for universities and research centers for the recruitment of new research personnel (FI-2019)”. V. Vilarrasa acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme through the Starting Grant GEoREST (www.georest.eu), Grant agreement no. 801809. V. Vilarrasa also acknowledges support by the Spanish Ministry of Science and Innovation (Project CEX2018-000794-S). S.D. Simone acknowledges financial support from the SAD2018 project funded by the Brittany Region and from ANR LabCom Project eLabo ANR-17-LCV2-0012. M. Saaltink acknowledges financial support from the “HEATSTORE” project (170153-44011), which has been subsidized through the ERANET Cofund GEOTHERMICA (Grant agreement no. 731117), from the European Commission and the Spanish Ministry of Science, Innovation and Universities (PCI2018-092933). F. Parisio acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—project number PA 3451/1-1. The authors thank Tomas Aquino for his advice on the integral solutions.

## References

Ahmadi, H., Lotfollahi-Yaghin, M. A., Yong-Bo, S., & Aminfar, M. H. (2012). Parametric study and formulation of outer-brace geometric stress concentration factors in internally ring-stiffened tubular KT-joints of offshore structures. Applied Ocean Research, 38, 74–91.

Bai, B., Li, X., Wu, H., Wang, Y., & Liu, M. (2017). A methodology for designing maximum allowable wellhead pressure for CO₂ injection: Application to the Shenhua CCS demonstration project, China. Greenhouse Gases: Science and Technology, 7(1), 158–181. https://doi.org/10.1002/ghg.1640

Bandilla, K. W., & Celia, M. A. (2017). Active pressure management through brine production for basin-wide deployment of geologic carbon sequestration. International Journal of Greenhouse Gas Control, 61, 155–167.

Bourne, S. J., & Oates, S. J. (2017). Extreme threshold failures within a heterogeneous elastic thin sheet and the spatial-temporal development of induced seismicity within the Groningen gas field. Journal of Geophysical Research: Solid Earth, 122(10), 299–310. https://doi.org/10.1002/2017JB014356

Buijze, L., Van den Bogert, P. A., Wassing, B., & Orlic, B. (2019). Nucleation and arrest of dynamic rupture induced by reservoir depletion. Journal of Geophysical Research: Solid Earth, 124, 3620–3645. https://doi.org/10.1029/2018JB016941

Buijze, L., van den Bogert, P. A., Wassing, B. B., Orlic, B., & ten Veen, J. (2017). Fault reactivation mechanisms and dynamic rupture modelling of depletion-induced seismic events in a Rotliegend gas reservoir. Netherlands Journal of Geosciences, 96(5), s131–s148.

Caine, J. S., Evans, J. P., & Forster, C. B. (1996). Fault zone architecture and permeability structure. Geology, 24(11), 1025–1028.

Castelletto, N., Gambolati, G., & Teatini, P. (2013). Geological CO₂ sequestration in multi-compartment reservoirs: Geomechanical challenges. Journal of Geophysical Research: Solid Earth, 118(5), 2417–2428. https://doi.org/10.1002/jgrb.50180

Chang, K. W., & Segall, P. (2016). Injection-induced seismicity on basement faults including poroelastic stressing. Journal of Geophysical Research: Solid Earth, 121(4), 2708–2726. https://doi.org/10.1002/2015JB012561

Chiaramonte, L., White, J. A., & Hao, Y. (2013). Probabilistic risk assessment of mechanical deformation due to CO₂ injection in a compartmentalized reservoir. Paper Presented at 47th US Rock Mechanics/Geomechanics Symposium, 23–26 June 2013. ARMA 13-577. American Rock Mechanics Association, San Francisco, CA.

Courant, R., & John, F. (1989). Introduction to calculus and analysis volume II. New York, NY: Springer-Verlag. https://doi.org/10.1007/978-1-4613-8958-3

De Simone, S., Carrera, J., & Vilarrasa, V. (2017). Superposition approach to understand triggering mechanisms of post-injection induced seismicity. Geothermics, 70, 85–97.

Deichmann, N., & Giardini, D. (2009). Earthquakes induced by the stimulation of an enhanced geothermal system below Basel (Switzerland). Seismological Research Letters, 80(5), 784–798.

Del Potro, R., & Diez, M. (2015). Induced seismicity in underground gas storage - the case of castor, offshore NE Spain. Paper Presented at 77th EAGE Conference and Exhibition - Workshops (Vol. 2015, pp. 1–3). Madrid: European Association of Geoscientists & Engineers.

DiBenedetto, E. (2016). Real analysis (2nd ed.). Boston, MA: Birkhäuser Advanced Texts: Basler Lehrbücher. https://doi.org/10.1007/978-1-4939-4005-9

Duverger, C., Godano, M., Bernard, P., Lyon-Caen, H., & Lambotte, S. (2015). The 2003–2004 seismic swarm in the western Corinth rift: Evidence for a multiscale pore pressure diffusion process along a permeable fault system. Geophysical Research Letters, 42(18), 7374–7382. https://doi.org/10.1002/2015GL065298

Ellsworth, W. L. (2013). Injection-induced earthquakes. Science, 341(6142), 1225942.

Ellsworth, W. L., Giardini, D., Townend, J., Ge, S., & Shimamoto, T. (2019). Triggering of the Pohang, Korea, earthquake (Mw 5.5) by enhanced geothermal system stimulation. Seismological Research Letters, 90(5), 1844–1858.

Ellsworth, D., Spiers, C. J., & Niemeijer, A. R. (2016). Understanding induced seismicity. Science, 354(6318), 1380–1381.

Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society London. Series A. Mathematical and Physical Sciences, 241, 376–396. http://dx.doi.org/10.1098/rspa.1957.0133

Foulger, G. R., Wilson, M. P., Gluyas, J. G., Julian, B. R., & Davies, R. J. (2018). Global review of human-induced earthquakes. Earth-Science Reviews, 178, 438–514.

Galis, M., Ampuero, J. P., Mai, P. M., & Cappa, F. (2017). Induced seismicity provides insight into why earthquake ruptures stop. Science Advances, 3(12) eaap7528.

Galis, M., Ampuero, J. P., Mai, P. M., & Kristek, J. (2019). Initiation and arrest of earthquake ruptures due to elongated overstressed regions. Geophysical Journal International, 217(3), 1783–1797.

Garagash, D. I., & Germanovich, L. N. (2012). Nucleation and arrest of dynamic slip on a pressurized fault. Journal of Geophysical Research, 117(B10310). https://doi.org/10.1029/2012JB009209

Géraud, Y., Diraison, M., & Orellana, N. (2006). Fault zone geometry of a mature active normal fault: A potential high permeability channel (Pirgaki fault, Corinth rift, Greece). Tectonophysics, 426(1–2), 61–76.
Goebel, T. H. W., Weingarten, M., Chen, X., Haffener, J., & Brodsky, E. E. (2017). The 2016 Mw5.1 Fairview, Oklahoma earthquakes: Evidence for long-range poroelastic triggering at > 40 km from fluid disposal wells. *Earth and Planetary Science Letters, 472*, 50–61.

Grigoli, F., Cesca, S., Rinaldi, A. P., Manconi, A., López-Comino, J. A., Clinton, J. F., et al. (2018). The November 2017 Mw 5.5 Pohang earthquake: A possible case of induced seismicity in South Korea. *Science, 360*(6392), 1003–1006.

Hansen, O., Gilding, D., Nazarian, B., Osdal, B., Ringrose, P., Kristoffersen, J. B., et al. (2013). Snøhvit: The history of injecting and storing 1 Mt CO₂ in the fluvial Tubåen Fm. *Energy Procedia*, 37, 3565–3573.

Häring, M. O., Schanz, U., Ladner, F., & Dyer, B. C. (2008). Characterisation of the Basel 1 enhanced geothermal system. *Geothermics*, 37(5), 469–495.

Haug, C., Nuchter, J.-A., & Henk, A. (2018). Assessment of geological factors potentially affecting production-induced seismicity in North German gas fields. *Geomechanics for Energy and the Environment, 16*, 15–31. https://doi.org/10.1016/j.gete.2018.04.002

Jansen, J. D., Singhal, P., & Vossepoel, F. C. (2019). Insights from closed-form expressions for injection-and production-induced stresses in displaced faults. *Journal of Geophysical Research: Solid Earth*, 124(7), 7193–7212. https://doi.org/10.1029/2019JB017932

Jha, B., & Juanes, R. (2014). Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering. *Water Resources Research*, 50(5), 3776–3808. https://doi.org/10.1002/2013WR015175

Juanes, R., Castiñeira, D., Fehler, M. C., Hager, B. H., Jha, B., Shaw, J. H., & Plesch, A. (2017). *Coupled flow and geomechanical modeling, and assessment of induced seismicity, at the Castor underground gas storage project*, Huelva: Mesa de la Ria.

Karvounis, D. C., Gischig, V. S., & Wiemer, S. (2014). *Towards a real-time forecast of induced seismicity for enhanced geothermal systems. Shale energy engineering 2014: Technical challenges, environmental issues, and public policy* (pp. 246–255). Pittsburgh, PA: American Society of Civil Engineers.

Kim, K. H., Ree, J. H., Kim, Y., Kim, S., Kang, S. Y., & Seo, W. (2018). Assessing whether the 2017 Mw 5.4 Pohang earthquake in South Korea was an induced event. *Science*, 360(6392), 1007–1009.

King, G. C., Stein, R. S., & Lin, J. (1994). Static stress changes and the triggering of earthquakes. *Bulletin of the Seismological Society of America*, 84(3), 935–953.

Langenbruch, C., & Zoback, M. D. (2016). How will induced seismicity in Oklahoma respond to decreased saltwater injection rates? *Science Advances*, 2(11) e1601542.

Lee, K. K., Ellsworth, W. L., Giardini, D., Townend, J., Ge, S., Shimamoto, T., et al. (2019). Managing injection-induced seismic risks. *Science*, 364(6442), 730–732.

Lehner, F. K. (2019). *An analysis of depletion-induced fault stressing report SR.19.01549: Shell Global Solutions International B.V.*

Lele, S., Hsu, S.-Y., Garzon, J., DeDontney, N., Searles, K., Gist, G., et al. (2016). *Geomechanical modeling to evaluate production-induced seismicity at Groningen field. Paper Presented at Abu Dhabi International Petroleum Exhibition and Conference*. Abu Dhabi: Society of Petroleum Engineers. https://doi.org/10.2118/183554-MS

Love, A. E. H. (1944). *A treatise on the mathematical theory of elasticity* (4th ed.). New York, NY: Dover Publications.

Makhnenko, R. Y., & Labuz, J. F. (2015). Dilatant hardening of fluid-saturated sandstone. *Journal of Geophysical Research: Solid Earth*, 120(2), 909–922. https://doi.org/10.1002/2014JB011287

Mathias, S. A., Hardisty, P. E., Trudell, M. R., & Zimmerman, R. W. (2009). Screening and selection of sites for CO₂ sequestration based on pressure buildup. *International Journal of Greenhouse Gas Control*, 3(5), 577–585.

Mura, T. (1987). *Micromechanics of defects in solids* (2nd ed.). Dordrecht: Martinus Nijho Publishers.

NAM (2016). *Winningsplan Groningen gasveld (Technical Report)*. Nederlands Aardolie Maatschappij. http://www.nam.nl/gas-en-oliewinning/groningen-gasveld/winningsplan-groningen-gasveld.html

Nordbotten, J. M., Celia, M. A., & Bachu, S. (2005). Injection and storage of CO₂ in deep saline aquifers: Analytical solution for CO₂ plume evolution during injection. *Transport in Porous Media*, 58, 339360.

Olivella, S., Carrera, J., Gens, A., & Alonso, E. E. (1994). Nonisothermal multiphase flow of brine and gas through saline media. *Transport in Porous Media*, 15(3), 271–293.

Olivella, S., Gens, A., Carrera, J., & Alonso, E. E. (1996). Numerical formulation for a simulator (CODE_BRIGHT) for the coupled analysis of saline media. *Engineering Computations*, 13(7), 87–112.

Orlic, B., & Wassing, B. (2013). A study of stress change and fault slip in producing gas reservoirs overlain by elastic and viscoelastic caprocks. *Rock Mechanics and Rock Engineering*, 46, 421–435. https://doi.org/10.1007/s00603-012-0347-6

Orlic, B., Wassing, B., & Geel, C. (2013). Field scale geomechanical modeling for prediction of fault stability during underground gas storage operations in a depleted gas field in the Netherlands. *Paper Presented at 47th US Rock Mechanics/Geomechanics Symposium*. Paper ARMA-13-300. 23-26 June. American Rock Mechanics Association, San Francisco, CA.

Rohmer, J., Pluymakers, A., & Renard, F. (2016). Mechano-chemical interactions in sedimentary rocks in the context of CO₂ storage: Weak acid, weak effects? *Earth-Science Reviews*, 157, 86–110.

Rudnicki, J. W. (1999). *Alteration of regional stress by reservoirs and other inhomogeneities: Stabilizing or destabilizing?* In G. Vouille, & P. Berest (Eds.). Proceedings of the ninth international congress on rock mechanics (Vol. 3, pp. 1629–1728). Paris: International Society for Rock Mechanics and Rock Engineering.

Rudnicki, J. W. (2002). *Eshelby transformations, pore pressure and fluid mass changes, and subsidence*. Paper presented at Second Biot conference on poromechanics (poromechanics II) (pp. 26–28). Grenoble: CRC Press.

Rudnicki, J. W. (2011). Eshelby's technique for analyzing inhomogeneities in geomechanics. In: Mechanics of crustal rocks (pp. 4372). Vienna: Springer.

Rutqvist, J., Birkholzer, J., Cappa, F., & Tsang, C. F. (2007). Estimating maximum sustainable injection pressure during geological sequestration of CO₂ using coupled fluid flow and geomechanical fault-slip analysis. *Energy Conversion and Management*, 48(6), 1798–1807.

Rutqvist, J., Rinaldi, A. P., Cappa, F., Jeanne, P., Mazzoldi, A., Urpi, L., & (2016). Fault activation and induced seismicity in geological carbon storage—Lessons learned from recent modeling studies. *Journal of Rock Mechanics and Geotechnical Engineering*, 8(6), 789–804.

Sacks, I. S., Linde, A. T., Suyehiro, S., & Snoke, J. A. (1978). Slow earthquakes and stress redistribution. *Nature*, 275(5681), 599–602.

Segall, P. (1985). Stress and subsidence resulting from subsurface fluid withdrawal in the epicentral region of the 1983 Coalinga earthquake. *Journal of Geophysical Research*, 90(B8), 6801–6816. https://doi.org/10.1029/JB090iB8p06801

Segall, P. (1989). Earthquakes triggered by fluid extraction. *Geology*, 17(10), 942–946. https://doi.org/10.1130/0091-7613

Segall, P. (1992). Induced stresses due to fluid extraction from axisymmetric reservoirs. *Pure and Applied Geophysics*, 139(3-4), 535–560. https://doi.org/10.1007/BF00879950

Segall, P., & Fitzgerald, S. (1998). A note on induced stress changes in hydrocarbon reservoirs and geothermal reservoirs. *Tectonophysics*, 289, 117–128. https://doi.org/10.1016/S0040-1951(97)00311-9
Segall, P., Grasso, J., & Mossop, A. (1994). Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France. *Journal of Geophysical Research*, 99(B8), 15423–15438. https://doi.org/10.1029/94JB00989

Shapiro, S. A., & Dinske, C. (2009). Fluid-induced seismicity: Pressure diffusion and hydraulic fracturing. *Geophysical Prospecting*, 57(2), 301–310.

Simpson, D. W., Leith, W. S., & Scholz, C. H. (1988). Two types of reservoir-induced seismicity. *Bulletin of the Seismological Society of America*, 78(6), 2025–2040.

Soltanzadeh, H. (2009). *Geomechanical analysis of caprock integrity*. Saskatoon: University of Saskatchewan.

Soltanzadeh, H., & Hawkes, C. (2008). Semi-analytical models for stress change and fault reactivation induced by reservoir production and injection. *Journal of Petroleum Science and Engineering*, 60, 71–85. https://doi.org/10.1016/j.petrol.2007.05.006

Soltanzadeh, H., & Hawkes, C. (2009). Assessing fault reactivation tendency within and surrounding porous reservoirs during fluid production or injection. *International Journal of Rock Mechanics and Mining Sciences*, 46, 1–7. https://doi.org/10.1016/j.ijrmms.2008.03.008

Terakawa, T., Miller, S. A., & Deichmann, N. (2012). High fluid pressure and triggered earthquakes in the enhanced geothermal system in Basel, Switzerland. *Journal of Geophysical Research*, 117(B07305). https://doi.org/10.1029/2011JB008980

Van den Bogert, P. (2015). Impact of various modelling options on the onset of fault slip and fault slip response using 2-dimensional finite-element modelling. Retrieved From Shell Global Solutions International, B.V. https://nam-feitenencijfers.data-app.nl/download/rapport/604ef59b-5ac5-4770-93b1-e501c87da3a9?open=true

Van Wees, J., Buijze, L., Van Thienen-Visser, K., Nepveu, M., Wassing, B., Orlic, B., & (2014). Geomechanics response and induced seismicity during gas field production in the Netherlands. *Geothermics*, 52, 206–219. https://doi.org/10.1016/j.geothermics.2014.05.004

Van Wees, J., Fokker, P., Van Thienen-Visser, K., Wassing, B., Osinga, S., Orlic, B., et al. (2017). Geomechanical models for induced seismicity in the Netherlands: Inferences from simplified analytical, finite element and rupture model approaches. *Netherlands Journal of Geosciences*, 96(5), S183–S202. https://doi.org/10.1017/njg.2017.38

Vilarrasa, V., Carrera, J., Olivella, S., Rutqvist, J., & Laloui, L. (2019). Induced seismicity in geologic carbon storage. *Solid Earth*, 10(3), 871–892.

Vilarrasa, V., Makhnenko, R., & Gheibi, S. (2016). Geomechanical analysis of the influence of CO₂ injection location on fault stability. *Journal of Rock Mechanics and Geotechnical Engineering*, 8(6), 805–818.

Villaseñor, A., Herrmann, R. B., Gaite, B., & Ugalde, A. (2020). Fault reactivation by gas injection at an underground gas storage off the east coast of Spain. *Solid Earth*, 11(1), 63–74.

Wang, L., Bai, B., Li, X., Liu, M., Wu, H., & Hu, S. (2016). An analytical model for assessing stability of pre-existing faults in caprock caused by fluid injection and extraction in a reservoir. *Rock Mechanics and Rock Engineering*, 49(7), 2845–2863.

Wu, H., Bai, B., & Li, X. (2018). An advanced analytical solution for pressure buildup during CO₂ injection into infinite saline aquifers: The role of compressibility. *Advances in Water Resources*, 112, 95–105. https://doi.org/10.1016/j.advwatres.2017.12.010

Wu, H., Bai, B., Li, X., Gao, S., Liu, M., & Wang, L. (2016). An explicit integral solution for pressure buildup during CO₂ injection into infinite saline aquifers. *Greenhouse Gas: Science and Technology*, 6(5), 633–647. https://doi.org/10.1002/ghg.1601

Zbinden, D., Pio Rinaldi, A., Urpi, L., & Wiemer, S. (2017). On the physics based processes behind production-induced seismicity in natural gas fields. *Journal of Geophysical Research: Solid Earth*, 122(5), 3792–3812. https://doi.org/10.1002/2017JB014003

Zhou, Q., Birkholzer, J. T., Tsang, C. F., & Rutqvist, J. (2008). A method for quick assessment of CO₂ storage capacity in closed and semi-closed saline formations. *International Journal of Greenhouse Gas Control*, 2(4), 626–639.