# Fluid flow migration, rock stress and deformation due to a crustal fault slip in a geothermal system: A poro-elasto-plastic perspective

Felipe Sáez-Leiva, Daniel Hurtado, Muriel Gerbault, Javiera Ruz-Ginouves, Pablo Iturrieta, José Cembrano

► To cite this version:

Felipe Sáez-Leiva, Daniel Hurtado, Muriel Gerbault, Javiera Ruz-Ginouves, Pablo Iturrieta, et al. Fluid flow migration, rock stress and deformation due to a crustal fault slip in a geothermal system: A poro-elasto-plastic perspective. Earth and Planetary Science Letters, 2023, 604, pp.117994. 10.1016/j.epsl.2023.117994. hal-03969333

HAL Id: hal-03969333

https://hal.science/hal-03969333

Submitted on 2 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

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

# Fluid flow migration, rock stress and deformation due to a crustal fault slip in a geothermal system: a poro-elasto-plastic perspective

Felipe Sáez-Leiva<sup>a,b,*</sup>, Daniel E. Hurtado<sup>a,f</sup>, Muriel Gerbault<sup>c</sup>, Javiera Ruz-Ginouves<sup>a,b,g</sup>, Pablo Iturrieta<sup>d,e</sup> and José Cembrano<sup>a,b</sup>

<sup>a</sup>Department of Structural and Geotechnical Engineering, School of Engineering, Pontificia Universidad Católica de Chile, Santiago, Chile

<sup>b</sup>Andean Geothermal Centre of Excellence (CEGA), Santiago, Chile

<sup>c</sup>GET/UMR5563 (UPS, CNRS, IRD, CNES); Obs. Midi-Pyrénées, Université Paul Sabatier, Toulouse, France

<sup>d</sup>Helmholtz Centre Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany

<sup>e</sup>University of Potsdam, Potsdam, Germany

<sup>f</sup>Schools of Engineering, Medicine and Biological Sciences, Institute for Biological and Medical Engineering, Pontificia Universidad Católica de Chile, Santiago, Chile

<sup>g</sup>Department of Geology, University of Otago, Dunedin, New Zealand

13

16

18

20

22

24

26

28

30

32

34

36

38

40

42

44

46

48

50

52

54

56

58

60

62

64

66

## ARTICLE INFO

### Keywords:

Finite Element Method

Poromechanics

Geothermal System

Fault Zone

Suction Pump

### Abstract

Geothermal systems are commonly genetically and spatially associated with volcanic complexes, which in turn, are located nearby crustal fault systems. Faults can alter fluid flow in their surroundings, potentially acting as barriers or conduits for fluids, depending on their architecture and slip-rate. However, this fundamental control on fluid migration is still poorly constrained. Most previous modeling efforts on volcanic and hydrothermal processes consider either only fluid flow in their formulations, or only a mechanical approach, and seldom a full, monolithic coupling between both. In this work, we present a poro-elasto-plastic Finite Element Method (FEM) to address the first-order, time-dependent control that a strike-slip crustal fault exerts on a nearby geothermal reservoir. For the model setting, we selected the Planchón-Peteroa geothermal system in the Southern Andes Volcanic Zone (SAVZ), for which the geometry and kinematics of a potentially seismogenic fault and fluid reservoir is constrained from previous geological and geophysical studies. We assess the emergence and diffusion of fluid pressure domains due to fault slip, as well as the development of tensile/dilational and compressive/contractional domains in the fault¹ surroundings. Mean stress and volumetric strain magnitudes in these domains range between ±1 [MPa] and ±10⁻⁴ [-], respectively. Our results show the appearance of negative and positive fluid pressure domains in these dilational and contractional regions, respectively. We also investigate the spatial and temporal evolution of such domains resulting from changes in fault permeability and shear modulus, fluid viscosity, and rock rheology. These variations in fluid pressure alter the trajectory of the reservoir fluids, increasing migration to the eastern half of the fault, reaching a maximum fluid flux of 8 to 70 times the stationary flux. Pressure-driven fluid diffusion over time causes fluid flow to return to the stationary state between weeks to months after fault slip. These results suggest that the mechanism that exerts a first-order control is similar to a suction pump, whose duration heavily depends on fault permeability and fluid viscosity. We also show how a von Mises plasticity criterion locally enhances fluid flow. The transient process analyzed in this work highlights the importance of addressing the solid-fluid coupling in numerical models for volcano-tectonic studies.

44

## 1. Introduction

45

47

49

51

53

55

57

Renewable energies such as geothermal power, have experienced exponential growth in their generation capacity over the last decades (Pandey et al., 2018). Apart from electricity generation from high-enthalpy reservoirs (> 150°C), low-temperature reservoirs (30° – 90°C) can also be used in a variety of ways such as air conditioning, water desalinization, among others (e.g. Muñoz et al., 2015; Pandey et al., 2018).

Amongst the most favorable areas for geothermal development, volcanic arcs in transpressional convergent margins such as the Andes, often host crustal fault systems that accommodate most of the strike-slip component of the

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

<footnote>*Corresponding author</footnote>

<footnote>fcsaez@uc.cl (F. Sáez-Leiva)</footnote>

<footnote>ORCID(s): 0000-0002-6687-737X (F. Sáez-Leiva)</footnote>
subduction vector (Teyssier et al., 1995). Depending on their architecture (Caine et al., 1996) and slip-rate (Sibson, 1990), faults can act either as barriers or conduits for fluids, thus controlling fluid flow migration in the upper crust. Moreover, as fluids reduce the normal stress on rocks, they control rock deformation, causing hydraulic fracturing and, consequently, increasing rock permeability. Even though previous studies have looked into fault-controlled geothermal systems (e.g. Bellani et al., 2004; Roquer et al., 2017; Rowland and Simmons, 2012), the spatial and temporal evolution of fluid migration, coupled to bedrock strain still lacks comprehensive understanding. This is relevant considering that fluid flow can induce seismicity due to fault reactivation (Moeck et al., 2009; Legrand et al., 2011), which is one of the main social obstacles to geothermal energy development (e.g. Trutnevyte and Ejderyan, 2017).

In this work, we analyze the time-dependent control that a strike-slip crustal fault exerts on a nearby geothermal reservoir, by defining a synthetic geometrical setting based on the Planchón-Peteroa geothermal system located in the Southern Andean Volcanic Zone (SAVZ). We simulate the evolution of fluid pressure, and bedrock strain and stress fields through poro-elasto-plastic numerical simulations. We aim to understand the role of fault systems in the spatial and temporal evolution of fluid migration and bedrock deformation at a time scale from hours to weeks after fault slip. To achieve this, we impose a left-lateral fault slip and solve for the coupled fields of fluid pressure and rock displacement, as well as for plastic strain. We deduce the Darcy velocity of the fluid, the volumetric strain and the rock stress field. Fluid flux evolution provides a first-order determination of the areas of possible hydraulic fracturing. Numerical verification was carried out with classical poroelastic and elasto-plastic benchmarks. The originality of our approach not only stands in the hydromechanical setting that we consider, but also in the implementation of a fully-coupled technique (details in Section 2 and 3). Finally, we discuss our results with respect to previous studies of seismo-tectonic transient processes (Dorbath et al., 2009; Legrand et al., 2011) and fluid migration mechanisms (Nur and Booker, 1972; Sibson, 1985, 1987, 2000).

## 2. State of the art and SAVZ geotectonic setting

Fluid migration in the upper crust can be altered due to fault slip, mainly through two mechanisms: the suction pump (Sibson, 1985, 1987, 2000) and the fault-valve (Sibson, 1990, 2000) mechanisms. The suction pump mechanism refers to fluid migration towards a slipping strike-slip fault system because of the emergence of dilational jogs, from pre-existing stepovers on the trace of the fault. These jogs experience a sudden dilation, causing a fluid pressure drop inducing suction on the surrounding fluid. Suction pumping can promote the migration of a volume of fluid between $10^6 – 10^7 m^3$ during crustal fault slip events (Sibson, 1987; Manga and Wang, 2015). The fault-valve mechanism, instead, refers to the migration of over-pressurized fluids, usually originated in deep regions, to shallower crustal domains where fluid pressure is lower, through a crustal fault which can cut through impermeable barriers. It can be driven by seismic fault slip, or by a reduction of the effective stress due to the over-pressurized fluids.

Field efforts have studied the influence of crustal faults over fluid migration according to the width ratio between their core and their damage zone (Caine et al., 1996). The higher the ratio, the more the fault acts as a barrier to fluid flow. Observations suggest that the damage zone is responsible for enhancing fault permeability. Furthermore, the fluid transmissibility of the fault surface is dependent on its seismic cycle. During a seismic event, transmissibility increases between 6-9 orders of magnitude, enhancing fluid flow through the fault zone (Gudmundsson, 2011). This increase, however, is transitory because of fracture sealing due to hydrothermal fluid precipitation (Dobson et al., 2003).

