Hydro-mechanical earthquake cycles in a poro-visco-elasto-plastic fluid-bearing fault structure

A major goal in earthquake physics is to derive a constitutive framework for fault slip that captures the dependence of shear strength on fault rheology, sliding velocity, and pore-fluid pressure. In this study, we present H-MEC (Hydro-Mechanical Earthquake Cycles), a newly-developed two-phase flow numerical code - which couples solid rock deformation and pervasive fluid flow - to simulate how crustal stress and fluid pressure evolve during the earthquake cycle on a fluid-bearing fault structure. This unified, continuum-based model, incorporates a staggered finite difference-marker-in-cell (SFD-MIC) method and accounts for full inertial (wave mediated) effects and fluid flow in poro-visco-elasto-plastic compressible medium. Global Picard-iterations and an adaptive time stepping allows the correct resolution of both long- and short-time scales, ranging from years during slow tectonic loading to milliseconds during the propagation of dynamic ruptures. We present a comprehensive in-plane strike-slip setup in which we test analytical poroelastic benchmarks of pore-fluid pressure diffusion from an injection point along a finite fault width. We then investigate how pore-fluid pressure evolution and solid-fluid compressibility control sequences of seismic and aseismic slip on geologic faults. While the onset of fluid-driven shear cracks is controlled by localized collapse of pores and dynamic self-pressurization of fluids inside the undrained fault zone, subsequent dynamic ruptures are driven by solitary pulse-like fluid pressure waves propagating at seismic speed. Furthermore, shear strength weakening associated with rapid self-pressurization of pore-fluid can account for the slip-fracture energy scaling observed for large earthquakes.


Introduction
There is a growing interest in understanding how geologic faults respond to transient sources of fluid. This issue has become a societal concern due to the increasing interest in geothermal energy (e.g., Grigoli et al., 2017;Terakawa et al., 2012), CO 2 storage (Zoback and Gorelick, 2012), or injection of waste waters associated with hydraulic fracking (Ellsworth, 2013). Natural and artificial sources of fluid can elevate pore-fluid which can destabilize active faults through either slow-slip or fast, seismic-wave-producing rupture events perceived as earthquakes. Examples of fluid-driven seismic and aseismic slip on natural faults are related to seismic swarms (e.g., Roland and McGuire, 2009;Ross et al., 2020;Chen et al., 2012;Hatch et al., 2020;Hainzl, 2004), aftershock sequences (e.g., Miller, 2020;Nur and Booker, 1972;Ross et al., 2017;Miller et al., 2004), and other natural phenomena such as slow-slip events and tectonic tremors (e.g., Jolivet and Frank, 2020;Bürgmann, 2018;Behr and Bürgmann, 2021). On the other hand, fluid injections associated with human activities are well known to lead to the reactivation of faults, slow-slip transients, and earthquakes, which occasionally are large enough to cause damage (e.g., Guglielmi et al., 2015;Deichmann and Giardini, 2009;Ellsworth et al., 2019;Keranen et al., 2013;Raleigh et al., 1976;Keranen et al., 2014;Bao and Eaton, 2016). Taken together, this body of work suggest that fluid-driven seismic and aseismic fault slip in nature and anthropogenic fluid injections are closely linked.
The propagation of fluid-driven shear cracks, the evolution of pore pressure, and the poroelastic response of faults, have been subjects of interest for decades in the earthquake modeling community (e.g., Rice and Cleary, 1976;Rice, 1992;Dunham and Rice, 2008;Rudnicki and Rice, 2006;Cruz-Atienza et al., 2018;Heimisson et al., 2019;Zhu et al., 2020;Petrini et al., 2020). Recent modeling efforts have mostly focused on the stability of frictional slip when a fluid is locally injected and diffuses within a fault (Bhattacharya and Viesca, 2019;Cappa et al., 2018;Garagash and Germanovich, 2012;Larochelle et al., 2021;Yang and Dunham, 2021), or when fluid injection initiates shear cracks on a strengthening rate-and-state frictional fault (Dublanchet, 2019). Despite the apparent relevance of fluid-driven slow and fast slip transients in a wide variety of natural and anthropogenic environments, the spatio-temporal evolution of sequences of seismic and aseismic slip in response to pore-fluid evolution remains poorly constrained. This is, in part, due to the challenge of solving a fully coupled solid-fluid system in which the dynamic evolution of fault slip, stress, pore-fluid pressure, and shear strength are unknown.
In this study, we present H-MEC (Hydro-Mechanical Earthquake Cycles), a continuum-based modeling approach to simulate fully dynamic sequences of seismic and aseismic slip (SEAS) in a poro-visco-elastoplastic compressible medium with rate-dependent strength, which considerably extends our previous numerical model developed for incompressible media with rate-independent strength (Petrini et al., 2020). The computational method is developed for a plane-strain problem on a simplified strike-slip fault. Our coupled formulation employs a Picard-type non-linear solver with a fully implicit, first order accurate time integrator that utilizes an adaptive time stepping to efficiently solve the system through multiple seismic cycles. We apply a staggered finite difference-marker-in-cell (SFD-MIC) method, as it is easy to program, efficient, and the marker-in-cell technique can be applied to simulate highly deforming materials, such as faults and geologic structures (Gerya and Yuen, 2007;Petrini et al., 2020).
The purpose of this study is to introduce a quantitative set of numerical experiments in which we test analytical benchmarks and illustrative examples of sequences of seismic and aseismic slip along with results from the solid-fluid coupling, fluid transport on-and off-fault, pore pressure evolution, and fault slip. Our focus is to demonstrate the validity of our methodology and highlight the mechanical processes and phenomena that arise from this solid-fluid coupling in a generic sense. We apply this methodology in the context of a simple, 2-D plane-strain shear model of a strike-slip fault in a uniform bulk domain, which is an idealization of a natural strike-slip fault. While parameter choices are chosen to be representative of a strike-slip fault, we are not attempting to model any specific fault or earthquake sequence. In Section 2, we present the conceptual and mathematical model and state the continuum problem solved in this work. In Section 3, we present results of two analytical benchmarks and a set of simulations of earthquake sequences. Finally, Section 4 we analyse and interpret the results and compare solutions for different confining pressure, solid and fluid compressibility, porosity, and fracture energy. Lastly, we discuss future avenues for research through the use of this methodology.

