Felipe Sáez-Leiva

Daniel E Hurtado

Muriel Gerbault

Javiera Ruz-Ginouves

Pablo Iturrieta

José Cembrano

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

Keywords:

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 -4 [-], 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.

Introduction

Renewable energies such as geothermal power, have experienced exponential growth in their generation capacity over the last decades [START_REF] Pandey | Geothermal reservoir modeling in a coupled thermo-hydro-mechanical-chemical approach: A review[END_REF]. 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. [START_REF] Muñoz | Estimating low-enthalpy geothermal energy potential for district heating in santiago basin-chile (33.5 • s)[END_REF][START_REF] Pandey | Geothermal reservoir modeling in a coupled thermo-hydro-mechanical-chemical approach: A review[END_REF].

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

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 [START_REF] Sibson | Stopping of earthquake ruptures at dilational fault jogs[END_REF][START_REF] Sibson | Earthquakes rupturing as a mineralizing agent in hydrothermal systems[END_REF][START_REF] Sibson | Fluid involvement in normal faulting[END_REF] and the fault-valve [START_REF] Sibson | Conditions for fault-valve behaviour[END_REF][START_REF] Sibson | Fluid involvement in normal faulting[END_REF] 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 [START_REF] Sibson | Earthquakes rupturing as a mineralizing agent in hydrothermal systems[END_REF][START_REF] Manga | Earthquake Hydrology[END_REF]. 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 [START_REF] Caine | Fault zone architecture and permeability structure[END_REF]. 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 [START_REF] Gudmundsson | Rock Fractures in Geological Processes[END_REF]. This increase, however, is transitory because of fracture sealing due to hydrothermal fluid precipitation [START_REF] Dobson | Porosity, permeability, and fluid flow in the Yellowstone geothermal system, Wyoming[END_REF].

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 [START_REF] Heiland | Permeability of triaxially compressed sandstone: Influence of deformation and strain-rate on permeability[END_REF][START_REF] Mitchell | Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones[END_REF]. In contrast, high-porosity rocks show a reduction in permeability at damage onset. [START_REF] Zhu | The transition from brittle faulting to cataclastic flow: Permeability evolution[END_REF] 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 shearenhanced compaction. This porosity-dependent behavior indicates that permeability in faults damage zones is controled by macrofractures [START_REF] Mitchell | Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones[END_REF]. 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 [START_REF] Sepúlveda | Thermomechanical behavior of a granodiorite from the liquiñe fractured geothermal system (39°s) in the southern volcanic zone of the andes[END_REF].