At sample scale, laboratory experiments have characterized permeability evolution upon hydromechanical loading. Triaxial drained compression tests show that low-porosity rocks can experience permeability reduction at the onset of loading due to elastic crack closure and compaction. Then, due to damage onset, rocks can undergo a permeability increase up to 3 orders of magnitude (Heiland, 2003; Mitchell and Faulkner, 2008). In contrast, high-porosity rocks show a reduction in permeability at damage onset. Zhu and Wong (1997) observed a reduction between 2-3 orders near compactive yield stress, which was attributed to pore collapse and an increase in pore tortuosity due to shear-enhanced compaction. This porosity-dependent behavior indicates that permeability in faults damage zones is controlled by macrofractures (Mitchell and Faulkner, 2008). Fracture networks can enhance permeability up to 10 orders of magnitude in intrusive rocks, such as granodiorites: from $10^{-19} – 10^{-18}$ to $10^{-10} – 10^{-9} m^2$ (Sepúlveda et al., 2020).
Numerical studies have addressed fluid flow in fault zones through different approaches, amongst which thermo-hydraulic coupling allows to evaluate convective flow directions and intermittence in both present-day and fossil geothermal systems (e.g. Kissling et al., 2018; Duwiquet et al., 2021). To assess the mechanical interaction of the host rock with fluid flow, purely mechanical models have used the plastic yield threshold as a proxy for fluid pathways at a regional or local scale, (e.g. Iturrieta et al., 2017; Kavanagh et al., 2018; Ruz-Ginouves et al., 2021). The next steps towards coupled hydromechanical and thermohydromechanical approaches are the so called one-way coupling models, in which the stress and strain fields from a mechanical model are used as input for a hydraulic or thermo-hydraulic model, where fluid flow evolution is assessed according to the (initial) mechanical conditions, or vice-versa. Such models allow to assess earthquakes control on fluid and heat flow evolution in fault-controlled geothermal systems, for example (e.g. Dempsey et al., 2013). However, these do not fully capture the poro-mechanical interaction between solid and fluid. Two-way coupling models in turn, address the coupling between the solid and fluid equations through a simultaneous integration, thus accounting for the full poro-mechanical effects (Prevost, 2013). Such models have been mostly used to constrain geothermal fluid flow distribution across a fault zone with spatially-variable permeability (e.g. Duwiquet et al., 2021), and crustal fault reactivation due to fluid migration (e.g. Alghannam and Juanes, 2020; Vilarrasa et al., 2021), CO₂ sequestration (e.g. Rutqvist et al., 2016) and nuclear waste storage (e.g. Rinaldi et al., 2022). However, to our knowledge, two-way coupling approaches remain to be comprehensively addressed in the specific setting of fluid flow responding to fault slip.

Two-way coupling models can be implemented sequentially or monolithically. Sequential coupling, usually achieved via a fixed-stress approach (Kim et al., 2011), requires less computational power and can solve both the solid and fluid/heat problems on dedicated softwares, such as the widely used TOUGH-FLAC software (e.g. Rinaldi et al., 2022). Unconditional numerical stability, however, is not always ensured as certain conditions have to be fulfilled when defining the sequential scheme (see Kim et al., 2011, for further details). Fully, or monolithically, coupled models, which solve both sets of equations simultaneously, are unconditionally stable. And even though such models are computationally more expensive than sequential models, the constant growth of computational power over the last decade make them increasingly achievable (e.g. Haagenson et al., 2020; Liu and Huang, 2021). Within this framework, plastic behavior has been accounted for in sequentially coupled approaches (e.g. Rinaldi et al., 2022), but seldom in fully-coupled approaches (Liu and Huang, 2021). In the present work we develop a fully-coupled hydromechanical method, considering an elasto-plastic solid, to address the spatial and temporal evolution of fluid flow in a high-enthalpy geothermal system.

High-enthalpy geothermal systems in transpressional convergent margins, such as the Andes, are genetically and spatially associated with volcanic complexes, their respective magmatic reservoirs, and regional crustal fault systems. The Southern Andes Volcanic Zone (SAVZ), located between 33°-46° S displays numerous examples of said geothermal systems, such as Mariposa (e.g. Hickson et al., 2011), Tolhuaca (e.g. Melosh et al., 2010), Villarrica (e.g. Pavez et al., 2020), and Planchón-Peteroa (e.g. Pearce et al., 2020). Even though further local studies are yet to be conducted, the Planchón-Peteroa inferred geothermal system displays the fundamental traits of an Andean high-enthalpy geothermal system: WNW-striking seismogenic faults, highly conductive rock volumes in between and close to these faults, and hotsprings. Magnetotelluric and seismic surveys in the region have constrained the geometry and location of these regional fault systems and the inferred fluid reservoirs and have determined that fluids are closely related to rock deformation (Pearce et al., 2020) (Figure 1). We focus on the influence of the seismically active WNW-striking Andean Transverse Fault (ATF) (seismic cluster Cls1) on the geothermal reservoir located NE of the Planchón volcano at 4-8 km depth, inferred from the C1 resistive anomaly.

## 3. Methodology

To address the nature of our setting, we implement a fully-coupled poro-elasto-plastic formulation using a finite element method. A small-strain formulation is considered as we aim at addressing first-order factors that influence fluid flow and bedrock strain in a geothermal system.

### 3.1. Poromechanics formulation and governing equations

Poromechanics studies the simultaneous deformation of a solid porous material and pore fluid flow (Coussy, 2004) through the definition of a Representative Volume Element (RVE) of a saturated porous medium. We adopt
<figure><img src="image_3.png" /><figcaption>**Figure 1:** (a) Regional scale map of a segment of the Southern Andes and the Southern Andes Volcanic Zone (SAVZ) (Cembrano and Lara, 2009). (b) Local scale map with the main geological structures such as volcanoes and faults, and the results of the geophysical surveys conducted by Pearce et al. (2020). MT anomalies and seismic hypocenters are located at 6 kilometers depth. The convergence vector is taken from Angermann et al. (1999); Métois et al. (2012), and the traces of the faults, as well as the focal mechanisms, are taken from Spagnotto et al. (2015).</figcaption></figure>

standard mechanical sign convention: tensile stresses and anticlockwise shear-induced rotations are positive. Bold fonts represent vectors and tensors and normal fonts denote scalars.

We consider solid and fluid as a continuum that satisfies two main equations, which we solve for solid displacement **u** and fluid pressure p. First, mass balance laws are imposed on the domain of analysis. Considering constant porosity and fluid density, and no fluid source/sink term, mass balance over the RVE is defined as:

$$ \alpha \frac{\partial \varepsilon_v(\mathbf{u})}{\partial t} + \frac{1}{M} \frac{\partial p}{\partial t} + \text{div}(\mathbf{q}(p)) = 0, \quad (1) $$

where α and M are Biot's coefficient and modulus respectively, ε_v(u) = div(u) is the volumetric strain of the solid skeleton, and q(p) is the Darcy velocity vector or the volume-averaged velocity vector. We define Darcy velocity in terms of matrix permeability κ, dynamic viscosity μ, and fluid pressure, and neglect the gravity term:
$$q(p) = -\frac{K}{\mu}\nabla p. \qquad (2)$$

Note that Eq. 2 assumes isotropic permeability. The RVE must satisfy the linear momentum balance equation which, after neglecting inertia and gravity terms, is defined as:

$$\mathrm{div}(\sigma(\mathbf{u}, p)) = 0, \qquad (3)$$

where $\sigma(\mathbf{u}, p)$ is the Cauchy stress tensor. Because of the solid-fluid interactions, $\sigma$ is defined by Eq. 4 in the elastic range:

$$\sigma(\mathbf{u}, p) = \left(K - \frac{2G}{3}\right) \mathrm{trace}(\varepsilon(\mathbf{u})) \mathbf{I} + 2G\varepsilon(\mathbf{u}) - \alpha p \mathbf{I}, \qquad (4)$$

where $K$ is the matrix bulk modulus, $G$ the shear modulus, and $\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla(\mathbf{u}) + \nabla^T(\mathbf{u}))$ the infinitesimal strain tensor. To model elasto-plasticity, we consider the following additive decomposition of the strain tensor:

$$\varepsilon = \varepsilon^e + \varepsilon^p, \qquad (5)$$

in which $\varepsilon^e$ and $\varepsilon^p$ are the elastic and plastic components of the strain tensor, respectively. Then, the following incremental constitutive function is used to define an increment on the stress tensor $\Delta\sigma$:

$$\Delta\sigma = \mathbf{D} : \Delta\varepsilon^e + \mathbf{D}^{\mathrm{plas}} : \Delta\varepsilon^p, \qquad (6)$$

