https://orcid.org/0000-0002-8864-4454
SUMMARY
Although pore fluid has received increased attention in the problems of earthquake rupture, its effects on coseismic rupture are not fully understood. In this paper, we numerically simulate dynamic rupture in saturated poroelastic media to study the interplay between poroelastic pressurization and rupture evolution. Particular attention is paid to the physical process of the supershear transition. The spontaneous rupture is assumed to be governed by slip weakening law and the fault is treated as a leaky interface where pore pressure equilibrium is delayed. Results show that the pore pressure induced by coseismic mechanical deformation can strengthen the fault on the extensional side but weaken the fault on the compressional side. As a result, the poroelastic pressurization on the compressional side can promote the occurrence of the transition to supershear rupture while its strengthening effect on the extensional side can suppress this transition even in the state of high-initial shear stress. Meanwhile, it is found that the permeability distribution along the fault plane also affects rupture evolution. Faults with nucleation zones located in low-permeability regions are more likely to undergo a supershear rupture.
1 INTRODUCTION
Non-steady shear rupture along a fault results in the generation of seismic radiation and strong ground motion (Dunham & Bhat 2008). Ruptures can often propagate at a speed that exceeds the transverse wave speed for several large strike-slip earthquakes, including the 2001 Kunlunshan earthquake (Bouchon & Vallée 2003), the 2018 Palu earthquake (Socquet et al. 2019), etc. The occurrence of such supershear rupture has been confirmed in laboratory studies (Xia et al. 2004). Earthquakes with fast rupture speeds can radiate larger energy and cause bigger damage (Hu et al. 2020). Meanwhile, more rapid rupture propagation implies less time is left for early warning of an earthquake. Exploring the underlying physical mechanisms of supershear transition is important in probing fault behaviour and performing seismic hazard analysis.
The speed and features of ruptures depend on many factors in earthquake mechanics. A well-recognized factor controlling the fault rupture is the initial stress on the fault plane. As the rupture spreads, a peak in shear stress will initiate a daughter crack ahead of the main rupture if the initial stress exceeds a critical value (Burridge 1973). The daughter crack is unstable, and the rupture eventually jumps to a supershear speed. Andrews (1976) proposed a non-dimensional seismic ratio S to characterize the initial stress and to quantify whether a daughter crack appears. |$S{\rm{ = (}}{\tau }_\mathrm{ p}{\rm{ - }}{\tau }_{\rm{0}}{\rm{)/(}}{\tau }_{\rm{0}}{\rm{ - }}{\tau }_\mathrm{ r}{\rm{)}}$|, in which |${\tau }_{\rm{0}}$| is the initial stress on the fault plane, |${\tau }_\mathrm{ p}$| is the peak strength and |${\tau }_r$| is the residual strength. Transition to supershear rupture would occur if S < 1.77 in two dimensions (Andrews 1985) and S < 1.19 in three dimensions (Dunham 2007). Another factor that can induce a transition to supershear speeds is the presence of the ground surface (Zhang & Chen 2006; Kaneko & Lapusta 2010; Xu et al. 2015). For earthquakes that rupture to the ground surface, the propagation of rupture will strongly be affected by the interaction with the traction-free mechanical boundary. According to Burridge (1973) and Andrews (1976), supershear transition related to non-uniformly distributed prestress and fault strength is also discussed by a large number of studies (Fukuyama & Olsen 2002; Lapusta & Liu 2009). In recent years, considerable advances have been made in experimental and theoretical studies of supershear transition. Especially, a motion equation for characterizing supershear frictional rupture fronts has been proposed by Kammer et al. (2018), upon which rupture propagation speeds in laboratory experiments can be well predicted and described.
Since rocks are typical saturated porous media, the interaction between solid frame and pore fluid would strongly affect fault behaviour (Jha & Juanes 2014; Bhattacharya & Viesca 2019). During earthquakes, strong coupling between fluid and rocks causes pore pressure changes, which in turn affects the propagation process of rupture (Rudnicki & Rice 2006; Jha & Juanes 2014; Yang & Juanes 2018; Heimisson et al. 2019; Pampillón et al. 2023). Early studies (Rudnicki & Rice 2006; Dunham & Rice 2008) on steady-state sliding between dissimilar fault planes showed that pore pressure changes caused by slipping can change effective compressive stress on the slip surface, and the fault can be either weakened or strengthened. Viesca et al. (2008) further performed numerical simulations of rupture in poroelastic blocks and found large increases in pore pressure on the compressional side of the fault. Researches on frictional ruptures or mode II shear sliding have identified that the induced pore pressure is a significant factor influencing earthquake rupture (Rudnicki & Koutsibelas 1991; Jha & Juanes 2014; Pampillón et al. 2018, 2023; Torberntsson et al. 2018; Yang & Juanes 2018; Heimisson et al. 2019, 2021). Both thermal pressurization and poroelastic pressurization can contribute to rupture-induced pore pressure. For thermal pressurization, pore pressure is generated by released heat due to friction between two planes of a fault. Many works have emphasized the importance of thermal effects (Rice 2006; Yao et al. 2023). In contrast, poroelastic pressurization can be generated once the fault undergoes coseismic mechanical deformation. In this study, we focus on the effects of poroelastic pressurization and also emphasize the influences of fault permeability on rupture evolution.
Studies on fault structures (Chester & Chester 1998; Wibberley & Shimamoto 2003; Delle et al. 2017) show that the majority of slip is localized within a core, which consists of a narrow gouge surrounded by a zone of foliated cataclasite. The permeability of the gouge can be much lower (typically 4 to 6 orders of magnitude) than the external media (Caine et al. 1996). The extremely low permeability of the fault zone core would influence the fault behaviour by preventing fluid from flowing across it. Thus, it is essential to investigate the impact of low-permeability fault zone core on the evolution of fault slipping. The effects of internal structures on nucleation were studied by Heimisson et al. (2019), who showed that slow slip events can nucleate under mildly rate-strengthening friction. In consideration of fault permeability structure, Fan et al. (2019) studied the stability of the fault during fluid injection. The model agrees primarily with realistic situations, but the simulation is time consuming in meshing and computing. Song & Rudnicki (2017) proposed a theoretical description of the thin core of the fault zone with arbitrary permeability in characterizing a plane-strain shear dislocation propagating in a poroelastic solid. Through a linear approximation of the pressure in the fault zone core, the fault is modelled as a leaky plane across which the pore pressure can be equilibrated over a finite duration. On this basis, Heimisson et al. (2021) further introduced a bi-linear approximation that permits significant changes in pore pressure within the layer and conducted a range of important research on fault stability with dilatancy coupled to the poroelastic bulk (Heimisson et al. 2021, 2022; Heimisson 2024).
Even though previous studies have extensively explored the role of poroelasticity in the propagation of ruptures, the effects of fluid-structure interaction on coseismic rupture evolution have not been fully understood. Under the assumption that the fault zone core is impermeable, Mu et al. (2023) performed a dynamic rupture simulation. They found that rupture would transition into a supershear process in particular circumstances. A more elaborate study on the role of pore fluids in rupture propagation was presented by Pampillón et al. (2023) under the undrained assumption. By treating the fault as an impermeable plane, they investigated and emphasized the impact of poroelastic constants of rocks (the effective stress coefficient and the storage modulus) on rupture propagation. Meanwhile, the condition for supershear transition on an impermeable plane in fully saturated rocks was also recharacterized by an effective seismic ratio |$S^{\prime}$|. Natural faults often possess different permeability for fluid flow across the fault zone core (Yao et al. 2023). If the fault is sealed by the fault gouge of low permeability (e.g. mature fault), it is reasonable to idealize the fault zone core as an impermeable surface (Rudnicki & Koutsibelas 1991; Heimisson et al. 2019; Mu et al. 2023; Pampillón et al. 2023). In contrast, for newly formed faults in intact rocks, pores of the two sides of the fault plane are connected and hence the fault zone core is permeable. To further understand the coupled rock-fluid effects on rupture propagation, the influence of permeability of the fault zone core should be investigated and detailed analyses of the rupture evolution related to poroelasticity need to be done.
Theoretical analysis by Heimisson et al. (2021, 2024) has stressed the importance of the permeability of the fault zone core on fault stability. There are three characteristic timescales that can be used to characterize the fluid diffusion during the frictional sliding of a fault, namely the timescale of frictional sliding, diffusion in poroelastic bulk, and flux across the fault zone core. The timescale of frictional sliding during rupture propagation can be defined by |${t}_{\mathrm{ rup}}{\rm{ = }}L{\rm{/}}{V}_{\mathrm{ rup}}$|, where L is the rupture length and |${V}_{\mathrm{ rup}}$| is the rupture speed. The diffusion timescale is |${t}_\mathrm{ b}{\rm{ = }}{L}^2{\rm{/(4}}{\pi }^2c{\rm{)}}$|, where c is the hydraulic diffusivity coefficient. The timescale of flux across the fault zone core is |${t}_\mathrm{ f}{\rm{ = (}}{\kappa }^{\rm{2}}{\varepsilon }^{\rm{2}}{\rm{)/(}}\kappa _{\textrm{layer}}^{\rm{2}}c{\rm{)}}$|, where |$\kappa $| is the Darcy's permeability in poroelastic bulk, |$\varepsilon $| and |${\kappa }_{\textrm{layer}}$| are the half thickness and Darcy's permeability of the fault zone core, respectively. Non-dimensional timescale |${t}_{\mathrm{ rup}}{\rm{/}}{t}_\mathrm{ b}{\rm{ = (}}4{\pi }^2c{\rm{)/(}}L{V}_{\mathrm{ rup}}{\rm{)}}$| is small due to the small value of diffusivity coefficient c and large value of the rupture length L and rupture speed |${V}_{\mathrm{ rup}}$|. This means that pore fluid has no time to penetrate the bulk through diffusion in the short time of rupture propagation. The second non-dimensional timescale characterizing flux through the shear zone core during rupture propagation is |${t}_{\mathrm{ rup}}{\rm{/}}{t}_\mathrm{ f}{\rm{ = (}}\kappa _{\textrm{layer}}^{\rm{2}}cL{\rm{)/(}}{\kappa }^{\rm{2}}{\varepsilon }^{\rm{2}}{V}_{\mathrm{ rup}}{\rm{)}}$|. If |${\kappa }_{\textrm{layer}}$| is much lower than |$\kappa $|, |${t}_{\mathrm{ rup}}{\rm{/}}{t}_\mathrm{ f}$| will be small compared to unity, which means the propagation of rupture occurs much faster than the fluid flux in the thin layer, and this corresponds to the limit where the fault is sealed and impermeable. On the other hand, if |${t}_{\mathrm{ rup}}{\rm{/}}{t}_\mathrm{ f}$| is large, the fluid flux can occur in a shorter time than the rupture process, and the fault will behave as permeable. Therefore, the influence of fault permeability on rupture propagation deserves further detailed study.
The place where pore pressure is used to determine the effective normal stress remains largely unsolved. It is a complex multiscale problem (Viesca et al. 2008; Heimisson et al. 2022) because the determination of the relevant pore pressures requires upscaling from the small scale of poromechanical processes inside the fault zone to the meter and kilometer scales. Heimisson et al. (2021, 2022, 2024) selected the average pore pressure across the shear zone to calculate the effective normal stress in fault stability analyses. Numerical simulations presented by Yang & Juanes (2018) and Jha & Juanes (2014) showed that the onset of slip occurs where pore pressure is higher and thus suggested that the maximum pressure on both sides controls the slip. Accordingly, Pampillón et al. (2023) adopted the criterion of maximum pressure of both fault sides in their dynamic rupture simulation. Field observations have suggested that the principal slip zone of faults often lies at the boundary of the fault zone core (Chester et al. 1993; Dor et al. 2006). Correspondingly, the relevant pore pressure is often taken to be the value at a boundary of the fault zone core (Rudnicki & Koutsibelas 1991; Rudnicki & Rice 2006; Dunham & Rice 2008; Heimisson et al. 2019). In this work, our key task is to discuss the effects of poroelasticity and permeability on rupture evolution instead of clarifying where pore pressure should be used for computing effective stress. We follow Heimisson et al. (2019) by assuming the slip is localized at the negative side of the layer, as depicted in Fig. 1(a), so that the pore pressure used to compute the effective normal stress is selected as p = |${p}^{\rm{ - }}$|. One advantage of this configuration is that it allows us to investigate different impacts of poroelasticity on fault rupture propagation when the fluid-saturated rocks suffer either compressive or extensive deformations.
(a) Schematic of the low-permeability fault zone core with slip occurring on one side of the layer. Compression/extension of the near-fault rocks causes pressure changes during bidirectional propagation of rupture. (b) 2-D fault model for spontaneous rupture simulation. The red segment in the centre represents the 3-km-long nucleation patch. The computational region is surrounded by absorbing boundaries. (c) Schematic representation of pore pressure distribution in the thin layer.
In this paper, we numerically study rupture propagation in fluid saturated poroelastic media. We begin by introducing the governing equations of poroelasticity, providing a model description, and presenting the theoretical treatment of the fault zone core. In this preliminary exploration, analyses mainly focus on the effects of poroelastic pressurization in the plane strain state and intentionally leave out the influence of the ground surface. Particular attention is paid to the impact of the poroelastic pressurization and permeability structure on supershear transition. Considering the spatial heterogeneity of the permeability, we also discuss scenarios where the rupture propagates along a piecewise permeable plane. Finally, conclusions are given.
2 FORMULATION AND MODEL DESCRIPTION
2.1 Representation of the poroelasticity
The coupling between pore fluid and solid frame in fluid saturated rocks can be macroscopically characterized by Biot's theory of poroelasticity (Biot 1941, 1962). For an isotropic poroelastic material, the constitutive equations are
where |${{\boldsymbol \sigma }}$| is the stress tensor and p is the pore pressure. The sign convention for normal stresses is that tensile stress is positive and compressive stress is negative. |$\lambda $| and G are the drained Lamé constants, and |${{\bf u}}$| is the displacement vector of the solid phase. As defined by Biot (1941), |$\zeta $| is the amount of fluid volume entering unit volume of the solid frame, and it is related to the filtration displacement w through |$\zeta {\rm{ = - }}\nabla \cdot {{\bf w}}$|. The filtration displacement is represented by |${{\bf w}}\,{\rm{ = }}\,\phi {\rm{(}}{{{\bf u}}}_\mathrm{ f}{\rm{ - }}{{\bf u}}{\rm{)}}$|, where |${{{\bf u}}}_\mathrm{ f}$| is the average pore fluid displacement and |$\phi $| is the rock porosity. M and |$\alpha $| are the Biot storage modulus and Biot–Willis coefficient, respectively.
in which |${K}_\mathrm{ s}$| and |${K}_\mathrm{ f}$| are the bulk moduli of solid grains and pore fluid, respectively, and
where |${K}_\mathrm{ d}$| is the drained bulk modulus of the rock skeleton. M represents the pressure change caused by fluid volume change at a constant volumetric strain.
The motion equation which relates stresses to the kinematic quantities is
where |$\rho {\rm{ = (1 - }}\phi {\rm{)}}{\rho }_\mathrm{ s}{\rm{ + }}\phi {\rho }_\mathrm{ f}$| is the overall density, |${\rho }_\mathrm{ f}$| is the pore fluid density, and |${\rho }_\mathrm{ s}$| is the solid density. The double overdots on field symbols u and w represent the second derivative with respect to time t.
The fluid movement in porous rocks is characterized by the generalized Darcy's law
where |$\kappa $| denotes Darcy's permeability, |$\eta $| the fluid viscosity. |$\rho ^{\prime} = {\rho }_\mathrm{ f}/\phi + {\rho }_\mathrm{ a}/{\phi }^2$| is a shorthand notation and |${\rho }_\mathrm{ a}\,{\rm{ = (}}{\tau }_\mathrm{ t}{\rm{ - 1)}}\phi {\rho }_\mathrm{ f}$| is a so-called added mass density with |${\tau }_\mathrm{ t}$| denoting the tortuosity.
2.2 Computational model and fault plane treatment
We study the impact of fault zone core permeability on rupture and explore its influences on the rupture speed transition in fully fluid-saturated poroelastic media. The numerical method developed by Mu et al. (2023) is employed to simulate the spontaneous propagation of rupture along a finite fault. They use a staggered-grid finite-difference method with governing equations of first-order velocity–stress formulations and split-node technique (Dalguer & Day 2007) to treat discontinuous components of physical quantities on the fault plane. As depicted in Fig. 1(b), rupture is allowed on a 40-km-long fault segment along the x-axis. This fault segment is embraced in a homogeneous poroelastic bulk. Spatial boundaries of the computational domain are described by absorbing boundaries (Martin et al. 2008). We assume that the fault slip occurs on an interface between the lower side of the layer and the host material. For the sake of computational efficiency, we take advantage of the antisymmetry of the mode-II problem (Song & Rudnicki 2017; Heimisson et al. 2019) so that only the poroelastic bulk at the negative side is modelled. Parameters of the poroelastic bulk are listed in Table 1. In numerical computation, a uniform grid interval of 50 m is adopted, the time interval is 1 ms, and the time duration is 9 s.
Fluid saturated porous medium parameters in spontaneous rupture simulation.
| Definition | Unit | Symbol | Value |
|---|---|---|---|
| Porosity | – | |$\phi $| | 0.14 |
| Permeability | Darcy | |$\kappa $| | 0.23 |
| Drained bulk modulus | GPa | |${K}_\mathrm{ d}$| | 21.09 |
| Frame shear modulus | GPa | G | 18.90 |
| Solid bulk modulus | GPa | |${K}_\mathrm{ s}$| | 35.70 |
| Pore fluid bulk modulus | GPa | |${K}_\mathrm{ f}$| | 2.25 |
| Solid density | kg m−3 | |${\rho }_\mathrm{ s}$| | 2650 |
| Pore fluid density | kg m−3 | |${\rho }_\mathrm{ f}$| | 1000 |
| Pore fluid viscosity | Pa s | |$\eta $| | 0.001 |
| Tortuosity | – | |${\tau }_\mathrm{ t}$| | 2.67 |
| Definition | Unit | Symbol | Value |
|---|---|---|---|
| Porosity | – | |$\phi $| | 0.14 |
| Permeability | Darcy | |$\kappa $| | 0.23 |
| Drained bulk modulus | GPa | |${K}_\mathrm{ d}$| | 21.09 |
| Frame shear modulus | GPa | G | 18.90 |
| Solid bulk modulus | GPa | |${K}_\mathrm{ s}$| | 35.70 |
| Pore fluid bulk modulus | GPa | |${K}_\mathrm{ f}$| | 2.25 |
| Solid density | kg m−3 | |${\rho }_\mathrm{ s}$| | 2650 |
| Pore fluid density | kg m−3 | |${\rho }_\mathrm{ f}$| | 1000 |
| Pore fluid viscosity | Pa s | |$\eta $| | 0.001 |
| Tortuosity | – | |${\tau }_\mathrm{ t}$| | 2.67 |
Fluid saturated porous medium parameters in spontaneous rupture simulation.
| Definition | Unit | Symbol | Value |
|---|---|---|---|
| Porosity | – | |$\phi $| | 0.14 |
| Permeability | Darcy | |$\kappa $| | 0.23 |
| Drained bulk modulus | GPa | |${K}_\mathrm{ d}$| | 21.09 |
| Frame shear modulus | GPa | G | 18.90 |
| Solid bulk modulus | GPa | |${K}_\mathrm{ s}$| | 35.70 |
| Pore fluid bulk modulus | GPa | |${K}_\mathrm{ f}$| | 2.25 |
| Solid density | kg m−3 | |${\rho }_\mathrm{ s}$| | 2650 |
| Pore fluid density | kg m−3 | |${\rho }_\mathrm{ f}$| | 1000 |
| Pore fluid viscosity | Pa s | |$\eta $| | 0.001 |
| Tortuosity | – | |${\tau }_\mathrm{ t}$| | 2.67 |
| Definition | Unit | Symbol | Value |
|---|---|---|---|
| Porosity | – | |$\phi $| | 0.14 |
| Permeability | Darcy | |$\kappa $| | 0.23 |
| Drained bulk modulus | GPa | |${K}_\mathrm{ d}$| | 21.09 |
| Frame shear modulus | GPa | G | 18.90 |
| Solid bulk modulus | GPa | |${K}_\mathrm{ s}$| | 35.70 |
| Pore fluid bulk modulus | GPa | |${K}_\mathrm{ f}$| | 2.25 |
| Solid density | kg m−3 | |${\rho }_\mathrm{ s}$| | 2650 |
| Pore fluid density | kg m−3 | |${\rho }_\mathrm{ f}$| | 1000 |
| Pore fluid viscosity | Pa s | |$\eta $| | 0.001 |
| Tortuosity | – | |${\tau }_\mathrm{ t}$| | 2.67 |
If the permeability of the fault zone core is relatively low, the pressure equilibrium across the fault zone core will be delayed. As shown in Fig. 1(c), positive pore pressure will be induced on the compressional side, while negative pore pressure will be induced on the extensional side. In the next section, we will show that the permeability of the fault slip zone is important for understanding how rupture propagates. Typically, the fault slip zone in mature faults consists of the fault gouge formed due to historical slips. The fault gouge is extremely thin and its permeability is much lower than that of the surrounding damage zone and intact host rocks (Wibberley & Shimamoto 2003), which makes it behave as a barrier for fluid flow. To characterize the impact of fault permeability on pore pressure distribution across the fault, we employ a boundary layer approach to model the fault gouge (Song & Rudnicki 2017; Heimisson et al. 2019). In this approach, the fault gouge with a thickness of |${\rm{2}}\varepsilon $|, as shown in Fig. 1(a), is treated as a mathematical zero-thickness interface for the outer region. The permeability of the fault gouge is denoted by |${\kappa }_{\textrm{layer}}$| while a different permeability of the host rock is presented by |$\kappa $|.
In the case of a thin shear zone layer embraced in a homogeneous poroelastic block, the pore pressure at the layer boundary can be formulated by a leaky interface condition (Song & Rudnicki 2017). Since the layer possesses a very small size, the effects of inertia can be disregarded (Rice 2006; Platt et al. 2014). In obtaining the leaky interface condition, quasi-static Darcy's law is employed to govern the fluid flow across the fault. Neglecting the acceleration term in eq. (6), we find the filtration velocity only relates to the pore pressure gradient. Consequently, the conservation of mass in fluid flux on the layer boundary leads to a leaky interface condition
where |${p}^{\rm{ - }}$| denotes pore pressure on the lower boundary of the thin layer, and |$R\,{\rm{ = }\,}\varepsilon \kappa {\rm{/}}{\kappa }_{\textrm{layer}}$| is a hydraulic resistance coefficient (the leaky interface condition can also be formulated by the pore pressure |${p}^{\rm{ + }}$| on the upper boundary and its gradient in the fault-normal direction through |${p}^{\rm{ + }}{\rm{ = }\,}R\partial {p}^{\rm{ + }}{\rm{/}}\partial y$|). Eq. (7) characterizes the fluid flow across the fault by relating pore pressure and its normal gradient in terms of hydraulic properties. Note that fluid flow in the direction parallel to the layer is neglected. In the case of a low-permeability fault zone core, the horizontal flow will be dominated by that in the host material. In this case, the horizontal flow in the fault zone core can be neglected. In the case of a high-permeability fault zone core, induced pore pressure is small so the horizontal flow will be weak. In the second case, the horizontal flow becomes insignificant. If the fault zone core has anisotropic permeability, horizontal flow might be important. Nevertheless, eq. (7) is a simple model for exploring the interaction of slip-induced pore pressure changes and earthquake rupture. The permeability effect of the fault zone core can be described through the hydraulic resistance coefficient R. For |$R{\rm{ = 0}}$|, the fault is fully permeable. In this case, pore pressure is equilibrated instantaneously, and the mechanical responses are equivalent to a fault embraced in elastic media but with undrained elastic constants of the host rocks. For |$R\,{\rm{ = }}\,\infty $|, the fault is impermeable for flux in the normal direction. R can be inferred from fault zone structure and permeability data. For the median tectonic line in southwest Japan reported by Wibberley & Shimamoto (2003), parameters can be selected that |${\kappa }_{\textrm{layer}}$| = 0.003 12 × 10−18 m2, |$\kappa $| = 3430 × 10−18 m2, ε = 0.1 m and thus we can compute R = 1.1 × 105 m. Although eq. (7) was derived given both sides of the layer consist of the same material, we emphasize that its applicability can be extended for bimaterial interfaces when the fault zone core is impermeable (|$R\,{\rm{ = }}\,\infty $|) or fluid pressure is instantaneously equilibrated (|$R\,{\rm{ = 0}\,}$|).
Slip-weakening friction law (Ida 1972; Palmer & Rice 1973) is used to govern the fault rupture. Friction coefficient |${\mu }_\mathrm{ f}$| is assumed to be linearly dependent on the slip distance l through
where |${\mu }_\mathrm{ s}$| and |${\mu }_\mathrm{ d}$| are static and dynamic coefficients of friction, respectively, and |${D}_\mathrm{ c}$| is the critical slip-weakening distance. Typical values |${\mu }_\mathrm{ s}$| = 0.677 and |${\mu }_\mathrm{ d}$| = 0.525 are applied on the whole fault plane and |${D}_\mathrm{ c}$| is chosen as a constant value of 0.4 m in the context of this paper.
Coseismic mechanical deformation induces non-zero pore pressure during rupture propagation. Here, Terzaghi effective stress |$\sigma _n^{\prime}$| = |${\sigma }_n{\rm{ - }}p$|, where p = |${p}^{\rm{ - }}$|, is employed to determine the effective peak strength |$\tau _\mathrm{ p}^{\prime}$| and effective residual strength |$\tau _\mathrm{ r}^{\prime}$| acting on the fault through |$\tau _\mathrm{ p}^{\prime}\,{\rm{ = }}\,{\mu }_\mathrm{ s}\sigma _\mathrm{ n}^{\prime}$| and |$\tau _\mathrm{ r}^{\prime}\,{\rm{ = }}\,{\mu }_\mathrm{ d}\sigma _\mathrm{ n}^{\prime}$|. |${\sigma }_\mathrm{ n}$| is the contact traction in the normal direction and is defined as positive when the fault edges are compressive. The fault begins to slip when the shear stress |${{\rm{\sigma }}}_{xy}$| acting on the fault reaches the effective peak strength |$\tau _\mathrm{ p}^{\prime}$|.
Tectonic stress on the fault plane is the driving force for rupture propagation. Currently, in-site direct measurements of the stress are difficult for deep crustal rocks (Fuchs & Müller 2001) and the understanding of tectonic stress field is insufficient (Shebalin & Narteau 2017). In our simulation, we specify the initial normal stress |${\sigma }_n$| = 120 MPa, and thus its product with friction coefficients determines the peak strength |${\tau }_\mathrm{ p}$| = 81.24 MPa and residual strength |${\tau }_\mathrm{ r}$| = 63 MPa. These values are consistent with theoretically predicted values at a typical seismogenic depth of about 10 km (Aochi & Madariaga 2003) and they have been used in simulations of spontaneous rupture (Dalguer & Day 2007). We assume the initial shear stress to be uniform along the fault and consider its values from 68 to 73 MPa (|${\tau }_{\rm{0}}$| = 68, 71, and 73 MPa), which correspond to earthquakes occurring at a depth of about 10 km with stress drops of 5–10 MPa (McGarr 1984). Besides, in our simulation, the earthquake rupture is nucleated at the center 3-km-long zone within which the shear stress is set to be |${\tau }_{nuc}$| = 81.6 MPa, which is a little higher than the peak strength, and then it propagates spontaneously over the rest of the fault plane.
Initial shear stresses of 68, 71, and 73 MPa correspond to different seismic ratios of 2.65, 1.28, and 0.82, respectively. According to the Burridge–Andrews mechanism (Burridge 1973; Andrews 1976, 1985), the supershear transition would occur on slip-weakening faults if the seismic ratio S is smaller than a threshold of 1.77. There exists a critical rupture length below which supershear propagation is impossible (Mello et al. 2010; Kammer et al. 2018). This critical length has been defined in both a fracture mechanics approach (Kammer et al. 2018) and an analytical-numerical approach given a linear slip-weakening friction law (Mello et al. 2010). Here, the rupture length required for supershear transition on slip-weakening faults can be estimated by the second approach (Mello et al. 2010) as |$L{\rm{ = 9}}{\rm{.8(1}}{\rm{.77 - }}S{{\rm{)}}}^{{\rm{ - 3}}}{L}_\mathrm{ c}$|, where S is the seismic ratio and |${L}_\mathrm{ c}$| is the critical (plane strain) crack half length. |${L}_\mathrm{ c}\,{\rm{ = }}\,[ {G( {{\tau }_\mathrm{ p}{\rm{ - }}{\tau }_\mathrm{ r}} ){D}_\mathrm{ c}} ]{\rm{/}}[ {{\rm{\pi (1 - }}\nu {\rm{)(}}{\tau }_{\rm{0}}{\rm{ - }}{\tau }_\mathrm{ r}{{\rm{)}}}^{\rm{2}}} ]$|, where G is the shear modulus of solid frame, |${D}_\mathrm{ c}$| is the critical slip weakening distance and |$\nu $| is Poisson's ratio (an undrained value is used in our poroelastic model). For a high-initial shear stress of 73 MPa, the rupture length for supershear transition is approximated as 6 km, and the rupture in a plane strain elastic space would transition to a supershear speed very soon. While for initial shear stress |${\tau }_{\rm{0}}$| = 71 MPa, the required rupture length is approximated as 70 km, which is longer than the length of our fault model. Thus, supershear rupture would not emerge in our computational domain at this initial shear stress in the case of a fully permeable fault, despite that the computed seismic ratio S is smaller than the threshold. In this paper, discussions are restricted in our simulation domain, but the possibility of supershear transition is not excluded if rupture propagates a larger length.
3 RESULTS
3.1 Effects of fault leakage on rupture evolution due to poroelastic pressurization
We begin by comparing two limiting cases: (1) the fault zone core is fully permeable (R = 0); (2) the fault zone core is sealed by the fault gouge of extremely low permeability (R = ∞). If the fault zone core is fully permeable, pore pressure is instantaneously equilibrated so that no pore pressure change is induced on the slipping surface during the dynamic propagation of rupture. As a result, the fully permeable fault is equivalent to a fault embraced in an elastic solid with the same undrained elastic constants as the fluid saturated rock. The impermeable fault model is a more common scenario as the fault zone core is usually sealed by fault gouge. Fig. 2 exhibits the wavefield snapshots of the x-component velocity for three different initial shear stresses. At a low stress |${\tau }_{\rm{0}}$| = 68 MPa (S = 2.65), the rupture evolution on an impermeable fault is almost the same as on a permeable fault. The rupture propagation maintains a sub-Rayleigh process on both sides of the nucleation zone. In contrast, when the initial stress increases to 71 MPa (S = 1.28), a Mach cone is observed on the compressional side of an impermeable fault while the rupture on a permeable fault maintains sub-Rayleigh on both sides. Since the required rupture length for S = 1.28 is estimated as 70 km, this result indicates that the transition to supershear rupture is promoted and it occurs much earlier on the compressional side of an impermeable fault. At a high stress |${\tau }_{\rm{0}}$| = 73 MPa (S = 0.82), the rupture on the impermeable fault also evolves to be a single-side supershear rupture. The transition to supershear rupture is suppressed and delayed on the extensional side of an impermeable fault. In contrast, at this very small seismic ratio S of 0.82, the earthquake on a permeable fault would transition into a double-side supershear rupture soon. These results clearly show that the occurrence of the supershear rupture is affected by both initial shear stress and hydraulic properties of the fault plane.
Wavefield snapshots of the x-component velocity for both impermeable and permeable faults with different initial shear stresses: (a) |${\tau }_{\rm{0}}$| = 68 MPa (S = 2.65), (b) |${\tau }_{\rm{0}}$| = 71 MPa (S = 1.28), and (c) |${\tau }_{\rm{0}}$| = 73 MPa (S = 0.82).
To gain further insights into the associated effects of initial shear stress and fault permeability on the supershear transition, we give in Fig. 3 the plot of the rupture evolution diagram for a broad range of the fault hydraulic resistance coefficient |$R$| and seismic ratio S. The region where there is no supershear rupture in our computational domain is shown in ivory colour, while the region where supershear rupture transition occurs is shown in light pink colour. The distance from the coordinate origin to the place where the supershear transition starts is marked out by circles with different colours and sizes at the boundary between these two regions. We can see that a large hydraulic resistance can promote the supershear rupture transition on the compressional side (Fig. 3a). For example, at a given initial shear stress level of S = 1.3, the supershear rupture will emerge in our computational domain if the hydraulic resistance coefficient is up to about R = 103. We note that although S = 1.3 < 1.77, the supershear rupture cannot appear on the fully permeable fault plane in our limited computational domain. This is because for S = 1.3, the required rupture length is approximated as 86 km, much larger than our computational domain. But, the supershear transition distance on the compressional side can be shortened by the increased hydraulic resistance coefficient as the induced positive pore pressure will decrease the frictional strength by reducing the effective normal stress, |$\sigma _{\rm{n}}^{\prime}$| = |${\sigma }_{\rm{n}}{\rm{ - }}p$|. In contrast, on the extensional side, coseismic stretching deformation induces a negative pore pressure, increasing the effective normal contact traction and strengthening the fault. As shown in Fig. 3(b), the transition to supershear rupture is suppressed and delayed, and the transition distance lengthens with increasing hydraulic resistance. If the hydraulic resistance is large enough, there will be no longer supershear transition in our computational domain.
Rupture evolution diagram in the computational domain with changes in fault permeability and seismic ratio S: (a) on the compressional side; (b) on the extensional side. The supershear transition distances from the coordinate origin are marked by circles at the boundary between the supershear rupture and sub-Rayleigh rupture regions. Bigger and brighter circle indicates a larger distance.
When rupture occurs, the fault undergoes either compressional or extensional deformation, which can cause positive or negative pore pressure. The fault zone core with low permeability delays pressure equilibration by impeding flow across it. Fig. 4(a) depicts pore pressure distribution on the slipping surface with initial shear stress |${\tau }_{\rm{0}}$| = 71 MPa at t = 4.2 s for various resistance coefficients. Fig. 4(b) shows the corresponding distribution of shear stresses at these instances. We can see that a higher pore pressure is induced for faults with a smaller permeability value, and the maximum change in pore pressure is localized at the rupture front. This is because for short-term response, pore fluids are trapped and they cannot penetrate the low-permeable fault zone core, leading to a higher pore pressure change in response to the mechanical solid deformation. On the compressional side, the positive pore pressure leads to a drop in effective normal contact traction, facilitating the spread of the rupture. A noteworthy result is that the rupture transitions into a supershear process on an impermeable plane (|$R = \infty $|). This supershear transition can be identified through the distribution of shear stresses, that is the shear stress peak preceding the rupture front disappears after supershear rupture occurs. The earlier occurrence of supershear rupture on an impermeable fault is attributed to poroelastic pressurization. In contrast, as shown in Figs 4(a) and (b), the rupture speed on the extensional side exhibits a little reduction since the arrival of pore pressure peak and rupture front delay slightly. This decrease in rupture speed stems from the fault strengthening induced by the pore pressure change in extensional deformation.
(a) Pore pressure distribution and (b) shear stress distribution along the fault for different hydraulic properties at t = 4.2 s.
We note that for R > 104, the simulation results are almost the same. Since the hydraulic resistance coefficient computed from an actual fault can be R = 1.1 × 105 (see Section 2.2), we will simplify faults as impermeable when discussing the physical mechanisms of rupture transition due to poroelastic pressurization. Unless otherwise specified, the fault zone core is treated as impermeable and the initial shear stress is set to |${\tau }_{\rm{0}}$| = 71 MPa in the following discussion.
3.2 Physical mechanism of supershear transition caused by coseismic poroelastic weakening
To understand the physical mechanism of supershear transition, we further explore how poroelastic pressurization accelerates the rupture. Fig. 5(a) depicts the detailed development of the cohesive zone on the compressional side of an impermeable fault (|$R = \infty $|). The cohesive zone is a small portion behind the rupture tip where shear stress decreases from the peak strength to the residual strength. Solid and dashed lines plot the leading and trailing edges of it. In our simulation, the dimensionless ratio of the cohesive zone width to the grid interval is greater than 3.3 (critical value derived by Day et al. 2005), above which the solution is numerically stable). Slope of the leading edge is related to the rupture speed. As can be seen in Fig. 5(b), the rupture speed gradually approaches the Rayleigh wave speed soon after the rupture initiates and then suddenly changes to a value higher than the transverse wave speed when the rupture propagates about 4 km away from the coordinate origin. From this locally enlarged view, we can see that the leading edge jumps ahead at about 4.0 km and the subsequent rupture speed is greater than the transverse wave speed. The transition is completed at about 6.0 km where the trailing edge also jumps ahead. The travelling leading edge bypasses the forbidden zone within which the speed is between the Rayleigh and transverse wave speeds (also illustrated in Fig. 5c).
(a) Spatio-temporal plot of cohesive zone on an impermeable fault plane (|$R = \infty $|). (b) Locally enlarged view of supershear transition region. (c) Rupture speed and slip rate. Black line marks the Rayleigh speed (Rayleigh) and pink line marks the transverse wave speed (S). (d) Pore pressure and peak strength changes in the same region.
Within the sub-Rayleigh region, the slip rate increases as the rupture propagates. This high-slip rate would cause severe compression because larger solid phase strain would be accumulated within the same time interval if the slip rate is larger. As shown in Fig. 5(d), larger pore pressure change causes stronger fault weakening (Fig. 5d). As soon as the supershear transition finishes, the slip rate drops sharply to a smaller value. Accordingly, the pore pressure also drops, and thus the peak strength rises as depicted in Fig. 5(d). We note that physical quantities in Figs 5(c) and (d) are computed at the rupture front, and there are sudden jumps when supershear transition occurs. Therefore, there is a gap for computed values of these physical quantities between 3.5 and 4 km in Figs 5(c) and (d).
Figs 6(a) and (b) depict time histories of shear stress at 2.5, 3.0, 3.5, 4.0, 5.0, and 6.0 km on two sides of the impermeable fault. The results for a permeable fault are shown in Fig. 6(c). It is shown that during sub-Rayleigh rupture, the shear stress peak on the compressional side of the impermeable fault is higher than results of the permeable fault. In comparison, the shear stress peak on the extensional side of the impermeable fault is lower. We note that the shear stress peak propagates at the speed of the transverse wave. This phenomenon is analogous to that in the elasticity model predicted by Burridge (1973) and observed by Svetlizky et al. (2016).
Time histories of shear stress at different locations (2.5, 3.0, 3.5, 4.0, 5.0, and 6.0 km away from the coordinate origin, respectively): (a) on the compressional side of an impermeable fault; (b) on the extensional side of an impermeable fault; (c) on a fully permeable fault.
To explore the effects of fault permeability on shear stress peak, we further discuss radiated and dissipated energy. In an earthquake rupture, suddenly released strain energy is transferred into energy radiated as elastic waves and energy dissipated to extend the rupture (Cocco et al. 2023; Kammer et al. 2024). In the linear poroelastic model, the dissipated energy is divided into energy density for creating a new rupture surface (fracture energy) and energy density through generating frictional heat on the fault (frictional work). Here, we discuss the shear stress evolution ahead of the rupture front by analysing the fracture energy (Gc) and the radiated energy (Er). The definition of fracture energy density is (Andrews 2005; Cocco et al. 2023)
where |${\sigma }_{xy}$| and |$\tau _\mathrm{ r}^{\prime}$| are the shear stress and the effective residual frictional stress on the fault plane, respectively; |$l$| is the accumulative slip. Eq. (9) prescribes the amount of energy to create a unit area of a rupture. Radiated energy can be estimated by the stress drop |$\Delta \tau = {\tau }_{\rm{0}} - \tau _\mathrm{ r}^{\prime}$| and accumulative slip l(x) at location x of the fault using the following definition (Rivera & Kanamori 2005; Khademian et al. 2020)
where |${L}_{\rm{1}}$| is the rupture length propagating to the left from the nucleation zone and |${L}_{\rm{2}}$| is that to the right. Thus, the total rupture length is |${L}_{\rm{1}}\,{\rm{ + }}\,{L}_{\rm{2}}$|. We aim to study the differences in radiated energy on two sides of the nucleation zone. To this end, we partition the radiated energy in a simple way that the radiated energy due to the rupture on the left side is represented by |${E}_\mathrm{ r}\,{\rm{ = }\,}\frac{{\rm{1}}}{{\rm{2}}}\mathop \smallint \nolimits_{{\rm{ - }}{L}_{\rm{1}}}^{\rm{0}} \Delta \tau \cdot l{\rm{(}}x{\rm{)d}}x$| while that on the right side is defined by |${E}_\mathrm{ r} = \frac{1}{2}\int_{0}^{{{L}_2}}{{\Delta \tau \cdot l(x){\rm{d}}x}}$|. We note that the fracture energy density is calculated at a specific point on the fault, while radiated energy can be estimated in terms of the whole rupture surface (Cocco et al. 2023). To exhibit energy variation when rupture propagates to different locations, we calculate the radiated energy every simulation time step and obtain the total radiated energy curve as a function of rupture lengths.
Fig. 7(a) compares fracture energy distribution between the permeable and impermeable fault planes. It is shown that the fracture energy on the compressional or extensional side of an impermeable fault is lower or higher than that on the permeable fault. On the impermeable fault, we find a sharp rise in fracture energy on the compressional side, and its location implies the one where supershear rupture emerges. This phenomenon indicates that more energy is required to create a new surface when the supershear transition begins. At the location where supershear rupture appears, we also plot the energy release rate predicted by the theory of Kammer et al. (2018) in the figure. The energy release rate |${G}_\mathrm{ r}$| can be expressed in terms of rupture speed |${V}_\mathrm{ r}$|, shear stress on the fault |${\tau }_{\rm{0}}$|, and the rupture half-length l through (Broberg 1994; Kammer et al. 2018)
where |$\mu $| is the shear modulus. |$B{\rm{(}}{V}_\mathrm{ r}{\rm{/}}{V}_\mathrm{ P}{\rm{)}}$| and |${\Gamma }_{\rm{D}}{\rm{(}}\gamma {\rm{)}}$| are known functions given in Broberg (1994) (eqs 69 and 70 in Broberg 1994), and they are not shown here for brevity). The required rupture speed for the calculation is extracted from numerical simulation. Kammer et al. (2018) proposed a crack propagation criterion by balancing the energy release rate and the fracture energy and applied it to compute the steady state supershear rupture speed. In Fig. 7(a), it can be seen that the obtained energy release rate is almost consistent with the fracture energy calculated by eq. (9). This means that for ruptures propagating dynamically, there is also a balance between the energy release rate and the fracture energy when supershear rupture initiates. A minor difference is attributed to different treatments of cohesive zone and numerical errors.
(a) Fracture energy distribution along the fault planes. Five-pointed star marks the energy release rate calculated through the theory of Kammer et al. (2018). (b) Accumulative energy radiated through elastic waves as rupture lengths. (c) Accumulative slip along permeable and impermeable faults with a time interval of 1.2 s between adjacent curves.
Fig. 7(b) shows the total radiated energy on two sides of the nucleation zone. The comparison with a permeable fault shows that radiated energy from an impermeable fault is larger on the compressional side but smaller on the extensional side. Rupture will become asymmetric if the fault zone core is impermeable. This asymmetry can be seen in the wavefield snapshots (Fig. 2b) and slip distribution (Fig. 7c). On the extensional side, the decreased slip and negative pore pressure, which increases the residual strength and thus reduces stress drop, can cause a decrease in radiated energy. In contrast, positive pore pressure on the compressional side can decrease the effective normal stress and increase the stress drop, which leads to an increase in radiated energy. Since the shear stress peak on the slipping surface is formed under the combination of radiated longitudinal and transverse waves (Dunham 2007), the increase in radiated energy accounts for the faster increase in the amplitude of shear stress peak (Fig. 6), which eventually promotes earlier transition to a supershear rupture process. Because of the combined effects of the variation in stress drop and the disturbance in cumulative slip, the radiated energy exhibits a sharp drop at the location where supershear rupture transition completes. Subsequently, the radiated energy increases again as the rupture continues to lengthen. These results of fracture energy density and radiated energy show that the coseismic poroelastic pressurization can change the energy budget of a rupture event. Specifically, on the compressional side of an impermeable fault, rupture requires lower fracture energy and more energy is radiated forward through elastic waves, increasing the shear stress peak ahead of rupture front and promoting an earlier occurrence of supershear rupture.
The permeability distribution can be non-uniform even on the same fault plane (Caine et al. 1996; Mizoguchi et al. 2008). Therefore, in the following, we also exhibit and discuss the effects of the heterogeneity of fault zone core permeability on rupture propagation. For simplicity, we confine the discussion to two typical piecewise permeable scenarios (Fig. 8a): (1) The nucleation zone is located within a zone with impermeable fault surface while the fault surface of the exterior process zone is permeable; (2) the nucleation zone is located within a zone with permeable fault surface while the fault surface of the exterior process zone is impermeable.
(a) Schematic diagram of permeability distribution for two piecewise permeable scenarios. (b) Wavefield snapshots of the x-component velocity for these two scenarios. (c) Distribution of fracture energy along piecewise permeable faults.
The wavefield snapshots for the two scenarios are shown in Fig. 8(b). For scenario (1), the supershear rupture occurs on the compressional side of the impermeable fault segment. After it intrudes the permeable region, the rupture keeps on propagating at a supershear speed. For scenario (2), however, no supershear transition is observed even if the rupture propagates into the exterior impermeable region. The supershear transition is delayed. In Fig. 8(c), we plot fracture energy density distribution along the fault for the two scenarios. The comparison between the two cases shows that the fracture energy density of the exterior impermeable plane of scenario (2) is even much smaller than that of the interior impermeable plane of scenario (1). However, the supershear transition is not observed in the impermeable process zone in scenario (2). At this point, a question immediately appears that why is the supershear rupture delayed given a lower value of fracture energy?
We further plot in Fig. 9 the shear stress time histories received at 5 and 19 km on the fault plane of the compressional side for scenario (2). Corresponding results on a uniformly permeable fault are also shown in dashed lines. It can be seen that the shear stress peak runs away from the rupture front after the rupture propagates for a long distance. For scenario (2), when rupture intrudes into the impermeable zone, fracture energy should decrease immediately due to coseismic poroelastic pressurization, increasing the radiated energy. However, since the shear stress peak has propagated for a long distance of 10 km, the new portion of forward radiated energy cannot catch up with the shear stress peak. Consequently, the amplitude of the shear stress peak cannot be increased efficiently. Instead, the extra radiated energy produces a sidelobe in the shear stress. Above results indicate that the heterogeneity of the permeability on the slipping surface can significantly affect the supershear transition. The supershear rupture is readier to occur if the fault zone core around the nucleation zone is impermeable.
Shear stress time histories received at 5 and 19 km on the compressional side for piecewise permeable scenario (2) (black solid lines) and results received at corresponding locations on a permeable fault are also shown (orange–yellow dashed lines).
4 CONCLUSIONS
The effects of fault zone core permeability on the rupture evolution in poroelastic rocks are explored in this paper and the physical mechanism of supershear transition induced by coseismic poroelastic pressurization is investigated. Coseismic pore pressure on the compressional side of a low-permeable plane can decrease the effective normal contact traction and weaken the fault by reducing the effective fault frictional strength. Correspondingly, the fracture energy required at the rupture front for forming a new rupture surface decreases and more energy is radiated forward through elastic waves. As a result, the amplitude of the shear stress peak increases, instigating the supershear transition. In contrast, on the extensional side, coseismic stretching deformation induces a negative pore pressure which strengthens the fault by increasing the effective normal contact traction. Rupture dissipates more energy so that radiated energy declines. This results in a decrease in the amplitude of shear stress peak preceding the rupture front, suppressing the supershear transition. Consequently, on this extensional side, supershear rupture would be delayed on a low-permeable plane even if the initial shear stress is very high. It is found that the permeability distribution on the fault plane also influences supershear transition. For piecewise permeable fault plane, supershear rupture seems readier to occur if permeability of the fault zone core around the nucleation zone is low.
ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China under Projects 42074057, 12272107, 42474074, and 12472087.
DATA AVAILABILITY
All data used have been given in this article. FORTRAN source code used for calculation in this study is available at https://doi.org/10.6084/m9.figshare.27291423.
