Fault friction and Mantle viscosity
The stages of the seismic cycle in faults can be described using constitutive relationships of rate – and state – friction that explain the evolution of sliding velocity over time. These laws had been used to model several ruptures styles in nature69 and to characterize the seismic cycle at different tectonic settings. This leads to the possibility of understanding fast and slow ruptures at different fault configurations. Here, we used a constitutive framework obtained from the assumption of a micro-physical model of the rate – and state – dependent friction under isothermal conditions. In this context, the sliding velocity depends on the density of the real area of contact as follows70:
$$V = {V}_{0}{\left(\frac{\tau }{\mathcal{A}\chi }\right)}^{\frac{{\mu }_{0}}{a}}exp\left[-\frac{Q}{R}\left(\frac{1}{T}-\frac{1}{{T}_{0}}\right)\right]$$
1
where \(V\) is the sliding velocity, \(\tau\) is the norm of the shear stress resolved on the fault plane, \(T\) is the absolute temperature, \({\mu }_{0}\) is the static coefficient of friction at \({V}_{0}\) and \({T}_{0}\), \(\mathcal{A}\) is the real area of contact density (real area of contact divided by the nominal surface area), \(\chi\) is the indentation hardness, \(a\) << 1 is the direct effect, \(Q\) is the activation energy and \(R\) is the gas constant. The real area of contact depends on the effective normal stress and changes as a function of the shape of the surrounding contact junctions71. This last feature is characterized by a state variable that represents the age of the grain which coincides with the age of contact defined by reference72. The relationship between the effective normal stress and the real area of contact is given by70:
$$\mathcal{A} = \frac{{\mu }_{0}{\sigma }_{eff}}{\chi }{\left(\frac{\theta {V}_{0}}{L}\right)}^{\frac{b}{{\mu }_{0}}}$$
2
where \({\sigma }_{eff}\) is the effective normal stress, \(\theta\) is the state variable, \(L\) is the characteristic weakening distance associated with the gouge thickness or fault roughness and \(b\) << 1 is the evolution effect. By combining Eqs. (1) and (2) the multiplicative form of rate – and state – friction at isothermal conditions is obtained70,73:
$$\tau ={\mu }_{0}{\sigma }_{eff}{\left(\frac{V}{{V}_{0}}\right)}^{\frac{a}{{\mu }_{0}}}{\left(\frac{\theta {V}_{0}}{L}\right)}^{\frac{b}{{\mu }_{0}}}$$
3
The state variable under isothermal conditions follows an evolution law which allows healing at the stationary contacts5,70:
$${\theta }^{\cdot }=exp\left[-\frac{H}{R}\left(\frac{1}{T}-\frac{1}{{T}_{0}}\right)\right]-\frac{\theta V}{L}$$
4
Fault slip dynamics is controlled by the spatial distribution of \({\mu }_{0}\), \({\sigma }_{eff}\), \(a\), \(b\) and \(L\), which are governed by the thermal state at the fault interface and by the mineral assembly of the respective stable lithological facie under that state. We assume that those parameters do not change over time and are not affected by heat and fluid transport, which means that the value that controls the frictional behavior on the fault at steady state will be \(a-b\). Depending on its temperature distribution, effective normal stress and rock composition, any given fault can experience a velocity weakening (\(a-b<0\)) or a velocity strengthening (\(a-b>0\)) behavior74. Velocity weakening areas will experience unstable slip, while velocity strengthening regions should slip aseismically75. Nonetheless, previous works had shown that velocity strengthening regions can propagate ruptures76.
As it can be anticipated from Eq. (3) and Eq. (4), there is a very wide spectrum of sliding mechanisms that will be generated at the fault over time. This introduces high uncertainties about the distribution of the parameters that govern rate – and state – friction relations. Nevertheless, from a dimensional analysis of the frictional equations, reference73 defined two dimensionless parameters that can be understood as two degrees of freedom for the frictional behavior of a given fault:
$${R}_{u}= \frac{(b-a){\sigma }_{eff}}{G}\frac{W}{L}$$
5
$${R}_{b}=\frac{b-a}{b}$$
6
Where \({R}_{u}\) is known as the Dieterich-Ruina-Rice number and depends on the geometry of the fault given by the width \(W\) of the velocity weakening patch, on the frictional behavior given by \(a-b\), on the effective normal stress and on the shear modulus \(G\) at plane strain conditions. For a velocity weakening segment \({R}_{u}\) takes positive values and represents the relation between the segment dimension and the nucleation size. Grater values of \({R}_{u}\) will leads to the nucleation of smaller instabilities with respect to \(W\) 73.
The other parameter is \({R}_{b}\), which controls the transient evolutionary effects and the strengthening behavior. Velocity strengthening patches will have \({R}_{b}\) < 0, velocity weakening domains are characterized by 0 < \({R}_{b}\) < 1 and velocity neutral regimes are defined when \({R}_{b}\) ≈ 0. \({R}_{b}\) also controls the rupture propagation in the same way as the ratio \(a/b\) 77. At velocity weakening domains and high sliding velocities, when \(a/b\) is close to the unity, then the rupture will tend to extend towards the whole segment.
The combination of \({R}_{u}\) and \({R}_{b}\) gives a two-dimensional space where different rupture styles can arise. By constraining the spectrum of frictional parameters, the actual slip behavior of a given fault can be approximated. We used this approach to feed the initial state of all the parameters in our model. We collected a wide range of datasets to constrain the pressure-temperature conditions and the slip mechanism at the plate interface. Then, we assumed an initial composition for the oceanic plate and for the sediments at the trench. From that point, we used previous stability fields of subduction related metamorphic facies to define the metamorphic state along dip. Finally, we defined the frictional parameters using previous gouge data from different lithologies that represent the metamorphic facies at the plate interface.
The viscoelastic behavior was calculated using a power-creep law that assumes pure dislocation creep over minerals with a constant water concentration31:
$$\epsilon =A{\sigma }^{n}{{C}_{OH}}^{r}{d}^{-m}exp\left(\frac{-Q+p\varOmega }{RT}\right)$$
7
Where 𝜀 is the strain rate, \(A\) is a pre-exponential factor, 𝜎 is the deviatoric stress, \({C}_{OH}\) is the water concentration, d is grain size, Q is the activation energy, p is the lithostatic pressure, 𝛺 is the activation volume, R is the constant for ideal gasses, T is the temperature and n, r and m are experimentally derived exponents. To yield an expression independent of the grain size, we assumed m = 0 and n > 131,78.
Assuming this configuration for the subduction fault, we made a meshed subduction system (Extended Data Fig. 1) to run a numerical model of the subduction seismic cycle. We discretized the system into a triangular mesh of 857 cells with a size of 20 km (Extended Data Fig. 1), while the points elements at the fault were separated by 200 m along the Slab2 geometry18 allowing numerical convergence. Deformation at each triangle was computed considering a Burgers assembly where the total anelastic strain rate is given by 𝜀T=𝜀M+𝜀K, in which 𝜀M and 𝜀K are the instantaneous strain rates in the Kelvin and Maxwell elements respectively32. Power-creep law parameters for the oceanic and continental mantle (Extended Data Table 2) were assumed following previous experimental31,32. We used the numerical code Unified Cycle of Earthquakes (UniCyclE) that employs the integral method5,79,80, including surface and volume elements. The changes in fault traction will depend on the surrounding fault slip and distributed strain. We assumed a background strain rate in the horizontal direction of ε = -1x10−14 s− 1 and predefined a set of synthetic GNSS stations at the surface of our mesh to compare those results with the actual data.
Temperature model
The subduction thermal model was constructed considering an oceanic slab subducting with a kinematically prescribed velocity (66 mm/year)17 beneath a fixed and layered continental plate. We used the continental lithospheric structure of reference26 considering an upper and lower crust, an oceanic crust, and a continental and oceanic mantle. In this framework the downgoing slab drives flow within the overlying viscous mantle wedge.
Velocity and pressure at the mantle wedge are found by solving the equations of mass and momentum conservation. Following reference81, we solve the equation of conservation of mass:
And the conservation of momentum:
$$\nabla \cdot \tau -\nabla P=0$$
9
Where v is velocity, 𝜏 is the deviatoric stress tensor, and P the dynamic pressure. The deviatoric stress tensor is given by:
$$\tau =2\eta \epsilon$$
10
In which 𝜂 is the effective dynamic viscosity and the components of the strain rate tensor ε are:
$${\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial {v}_{i}}{\partial {x}_{j}}+\frac{\partial {v}_{j}}{\partial {x}_{i}}\right)$$
11
We assume that viscosity is obtained by dislocation creep:
$$\eta =Aexp\left(\frac{Q}{nRT}\right){\epsilon }^{\frac{(1-n)}{n}}$$
12
Where rheological parameters are the same as Eq. (7). We considered a temperature and stress dependent wet olivine rheology for the mantle wedge31.
Boundary conditions for mantle wedge are (1) no slip below the rigid overriding plate and (2) constant velocity with a value equal to plate convergence along the top of the slab and below the maximum decoupling depth68.
The thermal field within the entire subduction zone was computed by solving the steady-state heat advection-diffusion equation, so temperature is obtained by solving:
$$\rho {c}_{p}(v\cdot \nabla )T=\nabla \cdot (k\nabla T)+f$$
13
In which 𝜌 is density, cp is the specific heat, k is thermal conductivity, and f is a parameter that describe possible heat sources.
To solve the system, we assumed a constant temperature at the surface of the model (0°C) and at the bottom of the oceanic lithosphere (1450°C). On the slab inflow boundary, we prescribed a geotherm obtained using a half-space cooling model82 considering an age of 14 Ma83. At the right edge of the model, we assumed a conductive geotherm calculated using the thermal properties of Extended Data Table 3, which is then applied along the landward boundary of the overriding plate. Along the inflow part of the wedge, we prescribed a geotherm which was obtained using an adiabatic gradient of 0.4°C/km and a mantle potential temperature of 1300°C. For the outflow boundary of the wedge, zero heat flux is assumed.
Finally, we solve equations (8), (9) and (11) using a stabilized finite element84 which allows us to use equal order of polynomial interpolation at each unknown point.
Pressure and temperature path
To compute the PTp, we used the lithospheric structure of reference26 and we take the ponderate density structure of reference85. We compute the vertical stress with:
$${\sigma }_{V}=g\rho h$$
14
Where g is the gravity acceleration constant, ⍴ is density and h is the thickness of the lithospheric layer (upper crust, lower crust and lithospheric mantle). By combining the temperature distribution at the plate interface and the \({\sigma }_{V}\) over it, we obtained the pressure-temperature path from Extended Data Fig. 1.
Coseismic slip model
We inverted the coseismic slip for the dip direction using the uplift displacements on the surface from the work of reference12 in between 42°S and 45°S.
To find the solution, we applied the Bayesian inversion approach of reference66, which allows us to obtain positive model parameters and associated uncertainties. This constraint allows us to obtain a more realistic slip solution as the slip direction is positive with respect to the dip direction, given the subduction stress regime. We impose a slip solution to be a multivariate folded normal distribution, and simultaneously we seek to find a solution for \(Gm=d\), where \(d\) are the deformation data, \(m\) is the backslip on each subfault in the plate contact, and \(G\) is the matrix which contains the Green’s functions. The Bayesian formulation is given by:
$$p\left(s|d,H\right)=\frac{p\left(d\right|s\left)p\right(s\vee H)}{p(d\vee H)}$$
15
with
$$p(s\vee H)={\left(2\pi \right)}^{\frac{-{N}_{m}}{2}}{\left|S\right|}^{\frac{-1}{2}}exp\left[\frac{1}{2}{({s}^{p}-s)}^{T}{S}^{-1}({s}^{p}-s)\right]$$
16
$$p(d\vee s)={\left(2\pi \right)}^{\frac{-{N}_{d}}{2}}{\left|D\right|}^{\frac{-1}{2}}exp\left[\frac{-1}{2}{(d-G|s\left|\right)}^{T}{D}^{-1}(d-G|s\left|\right)\right]$$
17
To enforce positivity constraints on the model parameters, we impose the following changes of variables:
$$m\left(s\right)={\left(\right|{s}_{1}|,|{s}_{2}|...,|{s}_{m}\left|\right)}^{T}$$
18
which guarantees us a positive solution and results in a Folded Normal Distribution in the posterior distribution.
We used the maximum evidence criteria to perform Bayesian model selection, namely, to search for the hyperparameters of a particular hypothesis \(H\). \(D\) and \(S\) are the covariance matrices of the likelihood and prior, respectively. The mean of the prior is \({s}^{p}\). We then obtain the correlation between the slip parameters, based on the information of the same data.
References
69. Lapusta, N. & Rice, J. R. Nucleation and early seismic propagation of small and large events in a crustal earthquake model. J Geophys Res Solid Earth 108, (2003).
70. Barbot, S. Modulation of fault strength during the seismic cycle by grain-size evolution around contact junctions. Tectonophysics 765, 129–145 (2019).
71. Dieterich, J. H. & Kilgore, B. Implications of fault constitutive properties for earthquake prediction. Proceedings of the National Academy of Sciences 93, 3787–3794 (1996).
72. Dieterich, J. H. Modeling of rock friction 1. Experimental results and constitutive equations. in Journal of Geophysical Research: Solid Earth vol. 84 2161–2168 (Blackwell Publishing Ltd, 1979).
73. Barbot, S. Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. Tectonophysics 768, (2019).
74. Marone, C. LABORATORY-DERIVED FRICTION LAWS AND THEIR APPLICATION TO SEISMIC FAULTING. Annu Rev Earth Planet Sci 26, 643–696 (1998).
75. Bürgmann, R. The geophysics, geology and mechanics of slow fault slip. Earth Planet Sci Lett 495, 112–134 (2018).
76. Kozdon, J. E. & Dunham, E. M. Rupture to the Trench: Dynamic rupture simulations of the 11 march 2011 Tohoku earthquake. Bulletin of the Seismological Society of America 103, 1275–1289 (2013).
77. Rubin, A. M. & Ampuero, J. P. Earthquake nucleation on (aging) rate and state faults. J Geophys Res Solid Earth 110, 1–24 (2005).
78. Hirth, G. & Kohlstedt, D. Rheology of the Upper Mantle and the Mantle Wedge: A View from the Experimentalists. Geophysical Monograph Series 138, 83–105 (2004).
79. Barbot, S., Moore, J. D. P. & Lambert, V. Displacement and stress associated with distributed anelastic deformation in a half-space. Bulletin of the Seismological Society of America 107, 821–855 (2017).
80. Muto, J. et al. Coupled afterslip and transient mantle flow after the 2011 Tohoku earthquake. Sci Adv 5, (2019).
81. van Keken, P. E. et al. A community benchmark for subduction zone modeling. Physics of the Earth and Planetary Interiors 171, 187–197 (2008).
82. Turcotte, D. & Schubert, G. Geodynamics. (Cambridge University Press, 2014). doi:10.1017/CBO9780511843877.
83. Tebbens, S. F. & Cande, S. C. Southeast Pacific tectonic evolution from Early Oligocene to Present. J Geophys Res Solid Earth 102, 12061–12084 (1997).
84. Araya, R., Cárcamo, C., Poza, A. & Poza, A. H. A stabilized finite element method for the Stokes-Temperature coupled problem. (2022).
85. Tassara, A. Control of forearc density structure on megathrust shear strength along the Chilean subduction zone. Tectonophysics 495, 34–47 (2010).
86. den Hartog, S. A. M., Niemeijer, A. R. & Spiers, C. J. New constraints on megathrust slip stability under subduction zone P-T conditions. Earth Planet Sci Lett 353–354, 240–252 (2012).
87. Phillips, N. J., Belzer, B., French, M. E., Rowe, C. D. & Ujiie, K. Frictional Strengths of Subduction Thrust Rocks in the Region of Shallow Slow Earthquakes. J Geophys Res Solid Earth 125, (2020).
88. Kirkpatrick, J. D., Fagereng, Å. & Shelly, D. R. Geological constraints on the mechanisms of slow earthquakes. Nature Reviews Earth and Environment vol. 2 285–301 Preprint at https://doi.org/10.1038/s43017-021-00148-w (2021).