where $\mathbf{D}$ and $\mathbf{D}^{\mathrm{plas}}$ are the fourth-order elastic and elastoplastic tangent operators, respectively. A von Mises plasticity model with isotropic linear hardening is implemented for the solid skeleton, as described in de Souza Neto et al. (2008). Plastic onset is defined by the yield surface $f(\sigma, \bar{\varepsilon}^p)$:

$$f(\sigma, \bar{\varepsilon}^p) = \sqrt{3J_2(s(\sigma))} - (\sigma_{y0} + H\bar{\varepsilon}^p) \le 0, \qquad (7)$$

where $\bar{\varepsilon}^p$ is the accumulated plastic strain, $H$ is the hardening modulus, $\sigma_{y0}$ is the initial yield strength, and $J_2(s(\sigma)) = \frac{1}{2}s(\sigma) : s(\sigma)$ is the second invariant of the deviatoric stress tensor $s = \sigma - \frac{1}{3}\mathrm{trace}(\sigma)\mathbf{I}$. Eq. 7 states that the von Mises plasticity model is pressure-insensitive, only involving the deviatoric stress tensor. Linear hardening is considered as this considerably simplifies the plastic return mapping procedure, reducing it to a single nonlinear equation. This rather standard approach from the engineering numerical community provides unconditional stability, and is suitable for the derivation of tangent operators - essential to our formulation (de Souza Neto et al., 2008). Whereas this behavior is less frequent in the geomechanical community when applied to rock masses, it is argued to mimic natural softening due to grain size reduction in developing gouge zones, a process which actually also occurs upon localisation in pressure-dependent yield materials (e.g. Gerbault et al., 1998; Rybacki et al., 2021).

Poro-elasto-plasticity is defined by Eqs. 1-3, considering the incremental constitutive function 6. Although pressure-dependent criteria, such as the Drucker-Prager model, have been identified to be more adequate for upper crustal rock behavior and are expected to affect pore fluid pressure, we decide to proceed here step by step and to start with a von Mises approach, for two main reasons: 1) a recent compilation study from borehole measurements in saturated felsic crustal domains indicates that, at depths greater than ca. 5 km depth, the pressure dependency of the yield stress tends to vanish in fluid saturated crusts (e.g. Suppe, 2014), and 2) the results presented below actually demonstrate how even a von Mises criterion influences the fluid flow. Hence, we postpone the study of a pressure-dependent yield criterion for future work. This point is further discussed in Section 5.
## 3.2. Numerical Scheme

The poro-elasto-plastic formulation is numerically solved using a Finite-Element scheme implemented in the Python opensource library FEniCS (a collaborative effort to solve complex differential equations, Alnæs et al. (2015)). We use the latest legacy version (2019.1.0) of FEniCS which documents previous poromechanical studies (e.g. Haagenson et al., 2020). Our implementation manages the coupling between the poro-elastic and the elasto-plastic equations and, because this is not commonly found in the literature (to our knowledge), we provide the key mathematical reasoning. We start from the strong formulation (S) described by Eqs. 1 and 3, previously time-discretized using a backward Euler scheme, with test functions φ and η for pressure and solid displacement, respectively. Multiplying (S) by the test functions results in the following weak formulation (W):

$$ 
\begin{aligned}
\mathbf{G}(\varphi, \mathbf{u}, \mathbf{p}) ={}& \alpha \int_{\Omega} (\mathrm{div} \, \mathbf{u}^{t+1} - \mathrm{div} \, \mathbf{u}^t) \varphi \, d\Omega + \frac{1}{M} \int_{\Omega} (\mathbf{p}^{t+1} - \mathbf{p}^t) \varphi \, d\Omega \\
& + \frac{\kappa}{\mu} \Delta t \left( \int_{\Omega} \nabla \mathbf{p}^{t+1} \cdot \nabla \varphi \, d\Omega - \int_{\partial \Gamma_q} (\nabla \mathbf{p} \cdot \mathbf{n}) \varphi \, d\Gamma \right) = 0 \quad \forall \varphi \in H^1(\Omega), \quad t \in [t, t+1],
\end{aligned}
\quad (8) $$

$$ \mathbf{G}(\eta, \mathbf{u}, \mathbf{p}) = \int_{\Omega} \sigma^{t+1} : \varepsilon(\eta) \, d\Omega - \int_{\partial \Gamma_t} \bar{\mathbf{t}} \cdot \mathbf{n} \, d\Gamma = 0 \quad \forall \eta \in H^1(\Omega), \quad t \in [t, t+1], \quad (9) $$

where $\bar{\mathbf{t}}$ is traction. Note that the boundary of the domain is divided as $\Gamma = \Gamma_p \cap \Gamma_q = \Gamma_u \cap \Gamma_t$, where $\Gamma_p$ and $\Gamma_q$ are the prescribed fluid pressure and Darcy velocity boundaries, and $\Gamma_u$ and $\Gamma_t$ are the prescribed solid displacement and traction boundaries. (W) is then linearized using Gâteaux's differential to obtain its residual form and the tangent operator. Then, these equations are spatially discretized using a finite element scheme. Due to the solid-fluid coupling, continuous Lagrange elements of second and first-order are used for $\mathbf{u}$ and $\mathbf{p}$, respectively. This combination of elements, usually referred as a Taylor-Hood element, fulfills the inf-sup condition and therefore, is numerically stable (Markert, 2008).

At every time step, we obtain spatial fields for the fluid pressure, solid displacement, and accumulated plastic strain. Then, using Eqs. 2 and 6 we calculate the Darcy velocity vector and the stress tensor. Our numerical implementation of coupled poroelastic and elasto-plastic behaviors is benchmarked with two standard tests, the Cryer's sphere problem (Cryer, 1963) and the thick-walled hollow cylinder loading test (Nadai, 1950). Further details on this implementation as well as the benchmarks are available in our GitHub repository: https://github.com/FNSL1996/PEP_FEniCS.

## 3.3. Models setup and tests

We consider a spatial domain that geometrically and kinematically resembles the Planchón-Peteroa geothermal system, according to Pearce et al. (2020). The main structures we refer to are shown in Figure 1. We consider the following domains:

a) A parallelepiped domain that encompasses the system's location. Its dimensions are $50 \times 50 \times 12 \, \text{km}^3$ so as to reduce boundary effects on the results obtained near the geothermal reservoir and the fault zone. The total number of elements reaches 61711.

b) A fault zone intended to kinematically mimic the strike-slip component of the active blind ATF, inferred from the WNW-aligned Cls1 seismicity cluster (Pearce et al., 2020). This is represented using a parallelepiped of dimensions $10 \times 5 \times 0.5 \, \text{km}^3$ (length, depth, width), meshed with a 250 m resolution. A 2-meters left-lateral slip is applied on the northern side of the fault to represent the ATF kinematics, once an initial fluid flow steady state is reached.

c) An ellipsoidal hollow cavity that represents the inferred geothermal reservoir from the C1 resistivity anomaly (Pearce et al., 2020). Its semiaxes dimensions are $4 \times 2 \times 2 \, \text{km}^3$. A 5 MPa fluid pressure is imposed on its surface and pointing outwards, to represent the presence of fluids in the reservoir. This boundary condition is applied as a ramp function until reaching stationary state.