Numerical studies have addressed fluid flow in fault zones through different approaches, amongst which thermohydraulic coupling allows to evaluate convective flow directions and intermittence in both present-day and fossil geothermal systems (e.g. [START_REF] Kissling | Modelling of hydrothermal fluid flow and structural architecture in an extensional basin, ngakuru graben, taupo rift, new zealand[END_REF][START_REF] Duwiquet | Crustal fault zones (cfz) as geothermal power systems: A preliminary 3d thm model constrained by a multidisciplinary approach[END_REF]. 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. [START_REF] Iturrieta | States of stress and slip partitioning in a continental scale strike-slip duplex: Tectonic and magmatic implications by means of finite element modeling[END_REF][START_REF] Kavanagh | A review of laboratory and numerical modelling in volcanology[END_REF][START_REF] Ruz-Ginouves | 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[END_REF]. 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. [START_REF] Dempsey | Hydrological effects of dip-slip fault rupture on a hydrothermal plume[END_REF]. 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 [START_REF] Prevost | One-way versus two-way coupling in reservoir-geomechanical models[END_REF]. Such models have been mostly used to constrain geothermal fluid flow distribution across a fault zone with spatially-variable permeability (e.g. [START_REF] Duwiquet | Crustal fault zones (cfz) as geothermal power systems: A preliminary 3d thm model constrained by a multidisciplinary approach[END_REF], and crustal fault reactivation due to fluid migration (e.g. [START_REF] Alghannam | Understanding rate effects in injection-induced earthquakes[END_REF][START_REF] Vilarrasa | Unraveling the causes of the seismicity induced by underground gas storage at castor, spain[END_REF], CO 2 sequestration (e.g. [START_REF] Rutqvist | Fault activation and induced seismicity in geological carbon storage -lessons learned from recent modeling studies[END_REF] and nuclear waste storage (e.g. [START_REF] Rinaldi | Tough3-flac3d: a modeling approach for parallel computing of fluid flow and geomechanics[END_REF]. 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 [START_REF] Kim | Stability, accuracy and efficiency of sequential methods for coupled flow and geomechanics[END_REF], 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. [START_REF] Rinaldi | Tough3-flac3d: a modeling approach for parallel computing of fluid flow and geomechanics[END_REF]. 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. [START_REF] Haagenson | A generalized poroelastic model using fenics with insights into the noordbergum effect[END_REF][START_REF] Liu | Finite element modeling of spherical indentation in a poro-elasto-plastic medium via step displacement loading[END_REF]. Within this framework, plastic behavior has been accounted for in sequentially coupled approaches (e.g. [START_REF] Rinaldi | Tough3-flac3d: a modeling approach for parallel computing of fluid flow and geomechanics[END_REF], but seldom in fully-coupled approaches [START_REF] Liu | Finite element modeling of spherical indentation in a poro-elasto-plastic medium via step displacement loading[END_REF]. 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. [START_REF] Hickson | The mariposa geothermal system, chile[END_REF], Tolhuaca (e.g. [START_REF] Melosh | Exploration results and resource conceptual model of the tolhuaca geothermal field, chile[END_REF], Villarrica (e.g. [START_REF] Pavez | Visualizing preferential magmatic and geothermal fluid pathways via electric conductivity at villarrica volcano, s-chile[END_REF], and Planchón-Peteroa (e.g. [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF]. Even though further local studies are yet to be conducted, the Planchón-Peteroa inferred geothermal system displays the fundamental traits of an Andean highenthalpy 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 [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF] (Figure 1). We focus on the influence of the seismically active WNWstriking 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.

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.

Poromechanics formulation and governing equations

Poromechanics studies the simultaneous deformation of a solid porous material and pore fluid flow [START_REF] Coussy | Poromechanics[END_REF] through the definition of a Representative Volume Element (RVE) of a saturated porous medium. We adopt Figure 1: (a) Regional scale map of a segment of the Southern Andes and the Southern Andes Volcanic Zone (SAVZ) [START_REF] Cembrano | The link between volcanism and tectonics in the southern volcanic zone of the chilean andes: A review[END_REF]. (b) Local scale map with the main geological structures such as volcanoes and faults, and the results of the geophysical surveys conducted by [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF]. MT anomalies and seismic hypocenters are located at 6 kilometers depth. The convergence vector is taken from [START_REF] Angermann | Space-geodetic estimation of the nazca-south america euler vector[END_REF]; [START_REF] Métois | Interseismic coupling, segmentation and mechanical behavior of the central chile subduction zone[END_REF], and the traces of the faults, as well as the focal mechanisms, are taken from [START_REF] Spagnotto | New evidences of rupture of crust and mantle in the subducted Nazca plate at intermediate-depth[END_REF].

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 𝒖 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:

𝛼 𝜕𝜀 𝑣 (𝒖) 𝜕𝑡 + 1 𝑀 𝜕p 𝜕𝑡 + div (𝒒 (p)) = 0, (1)
where 𝛼 and 𝑀 are Biot's coefficient and modulus respectively, 𝜀 𝑣 (𝒖) = div(𝒖) is the volumetric strain of the solid skeleton, and 𝒒 (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:

𝒒 (p) = - 𝜅 𝜇 ∇p.
(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:

div (𝝈 (𝒖, p)) = 0, (3)
where 𝝈 (𝒖, p) is the Cauchy stress tensor. Because of the solid-fluid interactions, 𝜎 is defined by Eq. 4 in the elastic range:

𝝈 (𝒖, p) = ( 𝐾 - 2𝐺 3
) trace (𝜺 (𝒖)) 𝑰 + 2𝐺𝜺 (𝒖) -𝛼p𝑰, ( 4
)
where 𝐾 is the matrix bulk modulus, 𝐺 the shear modulus, and

𝜺 (𝒖) = 1 2 ( ∇(𝒖) + ∇ 𝑇 (𝒖) )
the infinitesimal strain tensor. To model elasto-plasticity, we consider the following additive decomposition of the strain tensor:

𝜺 = 𝜺 𝑒 + 𝜺 𝑝 , ( 5
)
in which 𝜺 𝑒 and 𝜺 𝑝 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 Δ𝝈:

Δ𝝈 = 𝑫 ∶ Δ𝜺 𝑒 + 𝑫 𝑝𝑙𝑎𝑠 ∶ Δ𝜺 𝑝 , (6)
where 𝑫 and 𝑫 𝑝𝑙𝑎𝑠 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 [START_REF] De Souza Neto | Computational methods for plasticity[END_REF]. Plastic onset is defined by the yield surface 𝑓 (𝝈, ε𝑝 ):

𝑓 (𝝈, ε𝑝 ) = √ 3𝐽 2 (𝒔 (𝝈)) - ( 𝜎 𝑦0 + 𝐻 ε𝑝 ) ≤ 0, (7)
where ε𝑝 is the accumulated plastic strain, 𝐻 is the hardening modulus, 𝜎 𝑦0 is the initial yield strength, and

𝐽 2 (𝒔 (𝝈)) = 1 2 𝒔 (𝝈) ∶ 𝒔 (𝝈)
is the second invariant of the deviatoric stress tensor 𝒔 = 𝝈 -1 3 trace (𝝈) 𝑰. 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 [START_REF] De Souza Neto | Computational methods for plasticity[END_REF].

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. [START_REF] Gerbault | Prediction of faulting from the theories of elasticity and plasticity: what are the limits[END_REF][START_REF] Rybacki | Semi-brittle deformation of carrara marble: Hardening and twinning induced plasticity[END_REF].

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. [START_REF] Suppe | Fluid overpressures and strength of the sedimentary upper crust[END_REF], 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 pressuredependent yield criterion for future work. This point is further discussed in Section 5.

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, [START_REF] Alnaes | The fenics project version 1.5[END_REF]).

We use the latest legacy version (2019.1.0) of FEniCS which documents previous poromechanical studies (e.g. [START_REF] Haagenson | A generalized poroelastic model using fenics with insights into the noordbergum effect[END_REF]. 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):

𝐆 (𝜑, 𝒖, p) = 𝛼 ∫ Ω (div 𝒖 𝑡+1 -div 𝒖 𝑡 )𝜑 𝑑Ω + 1 𝑀 ∫ Ω (p 𝑡+1 -p 𝑡 )𝜑 𝑑Ω + 𝜅 𝜇 Δ𝑡 ( ∫ Ω ∇p 𝑡+1 ⋅ ∇𝜑 𝑑Ω -∫ 𝜕Γ 𝑞 (∇p ⋅ 𝑛)𝜑 𝑑Γ ) = 0 ∀𝜑 ∈ 𝐻 1 (Ω), 𝑡 ∈ [𝑡, 𝑡 + 1], (8)
𝐆 (𝜼, 𝒖, p) = ∫ Ω 𝝈 𝑡+1 ∶ 𝜺(𝜼) 𝑑Ω -∫ 𝜕Γ 𝑡 t ⋅ 𝑛 𝑑Γ = 0 ∀𝜂 ∈ 𝐻 1 (Ω), 𝑡 ∈ [𝑡, 𝑡 + 1], ( 9
)
where t is traction. Note that the boundary of the domain is divided as Γ = Γ 𝑝 ∩ Γ 𝑞 = Γ 𝑢 ∩ Γ 𝑡 , where Γ 𝑝 and Γ 𝑞 are the prescribed fluid pressure and Darcy velocity boundaries, and Γ 𝑢 and Γ 𝑡 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 𝒖 and p, respectively. This combination of elements, usually referred as a Taylor-Hood element, fulfills the inf-sup condition and therefore, is numerically stable [START_REF] Markert | A biphasic continuum approach for viscoelastic high-porosity foams: Comprehensive theory, numerics, and application[END_REF].

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 [START_REF] Cryer | A comparison of the three-dimensional consolidation theories of Biot and Terzaghi[END_REF] and the thick-walled hollow cylinder loading test [START_REF] Nadai | Theory of Flow and Fracture of Solid[END_REF]. Further details on this implementation as well as the benchmarks are available in our GitHub repository: https://github.com/FNSL1996/PEP_FEniCS.

Models setup and tests

We consider a spatial domain that geometrically and kinematically resembles the Planchón-Peteroa geothermal system, according to [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF]. 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 × 50 × 12 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 [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF]. This is represented using a parallelepiped of dimensions 10 × 5 × 0.5 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 [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF]. Its semiaxes dimensions are 4 × 2 × 2 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). We investigate the fluid flow, bedrock strain and stress at a kilometric scale, and as a first-order approximation, we 230 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 [START_REF] Cheng | Poroelasticity (Theory and Applications of Transport in Porous Media[END_REF], 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 (𝑆𝑅), the fault zone shear modulus (𝐺) and permeability (𝜅), the bedrock's yield strength (𝜎 𝑦0 ) 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 𝑀𝐴𝑋 , cf. Figure 5f).

G 𝜈 K K 𝑠 K 𝑓 𝛼 𝜙 𝜅 M [GPa] [-] [GPa] [GPa] [GPa] [-] [-] [m 2 ] [GPa] 6 0.2 8 36 2.25 0.778 0.19 1.90⋅10 -13 9.92
Table 2 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 𝑓 its shear modulus, 𝜅 𝑓 its permeability. 𝜎 𝑦0 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 𝑏 = 6 [GPa], and 𝜅 𝑏 = 1.9⋅10 -13 [𝑚 2 ], respectively. The last column displays a key result: NF 𝑀𝐴𝑋 is the maximum normalized flux measured in the area defined in Figure 5f)). See text for details. 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 [START_REF] Etheridge | Differential stress magnitudes during regional deformation and metamorphism: Upper bound imposed by tensile fracturing[END_REF], and it corresponds to experimental estimates [START_REF] Duboeuf | Aseismic deformations perturb the stress state and trigger induced seismicity during injection experiments[END_REF], 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. [START_REF] Kanamori | The diversity of the physics of earthquakes[END_REF]. Note that, because of the fault's thickness, slip-rates are not expected to have a considerable impact on bedrock stress.

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.

Solid deformation and stress state

We compute the mean stress σ and the stress shape ratio 𝜙 (Eq. 10), which give insights on the rock's confinement and on the relative magnitude between principal stresses, respectively.

σ = 1 3 trace (𝝈 (𝒖, p)) , and 𝜙 = 𝜎 2 -𝜎 3 𝜎 1 -𝜎 3 , ( 10
)
where 𝜎 1 , 𝜎 2 , and 𝜎 3 are the most compressive, intermediate, and least compressive principal stresses, respectively.

𝜙 is defined between 0 and 1, and addresses the stress tensor "shape" by comparing the magnitude of the principal stresses. A value of 𝜙 = 0.0 implies that 𝜎 2 = 𝜎 3 , therefore a preferential compressive direction, similar to a uniaxial compression test. If 𝜙 = 1.0 then 𝜎 1 = 𝜎 2 , which establishes a preferential tensile or least compressive direction, similar to a uniaxial tensile test. Hereafter we use 𝜙 to refer to the following stress states: 𝜙 = 0: uniaxial compression, 𝜙 ∈ (0, 0.5): transpression, 𝜙 = 0.5: strike-slip regime, 𝜙 ∈ (0.5, 1): transtension, and 𝜙 = 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 ± 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 ∼ ± 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 𝜎 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.

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 ∼ 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.

Fluid flux variation over time as a function of varying parameters

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.

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.

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.

Influence of Fault permeability

An increased fault permeability, mimicking fractured fault zones as seen in the SAVZ [START_REF] Sepúlveda | Thermomechanical behavior of a granodiorite from the liquiñe fractured geothermal system (39°s) in the southern volcanic zone of the andes[END_REF], 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.

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 𝑚 3 /s vs 0.202 𝑚 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.

Discussion

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. [START_REF] Mcclure | Discrete Fracture Network Modeling of Hidraulic Stimulation[END_REF][START_REF] Wang | Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using xfem with cohesive zone method[END_REF]. Because fluid flow in geothermal systems is heavily controlled by the presence of fracture networks [START_REF] Sanderson | The use of topology in fracture network characterization[END_REF], advantages and limitations of our modelling approach exist when comparing it against those approaches. XFEM is able to represent crack propagation and solve a fully-coupled solid-fluid medium [START_REF] Wang | Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using xfem with cohesive zone method[END_REF]. 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 [START_REF] Mcclure | Discrete Fracture Network Modeling of Hidraulic Stimulation[END_REF]. 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. [START_REF] Tapponnier | Slip-line field theory and large-scale continental tectonics[END_REF]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. [START_REF] Ord | The strength of the continental crust, detachment zones and the development of plastic instabilities[END_REF][START_REF] Suppe | Fluid overpressures and strength of the sedimentary upper crust[END_REF]. 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.

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 [START_REF] Nur | Aftershocks caused by pore fluid flow?[END_REF]: 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. [START_REF] Schultz | Growth of geologic fractures into large-strain populations: review of nomenclature, subcritical crack growth, and some implications for rock engineering[END_REF]. Such wing-like fractures are observed in many field areas (e.g. [START_REF] Kim | Structural similarity and variety at the tips in a wide range of strike-slip faults: a review[END_REF]. 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 [START_REF] Sielfeld | Oblique-slip tectonics in an active volcanic chain: A case study from the southern andes[END_REF] 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 [START_REF] Sibson | Stopping of earthquake ruptures at dilational fault jogs[END_REF][START_REF] Sibson | Earthquakes rupturing as a mineralizing agent in hydrothermal systems[END_REF][START_REF] Sibson | Fluid involvement in normal faulting[END_REF]. 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 [START_REF] Dobson | Porosity, permeability, and fluid flow in the Yellowstone geothermal system, Wyoming[END_REF].

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 𝜔:

𝜔 = 𝜅 𝜙𝜇𝛽 , ( 11
)
where 𝜙 is the rock porosity, and 𝛽 is the matrix compressibility (e.g. [START_REF] Reid | Massive collapse of volcano edifices triggered by hydrothermal pressurization[END_REF], and references therein). The higher 𝜔, 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 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. [START_REF] Nur | Aftershocks caused by pore fluid flow?[END_REF][START_REF] Bosl | Aftershocks and pore fluid diffusion following the 1992 landers earthquake[END_REF]. Aftershock seismicity due to fluid migration has been inferred in settings with fluids of different viscosity, ranging from CO 2 [START_REF] Miller | Aftershocks driven by a high-pressure co2 source at depth[END_REF], water, magma, or a mixture of both [START_REF] Legrand | The fluid-driven tectonic swarm of aysen fjord, chile (2007) associated with two earthquakes (mw=6.1 and mw=6.2) within the liquińe-ofqui fault zone[END_REF][START_REF] White | Volcano-tectonic earthquakes: A new tool for estimating intrusive volumes and forecasting eruptions[END_REF]. Its duration ranges from weeks to a month with water-viscosity fluids [START_REF] Miller | Aftershocks driven by a high-pressure co2 source at depth[END_REF] and up to six months with water-magma mixtures [START_REF] Legrand | The fluid-driven tectonic swarm of aysen fjord, chile (2007) associated with two earthquakes (mw=6.1 and mw=6.2) within the liquińe-ofqui fault zone[END_REF]. 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. [START_REF] Dorbath | Seismic response of the fractured and faulted granite of soultz-sous-forets (france) to 5 km deep massive water injections[END_REF][START_REF] Moeck | Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir[END_REF] 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 [START_REF] Pearce | Reactivation of fault systems by compartmentalized hydrothermal fluids in the southern andes revealed by magnetotelluric and seismic data[END_REF]. 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 [START_REF] Dorbath | Seismic response of the fractured and faulted granite of soultz-sous-forets (france) to 5 km deep massive water injections[END_REF].

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). 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.

Acknowledgements

F

  . Sáez-Leiva et al.: Preprint submitted to Elsevier

Figure 2 :

 2 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.

Figure 3 :

 3 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 (𝜙) in planar view at 5 km depth. Plastic envelope is represented as a yellow contour in the middle column, and 𝝈 𝟏 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. 𝜎 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.

Figure 4a :

 4a 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).

Figure 4b :

 4b 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).

  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 2 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 F. Sáez-Leiva et al.: Preprint submitted to Elsevier computed over time (Figure 5a-e). 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 5 :

 5 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 2 and is used for the flux calculation.

Figure 6 :

 6 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.

Figure 7 :

 7 Figure 7: (a) Suction pump cyclic mechanism at a dilational jog, showing a coseismic reduction in fluid pressure (P 𝑓 ) below hydrostatic pressure (P ℎ ). 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.

Figure 8 :

 8 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.

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.

  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.

Table 1

 1 Mechanical and hydraulic parameters (after Cheng, 2016). From left to right: G is the shear modulus, 𝜈 the Poisson's modulus, K the rock bulk modulus, K 𝑠 the solid bulk modulus, K 𝑓 the fluid bulk modulus, 𝛼 Biot's constant, 𝜙 rock porosity, 𝜅 rock permeability, and 𝑀 Biot's modulus.