Methodology
In this section we present the conceptual, mathematical, and physical foundations of the model before presenting the numerical simulations that will follow.

Governing equations
The mathematical description of a deforming continuum in a stationary gravity field is given by the following set of conservation equations for total momentum (solid matrix and fluid; Eq. 1), fluid momentum (Eq. 2), fully compressible solid mass (Eq. 3), and fully compressible fluid mass (Eq. 4):  [s] )), i.e., the volume flux relative to the moving solid matrix, k [φ] is the porosity-dependent permeability, η [f ] is the fluid viscosity, p [t] and p [f ] are the solid and fluid pressure, respectively, K [d] is the drained bulk modulus, η [φ] is the effective compaction viscosity, v [f ] is the fluid velocity vector, α is the Biot-Willis coefficient (Biot and Willis, 1957), B is the Skempton coefficient (Skempton, 1960;Bishop, 1973), and g is the gravity vector. The total momentum (Eq. 1) neglects differences in the solid and fluid accelerations (Gerya and Yuen, 2007;Petrini et al., 2020).
The porosity-dependent permeability is computed as where k * and φ * are the reference permeability and reference porosity (Table 1), respectively, and n is an porosity-dependent exponent, which for natural pores is assumed to be 3 (Rice, 1992;Connolly and Podladchikov, 2000), and implies that permeability changes as a cube of increasing porosity. While different k * are tested in this study, φ * is set to be 0.01 (i.e. 1%, Peacock et al., 2011) and does not evolve with time.
The drained bulk modulus (K [d] ), the Biot-Willis coefficient (α), and the Skempton coefficient (B) are defined as where K [φ] the effective bulk modulus of pores, K [s] is the solid bulk modulus, β [d] is the compressibility of the solid skeleton (i.e., drained compressibility of the porous medium), β [s] is the compressibility of the solid phase, and β [f ] is the compressibility of the fluid phase.
The stress is decomposed into a deviatoric and volumetric component: where τ is the deviatoric stress tensor, and I is the identity tensor. The total pressure (p [t] ) and density (ρ [t] ) are coupled with the porosity and computed from their respective fluid and solid matrix quantities The effective bulk modulus of pores (K [φ] ), as well as the effective compaction viscosity (η [φ] ) are computed following Yarushina and Podladchikov (2015) and Petrini et al. (2020) as where m is a geometrical factor, µ is the shear modulus, and η s is the shear viscosity of the solid matrix. In this study m takes a value of 1 (i.e., cylindrical pores), which reduces Eqs. (9a) and (9b) to µ * /φ and η s /φ. For in-plane (mode-II) slip, the effective shear modulus, which is simply the effective elastic modulus for quasi-static plane-strain deformation (Rubin and Ampuero, 2007), is equal to Where µ is the shear modulus, ν is Poisson's ratio.
The continuity equations of both solid (Eq. 3) and fluid (Eq. 4) contain poroelasticity terms consistent with Biot's theory (Biot, 1941;Gassmann, 1951;Yarushina and Podladchikov, 2015), which allow for compressibility of the solid matrix and fluid as well as viscous and elastic (de)compaction of the interconnected porous space. Being Eqs.
(3) and (4) fully coupled, convergence/divergence of the solid matrix occurs in response to local compaction/decompaction processes, forming a fully coupled hydro-mechanical system.