The simulations setup is shown in Figure 2. The mesh is built using Gmsh (https://www.gmsh.info). Visualization of the results is achieved using Paraview 5.9 (https://www.paraview.org).
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Model setup in (a) 3D view, (b) plane view at 5 km depth, and (c) vertical section. The geothermal reservoir is represented by the red ellipsoid and the crustal fault by the light brown parallelepiped. All sides of the domain are free-slip with zero fluid flux conditions, except for the top face, which is set free. A 5 MPa fluid pressure is imposed on the inner surface of the geothermal reservoir. After the stationary state is reached, a 2-meters left-lateral strike-slip motion is applied on the northern face of the fault (indicated with a blue arrow). This slip is applied on a so-called slip patch (blue ellipse shown in a). Outside of this patch, slip motion is accommodated according to the mechanical properties of the fault and bedrock. Additional figures for the slip profile are available in the repository.</figcaption></figure>

We investigate the fluid flow, bedrock strain and stress at a kilometric scale, and as a first-order approximation, we define our medium as a homogeneous domain. Since our crustal rock mass is made of highly fractured rocks of various
sedimentary and volcanic origin, we assumed the mechanical and hydraulic properties of a Berea Sandstone, which has been thoroughly characterized experimentally (Cheng, 2016, Table 1). While the importance of heterogeneous bedrock properties on fluid flow in fractured rock media has long been recognised, we aim first at quantifying the primary role of basic parameters within the specific setting that we designed. For the same reason, depth-dependent variations in permeability and fluid viscosity were also neglected here. These basic parameters are the fault slip-rate (SR), the fault zone shear modulus (G) and permeability (κ), the bedrock’s yield strength (σ<sub>y0</sub>) and the fluid viscosity (μ). We selected 8 model cases that are exposed below. The mechanical and hydraulic properties are listed in Table 2, as well as the resulting maximum normalized fluid flux (NF<sup>MAX</sup>, cf. Figure 5f).

240

<caption>**Table 1**</caption>

Mechanical and hydraulic parameters (after Cheng, 2016). From left to right: G is the shear modulus, v the Poisson's modulus, K the rock bulk modulus, K<sub>s</sub> the solid bulk modulus, K<sub>f</sub> the fluid bulk modulus, α Biot's constant, φ rock porosity, κ rock permeability, and M Biot's modulus.

<table><caption>**Table 2**</caption><thead><tr><th>G [GPa]</th><th>ν [-]</th><th>K [GPa]</th><th>K<sub>s</sub> [GPa]</th><th>K<sub>f</sub> [GPa]</th><th>α [-]</th><th>φ [-]</th><th>κ [m<sup>2</sup>]</th><th>M [GPa]</th></tr></thead><tbody><tr><td>6</td><td>0.2</td><td>8</td><td>36</td><td>2.25</td><td>0.778</td><td>0.19</td><td>1.90·10<sup>-13</sup></td><td>9.92</td></tr></tbody></table>

Main properties and results obtained from the model runs. P-Elastic and P-EPlas correspond to the poro-elastic and the poro-elasto-plastic rheologies, respectively. "Homogeneous" means that identical properties are used for the bedrock and the fault zone. Bold parameters are those that are modified with respect to the reference model case SR_1. The next columns display the chosen slip-rate imposed along the fault, G<sub>f</sub> its shear modulus, κ<sub>f</sub> its permeability. σ<sub>y0</sub> is the yield strength of both the bedrock and fault (it does not apply in case of Elastic rheology, marked with symbol '*'), μ is the fluid viscosity. The shear modulus and the permeability of the bedrock are the same in all simulations: G<sub>b</sub> = 6 [GPa], and κ<sub>b</sub> = 1.9·10<sup>-13</sup> [m<sup>2</sup>], respectively. The last column displays a key result: NF<sup>MAX</sup> is the maximum normalized flux measured in the area defined in Figure 5f). See text for details.

<table><thead><tr><th>Model Name</th><th>Rheology and Heterogeneity</th><th>Slip-rate [m/s]</th><th>G<sub>f</sub> [GPa]</th><th>κ<sub>f</sub> [m<sup>2</sup>]</th><th>σ<sub>y0</sub> [MPa]</th><th>μ [Pa·s]</th><th>NF<sup>MAX</sup> [-]</th></tr></thead><tbody><tr><td>Elastic_SR_1</td><td><b>P-Elastic</b>, Homogeneous</td><td><b>0.1</b></td><td>6</td><td>1.90·10<sup>-13</sup></td><td>*</td><td>10<sup>-3</sup></td><td>~10</td></tr><tr><td>SR_01</td><td>P-EPlas, Homogeneous</td><td><b>0.1</b></td><td>6</td><td>1.90·10<sup>-13</sup></td><td>5</td><td>10<sup>-3</sup></td><td>~10</td></tr><tr><td><b>SR_1</b></td><td>P-EPlas, Homogeneous</td><td>1</td><td>6</td><td>1.90·10<sup>-13</sup></td><td>5</td><td>10<sup>-3</sup></td><td>~12</td></tr><tr><td>SR_10</td><td>P-EPlas, Homogeneous</td><td><b>10</b></td><td>6</td><td>1.90·10<sup>-13</sup></td><td>5</td><td>10<sup>-3</sup></td><td>~12</td></tr><tr><td>Yield_SR_1</td><td>P-EPlas, Homogeneous</td><td>1</td><td>6</td><td>1.90·10<sup>-13</sup></td><td><b>2</b></td><td>10<sup>-3</sup></td><td>~11</td></tr><tr><td>Shear_SR_1</td><td>P-EPlas, <b>Compliant Fault</b></td><td>1</td><td><b>3</b></td><td>1.90·10<sup>-13</sup></td><td>5</td><td>10<sup>-3</sup></td><td>~8</td></tr><tr><td>Permeability_SR_1</td><td>P-EPlas, <b>Permeable Fault</b></td><td>1</td><td>6</td><td><b>4.86·10<sup>-10</sup></b></td><td>5</td><td>10<sup>-3</sup></td><td>~15</td></tr><tr><td>Viscosity_SR_1</td><td>P-EPlas, Homogeneous</td><td>1</td><td>6</td><td>1.90·10<sup>-13</sup></td><td>5</td><td><b>10<sup>-2</sup></b></td><td>~75</td></tr></tbody></table>

241

242

243

244

245

246

247

In all models, a fluid pressure of 5 MPa was imposed on the inner surface of the reservoir to simulate fluid flow coming out from this structure. This pressure magnitude was chosen because it stands in the range of rock masses tensile strength between 1-10 MPa (Etheridge, 1983), and it corresponds to experimental estimates (Duboeuf et al., 2021), consistent with the idea that fluid overpressure drives opening pathways. Once the stationary state is reached, a left-lateral slip of 2 meters was imposed on the northern side of the fault, at rates between 0.1 - 10 m/s, consistent with seismic slip-rates (e.g. Kanamori, 2004). Note that, because of the fault's thickness, slip-rates are not expected to have a considerable impact on bedrock stress.

248

## 4. Results

The following results describe the models evolution from the stationary state up to 31 days after fault slip was imposed. Firstly, we describe the solid deformation and fluid flow behavior for the reference model SR_1 (Sections 4.1 and 4.2). Then we compare the other simulations by quantifying the influence of varying each tested parameter.

249

250

251
## 4.1. Solid deformation and stress state

We compute the mean stress $\bar{\sigma}$ and the stress shape ratio $\phi$ (Eq. 10), which give insights on the rock's confinement and on the relative magnitude between principal stresses, respectively.

$$ \bar{\sigma} = \frac{1}{3} \text{trace}(\sigma(u,p)), \quad \text{and} \quad \phi = \frac{\sigma_2 - \sigma_3}{\sigma_1 - \sigma_3}, \qquad (10) $$

where $\sigma_1$, $\sigma_2$, and $\sigma_3$ are the most compressive, intermediate, and least compressive principal stresses, respectively. $\phi$ is defined between 0 and 1, and addresses the stress tensor "shape" by comparing the magnitude of the principal stresses. A value of $\phi = 0.0$ implies that $\sigma_2 = \sigma_3$, therefore a preferential compressive direction, similar to a uniaxial compression test. If $\phi = 1.0$ then $\sigma_1 = \sigma_2$, which establishes a preferential tensile or least compressive direction, similar to a uniaxial tensile test. Hereafter we use $\phi$ to refer to the following stress states: $\phi = 0$: uniaxial compression, $\phi \in (0, 0.5)$: transpression, $\phi = 0.5$: strike-slip regime, $\phi \in (0.5, 1)$: transtension, and $\phi = 1$: uniaxial tension.

Figure 3 shows the evolution of the solid's volumetric strain, mean stress and stress shape ratio, for simulation SR_1. Overall, the 8 selected models show results similar to those illustrated in this reference model:

1. Regions of compressive and tensile mean stress of approximately $\pm 1$ MPa appear after fault motion at the western and eastern half of the northern side of the fault, respectively. This magnitude depends on our choice of parameters (mainly the shear modulus and the amount of applied slip motion).

2. Contractional and dilational domains of volumetric strain $\sim \pm 0.0001$ mostly coincide with the compressive and tensile regions, respectively. A smaller dilational domain is observed in the western tip of the fault. This asymmetry at each side of the fault responds to the fault slip being applied on the northern side of the fault, where largest intensities occur. This affects fluid flow patterns, as discussed below.

3. Stress rotation occurs around the fault after slip: this is represented with orientations of $\sigma_1$ in Figure 3c, which go from radial relative to the reservoir before fault slip, to fault-parallel and fault-normal in the western and eastern half of the northern side of the fault, respectively. Transtensive domains, represented by dark grey colours develop along the north-eastern half of the fault zone and at its western termination. A smaller transpressive domain, represented by light grey colors, appears at the western half of the fault.

4. Mean stress, stress shape ratio, and volumetric strain do not show significant spatial variation after full fault slip was reached, whatever the applied slip-rate.

## 4.2. Fluid flow and pressure

Figure 4 shows the spatial and temporal evolution of fluid flow for reference model SR_1 through the fluid pressure field and Darcy velocity, represented using streamlines. Similarities are observed throughout the models:

1. A stationary fluid flow from the reservoir is achieved before the application of fault slip (Figure 4a, steady state fluid velocities reach $10^{-12}$ m/s near the fault in SR_1). Fluid migrates to the top surface due to its free flow condition.

2. After fault slip, positive and negative fluid pressure domains appear at the western and eastern half of the northern side of the fault, respectively. A smaller domain of negative fluid pressure is seen at the western tip of the fault, occupying a similar position to the contractional and dilational regions (Figure 3), showing the same asymmetry.

3. Fluid flow from the cavity was immediately altered after fault slip, migrating towards the eastern negative fluid pressure region and, to a lesser extent, towards the negative pressure domain in the western tip. The fluid's Darcy velocity increased to $\sim 10^{-6} - 10^{-7}$ m/s, at the eastern and western negative fluid pressure domain (Figure 4b-c).

4. Fluid pressure regions caused by fault slip start to dissipate one day after the applied slip, returning fluid flow to its stationary state.
Model SR_1

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Temporal sequence of solid stress and strain in reference model SR_1. Columns from the left to right display volumetric strain, mean stress, and stress shape ratio ($\phi$) in planar view at 5 km depth. Plastic envelope is represented as a yellow contour in the middle column, and $\sigma_1$ directions are shown in the right column as light blue segments. The grid displays 5 km-wide segments. At the stationary state before fault slip, an overall dilational stress state is seen around the reservoir. After fault slip, dilational and contractional domains develop around the fault, cf. panel (b). These domains remain almost unchanged throughout the rest of the simulation. Transtensive domains appear at the eastern half of the northern side of the fault, and NW to the western tip of the fault. A transpressive domain appears along the western half of the fault. $\sigma_1$ exhibits rotation after fault slip, especially at the northern side of the fault, from radial relative to the reservoir before fault slip to fault-parallel and fault-normal in the western and eastern half of the fault, respectively.</figcaption></figure>
Model SR_1

<figure><img src="image_6.png" /><figcaption>**Figure 4a:** Temporal sequence of fluid pressure and Darcy velocity (model SR_1) in 3D view (left column), and in a planar view at 4 km depth (right column). The grid displays 5 km-side squares. Thresholds have been applied in the left column to the fluid pressure magnitudes. Stationary state (a), where fluid migrates directly to the free surface. After fault slip, fluid flow migrate into the negative fluid pressure domains (b). As time passes, fluid pressure dissipates (c), reducing fluid migration to the westernmost negative pressure domain (continues on the next page).</figcaption></figure>
Model SR_1

<figure><img src="image_7.png" /><figcaption>**Figure 4b:** (continued): Fluid migration reduces to the westernmost negative pressure domain a day after fault slip (d), then to the easternmost negative pressure domain (e), then returning to the stationary state 30 days after fault slip (f).</figcaption></figure>

**4.3. Fluid flux variation over time as a function of varying parameters**

Despite the similarities, fluid flow from the reservoir varies depending on mechanical and hydraulic properties of the rock and fluid. To quantify them in terms of time and amount of mobilized fluid, a fluid flux calculation was conducted. Fluid flux was measured in an area of 2 × 3 km² located between the cavity and the eastern negative fluid pressure region (shown in blue in Figure 5f). A normalized fluid flux (NF), taking the stationary flux as unit, was then
295

computed over time (Figure 5a-e).

296

297

298

299

300

All models show four typical stages: (I) a continuous increase in NF by the end of fault slip, (II) a *plateau* stage where the NF remains relatively constant over time, (III) a local or global maximum (depending on the model's parameters), and (IV) an asymptotic return to the stationary flux. A summary of these stages is shown Figure 7c.

<figure><img src="image_8.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>f</figcaption></figure>

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

**Figure 5:** (a-e) Normalized flux variation over time for all simulations, to highlight the influence of specific parameters (listed in Table 2). (a) SR_1, SR_01, and SR_10: Slip-rate variation, (b) Elastic_SR_1, SR_1, and Yield_SR_1: influence of Plasticity. (c) SR_1 and Shear_SR_1: fault shear modulus contrast; (d) SR_1 and Permeability_SR_1: fault permeability contrast; (e) SR_1 and Viscosity_SR_1: Fluid viscosity variation. Note that the x-axes use logarithmic scale. (f) Area of maximum Darcy velocity in the normal direction to the fault plane. This region is 2x3 km² and is used for the flux calculation.

301

302

303

304

305

306

307

308

309

### 4.3.1. Effect of fault slip-rate

Figure 5a illustrates the effect of fault slip-rate on the NF, over a range between 0.1 - 10 m/s. A high fault slip-rate triggers a faster response and an increased NF, consistent with the slip duration (20, 2, and 0.2 seconds for the increasing slip-rate). However, the reduction rate is of the order of ~ 6 seconds independently of the slip-rate.

In terms of maximum NF, both SR_1 and SR_10 achieve 12 times the steady state flux at the end of stage I, whereas SR_01 reaches only 10 times the steady flux. This lower NF results from a lower volumetric strain-rate; because of the coupling on Eq. 1, the fluid pressure gradient is also lower, reducing fluid flux. All three models show the same behavior after entering stage II.
### 4.3.2. Influence of plasticity

Figure 5b displays the fluid flux evolution for the poroelastic model Elastic_SR_1 and the poro-elasto-plastic models SR_1 and Yield_SR_1. Their main difference is that Elastic_SR_1 does not display a maximum flux at the end of fault slip. We interpret this as a result from the upper bound effect of plasticity on the stress field. An elasto-plastic domain cannot increase the stress magnitudes infinitely upon fault slip, thus inducing sharp local deformation patterns (dilation), resulting in higher local pressure gradients and enhanced NF, compared to an elastic domain.

The NF is higher in model SR_1 than in Yield_1. This counter-intuitive effect results from less stress contrast around the fault in lower yield strength models, resulting in smaller pressure gradients (because of the coupling between volumetric strain and fluid pressure, Eqs. 1, 2, and 3). Additional tests are available on the repository.

### 4.3.3. Influence of Fault Shear modulus

A decreased fault shear modulus causes generally smaller NF, as shown in Figure 5c, compared to an homogeneous medium: maximum NF is 8 and 12, respectively. We interpret this as a concentration of strain on the northern wall of the fault, where slip is applied, because of its more compliant nature. Therefore, the negative pressure domain is skewed towards that wall, compared to the reference model, causing the pressure gradient on the southern side of the fault to be smaller, reducing the flux.

Shear_SR_1 shows a global maximum at stage III, unlike SR_1 that displays it at stage I. We propose that, as the negative fluid pressure dissipates faster inside the fault zone ~ 15 mins after slip, it creates a temporary higher pressure gradient also in the surrounding bedrock, which increases fluid flux until ~ 6 hours after slip. This gradient then decreases, and the NF returns to stationary state.

### 4.3.4. Influence of Fault permeability

An increased fault permeability, mimicking fractured fault zones as seen in the SAVZ (Sepúlveda et al., 2020), produces a higher NF at the end of fault slip compared to the reference model (14.5 vs. 12, Figure 5d). Unlike other models, Permeability_SR_1 does not display a plateau stage. Instead it decreases slightly after fault slip, then quickly re-increases to a global maximum. The fault's high permeability causes faster fluid pressure dissipation after fault slip.

The revealing aspect of this model is the flow inversion occurring about 1 hour after fault slip. Because of the fault's higher permeability, fluid pressure dissipates more quickly inside the fault zone than in the bedrock. This implies that the fluid pressure inside the fault remains higher in its southern part than in the surrounding bedrock, between 1-10 hours after slip, inducing fluid to flow from inside the fault to the outside. Afterwards, fluid pressure dissipates in the bedrock, and fluid flow recovers its original direction. Figure 6 illustrates this process. Note that this reverse fluid flow is local, depending on the area where the flux calculation is done. A more detailed study may establish the spatial distribution of such flow reversal.

### 4.3.5. Influence of Fluid viscosity

A fluid viscosity an order of magnitude greater than in our reference model shows a similar behavior, but much higher NF at all stages (Figure 5e). However, we note that the steady state flux prior to the application of fault slip motion is almost two orders of magnitude lower than in the reference model (0.003 $m^3/s$ vs 0.202 $m^3/s$). Therefore, we suggest that rather than implying that fault slip causes higher fluxes when dealing with high viscosity fluids, fault slip causes a bigger disturbance on the fluid flux, relative to its stationary state. This effect lasts longer in model Viscosity_SR_1, since even 31 days after fault slip, NF is still ~ 12 times that of the stationary level. We interpret that, because of their higher viscosity, fluids take longer to return to the stationary state.

## 5. Discussion

### 5.1. Numerical and Mechanical validation

To assess fluid flow and bedrock strain evolution, we use a continuous approach instead of Extended Finite Element Method (XFEM) or Discrete Fracture Network (DFN) formulations in which fractures are modeled (e.g. McClure and Horne, 2013; Wang, 2015). Because fluid flow in geothermal systems is heavily controlled by the presence of fracture networks (Sanderson and Nixon, 2015), advantages and limitations of our modelling approach exist when comparing
<figure><img src="image_14.png" /><figcaption>**Figure 6:** Comparison of the evolution of fluid pressure in model SR _1 (left column) and Permeability _SR _1 (right column), at different moments.**
a) Pressure isocontours indicate a bigger gradient in the reference model. Both gradients go from the northern side of the fault towards the reservoir, producing fluid flow onto the fault. b) The right column shows an increase in fluid pressure inside the fault, due to faster dissipation. This inner-fault fluid pressure exceeds that in the bedrock (highlighted area), provoking a reversal in the pressure gradient and, therefore, in fluid flow. c) The negative pressure domain imposed by the fault slip starts to dissipate as well, reestablishing the original gradient.</figcaption></figure>