Visco-elasto-plastic rheology
The deviatoric deformation occurs as a combination of elastic, viscous, and plastic deformation. Components of the deviatoric strain-rate tensor are defined via spatial grandients of the velocity under the assumption of infinitesimal strain as follows: where i and j are coordinate indexes and x i and x j are spatial coordinates such that in 2-D we can define four tensor components. We assume that strain rate can be additively decomposed into its elastic, viscous, and plastic components aṡ The elastic part is defined as  where µ is the shear modulus, while the differential operatorD/Dt denotes the co-rotational time derivative -that is, the local rates of change in the deviatoric stress with respect to translations, rotations, and/or deformation.
Viscous creep obeys Newton's law of viscous friction and is defined here via a single effective shear viscosity according to where the τ II andε II(v) denote the square root of the second invariant of the deviatoric stress and the viscous strain-rate tensors, respectively. The effective shear viscosity η s is then used to define the viscous contribution to the strain-rate tensor: Plastic (irreversible) deformation and frictional sliding is governed by a rate-dependent strength of a non-associated plasticity model (Yi et al., 2018) whereε II [p] is the square root of the second invariant of the deviatoric plastic strain ratė ε 0 is the reference strain rate, γ is the rate-strengthening exponent, and τ 0 represents the inclusion of the Drucker-Prager yield function defining the maximum allowable shear stress a rock can withstand (Prager and Drucker, 1952) where c is the cohesional strength and f is the reference, static friction coefficient. The slip rate V (t) is then computed as the magnitude of the second invariant of deviatoric plastic strain rateε II(p) over the thickness of the fault zone h In our continuum model, plastic deformation is computed as volumetric strain over a finite thickness, and it is represented by a tensor. This means that plastic deformation can spontaneously localize anywhere (Dal Zilio et al., 2022). Note that the rate-dependent strength (Eq. 16) relies on the presence of γ, which mimic the positive direct effect in the rate-and state-dependent friction formulation (e.g., Scholz, 1998), a feature that has ample laboratory confirmation (Dieterich, 1979(Dieterich, , 1981Ruina, 1983;Marone, 1998). On the other hand, our rate-dependent strength formulation does not require the presence of an evolutionary effect involving a decrease in friction. This is because the weakening mechanism in our fully coupled system is controlled by the evolution of pore-fluid pressure (p [f ] ), which is encapsulated in τ 0 (Eq. 18).
The non-associated plastic flow law is defined through the plastic flow potential (Q), which reflects the amount of mechanical energy per unit volume that supports plastic deformation (Vermeer, 1998) In this study, we assume the dilation angle ψ to be zero, which implies that the volume does not change during plastic yielding (Q = τ II − c). Accordingly, the deviatoric plastic strain is defined as with χ indicating the plastic multiplier, which serves to ensure that, if yielding occurs, the deviatoric stress will always satisfy the yield criteria (τ II = τ y ).
The deviatoric components (τ ij ) of the stress tensor σ in Eq. (1) are formulated from the visco-elastoplastic constitutive relationships (Eq. 12) by using an implicit first-order finite-difference scheme in time in order to represent objective time derivatives of visco-elastic stresses (e.g., Moresi et al., 2003) where Z is the visco-elasticity factor (Schmalholz et al., 2001) and η vp is the effective visco-plastic viscosity that characterizes the intensity of the plastic deformation whereε II [p] is the square root of the second invariant of the deviatoric plastic strain rate, and η m is the matrix viscosity (Katz et al., 2006) in which λ = −29 defines an experimentally derived porosity-weakening coefficient (Katz et al., 2006).

Global visco-elasto-plastic Picard iterations
We accurately satisfy the plastic yielding condition at the Eulerian nodal points of the staggered grid by using global Picard iteration (Gerya, 2019). As a result, we discretize the rate-strengthening yielding condition (Eq. 16) to satisfy the condition τ yield = τ II : The Picard iteration is performed by repeating solution of Eqs (1−4) and evaluating [ τ II ] elastic at each basic node as Where η vp is the current local value of the effective visco-plastic viscosity of the solid matrix, which defines the relationship between the deviatoric stress and the deviatoric strain rate for a viscous rheology. When [ τ II ] elastic > τ yield , we compute the new visco-plastic viscosity at the node as If η vp < η s , the new visco-plastic viscosity from Eq. (29) is used for the next iteration. The global plastic yielding condition error is then computed on the nodal points where N is the cumulative number of nodal points at which either the previously computed or new value of η vp satisfies the condition η vp < η s . The stopping criteria of the Picard iteration is reached when ∆τ yield has decreased below a desirable level of tolerance (δ err ) where δ err = 10 Pa.
An important mechanism for fluid pressurization is the feedback between visco-plastic viscosity and compaction viscosity (often called bulk viscosity), which characterizes solid matrix resistance to reversible (i.e., elastic, β [d] ) and irreversible (i.e., visco-plastic, η [φ] ) pore compaction/decompaction. When plastic yielding occurs, we first compute the plastic viscosity (η p ) at the basic nodes (b) as We then compute the new plastic viscosity at the center of cells (i.e., pressure nodes, p) by using the harmonic average of the plastic viscosity values from the four surrounding basic nodes (b) (Gerya, 2019) Finally we compute the new compaction viscosity at the pressure nodes as and the new visco-plastic viscosity at the pressure nodes as where η 0 is the reference solid viscosity on the pressure nodes (Table 1).

Adaptive time-stepping
To resolve different time-scales -ranging from long-term (slow) deformation to periodic (fast) slip events -a variable time stepping is required during computation. To accomplish this, we require that the time step is the minimum of the time steps needed to resolve (1) The first time step requires that the slip acceleration per time step is limited to a fraction of the grid size as where δ d = 10 −5 defines the maximum grid fraction. Similarly, we require that the displacement per time step is limited by the same fraction of the grid size Based on a previous study (Herrendörfer, 2018), we apply a visco-elasto-plastic time step to capture the relaxation time scale η vp /µ * by a fraction ξ where ξ = 0.2. Notably, the visco-elasto-plastic relaxation time step combines both the visco-elastic relaxation η vep /µ * and the elasto-plastic relaxation time scale As a result, ∆t vep provides a constraint on the final time step when viscous deformation is dominant, thus ensuring the resolution of the relaxation time. Importantly, we do not apply any minimum time step cutoff. This is in contrast to previous studies (e.g., Lapusta and Liu, 2009), in which the time step is only adapted until it reaches ∆t min = (1/3) ∆ x / c s , where c s is the shear wave speed (see Table 1). We do not adopt a minimum time step cutoff because during the propagation of fully dynamic ruptures the actual time step can decreases below the ∆t min threshold. As a result, applying a cutoff would reduces the convergence rate during Picard iterations (Section 2.3).