it against those approaches. XFEM is able to represent crack propagation and solve a fully-coupled solid-fluid medium (Wang, 2015). However, their development is still, to our knowledge, limited to small domains, with high numerical cost. Similarly, DFN three-dimensional models still need to be developed (McClure and Horne, 2013). Even though traditional FEM is not able to represent fractures, its robustness allow to tackle large and complex 3D geometries, at a lower computational cost. Furthermore, our results show that a FEM-based poro-elasto-plastic model allows to assess
bulk fluid flow and strain at a kilometric scale.

Because of the explicit solid-fluid coupling on Eq. 1, plastic strain affects fluid flow, which is confirmed by the NF evolution differences between a poroelastic and poro-elasto-plastic model in Figure 5b. Although the von Mises criterion simplifies greatly the yield behavior of crustal rocks, it has been widely used in mechanical models (e.g. Tapponnier and Molnar, 1976; Iturrieta et al., 2017, and references therein). We also recall studies indicating that, at the transition from pressure-dependent behavior to ductile behavior from 5 to 15 km depth, yield strength displays a pressure-independent trend (e.g. Ord and Hobbs, 1989; Suppe, 2014). A key result from our study is that this simple yield criterion already impacts fluid flow in stages I and II due to local deformation patterns, generating greater fluid flow. This impact should be enhanced when considering pressure-sensitive plasticity. Future models accounting for pressure-dependent yield will further unravel the interplay between brittle rock failure and fluid flow.

## 5.2. Suction pump-like mechanism

When contrasting our results against geological studies on fluid migration, we find common key points. For instance, when considering fluid migration from the geothermal reservoir into the bedrock, fluid migrates to regions that experienced an increase in tensile mean stress, as previously postulated by Nur and Booker (1972): here, the emergence of a tensile domain after slip is seen on the eastern half of the fault (Figure 3b), where the largest influx of fluids from the reservoir occurs (Figure 4c). This domain coincides with a dilational region, which can be a proxy of dilational wing cracks, as seen in nature (e.g. Schultz, 2000). Such wing-like fractures are observed in many field areas (e.g. Kim and Sanderson, 2006). For example, in the Tatara San Pedro area, Ruz-Ginouves et al. (2021) identified with numerical simulations a dilational region to the NE of the Melado fault that can be interpreted as a wing crack termination. This supports structural data from Sielfeld et al. (2019) who also identified a graben structure.

Furthermore, our results lead us to postulate that the first-order control on fluid flow, rock stress and strain is that of a suction pump-like mechanism (Sibson, 1985, 1987, 2000). Figure 7a shows the fluid pressure evolution at dilational jogs subject to suction pumping. A rapid decrease due to fault slip is followed by a progressive recovery to pre-slip levels. Figure 7b displays the evolution of fluid pressure, normalized to the stationary state, at the eastern negative fluid pressure domain indicated in Figure 4c. Our modeled fluid pressure follows the same behavior as the suction pump mechanism: a sudden decrease upon fault slip followed by a progressive return to pre-slip levels. Furthermore, our results place quantitative time constraints attesting that this mechanism is active at time scales ranging from a day up to a month (considering hydrothermal fluid). An illustration of the key stages of this mechanism is shown in Figure 8.

Note also that the asymmetry of the fluid pressure domains (cf. Figure 8), depends on the wall of the fault onto which slip is applied. For instance, the northern side domains are expected to shift to the southern side if slip was applied on the southern wall. This effect will have to be further explored in future studies. Other factors not taken into account here such as permeability evolution, and lithostatic pressure, among others, will also control the duration of this mechanism, and require future work. Variations in permeability with increasing depth, temperature and composition are known to influence the fluid flow via the Darcy relationship. Fracture sealing in turn, has been associated with a reduction of up to three orders of magnitude in matrix bulk permeability due to hydrothermal silica precipitation, thus reducing fluid flow magnitudes and durations (Dobson et al., 2003).

## 5.3. Fault reactivation, aftershock seismicity, and implications for geothermics

Our results show that a suction pump-like mechanism causes fluid influx from the geothermal cavity into the easternmost dilational domain after fault slip. Flux variation in the selected region of interest (Figure 7c) shows that a NF 8-15 times the stationary flux sustains up to 31 days after fault slip. Two of our models, however behave distinctively, and can be explained with the concept of hydraulic diffusivity $\omega$:

$$ \omega = \frac{\kappa}{\phi\mu\beta}, \qquad (11) $$

where $\phi$ is the rock porosity, and $\beta$ is the matrix compressibility (e.g. Reid, 2004, and references therein). The higher $\omega$, the faster fluid pressure dissipates. Hence, as fluid viscosity increases, greater fluid pressure (hence fluid flow) are sustained over longer time, as seen in Viscosity_SR_1. The opposite occurs in Permeability_SR_1, where
<figure><img src="image_15.png" /></figure>

**Figure 7:** (a) Suction pump cyclic mechanism at a dilational jog, showing a coseismic reduction in fluid pressure ($P_f$) below hydrostatic pressure ($P_h$). A progressive increase to hydrostatic pressure occurs throughout the interseismic period (modified from Sibson, 2000). (b) Fluid pressure evolution at the eastern negative fluid pressure domain in our models (Figure 4). Fluid pressure is normalized to the stationary state. A rapid decrease in fluid pressure is seen during fault slip, followed by a progressive increase back towards stationary level, 30 days after fault slip, showing remarkable resemblance to the suction pump mechanism. (c) Stages of fluid flux temporal evolution at the area of interest (Figure 5f). Stage I consists in a continuous increase of the normalized flux from the onset to the end of fault slip. Stage II shows its subsequent decrease to a constant, plateau value. Stage III shows a second increase of the normalized flux, slower than Stage I. Last, Stage IV shows an asymptotic decrease back to the stationary flux. Flux duration and peak vary according to the model parameters. Note the logarithmic scale used for the x-axis.

fluid flow reduces to levels below the stationary flux just an hour after fault slip.

Due to the time period over which fluid flux is sustained, we infer that, in nature, hydraulic fracturing would enable the opening of new pathways for fluid migration, causing a local increase in rock permeability, and generating aftershock seismicity, as has already been proposed (e.g. Nur and Booker, 1972; Bosl and Nur, 2002). Aftershock seismicity due to fluid migration has been inferred in settings with fluids of different viscosity, ranging from CO₂
413

414

415

416

417

418

419

420

421

422

423

424

425

426

(Miller et al., 2004), water, magma, or a mixture of both (Legrand et al., 2011; White and McCausland, 2016). Its duration ranges from weeks to a month with water-viscosity fluids (Miller et al., 2004) and up to six months with water-magma mixtures (Legrand et al., 2011). Our results are consistent with this timing, with a return to stationary levels 31 days after fault slip for water-viscosity fluids. As for higher viscosity fluids, Viscosity_SR_1 shows that a viscosity one order of magnitude higher increases NF levels for longer than a month, which is consistent with an increase in the duration of aftershock seismicity with increasing fluid viscosity.

A sustained fluid flux as we modeled here may also induce seismicity in surrounding faults. Injection monitoring at Enhanced Geothermal Systems (EGS) (e.g. Dorbath et al., 2009; Moeck et al., 2009) has shown that fluid flux can cause seismicity for days up to a month. If sustained fluid flux migrates close to an existent crustal fault system, fault reactivation can occur. This has direct implications for our setting, which displays a high potential for energy generation (Pearce et al., 2020). Apart from the ATF structure that we modeled, this setting is located within the EFFS, which is intersected by additional ATF structures. Fluid migration, both due to AFT slip or geothermal exploitation, could cause fault reactivation on these structures, which may be found difficult to anticipate and control (Dorbath et al., 2009).

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

**Figure 8:** Reservoir fluids trajectories variation over time due to the suction pump mechanism. (1) Starting at a stationary pre-slip state, fluid migrates towards the surface. (2) After the fault slips, negative (blue prisms) and positive (pink prisms) fluid pressure domains emerge because of the motion. Such domains cause the fluids to modify their trajectories, migrating into the negative fluid pressure domains, particularly, into the easternmost one. This state lasts for a few hours after fault slip. (3) As time passes, fluid pressure starts to dissipate, causing the reduction of the negative and positive fluid pressure regions generated due to fault motion. This reduces fluid migration into such domains days after fault slip. Fluid flow trajectories progressively return to the stationary state months after fault slip.
# 6. Conclusions

In this work, we approached the control exerted by a crustal fault slip on a nearby geothermal reservoir through poro-elasto-plastic simulations. Due to the fully-coupled nature of this formulation, fluid flow, bedrock strain and stress were analyzed simultaneously over a period of 31 days after fault slip. Parametric tests allowed us to compare and contrast our results with respect to key parameters. The main conclusions are:

1. Fault slip induces contractional and dilational domains at the eastern and western sides of the fault after it slips, respectively. As a consequence, positive and negative fluid pressure develop in these domains, respectively. Fluid pressure dissipates over time, making those positive and negative pressure areas disappear almost completely after one month (cf. Figure 8).

2. Maximum normalized fluid flux (NF) occurs along the fault's eastern half. NF temporal evolution shows a behavior in 4 stages, displayed Figure 7c. Note that the first three stages occur within approximately the first hour after fault slip.

3. A higher fluid viscosity led to fluid flow lasting longer than a fluid with water-type viscosity. Additionally, an increased fault permeability induced greater NF, returning faster to the stationary state, an hour after fault slip. These results are expected when considering the system's hydraulic diffusivity: a higher viscosity fluid reduces the hydraulic diffusivity, thus promoting long lasting fluid pressure gradients. The opposite effect occurs with a higher permeability domain, as hydraulic diffusivity increases.

4. We show that accounting for plastic yielding locally enhances the fluid flow, even without a pressure-dependent formulation. We also show that it enhances it in-as-much as this yield strength is high enough to allow for sufficient pressure gradients in the bedrock mass, that drive and focus fluid flow.

5. We postulate that the first-order control over fluid flow for a typical fault and geothermal system in the SAVZ is that of a suction pump-like mechanism. The emergence of a dilational domain is consistent with the development of dilational wing cracks in natural examples. EGS injection studies showed that an increase in fluid flux could trigger induced seismicity for ~ 10 days, which is consistent with our results. This good correlation suggests that a suction pump-like mechanism controls both aftershock and induced seismicity in such geothermal settings.

6. To further address fluid flow and bedrock strain within our geotectonic setting, factors such as lithostatic pressure, fracture sealing, depth and temperature-dependent, anisotropic permeability will be accounted for in the next stages of our numerical developments.

# 7. Acknowledgements