Spatial discretization
We adopt a two-dimensional (2-D), staggered-grid finite difference-marker-in-cell (SFD-MIC) method, in which the grid is comprised of quadrilaterals cells parallel to the underlying coordinate basis vectors (Gerya and Yuen, 2007). Discretization of governing equations is given in Appendix A. The numerical method is developed for an in-plane problem using a fully-staggered grid, which means that physical variables are defined at different geometric points (Fig. A.11). In order to model multiple rock types and permit these rocks to undergo large displacement, as it often occurs on geologic faults, the Eulerian spatial discretization is combined with a Lagrangian marker-in-cell advection scheme. A set of Lagrangian particles (or material points) advects through the Eulerian domain used to discretize the velocity, pressure, and stress field. The Lagrangian particles are used to assign rock properties (e.g., density, viscosity, shear modulus, and hydro-mechanical properties), and, during the simulation, they store deviatoric stresses, pressure, and velocities for evaluating the respective time derivatives. A 4-th Runge-Kutta scheme is then used to advect markers to their new position according to the velocity field, whereas stresses are rotated according to the vorticity field. Lastly, a standard distance-dependent bi-linear interpolation scheme is used to interpolate physical parameters from the Lagrangian markers to the staggered Eulerian nodes at the beginning of each time step and update marker values at the end of each time step.

Model setup and boundary conditions
The 2-D model setup consists of a dextral strike-slip fault (Fig. 1). We consider a mode-II plane-strain shear motion along a planar interface embedded in a poro-visco-elasto-plastic medium (Fig. 1). We utilize an x − y Cartesian coordinate system with y = 0 being the sliding interface. This simple setup is designed to enable comparison to both analytical benchmarks for numerical solution and classical earthquake cycle simulations and theoretical estimates of the critical length scales during nucleation and propagation of dynamic ruptures (e.g., Lapusta et al., 2000). Parameters for the reference model setup are given in Table 1.
The system of equations over the domain Ω = [0, L x ] × [0, L y ] has a unique solution given appropriate boundary conditions for the velocity and pressure (Fig. 1). For the specifics of the model setup presented in this study, the conditions on the left (x = 0) and right (x = L x ) boundaries are defined free slip as follow: A free slip condition requires that the normal velocity component on the boundary is zero and the other component do not change across the boundary. On the other hand, the fault is loaded from the upper (y = 0) and lower boundaries (y = L y ) by applying Dirichlet boundary conditions for horizontal velocity v x and zero velocity on the vertical component v y (Fig. 1): where V pl = 6.3 cm yr −1 corresponds to the long-term loading rate.
Pressure boundary conditions have to be defined for both total (p [t] ) and fluid (p [f ] ) pressure. Since the solid mass (Eq. 3) and fluid mass (Eq. 4) are intrinsically coupled, we first initialize the fluid pressure with a constant confining pressure (p conf ) and then we apply a constant effective pressure (p eff ) on the top and lower boundary: Fluid momentum (Eq. 2) is solved assuming free slip condition for the left and right boundaries ∂v y whereas an inward (ṗ + ) and outward (ṗ − ) fluid flux is defined on the top and lower boundary in the form of Neumann boundary conditions ∂v y

Fluid injection benchmarks
To test the robustness of the code in a broad range of applications relevant to solid-fluid interaction, in this section we compare our modeling results against analytical benchmarks of pore-fluid pressure diffusion. In particular, we consider two benchmark exercises (BP1 and BP2) from a set of quasi-static problems in a 2-D domain with a 1-D fault subjected to perturbations of fluid injection and along-fault pore-fluid diffusion. In order to verify these two benchmarks, we modify the model setup given in Section 2.6 by including an injection point on the left boundary (Fig. 2a). For the first benchmark problem (BP1), injection is modeled as a source of constant fluid pressure (∆p [f ] ) at x inj = x (0) . The resulting pore-fluid pressure diffusion along the interface can be solved analytically as (e.g., Viesca, 2021) where erfc is the complementary error function, l d (t) = √ 4 α t is the diffusion lengthscale, and α is the constant hydraulic diffusivity where η [f ] is the viscosity of the permeating fluid and β * is a storage coefficient reflecting the compressibility of the fluid and porous matrix. On the other hand, the second benchmark problem (BP2) assumes a constant volume rate injection from a point source (x inj = x (0) ), and the flow solution is given by Carslaw and Jaeger (1959): where Q w is the constant injection volume rate, h is the fault thickness, and E1 is the expo-nential integral function. Note that the fault hydraulic transmissivity represents the product k [φ] h.
For both benchmarks BP1 and BP2, fluid is injected into the fault and it is constrained to flow only within its hydraulic aperture, a scenario that would occur when the permeability of the host rock is negligible relative to the fault itself (i.e. no fluid leak-off). Also, the fault is assumed marginally pressurized -that is, ∆p [f ] is always below the total pressure (in order to avoid hydraulic fracturing). Figure 2 shows the results from the BP1 (Fig. 2b) and BP2 (Fig. 2c) and demonstrates the high accuracy of the numerical solution, which overlaps with the analytical solutions at different time stages. It is worth noting that the exponential integral function in Eq. (48) is singular at the injection point x = x inj (Fig. 2c). This means that, for the case of injection experiments at constant volume rate, the fluid pressure at the origin is unbounded.