This research is supported by Chilean FONDAP Project 15090013 CEGA (Centro de Excelencia en Geotermia de los Andes), and Chilean-French scientific cooperation ECOS-ANID Project 180027/PC18U08. JC acknowledges financial support from FONDECYT Project 1210591, which is funding research on Andean volcanotectonics in the SAVZ. FS thanks IRD (Institut de Recherche pour le Developpement France) for additional travel funding to France. The models were obtained thanks to the High-Performance Cluster of the School of Engineering PUC (https://deg.ing.uc.cl/informatica/cluster/). We acknowledge the significant help by Javier Bisbal with the computation of fluid fluxes, Javier Espinosa with the geological maps figures, and Fernanda Gallegos with the summary figures. Thorough reviews and insightful comments by Susan Ellis, an anonymous reviewer and associate editor Rebecca Bendick helped improving the clarity and quality of the manuscript.
465 References

466 Alghannam, M., Juanes, R., 2020. Understanding rate effects in injection-induced earthquakes. Nature Communications 11(1).

467 Alnæs, M.S., Blechta, J., Hake, J., Johansson, A., Kehlet, B. Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N., 2015. The fenics project version 1.5. Archive of Numerical Software 3(100), 9–23.

469 Angermann, D., Klotz, J., Reigber, C., 1999. Space-geodetic estimation of the nazca-south america euler vector. Earth and Planetary Science Letters 171(3), 329–334.

471 Bellani, S., Brogi, A., Lazzarotto, A., Liotta, D., Ranalli, G., 2004. Heat flow, deep temperatures and extensional structures in the larderello geothermal field (italy): constraints on geothermal fluid flow. Journal of Volcanology and Geothermal Research 132(1), 15–29.

473 Bosl, W.J., Nur, A., 2002. Aftershocks and pore fluid diffusion following the 1992 landers earthquake. Journal of Geophysical Research: Solid Earth 107(B12).

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

476 Cembrano, J., Lara, L., 2009. The link between volcanism and tectonics in the southern volcanic zone of the chilean andes: A review. Tectonophysics 471(1-2), 96–113.

478 Cheng, A.H., 2016. Poroelasticity (Theory and Applications of Transport in Porous Media: Volume 27). Springer.

479 Coussy, O., 2004. Poromechanics. Wiley.

480 Cryer, C.W., 1963. A comparison of the three-dimensional consolidation theories of Biot and Terzaghi. The Quarterly Journal of Mechanics and Applied Mathematics 16, 401–412.

482 Dempsey, D.E., Archer, R.A., Ellis, S.M., Rowland, J.V., 2013. Hydrological effects of dip-slip fault rupture on a hydrothermal plume. Journal of Geophysical Research: Solid Earth 118(1).

484 Dobson, P.F., Kneafsey, T.J., Hulen, J., Simmons, A., 2003. Porosity, permeability, and fluid flow in the Yellowstone geothermal system, Wyoming. Journal of Volcanology and Geothermal Research 123, 313–324.

486 Dorbath, L., Cuenot, N., Genter, A., Frogneux, M., 2009. Seismic response of the fractured and faulted granite of soultz-sous-forets (france) to 5 km deep massive water injections. Geophysical Journal International 177(2), 653–675.

488 Duboeuf, L., De Barros, L., Kakurina, M., Guglielmi, Y., Cappa, F., Valley, B., 2021. Aseismic deformations perturb the stress state and trigger induced seismicity during injection experiments. Geophysical Journal International 224, 1464–1475.

490 Duwiquet, H., Guillou-Frottier, L., Arbaret, L., Bellanger, M., Guillon, T., Heap, M.J., 2021. Crustal fault zones (cfz) as geothermal power systems: A preliminary 3d thm model constrained by a multidisciplinary approach. Geofluids 2021.

492 Etheridge, M.A., 1983. Differential stress magnitudes during regional deformation and metamorphism: Upper bound imposed by tensile fracturing. Geology 11(4).

494 Gerbault, M., Poliakov, A.N., Daignieres, M., 1998. Prediction of faulting from the theories of elasticity and plasticity: what are the limits? Journal of Structural Geology 20, 301–320.

496 Gudmundsson, A., 2011. Rock Fractures in Geological Processes. Cambridge University Press.

497 Haagenson, R., Rajaram, H., Allen, J., 2020. A generalized poroelastic model using fenics with insights into the noordbergum effect. Computers & Geosciences 135.

499 Heiland, J., 2003. Permeability of triaxially compressed sandstone: Influence of deformation and strain-rate on permeability. Pure and Applied Geophysics 160(5-6), 889–908.

501 Hickson, C.J., Ferraris, F., Rodriquez, C., Sielfeld, G., Henriquez, R., Gislason, T., Selters, J., Benoit, D., White, P., Southon, J., Ussher, G., Charroy, J., Smith, A., Lovelock, B., Lawless, J., Quinlivan, P., Smith, L., Yehia, R., 2011. The mariposa geothermal system, chile. Geothermal Resources Council Transactions (GRC) 35, 817–825.

504 Iturrieta, P.C., Hurtado, D.E., Cembrano, J., Stanton-Yonge, A., 2017. States of stress and slip partitioning in a continental scale strike-slip duplex: Tectonic and magmatic implications by means of finite element modeling. Earth and Planetary Science Letters 473, 71–82.

506 Kanamori, H., 2004. The diversity of the physics of earthquakes. Proceedings of the Japan Academy, Series B 80, 297–316.

507 Kavanagh, J.L., Engwell, S.L., Martin, S.A., 2018. A review of laboratory and numerical modelling in volcanology. Solid Earth 9(2), 531–571.

508 Kim, J., Tchelepi, H.A., Juanes, R., 2011. Stability, accuracy and efficiency of sequential methods for coupled flow and geomechanics. SPE Journal 16(2), 249–262.

510 Kim, Y.S., Sanderson, D.J., 2006. Structural similarity and variety at the tips in a wide range of strike-slip faults: a review. Terra Nova 18(5), 330–344.

512 Kissling, W., Villamor, P., Ellis, S., Rae, A., 2018. Modelling of hydrothermal fluid flow and structural architecture in an extensional basin, ngakuru graben, taupo rift, new zealand. Journal of Volcanology and Geothermal Research 357, 134–151.

514 Legrand, D., Barrientos, S., Bataille, K., Cembrano, J., Pavez, A., 2011. The fluid-driven tectonic swarm of ayson fjord, chile (2007) associated with two earthquakes (mw=6.1 and mw=6.2) within the liquiñe-ofqui fault zone. Continental Shelf Research 31(3-4), 154–161.

515 Liu, M., Huang, H., 2021. Finite element modeling of spherical indentation in a poro-elasto-plastic medium via step displacement loading. International Journal for Numerical and Analytical Methods in Geomechanics 45(10).

518 Manga, M., Wang, C.Y., 2015. Earthquake Hydrology. Treatise on Geophysics, 305–328.

519 Markert, B., 2008. A biphasic continuum approach for viscoelastic high-porosity foams: Comprehensive theory, numerics, and application. Archives of Computational Methods in Engineering 15(4), 371–446.

521 McClure, M.W., Horne, R.N., 2013. Discrete Fracture Network Modeling of Hidraulic Stimulation. Springer Publishing.

522 Melosh, G., Cumming, W., Benoit, D., Wilmarth, M., Colvin, A., Winick, J., Soto-Neira, E., Sussman, D., Urzúa-Monsalve, L., Powell, T., Peretz, A., 2010. Exploration results and resource conceptual model of the tolhuaca geothermal field, chile. Proceedings World Geothermal Congress 2010.

525 Métois, M., Socquet, A., Vigny, C., 2012. Interseismic coupling, segmentation and mechanical behavior of the central chile subduction zone. Geophysical Research: Solid Earth 117(B3).

526
527 Miller, S.A., Collettini, C., Chiaraluce, L., Cocco, M., Barchi, M., Kaus, B.J.P., 2004. Aftershocks driven by a high-pressure co2 source at depth. Nature 427(6976), 724-727.

528 Mitchell, T.M., Faulkner, D.R., 2008. Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones. Journal of Geophysical Research 113(B11).

531 Moeck, I., Kwiatek, G., Zimmermann, G., 2009. Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir. Journal of Structural Geology 31(10), 1174-1182.

533 Muñoz, M., Garat, P., Flores-Aqueveque, V., Vargas, G., Rebolledo, S., Sepúlveda, S., Daniele, L., Morata, D., Parada, M.N., 2015. Estimating low-enthalpy geothermal energy potential for district heating in santiago basin-chile (33.5 °s). Renewable Energy 76, 186-195.

535 Nadai, A., 1950. Theory of Flow and Fracture of Solid. McGraw-Hill. New York.

536 Nur, A., Booker, J.R., 1972. Aftershocks caused by pore fluid flow? Science 175(4024), 885-887.

537 Ord, A., Hobbs, B., 1989. The strength of the continental crust, detachment zones and the development of plastic instabilities. Tectonophysics 158, 269-289.

539 Pandey, S., Vishal, V., Chaudhuri, A., 2018. Geothermal reservoir modeling in a coupled thermo-hydro-mechanical-chemical approach: A review. Earth-Science Reviews 185, 1157-1169.

541 Pavez, M., Schill, E., Held, S., Díaz, D., Kohl, T., 2020. Visualizing preferential magmatic and geothermal fluid pathways via electric conductivity at villarrica volcano, s-chile. Journal of Volcanology and Geothermal Research 400, 106913.

543 Pearce, R.K., Sánchez De La Muela, A., Moorkamp, M., Hammond, J., Mitchell, T.M., Cembrano, J., Araya Vargas, J., Meredith, P.G., Iturrieta, P., Pérez-Estay, N., Marshall, N.R., Smith, J. Yañez, G., Ashley Griffith, W., Marquardt, C., Stanton-Yonge, A., Nuñez, R., 2020. Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data. Tectonics 39(12).

546 Prevost, J.H., 2013. One-way versus two-way coupling in reservoir-geomechanical models. Fifth Biot Conference on Poromechanics.

547 Reid, M.E., 2004. Massive collapse of volcano edifices triggered by hydrothermal pressurization. Geology 32(5).

548 Rinaldi, A.P., Rutqvist, J., Luu, K., Blanco-Martín, L., Hu, M., Sentís, M.L., Eberle, L., Kaestli, P., 2022. Tough3-flac3d: a modeling approach for parallel computing of fluid flow and geomechanics. Computational Geosciences, 1-18.

550 Roquer, T., Arancibia, G., Rowland, J., Iturrieta, P., Morata, D., Cembrano, J., 2017. Fault-controlled development of shallow hydrothermal systems: Structural and mineralogical insights from the southern andes. Geothermics 66, 156-173.

551 Rowland, J.V., Simmons, S.F., 2012. Hydrologic, magmatic, and tectonic controls on hydrothermal flow, taupo volcanic zone, new zealand: Implications for the formation of epithermal vein deposits. Economic Geologic 107(3), 427-457.

554 Rutqvist, J., Rinaldi, A.P., Cappa, F., Jeanne, P., Mazzoldi, A., Urpi, L., Guglielmi Y., Vilarrasa, V., 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.

557 Ruz-Ginouves, J., Gerbault, M., Cembrano, J., Iturrieta, P., Sáez Leiva, F., Novoa, C., Hassani, R., 2021. The interplay of a fault zone and a volcanic reservoir from 3d elasto-plastic models: Rheological conditions for mutual trigger based on a field case from the andean southern volcanic zone. Journal of Volcanology and Geothermal Research 418.

560 Rybacki, E., Niu, L., Evans, B., 2021. Semi-brittle deformation of carrara marble: Hardening and twinning induced plasticity. Journal of Geophysical Research: Solid Earth 126, e2021JB022573.

561 Sanderson, D.J., Nixon, C.W., 2015. The use of topology in fracture network characterization. Journal of Structural Geology 72, 55-66.

563 Schultz, R.A., 2000. Growth of geologic fractures into large-strain populations: review of nomenclature, subcritical crack growth, and some implications for rock engineering. International Journal of Rock Mechanics and Mining Sciences 37(1), 403-411.

565 Sepúlveda, J., Arancibia, G., Molina, E., Gilbert, J.P., Duda, M., Browning, J., Roquer, T., Morata, D., Ahrens, B., Bracke, R., 2020. Thermomechanical behavior of a granodiorite from the liquiñe fractured geothermal system (39°s) in the southern volcanic zone of the andes. Geothermics 87, 101828.

568 Sibson, R.H., 1985. Stopping of earthquake ruptures at dilational fault jogs. Nature 316(6025), 248-251.

569 Sibson, R.H., 1987. Earthquakes rupturing as a mineralizing agent in hydrothermal systems. Geology 15(8).

570 Sibson, R.H., 1990. Conditions for fault-valve behaviour. Geological Society, London, Special Publications 54(1), 15-28.

571 Sibson, R.H., 2000. Fluid involvement in normal faulting. Journal of Geodynamics 29(3-5), 469-499.

573 Sielfeld, G., Ruz, J., Brogi, A., Cembrano, J., Stanton-Yonge, A., Pérez-Flores, P., Iturrieta, P., 2019. Oblique-slip tectonics in an active volcanic chain: A case study from the southern andes. Tectonophysics 770, 228221.

574 de Souza Neto, E.A., Perić, D., Owen, D.R.J., 2008. Computational methods for plasticity. Theory and applications. Wiley.

575 Spagnotto, S.L., Triep, E.G., Giambiagi, L.B., Nacif, S.V., Álvarez, O., 2015. New evidences of rupture of crust and mantle in the subducted Nazca plate at intermediate-depth. Journal of South American Earth Sciences 58, 141-147.

577 Suppe, J., 2014. Fluid overpressures and strength of the sedimentary upper crust. Journal of Structural Geology 69, 481-492.

578 Tapponnier, P., Molnar, P., 1976. Slip-line field theory and large-scale continental tectonics. Nature 264, 319-324.

579 Teyssier, C., Tikoff, B., Markley, M., 1995. Oblique plate motion and continental tectonics. Geology 23(5), 447-450.

580 Trutnevyte, E., Ejderyan, O., 2017. Managing geoenergy-induced seismicity with society. Journal of Risk Research 21(10), 1287-1294.

581 Vilarrasa, V., de Simone, S., Carrera, J., Villaseñor, A., 2021. Unraveling the causes of the seismicity induced by underground gas storage at castor, spain. Geophysical Research Letters 48(7).

583 Wang, H., 2015. Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using xfen with cohesive zone method. Journal of Petroleum Science and Engineering 135, 127-140.

585 White, R., McCausland, W., 2016. Volcano-tectonic earthquakes: A new tool for estimating intrusive volumes and forecasting eruptions. Journal of Volcanology and Geothermal Research 309, 139-155.

587 Zhu, W., Wong, T.F., 1997. The transition from brittle faulting to cataclastic flow: Permeability evolution. Journal of Geophysical Research 102(B2), 3027-3041.