Earthquake cycles
The response of the fault to slow tectonic loading is characterized by long periods of quasi-static deformation followed by short, and fast, slip events ( Fig. 3a-c), in which the slip rate on the fault accelerates from ∼cm yr −1 to ∼m s −1 (Fig. 3d). Since the time step is set to be inversely proportional to slip velocity, relatively large time steps -of the order of a fraction of a year -are used in the inter-seismic period, while small time steps of the order of milliseconds are used to simulate the evolution of each dynamic rupture (Fig. 3e). Similarly, the effective pressure (p [t] − p [f ] ) drops by ∼20 MPa as a result of the poroelastic response of the fault zone during fast shearing (Fig. 3f).
Slip history on the fault indicate that, when assuming homogeneous properties, regular cycles of complete ruptures arises from the numerical experiments (Fig. 4). In the early stages of the simulations, the imposed loading rate increases stress linearly with time, whereas slip rate increases exponentially with time on the left and right side of the fault, i.e., where γ is high (Fig. 1b). As a result, while the two lateral segments of the fault steadily creep at rates comparable to the imposed loading rate (Fig. 4a), the seismogenic zone remains locked during inter-seismic periods and stress concentrates at the transition between high and low γ. In the later stages of the inter-seismic periods, the rate-strengthening behaviour allows creep to penetrate from both the left and the right segments towards the center of the fault, thus leading to a mechanical erosion of the locked patch. Due to a relatively large critical nucleation size (L c ), seismic events nucleate at the center of the seismogenic zone and the resulting dynamic ruptures typically propagate in a bilateral and symmetric fashion (Fig. 4a). For the reference model, the average slip rate is 0.82 m s −1 , whereas the average rupture speed is ∼2.35 km s −1 .
Evolution of the cumulative slip over multiple earthquake cycles is displayed in Figure 4b, with the accumulation of slip on the fault during inter-seismic periods represented by gray lines plotted every 4 years, whereas red lines display cumulative slip every 0.4 s when the maximum slip rate on the fault exceeds the following velocity threshold (V th ) which, according to our parameters, yields to a threshold of 8.1 cm s −1 . When the fault is experiencing a seismic event, the cumulative slip on the fault indicates that a complete fault rupture produces ∼5 m Figure 3: Spatial horizontal velocity and temporal evolution of earthquake cycles. Snapshots of the spatial horizontal velocity (vx) distribution during the phases of (a) nucleation, and (b-c) dynamic propagation. The response of the fault to tectonic loading is characterized by long inter-seismic periods of quasi-static deformation, in which the seismogenic zone is practically locked, whereas lateral segments of the fault zone steadily creeps at rates comparable to the plate convergence rate. Note that the color scale of vx velocities in (a-c) are not the same. Temporal evolution of (d) maximum slip velocity (e) time-step, and (f ) effective pressure on the fault. The time step is set to be inversely proportional to slip velocity, which results in large time steps in the inter-seismic periods, and small time steps (of the order of milliseconds) to simulate the evolution of each dynamic rupture.
of slip, which is followed by a transient post-seismic deformation of ∼0.8-1 m (blue lines; Fig. 4b). Our models show that dynamic ruptures occur due to an abrupt increase of pore-fluid pressure within the fault zone (Fig. 5a). The increase in pore-fluid pressure is simultaneously coupled to a rapid decrease in shear strength (Fig. 5b) and both shear and compaction viscosity (Fig. 5c), which are self-sustained by a localized strain rate. When an event begins to propagate dynamically, it produces a shear stress breakdown from its maximum static value to a low, dynamic value (Fig. 5b). Notably, at co-seismic slip rates, the effective shear viscosity drops by roughly 10 orders of magnitude (Fig. 5c), whereas static stress change immediately after the event results in a stress concentration at the transition between the lateral creeping segments and the seismogenic zone.
To investigate the influence that pore-fluid pressure level has on the spontaneous activation of seismic and aseismic slip on the fault, we perform further numerical experiments in which we decrease (model M2) and increase the effective pressure (p eff ) (model M3) compared to the value used for the reference model (Fig. 6). The two models illuminate the effect of pore-fluid pressure throughout the earthquake cycle, from the recurrence period to the nucleation phase of seismic events, as well as the relation between seismic and aseismic slip. In model M2 (Fig. 6b), due to a lower pore-fluid pressure, the effective pressure on the fault is higher (p eff = 60 MPa), and consequently the critical nucleation size decreases. As a result, seismic events nucleate earlier than in the reference model (Fig. 3d). Furthermore, since we assume a slightly asymmetric fault structure (Fig. 1b), seismic ruptures nucleate on one side of the seismogenic fault and propagate unilaterally on the other side of the fault segment in an asymmetric fashion. On the other hand, model M3 (Fig. 6c) assumes a higher pore-fluid pressure, which results in a lower effective pressure (p eff = 20 MPa). In this case, the average shear strength is significantly lower than in model M2 and plastic yielding occurs more easily. Consequently, the seismogenic zone promotes unstable slip, either as seismic partial ruptures or as slow-slip transients. The long-term slip history indicates a mixture of fast and slow events, with a recurrence interval that vary over time. These results are important, as they suggest that pore-fluid pressure evolution in the same fault segment can release tectonic stress through both slow and fast ruptures, without incorporating any rate-and state-dependent friction.
According to our results, it appears clear that visco-plastic compaction of fluid-filled porous media can lead to an abrupt self-pressurization of the fault zone, particularly under undrained conditions or when the permeability on the fault is sufficiently low to prevent fluid flux on the boundaries of the fault zone on the timescale of seismic ruptures. As a result, the self-pressurization of the fault zone can trigger the propagation of solitary pressure waves, which can travel at co-seismic speed. To analyze the effect of solid-fluid (de)compaction, we quantify the separate contributions of visco-plastic compaction of the solid skeleton (ζ [vp] ) and the decompaction of the fluid phase (ζ [e] ) as follow: (De)compaction along the fault during a representative major event is shown in Figure 7. Temporal evolution of the slip velocity displays a large event propagating unilaterally from the left to the right side of the seismogenic zone (Fig. 7a). An accurate quantification of visco-plastic compaction from Eq. (50) and elastic decompaction from Eq. (51) indicates that plastic yielding during the gradual nucleation phase is associated with a rapid increase in solid compaction (Fig. 7b), which is followed by a negative increase of fluid decompaction (Fig. 7c). As a result, pore-fluid pressure on the fault increases while shear strength decreases, thus causing a shear instability. Such shear instability grows on a relatively short time scale, causing the propagation of a dynamic rupture in the form of a pulse-like pressure wave (Fig. 7d). Notably, the strain difference between visco-plastic compaction and elastic decompaction represents the finite pore volume change. The combination of these two interconnected processes results in a small change in porosity and, in turn, permeability (Eq. 5). Furthermore, a dynamic rupture causes an increased in pore-fluid pressure due to elastically stressed pores, which is better resolved by our fully compressible model with rate dependent strength compared to the previous simplified incompressible formulation (Petrini et al., 2020). Such increases in pore-fluid pressure are released over a relatively long diffusion period, which essentially depends on the permeability of both the host rock and the fault zone.
It is important to note that the magnitude of visco-plastic and elastic (de)compaction depend on the compressibility of both the solid skeleton (β [s] ) and the fluid phase (β [f ] ), as they have a direct impact on the Biot-Willis coefficient (α ; Eq. 6b), and the Skempton coefficient (B ; Eq. 6c). To quantify the role of solid and fluid compressibility, we execute two parameter studies of solid-fluid compressibility assuming two end-member porosities (i.e., 1% and 10%; Tewksbury-Christle and , knowing that in nature the compressibility of porous rocks, and in particular fluids, can vary up to two orders of magnitude (e.g., Zimmerman et al., 1986;Span and Wagner, 1996;Mitchell and Soga, 2005). Our models produce distinctly different slip patterns within the range of solid-fluid compressibility (Fig. 8).
As illustrated in Figure 8a, an increase in fluid compressibility leads to a transition from regular seismic events, to slow-slip events and eventually stable creep. In contrast, solid compressibility plays only a subordinate role. Furthermore, when we increase the reference porosity from 1% to 10% (Fig. 8b), the parameter space in which we observe slow-slip events broaden. While high porosity serves as a damping factor to suppress the self-pressurization of pore-fluid and thus the occurrence of fast events, these results indicate that fluid compressibility can significantly affect the slip response of a fault loaded by tectonic stresses. As such, the higher the fluid compressibility, the lower the self-pressurization of the fault zone. Our numerical results can thus reproduce the full spectrum of fault slip behaviors under geophysically relevant conditions of pore-fluid pressure and fault composition, from slow to fast events, as observed for tectonic faults (e.g., Bürgmann, 2018;Obara and Kato, 2016;Jolivet and Frank, 2020).

Nucleation length, cohesive zone size, and fracture energy
The computational framework presented here is capable of resolving long-term deformation histories with continuous aseismic creep in the stable rate-strengthening fault regions throughout the loading period, nucleation of shear instabilities, co-seismic propagation of seismic ruptures, and post-seismic relaxation.
As the considered example shown in Figure 9, the algorithm is formulated with a rate-strengthening yield strength (Eq. 16), in which the rate-strengthening exponent (γ) is decisive for its success to reproduce long-term inter-seismic periods with essentially quasi-static deformation, aseismic slip, and most importantly, gradual nucleation of dynamic ruptures (Fig. 9a) (e.g., Dieterich, 1992;Lapusta et al., 2000). However, compared to the formulation of the classical rate-and state-dependent friction laws (e.g., Scholz, 1998), the evolutionary effect during the weakening phase is controlled by self-pressurization of pore-fluid pressure.
Our models indicate that shear cracks would become unstable when they reach the critical nucleation length (L c ) proposed by Andrews (1976), in which shear strength linearly decreases from static shear strength (τ s ) to a relatively low dynamic shear strength (τ d ) over a characteristic slip weakening distance (d c ) (Ida, 1972) where τ 0 represents the initial shear stress (Eq. 18). According to our model (Fig. 9a), Eq (52) predicts a critical nucleation length of ∼7.3 km. Comparing Fig. 4a and Fig. 6a with Fig. 9a, we notice that slip during and right after the nucleation phase are very similar despite the different critical nucleation length. This means that observing signals from the nucleation of fluid-driven shear cracks do not provide any evidence on the final size of the event. This suggests that the final size of seismic events is determined by the conditions on the fault, rather than by the nucleation process itself.  Figure 9b illustrates the weakening phases that buffer the pore pressure and slip velocity rise, leading to dynamic rupture. During Phase I, upon rapid loading, shear stress rises abruptly with strain as the fault is elastically loaded previous to the onset of slip. Phase II is characterized by a strain-hardening shear strength, which produces an early phase of slip (<1 cm of slip), and corresponds roughly to slip rates of the order of 1-10 cm s −1 . During Phase III, an abrupt weakening is initiated and shear strength drops while slip rate quickly accelerates up to ∼m/s. Peak in slip rate is reached when strength drop is completed, which is follow by a quasi steady-state phase where low shear strength is maintained with minor fluctuations.
When an event nucleates, shear stress drops while slip velocity increases rapidly. Since the rupture speed at this stage is still close to zero, the length along which the stress drop occurs represents the quasi static cohesive zone length (Λ). The cohesive zone size is an important resolution criterion in dynamic rupture because it gives the spatial length scale over which the shear stress drops from its peak to residual values at the propagating rupture front (Palmer and Rice, 1973;Day et al., 2005). Since fluid-driven shear cracks essentially mimic a linear slip-weakening law (Fig. 9c), the cohesive zone sizes at rupture speed c = 0 + (Λ 0 ) observed in our simulations correspond quite well to the estimate proposed by Day et al. (2005): where 9π/32 is a constant if the stress traction distribution within the cohesive zone is linear in space (Palmer and Rice, 1973). This latter study establishes that Λ 0 /∆x of 3 to 5 is required to resolve dynamic rupture. The ratio of nucleation zone size and cohesive zone size (Λ 0 /L c ) is between 0.3 and 0.5, as reported in previous studies (e.g., Lapusta et al., 2000). Therefore, resolving the cohesive zone is the more stringent numerical criterion here. In our models, we choose the cell size small enough to resolve the cohesive zone of 2.5 km with at least 10 cells. Hence the spatial discretization in the simulations, with the cell size of ∆x=250 m, is small enough to resolve the evolution of stress and slip rate.
Seismic ruptures are controlled by an energy balance involving elastic work, wave radiation, and dissipation by anelastic processes, including friction and plastic strain. As illustrated in Figure 9c, the frictional work during the weakening process -as the product of shear stress and slip -equates to a simple form of fracture energy (G ), which assumes a linear slip weakening (Palmer and Rice, 1973;Andrews, 1976;Okubo and Dieterich, 1984) To compare our modeling results with existing measurements of fracture energy from earthquakes globally, we compile fracture energy estimates made previously with events spanning five orders of magnitude in size, from borehole microseismicity to great earthquakes (Abercrombie and Rice, 2005;Rice, 2006;Malagnini et al., 2014;Viesca and Garagash, 2015;Nielsen et al., 2016a;Tinti et al., 2005) (Fig. 10). Previous inferences have shown a nonlinear scaling of fracture energy with slip (Abercrombie and Rice, 2005), from G ∝ δ 2 for small earthquakes to G ∝ δ 2/3 for large earthquakes (Fig. 10), which has been attributed to thermal pressurization of pore-fluid by the rapid shear heating of fault gouge (Viesca and Garagash, 2015). For the set of parameters considered in this study, our models produce both relatively large stress drop (τ s − τ d ) and large characteristic slip weakening distance (d c ), which results in high values of fracture energy (Fig. 9c). Remarkably, the range of fracture energy and slip from our dynamic analysis supports the observed G ∝ δ 2/3 scaling of fracture energy for large earthquakes, and suggests that self-pressurization of fluids is a viable mechanism for explaining widespread and prominent process of fault weakening.
It is important to note that in our models d c is not a constant, as it is treat as variable and it is dynamically determined during rupture itself (Fig. 9c). The weakening distance d c has been widely used in fault studies and often imposed as a constant (e.g., Tse and Rice, 1986;Ben-Zion and Rice, 1995;Lapusta et al., 2000;Rubin and Ampuero, 2005;Lapusta and Liu, 2009). However, ample observations suggest that d c should be treated as a variable (Cocco and Bizzarri, 2002;Nielsen et al., 2010), since the effective weakening distance depends on slip history and loading conditions (e.g., Guatteri et al., 2001; Figure 10: Compilation of fracture energy (G ) and slip (δ) from different earthquakes worldwide, from borehole microseismicity to great earthquakes (Abercrombie and Rice, 2005;Rice, 2006;Malagnini et al., 2014;Viesca and Garagash, 2015;Nielsen et al., 2016a;Tinti et al., 2005). Black dashed lines indicate the different scaling for small (G ∝ δ 2 ) and large (G ∝ δ 2/3 ) earthquakes. Tinti et al., 2004). Furthermore, laboratory measures of d c using rotary friction experiments performed at seismic slip rates (∼ 1 m s −1 ) yield weakening distances of the order of meters (Nielsen et al., 2016a, and references therein), substantial friction drops (Di Toro et al., 2011), and fracture energies arguably in the same order as the seismological estimates (Nielsen et al., 2016b).

Future avenues for research
Future studies will make use of this methodology to investigate a number of different phenomena, either as a natural process or induced by human activities. The proposed setup is rather simple, but it can be applied to investigate, e.g., fluid-driven swarms and foreshock sequences (e.g., Ross et al., 2020;Miller, 2020). The developed methodology presents a promising avenue for investigating slip that occurs in response to fluid injection into a permeable fault governed by either rate-strengthening yield strength or rate-and-state friction (e.g., Yang and Dunham, 2021). Future work will also focus on more realistic scenarios of geologic faults, which will allow to investigate the full spectrum of slip behaviour observed in subduction zones and strike-slip faults (e.g., Jolivet and Frank, 2020;Behr and Bürgmann, 2021).
Over the last decades, seismic and geodetic observations have revealed complex interactions of earthquakes, slow-slip transients, and stable creep (Bürgmann, 2018). Investigating whether pore-fluid diffusion affects variations in earthquake recurrence intervals, or fault tremor sequences with alternating longer and shorter repeat times (so-called period-doubling events) (Shelly, 2010), represents an important topic in earthquake physics. Also, systematic investigation of the conditions under which slow-slip events are mostly controlled by visco-plastic bulk properties (e.g., Viesca and Dublanchet, 2018;Lambert and Barbot, 2016;Gao and Wang, 2017;Fagereng and Beall, 2021;Dal Zilio et al., 2022), rather than fault frictional properties (e.g., Liu and Rice, 2007;Dal Zilio et al., 2020), is an important frontier in earthquake source modeling. Further developments can also incorporate fault dilation, temperature, shear heating, thermal runaway, and grain size evolution (e.g., Austin and Evans, 2007), which is needed to investigate the onset of highly localized viscous creep in pre-existing, fine-grained shear zones. As temperature may vary markedly with shear heating (Viesca and Garagash, 2015), the temperature-weakening behavior of steady-state friction may offer the conditions for episodic slow-slip transients, even under rate-strengthening fault behavior.

Conclusions
We have developed a finite difference method to account for fully coupled solid-fluid interaction over many earthquake cycles. The computational framework can model rate-strengthening plasticity in a poro-viscoelasto-plastic rheology. Numerical results are verified through convergence tests and comparisons with analytical benchmarks of pore-fluid pressure diffusion from an injection point along a finite fault width. Future work will include a deeper exploration of parameter space, including the effects of permeability, rate-strengthening, and viscosity, as our set of parameters were chosen primarily for efficiency of computation. For the parameter study in this work, we found that pore-fluid pressure, porosity, and solid-fluid compressibility, influence the occurrence of seismic and aseismic slip.
We have investigated spontaneously occurring shear instabilities on a mildly rate-strengthening fault zone with solid-fluid coupling. Remarkably, these instabilities are fundamentally different from standard instabilities with rate-and-state friction in that they are controlled by localized (de)compaction of pores and dynamic self-pressurization of fluids inside the undrained fault zone. Simulations show that these fluid-driven dynamic ruptures are controlled by solitary pulse-like pressure waves propagating at seismic speed. We have proposed a conceptual model for how this type of instability manifest in a 2-D planestrain shear model (mode-II). However, we recognize that in order to make a full comparison to geological settings we need to understand how the reported instabilities manifest in three dimensions. Despite this limitation, our work demonstrates how pore-fluid pressure, poroelastic effects, and inelastic deformation, can destabilize active faults and produce the full spectrum of slip, from fast-to slow-slip, consistent with observations.
Although multiple weakening mechanisms may operate on active faults, our results suggest that selfpressurization of fluid-filled porous rocks may be the dominant contributor to fault fracture energy, particularly for large earthquakes. Understanding how faults lose their shear strength during fluidinduced dynamic rupture is critical, as it may help to constrain the minimum level of shear stress a fault requires to become unstable. In a broader context, this study shows the importance of incorporating the realistic hydro-mechanical structure of faults to investigate sequences of seismic and aseismic slip. In particular, this study indicates that pore-pressure evolution can completely change both the failure process on the interface and the long-term slip history of geologic faults.

Declaration of Competing Interest
None.

Acknowledgments
We thank Nadia Lapusta, Massimo Cocco, Paul Selvadurai, Elias Heimisson, Federico Ciardo, Domenico Giardini, and Casper Pranger for discussions. This research was supported by the European Research Council (ERC) Synergy Grant FEAR (ID: 856559). Numerical simulations were performed on ETH cluster Euler.

Data availability
The data underlying this article will be shared upon request to the corresponding author.

Code availability
Computer code used within the manuscript is still under development and is available from the corresponding author upon request.