Rheological segmentation of the Cocos slab and its relation with the 2017 Mw8.2 Tehuantepec earthquake

Tectonic plates bend and deform when approaching a subduction zone, creating intense faulting and highly variable stress and strain fields across short distances inside the slab. In September 2017, a large M_w8.2 intraplate normal fault occurred in southern Mexico, with an epicentral area located within a seismic gap where no megathrust had struck in more than a century. Despite the relatively young and hot Cocos plate, this seismic episode ruptured almost the entire slab below the brittle–ductile transition zone that normally limits the depth extent of such events. Here, we present a high-resolution thermomechanical model of spontaneous subduction for this area, where bending-induced brittle and ductile deformation and grain plate damage are considered. Modeling results show that the 2017 M_w8.2 Tehuantepec normal fault earthquake occurred due to the reactivation of one of the outer-rise-formed abyssal faults. In addition, the hypocenter was located in a stable, hydrated region of the lithospheric mantle at the transition limit between the elastic and ductile regimes. We found that earthquake rupture orientation is consistent with a region where a clear localized shear band of reduced effective viscosity is predicted. We propose that the rupture of this large intraslab event propagated in the ductile portion of the slab initially by a transformational faulting process, followed by a thermal runaway mechanism at greater depths and higher temperature.

Introduction

In the forearc segment of subducting plates, the occurrence of large intraslab seismic events is relatively a rare phenomenon, especially in young oceanic plates. Bending of oceanic plates before entering into subduction generates stresses in the oceanic plate that can produce the so-called trench-outer-rise events, which are characterized by horizontal T-axes in the shallow part and P- axes in the deeper region (Seno and Yamanaka, 1996). This state of stress is not usually preserved after subduction when the slabs unbend at intermediate depths (Kawakatsu, 1986), and the occurrence of large intraslab events in the first tents of kilometers after subduction is expected to be uncommon. Although such large seismic events are sporadic, they have a real potential to cause severe local damage (i.e., the El Salvador earthquake, M_w 7.6, and the Geiyo earthquake, M_w 6.8, in Southwest Japan). Compared with megathrusts, intraslab earthquakes occurred in the past have shown a rather limited magnitude. Several relevant examples include the 1931 Oaxaca earthquake (M_w 7.7) in Mexico (Singh et al., 1985) and the 1970 Peru earthquake (M_w 7.9) (Abe, 1972). However, in 2017, a large Mw 8.2 intraslab normal fault earthquake unexpectedly struck southern Mexico and generated even a tsunami with peak waves reaching up to 1.75 m in amplitude (Melgar et al., 2018). Its seismological features, including the rupture process, aftershock sequence, and tsunami generation, have been extensively studied, providing valuable insights into the tectonic processes of the region. Meng et al. (2019) reported that the mainshock is characterized by a fast rupture velocity of 3.6 km/s, approaching 85% of the local shear wave velocity. Using relocated aftershocks, Suárez et al. (2019) determined a 160-km-long fault, sub-parallel to the oceanic trench, located beneath the interplate contact. In addition, the aftershock sequence revealed a secondary intraslab fault located ∼50 km down-dip from the mainshock rupture. This suggests that the Cocos slab in the Tehuantepec region is characterized by tensional stresses, possibly induced by strong slab pull forces (Suárez et al., 2019). Song and Ge (2019) revealed that the rupture of the 2017 Tehuantepec earthquake can be explained by a two-fault model, based on the intersection with the Tehuantepec fracture zone, which transversely intersects the Middle America Trench (MAT) (Figure 1) (Manea et al., 2005). In this model, the subducted Tehuantepec fracture zone acts as a temporarily barrier, where the first coseismic rupture is temporarily blocked and decelerates but crosses over and continues due to the accumulated stress. Despite the large magnitude, the intraslab earthquake generated a tsunami with a height of ∼3.5 m, as recorded in Puerto Chiapas, Mexico (Gusman et al., 2018; Adriano et al., 2018). Moreover, a detailed seafloor survey carried out in the epicentral area relatively shortly (in 2019) after the 2017 Tehuantepec earthquake reveals that no submarine landslides have occurred on the gulf of Tehuantepec continental slope, probably due to limited seafloor uplift during the intraslab event (Aguilar-Anaya et al., 2025). Apart from its unusually large magnitude (it is one of the largest instrumentally recorded earthquakes in Mexico) for this type of earthquake, there are a series of characteristics that make it more peculiar. It occurred in a region known as the Tehuantepec seismic gap (Figure 1), where no large megathrust earthquake had occurred in the last century or more (Suárez, 2021); its rupture area cut almost the entire Cocos slab, reaching regions with predicted temperature exceeding 1,000°C (Manea and Manea, 2008), and finally, the rupture area stopped in the region where the Tehuantepec fracture zone enters into subduction (Okuwaki and Yagi, 2017; Ye et al., 2017; Melgar et al., 2018). Moreover, its focal mechanism indicates an intraslab normal fault type, steeply dipping at ∼75° (Ye et al., 2017), and the epicenter is located only 55 km from the MAT (Figure 1A), which represents an uncommon location for large intraplate normal-faulting earthquakes (Astiz et al., 1988). Meng et al. (2019) reported that the mainshock is characterized by a fast rupture velocity of 3.6 km/s, approaching 85% of the local shear wave velocity. Using relocated aftershocks, Suárez et al. (2019) determined a 160-km-long fault, sub-parallel to the oceanic trench, located beneath the interplate contact. In addition, the aftershock sequence reveals a secondary intraslab fault located ∼50 km down-dip from the mainshock rupture. This suggests that the Cocos slab in the Tehuantepec region is characterized by tensional stresses, possibly induced by strong slab pull forces (Suárez et al., 2019). Therefore, these rather unique characteristics of the 2017 Mw8.2 Tehuantepec earthquake offer us a good opportunity to investigate why these types of events occur and what is the possible mechanism. The hypocenter of the 2017 M_w8.2 Tehuantepec earthquake is located at 50 km in depth (Figure 1B), where the lithostatic pressure approaches approximately 2 GPa and the temperature can reach 700°C. At such large temperature and pressure rates, laboratory experiments show that rocks deform by ductile creep rather than brittle deformation (Kirby, 1980; Kirby and Kronenberg, 1987). Nevertheless, the 2017 M_w8.2 Tehuantepec earthquake shows that this is rather possible, and therefore, some special mechanisms should be responsible for rupture propagation in regions where temperature can reach 1,200°C or even more.

Figure 1

Figure 1. (A) Topography/bathymetry map of the Mexican subduction zone showing the location of the 2017 M_w8.2 Tehuantepec intraslab normal fault earthquake. Red dots represent active volcanoes, and magenta diamonds show the location of intraslab earthquakes (Singh et al., 2020). Slab surface isodepths are derived from Slab 2.0 model (Hayes et al. 2018). EPR, East Pacific Rise; MAT, Middle America Trench; TFz, Tehuantepec fracture zone; NAM, North America Plate; CP, Caribbean Plate; RP, Rivera Plate. The semitransparent arrow shows the Cocos plate convergence rate (67 mm/yr) in the region of the 2017 M_w8.2 Tehuantepec earthquake. The red colored round region depicts the location of the Tehuantepec seismic gap bounded by two past rupture areas (yellow colored round regions), the 1965 M7.4 and the 1902 M7.7 megathrusts. (B) Slab geometry along the cross-section A-A′ shown in (A). The slab geometry, location of the 2017 M_w8.2 Tehuantepec earthquake hypocenter, and the strike-averaged amount of slip of the mainshock are derived from the study by Zhang and Brudzinski (2019). Maps were created based on ETOPO1 Global Relief Model dataset from the study by Amante and Eakins (2009) and generated with the open-source software ParaView (http://www.paraview.org) version 5.0.1, licensed under the CC BY 4.0 license (https://creativecommons.org/licenses/by/4.0/).

Numerical modeling became an essential instrument for understanding the mechanical behavior of subduction systems. Advances in computational power enabled higher-resolution models that capture small-scale processes in subduction zones. Currently, ongoing research focuses on improving computational efficiency, material property representation, and data integration. These advancements are paving the way for more accurate and predictive models of subduction zone dynamics (Gerya, 2022). One common approach in numerical modeling is to integrate thermal and mechanical processes to simulate the heat transfer and rheological changes in subducting slabs and also the surrounding mantle. Coupled thermomechanical models are used to study the thermal evolution of subduction zones and its impact on slab dehydration and mantle wedge melting (van Keken et al., 2008; Gerya and Meilick, 2011; Crameri et al., 2020; van Keken and Wilson, 2023). The mechanical behavior of subducting plates is strongly linked with the occurrence of intraslab earthquakes. As they travel long distances with almost no internal deformation, oceanic plates are considered mechanically strong before subduction. After subduction, slabs can preserve a cold core to greater depths, and numerical modeling of subduction predict rather high viscosity (i.e., 10²⁴ Pa s) specific for strong slabs (Zhong and Davies, 1999; Faccenna et al., 2007; Billen, 2008). As strong plates are difficult to segment along fault planes under the elevated lithostatic pressure during subduction, it is challenging to understand what the actual mechanism for intraslab earthquakes is.

Observations indicate that prior to subduction, plates undergo upward bending, creating tensional stress in the upper portion of the oceanic plate and compressional stress at deeper levels, below the neutral plane (Watts, 2001). Oceanic plate morphology off trenches (i.e., distances less than 80 km) (Nakanishi, 2017) indicates that plate bending is not mechanically an elastic process, but rather it reactivates in the form of normal faulting along pre-existing seafloor spreading fabric and bathymetric structures (Ranero et al., 2003). Therefore, slabs are rendered mechanically weak just before subduction, and weakening of the oceanic lithosphere is key for understanding how subduction is initiated (Bercovici and Richard, 2012).

As tectonic plates advance from oceanic ridges toward plate boundaries, where they eventually enter subduction, the grain size of minerals (as olivine and pyroxene) specific to the lithospheric mantle also evolves. Observations show the existence of localized shear bands coupled with grain-size reduction in mantle peridotites (Furusho and Kanagawa, 1999). Inspired by metallurgical research, Bercovici and Richard (2012) proposed a novel mechanism of lithospheric weakening by grain evolution and Zener pinning in two-phase lithospheric rocks. Lithospheric weakening by grain evolution refers to the processes by which changes in the size, shape, and arrangement of mineral grains within the oceanic lithosphere reduce its overall strength and rigidity. This weakening can significantly influence the mechanical behavior of the lithosphere, including its ability to deform and transmit stress. Grain evolution mechanically weakens the lithosphere through several processes, such as the grain-size reduction, where smaller grains increase the surface area of grain boundaries, which are weaker than the grains themselves. This process enhances considerably the deformation mechanisms like grain boundary sliding and diffusion creep, which are more efficient at smaller grain sizes, and as a result, the lithosphere becomes more ductile and less rigid. Zener pinning is specific for polycrystalline mediums and shows that grain size is significantly reduced in damaged regions, switching the rheology into the grain-size-dependent diffusion creep regime. In a two-phase medium (pyroxene + olivine) subjected to a strong deformation process (Supplementary Figure S1), the percolation of the olivine phase along the grain boundaries of the pyroxene phase has several effects: producing a rougher interface, obstruction of grain growth in a region near the interface between the two phases (known as Zener pinning), and simultaneously, splitting of pyroxene grains into new smaller grains. The severe splitting of pyroxene grains into new smaller grains leads to a reduction in the average grain size and produces a grain damage zone where shear is localized in narrow bands.

The physics behind rock deformation at high pressure and temperature is still not well understood, and the Zener pinning process in the lithosphere appears to play a key role in the formation of tectonic plates (Mulyukova and Bercovici, 2019). Gerya et al. (2021) numerically investigated the consequences of off trench bending induced for oceanic plates and found that grain-size reduction by Zener pinning in the outer rise focuses plate damage and mechanically segments the slab into relatively rigid blocks. This implementation of Zener pinning processes into numerical geodynamic models is a novel approach that reconciles the apparent inconsistency between strong tectonic plates at the surface and weak slabs, and is consistent with seismic tomography images beneath Japan, where the Pacific slab appears segmented into blocks with coherent v_p velocity anomalies (Tao et al., 2018). Motivated by the occurrence of the large 2017 M_w8.2 Tehuantepec intraslab earthquake, and based on the new numerical model of Gerya et al. (2021), here we present results from 2D high-resolution thermomechanical modeling of spontaneous subduction for southern Mexico, where bending-induced brittle–ductile deformation and grain damage [the race between normal grain growth and dynamic recrystallization changes grain size at time scales comparable to earthquake cycles (Mulyukova and Bercovici, 2019)] are considered. We aim to investigate the thermomechanical conditions inside the Cocos slab that led to the occurrence of the 2017 M_w8.2 Tehuantepec intraslab earthquake.

Materials and methods

Previous generations of 2D thermal models applied to the southern Mexico subduction zone (Manea and Manea, 2008) used a system of Stokes equations coupled with the steady-state heat transfer equation. This system is solved using the finite element solver PDE2D (http://pde2d.com/). One main limitation of these previous models is the shape of the slab that is considered fixed. Although this approach allows for the incorporation of complex geometries, material heterogeneities and nonlinear rheologies (e.g., viscoelasticity and plasticity) are not incorporated. Recent works emphasize the importance of integrating realistic rheologies, including temperature- and pressure-dependent behaviors (Billen, 2008; Čížková and Bina, 2013; Behr et al., 2022; Hummel et al., 2024) and, more recently, grain-size evolution due to brittle–plastic and ductile deformation (Gerya et al., 2021) (Supplementary Table S1). This novel approach explains several key features related with tectonic plates, such as the apparent inconsistency between strong plates and weak slabs, the formation of large-offset normal faults in the vicinity of trenches, and the segmented nature of seismic velocities within slabs.

Plasticity is implemented using a yield criterion which limits the creep viscosity, which is represented as a function of temperature and stress in terms of deformation invariants through experimentally determined flow laws. Constant grain size is assumed for the oceanic and continental crusts, and a ductile creep is used for the mantle considering grain-size reduction and growth processes assisted by Zener pinning for a bi-material mixture (olivine 60% and pyroxene 40%). The grain size is controlled by the roughness of the interface between the olivine and pyroxene phases and is computed based on several factors, including grain-growth rate and the fraction of mechanical work converted to interface damage (see Supplementary Material). We opt for free bending to initiate subduction, and the initial model setup (Supplementary Figure S2) incorporates a mechanically weak zone (see parameters in Supplementary Table S2) between the oceanic and continental domain. We also include a small vertical strip (5 km $\times$ 50 km) weak zone at the left-hand boundary to ensure oceanic lithosphere decoupling from the boundary. The model domain is 3,000 km wide, 1,000 km deep, and has a high-resolution grid (1 km × 1 km) centered in the region (1,000 km × 200 km), where subduction takes place. This high resolution is sufficient to allow the simulation to resolve shear localization. Outside the region of interest, the resolution is gradually decreased to 10 km $\times$ 10 km along the model domain boundaries. The grid is composed of 1,261 × 351 nodes, and we use 13 million randomly distributed markers used for rock composition (Supplementary Table S3). Free-slip mechanical boundary conditions are considered for all four sides. In terms of initial temperature distribution, for the oceanic plate, we use a vertical thermal structure as a function of plate age. The oceanic plate age varies from 0.001 Ma to the left-hand limit to a value between 25 and 35 Ma at the right-hand limit (Supplementary Table S2; Supplementary Figure S2).

Temperature at the surface is 0°C and mantle potential temperature (T_m) at the base of lithosphere, and for the mantle, we use an adiabatic gradient of 0.5°C/km. At the lateral sides, we consider zero heat flux as the thermal boundary condition. The continental plate consists of uniform sediments, and the upper and lower crusts with a total thickness of 30 km, followed by a 60-km-thick lithosphere (Supplementary Figure S2). To avoid surface oscillations and to ensure heat transfer from the surface of the crust, on top of the oceanic plate, we insert a weak zone and continental plate called “sticky” air/water with a constant viscosity of 10¹⁷ Pa s and a constant temperature of 0°C (Supplementary Figure S2). We performed a series of models where we varied the initial oceanic plate at the trench (25, 30, and 35 Ma), the overriding plate age (60, 70, and 80 Ma), mantle potential temperature T_m (1,623, 1,648, and 1,673 K), and the initial width of the weak zone between the oceanic and continental domains (80, 90, and 100 km).

Modeling results

For a young plate such as the Cocos plate in southern Mexico (Sdrolias and Muller, 2006), grain-size reduction and ductile damage in the outer rise region facilitate spontaneous subduction initiation. The numerical model incorporates the damaged fault zones, owing to the reactivation at the outer rise. We varied the initial oceanic plate age at the trench from 25 to 35 Ma, with 5 Myr increments, and other parameters were set as follows: overriding plate age 70 Ma, potential temperature T_m = 1623 K, and the initial width of the weak zone between the oceanic and continental domains set at 100 km. Modeling results show that for oceanic plate ages of 30 Ma (Figure 2; Supplementary Figures S3–S6) or 35 Ma (Supplementary Figure S7), predicted plate geometries are similar with the present-day slab shapes (Zhang and Brudzinski, 2019). However, for the model with a plate age of 25 Ma, the slab bends easier into the upper mantle, and slab geometry misfits the slab shape. At the same time, we observe that irrespective of oceanic plate age, model simulations predict slab segmentation through localized shear bands. For the model with an oceanic plate age of 35 Ma, we observe the presence of two V-shaped concentrated shear band regions in the area of the 2017 Tehuantepec intraslab earthquake. They merge into a single region but at higher depths (Supplementary Figures S8–S10). On the other hand, the influence of the overriding plate thermal structure is not significant. A thermal gradient specific for plate ages of 60 Ma and 80 Ma produces slab segmentation and geometries that fit the observations (Supplementary Figures S11–S14). Increasing the mantle potential temperature with 25 K (T_m = 1648 K) and 50 K (T_m = 1673 K), respectively (Supplementary Figures S15–S18), leads to increased slab dips, with the clear formation of slab segments delimited by focused damaged zones (Supplementary Figure S16), as in the previous simulations.

Figure 2

Figure 2. Modeling results: (A) effective viscosity, (B) temperature, (C) rock composition. Location of the 2017 M_w8.2 Tehuantepec earthquake hypocenter and the strike-averaged amount of slip of the mainshock are derived from the study by Zhang and Brudzinski (2019), and (D) grain size. White dashed curve represents the present-day Cocos slab geometry along profile A-A′ shown in Figure 1A. The black thin curve at the base of oceanic lithosphere in (A, D) represents the 1,300°C isotherm. The horizontal axis corresponds to the localization of the real MAT relative to the position of the 2017 M_w8.2 Tehuantepec earthquake.

The width of the weak zone plays a crucial role in subduction initiation. Models with an initial width of 80 km or less fail to subduct (Supplementary Figures S19–S22), whereas increasing the width to 90 km or more facilitates subduction. In the following sections, we discuss modeling results for one simulation (T_m = 1625 K, oceanic plate age = 30 Ma, overriding plate age = 70 Ma, and the width of the weak zone = 100 km) that best fit the present-day geometry and other observations. Owing to its weak mechanical properties, the initial stage of the modeling is marked by the weak zone dripping into the mantle, and slab starts to bend downward gently only after approximately 0.25 Myr of evolution (Supplementary Figure S3). As the subduction proceeds, the overriding portion initially located above the weak zone is exposed to extension that even creates a small short-lived transitory region of partial melt (Supplementary Movies SM1, SM2). Once the highly viscous Cocos plate starts to bend, we observe the formation of a wide region of bending-induced damage in front of the trench at the outer rise (Supplementary Figure S4). The viscosity in the damaged area characterized by a sequence of parallel shear bands decreases several orders of magnitude compared with the undamaged plate. Therefore, we can consider these low-viscosity shear bands the equivalent of plate faulting. In the outer rise, the faulting network has a lenticular shape controlled by the amount of slab bending in the region. It penetrates the first half (∼25 km) of the Cocos plate in the brittle region, and the lower part of the damaged region is limited to the ∼700°C isotherm.

As depicted in previous studies of outer rise faulting (Faccenda et al., 2008), our model also shows the formation of a dense arrangement of shear bands disposed symmetrically at ∼30° from the vertical (Figure 3A). Below the ∼700°C isotherm in the outer rise region, no brittle faulting is observed, indicating that the slab entered into ductile regime at that depth. The thin grain reduction band system in the outer rise starts forming as early as 0.25 Myr after model initiation and continues to grow in length and thickness as the slab bending further progresses (Supplementary Figure S5). In terms of Cocos plate damage before subduction in southern Mexico, the oceanic plate is clearly affected by abyssal faults disposed in a rather regular network subparallel with the MAT orientation (Figure 3).

Figure 3

Figure 3. (A) Effective viscosity distribution and the location of the 2017 M_w8.2 Tehuantepec earthquake hypocenter. (B) 3D view of the Cocos plate bathymetry [Global Multi-Resolution Topography (GMRT) Data Synthesis, distributed by OpenTopography: https://doi.org/10.5069/G9BG2M6R] and the forearc in the vicinity of MAT showing the distribution of reactivated seafloor faulting fabric. MAT, Middle America Trench; TFz, Tehuantepec fracture zone. (C) The observed bathymetry across the B-B′ cross-section shown in (B). Based on ocean floor bathymetry, dashed black lines mark our interpretation of depth distribution of abyssal faulting in the vicinity of MAT.

Analyzing the faulting patterns using effective viscosity along the downgoing plate (Figure 2A), we observed that brittle deformation is concentrated in narrow low-viscosity shear bands, or faults, which cut the oceanic plate up to the surface. In the outer rise region at the surface, we note values of faulting spacing that vary from ∼5 km off trench to ∼8 km near the trench (Figure 3A). The outer rise population of several abyssal faults is clearly visible in the bathymetry map (Figure 3B-inset), and this is partially due to a quite thin (<200 m) sedimentary layer of the Cocos plate in that area (Manea et al., 2003). As the subduction model evolves and the slab sinks further into the mantle, the brittle–plastic deformation and grain-size reduction wider regions shrink at depths until a localized narrow region (∼25 km in length) composed of only few shear bands is observed (Figure 2). We also observe that the cold core of the slab reaches 150 km depth or more (Figure 2B), and it is consistent with the intraslab seismicity in this area (Rebollar et al., 1999). Following the plate curvature deeper into subduction, faulting fabric gradually ceases in the upper part of the downgoing slab, owing to the increase in yielding stress due to the continuous increase in pressure (Figure 2A). However, due to the feedback between grain-size reduction and brittle and ductile deformation, slab damage merges at depth into a quite narrow area (Figure 2A; Supplementary Figure S4). This process represents the slab segmentation process that creates integral high-viscosity slab blocks loosely linked in a chain-like manner by a quite narrow faulted zone (Gerya et al., 2021). In our simulation, the location of the 2017 M_w8.2 Tehuantepec earthquake hypocenter corresponds to the 700°C isotherm and to one of the localized plate damage zones. The high temperature observed above the slab (Figure 2B) represents a remnant of the initial hot temperature in the weak zone, necessary to allow subduction initiation (Supplementary Figure S6; Supplementary Movie SM2).

Grain-size reduction inside the Cocos slab is unevenly distributed due to normal fault weakening. Regions with increased plate damage by fault weakening show concentration of the grain-size reduction area (Figure 2D). Our model prediction suggests that the 2017 M_w8.2 Tehuantepec earthquake might have occurred in a region inside the slab marked by the transition from a normal grain slab to grain-size reduction. At the same time, the dip of the mainshock slip (Figure 2C) is in a fairly good agreement with the orientation of one of the main shear bands in the localized damaged region (Figure 2A).

Discussion and conclusions

In this work, we present an updated thermomechanical model of subduction in southern Mexico in the region where a large M_w8.2 Tehuantepec intraslab earthquake occurred in 2017. The main goal is to better understand the thermomechanical conditions inside the Cocos slab that control the occurrence of such unusual event where the rupture area propagated extremely deep and cut through almost the entire oceanic lithosphere. Currently, based exclusively on seismological evidence, it is difficult to identify the actual mechanism that could be the main cause of the large 2017 M_w8.2 Tehuantepec intraslab earthquake. However, using our numerical predictions, we can explore different hypotheses and propose a model that explains how it was possible to accommodate such an extensive rupture area through the Cocos slab (Melgar et al., 2018). Compared with the previous version of the thermal model (Manea and Manea, 2008), this updated simulation includes spontaneous bending, a realistic visco-plastic rheology, grain-size reduction, and growth processes assisted by Zener pinning for the ocean lithospheric mantle (see Methods). Figure 3A shows the viscosity distribution of the subducting Cocos plate, where ∼100 km-wide slab segments of high viscosity (∼10²⁵ Pa s) are formed, separated by narrow low-viscosity zones (∼10²⁰ Pa s) concentrated in shear bands. These slab segments are strong (i.e., have high viscosity), whereas the joints between these segments are mechanically weak. This specific viscosity pattern forms as a consequence of localized brittle–ductile deformation that is linked with grain-size reduction initially formed at the outer rise (Figure 2D). A dynamically self-consistent process is responsible for such Cocos slab segmentation, where ductile damage (through grain-size reduction) at the outer rise creates the conditions, which, at depth, when combined with fault weakening result in a localized slab segmentation. The feedback between faulting and grain-size reduction has to be strong in order to obtain such sharp slab segmentation. Supplementary Movie SM3 shows how an initial reduction of grain size is formed at the outer rise, and then, once formed, it propagates with the slab. Then, at depth, due to strong slab bending coupled with the grain-size reduction zone, a ductile strain-localization deformation (where slab viscosity is reduced along high-strain-rate shear bands) pattern occurs. The relationship between trench faulting at the outer rise and intraslab seismicity is well known along MAT, where there is a good correlation between bend faulting with sectors of the slabs with higher intermediate-depth seismicity (Ranero et al., 2005).

We explore the evolution of Cocos plate subduction in terms of slab shape, temperature, and viscosity and compare with observations, including the location and the mainshock slip distribution of the 2017 Mw8.2 Tehuantepec intraslab earthquake off trench bathymetry in southern Mexico, showing the existence of an extensive system of faulting (Figures 3A, B), which is the plate reaction to slab bending as the Cocos plate organizes the state of stress before entering subduction. In southern Mexico, south of the Tehuantepec fracture zone (Figure 1A), the fault pattern makes an angle with MAT of ∼30°, and this suggest that they are not new faults but rather a reactivated abyssal hill fabric formed at East Pacific Rise (Masson, 1991). When the system of abyssal faults near trench maintains a systematic angle (e.g., 30°), this strongly suggests that they are relict from mid-ocean ridge spreading and reactivated near trench. This contrasts with newly formed faults, which would align with the local stress field characterized by thrust faults aligned at a lower angle (e.g., 10°–30°) with the trench. Our numerical simulations satisfactorily capture the dynamics of outer rise deformation. We predict a faulting pattern where the slab is damaged upon entering subduction, with a system of localized shear zones spaced between 5 and 8 km, similar with the observations (Figure 3B, C).

The maximum depth extent of this fault system is approximately 20 km, and it is temperature controlled by the 700°C isotherm (Figure 4). This fault length is comparable with the maximum observed depth of abyssal (normal) faults along MAT, where they typically extend down to 15–20 km within the oceanic lithosphere (Ranero et al., 2003). In addition, Craig et al. (2014) revealed that several focal mechanisms recorded in the outer rise in southern Mexico along MAT are located based on their focal mechanism. The normal-fault events are concentrated in the upper part of the Cocos plate at depths <10 km, whereas the thrust events are located deeper (i.e., 26 km depth). This observation suggests that the transition from extension to compression in the oceanic plate takes place at ∼15 km depth below the seafloor or even more. As the oceanic plate proceeds into subduction due to the strong interaction between brittle and ductile damage localization, the ∼150-km-wide outer-rise area affected by the reactivated hill fabric merges into a fairly compact “V”-shaped deformation zone (Figure 4A). The depth orientation of one of our main deep shear bands in this is in good agreement with the orientation of the main shock slip of the 2017 M_w8.2 Tehuantepec intraslab earthquake. This suggests a close relationship between the two, at least for the upper half of the Cocos slab, where most of the slip is localized. In terms of temperature, we observe that the 2017 Tehuantepec intraslab earthquake is localized in the near proximity of the 700°C isotherm, which marks the transition from brittle to ductile regimes. Although the predicted slab temperature in the hypocenter area is high, the estimated amount of fluids potentially contained in the oceanic lithosphere show that the M_w8.2 Tehuantepec intraslab earthquake occurred in the stability field of chlorite (2.4 H₂O wt% - Figures 5A, B). As the intraslab earthquake occurred in the lithosphere, our models do not predict breakdown of chlorite lherzolite in the hypocenter region, and it is possible that a mechanical origin (that reflects strong internal deformation to accommodate downward bending of the Cocos slab) prevails. What is quite unusual is that slip also propagated downward below the hypocenter, however in a smaller amount (<1.5 m) (Melgar et al., 2018; Okuwaki and Yagi, 2017; Ye et al., 2017; Melgar et al., 2018; Meng et al., 2019). Here, the slab temperature exceeds even 1,200°C and extends in the region where the slab gradually returns to the larger grain size (Figure 5C), and weakening of the lithosphere is essential to permit slip propagation in this area. In subduction zones, under upper mantle conditions, smaller grain size reduces viscosity. Fine-grained rocks formed due to strong deformation (i.e., mylonites) tend to facilitate strain localization and weakening. It is known that the nature of deep intraslab earthquakes is still difficult to be explained due to high-temperature and high-pressure regime. In general, several individual hypotheses are proposed to explain what causes intraslab earthquakes within the subducting slabs under extreme P-T conditions: dehydration embrittlement (Green and Houston, 1995; Hacker et al., 2003; Frohlich, 2006), transformational instability or faulting (Green and Houston, 1995; Seno, 2009; Ferrand et al., 2017; Hasegawa and Nakajima, 2017), and thermal runaway (Hobbs et al., 1986; Ogawa, 1987; Kelemen and Hirth, 2007; Kita and Katsumata, 2015; Thielmann, 2018). In Figure 5C, we observe that the hypocenter of the 2017 M_w8.2 Tehuantepec intraslab earthquake is located in a region characterized by a reduced grain size. Actually, most of the reduced grain size area is located below the hypocenter, at depth below the brittle–ductile boundary. It is proposed that if the grain size is sufficiently reduced, then lithospheric rocks can easily deform and facilitate rock failure by transformational faulting. Grain-size reduction and rock failure by transformational faulting can be interconnected processes and play a meaningful role in the deformation and mechanical behavior of rocks, particularly in tectonic settings like subduction zones. Grain-size reduction indicates the decrease in the size of mineral grains within a rock; transformational faulting refers to the formation of faults or narrow shear zones, where deformation is accompanied by significant changes in mineralogy and microstructure (Hirth and Kohlstedt, 2003).

Figure 4

Figure 4. (A) Effective viscosity distribution. Dotted black curve represents the 700°C isotherm that delimits the brittle and ductile regimes inside the Cocos slab. The white thin continuous curve at the base of oceanic lithosphere represents the 1,300°C isotherm. (B) Temperature distribution. In both figures, we show the location of the 2017 M_w8.2 Tehuantepec earthquake hypocenter and the strike-averaged amount of slip of the mainshock (Zhang and Brudzinski, 2019). White dashed curve represents the present-day Cocos slab geometry (Zhang and Brudzinski, 2019) along the profile A-A′ shown in Figure 1A. Blue beach ball focal mechanisms represent historical seismicity from SSN (2017). Yellow diamonds represent aftershocks from the study by Zhang and Brudzinski (2019).

Figure 5

Figure 5. (A) Calculated temperature isocurves. The square inset depicts the hydrous phase distribution inside the Cocos slab, according to the phase diagram of hydrated lithospheric mantle (lower left inset) (Suenaga et al., 2021; Lee and Kim, 2021). The semitransparent area represents the oceanic crust. (B) A magnified representation of the square inset shown in (A). (C) Effective viscosity distribution for the inset region shown in (A). The semitransparent area is the predicted grain reduction area from Figure 1D. In all figures, we show the location of the 2017 M_w8.2 Tehuantepec earthquake hypocenter and the strike-averaged amount of slip of the mainshock (Zhang and Brudzinski, 2019). White dashed curve represents the present-day Cocos slab geometry (Zhang and Brudzinski, 2019) along profile A-A′ shown in Figure 1A.

In terms of mechanical strength, grain-size reduction weakens the rock, making it more susceptible to localized deformation and faulting. At the same time, fine-grained rocks are more prone to ductile deformation mechanisms, which can cause transition to faulting (Bercovici and Richard, 2014). On the other hand, transformational faulting is a driver of grain-size reduction, where intense deformation along faults can cause cataclasis and dynamic recrystallization, further reducing grain size. There is a feedback loop between the two processes: grain-size reduction weakens the rock, promoting faulting, and faulting generates further grain-size reduction through deformation. This feedback loop can lead to the development of highly localized shear zones or fault systems, including in regions where temperature is above the brittle–ductile limit.

Although we do not have incorporated transformational faulting mechanisms in our numerical simulations, we suggest that the hypothesis is worth future investigation for the case of the 2017 M_w8.2 Tehuantepec intraslab earthquake, where, due to the interplay between the bending-induced plate damage and grain reduction, our models predict an increased volume of the grain reduction zone just below the hypocenter (Figures 2, 5C). Grain-size reduction in this area is significant as it decreases from ∼3 mm inside the slab in the brittle region to ∼0.05 mm (Figure 2D). The occurrence of transformational faulting is hypothesized based on the presence, at the base of the brittle upper plate, of a zone of reduced grain size, as also predicted in the models of Gerya et al. (2021). The fine grain size of the parent rock is a prerequisite to allow garnet exsolution in orthopyroxene to have a significant rheological weakening effect. Garnet exsolution in orthopyroxene can separate from the host mineral orthopyroxene due to changes in temperature, pressure, or chemical composition. Fine-grained rocks, as predicted by our model (Figure 2D), tend to have a higher density of grain boundaries, which act as nucleation sites for garnet exsolution. Smaller grains provide more surface area for chemical reactions and phase transformations, facilitating the exsolution process. At the same time, garnet exsolution along grain boundaries can enhance grain boundary sliding, a deformation mechanism that weakens the rock (Faryad et al., 2009). In subducting slabs such as the Cocos slab, fine-grained rocks undergoing garnet exsolution have the potential to localize deformation, facilitating slab bending and segmentation. In our model, we observe four orders of magnitude reduction in slab viscosity (from 10²⁵ Pa s to 10²¹ Pa s) in the region below the brittle–ductile limit and focused toward the hypocentral area of the 2017 Tehuantepec earthquake (Figure 4A).

The abovementioned assumptions are based on the experimental results of Shi et al. (2022), who recorded acoustic emissions in peridotite samples deformed under conditions that led to the formation of ultrafine (sub-micrometric) garnet exsolution in orthopyroxene, which can further decrease the fine grain size of a shear zone. According to Shi et al. (2022), for the transformational faulting mechanism to be effective, the mantle peridotite should be in a temperature range, where localized shearing can occur, of 1,000–1,100 K (727°C–827°C). This is relatively consistent with our model prediction that shows a temperature range below the hypocenter of 700°C–900°C (Figure 5). As extrapolating the experimental results of Shi et al. (2022) to natural tectonic systems is not straightforward, another alternative explanation for the downward propagation of the 2017 Tehuantepec intraslab seismic rupture, independent of a specific occurrence of a metamorphic reaction, might be the sudden increase in strain rates (Ellis and Stöckhert, 2004).

Our models show that the rupture plane propagated further down under even higher P-T conditions and beyond the predicted grain reduction zone. As transformational faulting is limited by P-T conditions inside the slab, the rupture propagation at extreme temperature levels points toward a different rupture mechanism. We propose the conversion of the transformational faulting mechanism into a thermal runaway process. Yet, how this change actually takes place remains an interesting topic for future research. However, there are few lines of evidence that might help shed some light on this topic. Thermal runaway occurs when a small increase in temperature triggers processes that further increase temperature, creating a self-sustaining cycle. The common drivers for the subducted oceanic lithosphere include dehydration reactions and shear heating. Most hydrous minerals (e.g., serpentine and chlorite, Figure 5) would have already broken down at lower temperatures (T < 800°C). Therefore, the potential for thermal runaway driven by dehydration reactions is reduced. Frictional heating along the slab–mantle interface, or within the base of the slab, could still contribute to localized temperature increases, but the efficiency of shear heating depends on the strain rate and the rheology of the slab in that region. Our model predicts a temperature in the region of the thermal runaway process exceeding 1,000°C. We suggest that if thermal energy is generated faster enough than the lateral conductive transport, then mantle rocks are eventually locally destabilized, allowing the rupture area to propagate in high-temperature areas. At this high temperature and confining pressure, the presence of pore fluid is unlikely, and according to the experiments of John et al. (2009), the rocks could fail by the mechanism where ductile deformation in shear zones leads to heating, thermal softening, and weakening of rocks.

Information from seismological analysis of this large intraslab event helps place some constraints and limitations on our numerical results and interpretations. For example, Singh et al. (2014) reported a decrease in stress by ∼37–60 MPa for several large intraslab earthquakes beneath Mexico, and García et al. (2004) proposed that a median value of 30 MPa for stress decrease is a general characteristic for the intraslab Mexican earthquakes. However, Ye et al. (2017) calculated stress decrease of only 18 MPa for the M8.2 2017 Tehuantepec earthquake. This might suggest a reduced fault strength due to warmer or more ductile conditions inside the slab in the rupture zone. In addition, such relatively reduced amount of stress decrease can indicate that the earthquake, or part of it, may have occurred in a region of the slab where stress is more diffuse, rather than being concentrated on a single fault. Small stress decrease could indicate the presence of fluids from dehydration reactions within the subducting slab, which can reduce further effective stress along the rupture area. In terms of shear wave velocity for large intraslab earthquakes produced in different subduction environments, this ranges from 4.5 to 5.0 km/s, depending on the depth and composition of the subducting slab (Zhao et al., 1994; Kirby et al., 1996; Hacker et al., 2003; Abers, 2005). These velocities are higher than those in the surrounding mantle due to the slab’s cooler and more rigid nature. For the M8.2 2017 Tehuantepec earthquake, Ye et al. (2017) fitted the teleseismic data with a maximum rupture expansion velocity of only 3–4 km/s (the average source spectrum fits for a shear wave velocity of 3.75 km/s) and total rupture lengths of 160 km or longer. Such a low shear wave velocity is anomalously slow for a cold and rigid oceanic lithosphere, indicating some specific or additional physical conditions within the slab. In other subduction zones, such as Japan, Cascadia, and South America (Kawakatsu and Watada, 2007; Nakajima et al., 2009; Bostock et al., 2002; Audet et al., 2009; Sippl et al., 2022; Haberland et al., 2009), the results of seismological data analyses have identified regions within the slab with S-wave velocities lower than 4 km/s. These anomalies are often linked to hydration and serpentinization, partial melting, thermal anomalies, or areas of intense deformation. Our model shows that the P-T condition of the upper part of the rupture area (above the brittle–ductile limit, Figure 5A) allows for partial serpentinization of the slab. For the deeper part, however, dry slab conditions are expected, but partial melting of dry peridotite unlikely occurs due to the low P-T conditions predicted from our model along the rupture area.

In terms of a localized thermal anomaly, Calò (2021) identified several V_p/V_s seismic wave velocity anomalies within the Cocos slab in southern Mexico. However, these anomalies are located deeper (i.e., at 100 km depth) within the slab than the rupture area of the M_w8.2 2017 Tehuantepec earthquake. Our model predicts the existence of an overthickened area of grain-size reduction below the brittle–ductile region (Figure 5C). It is known that in subducting slabs, grain-size reduction represents a key process that influences the mechanical properties of the slab, particularly its rigidity (Karato and Wu, 1993). A reduction in rigidity due to grain-size reduction can make the slab more ductile and less capable of transmitting seismic waves efficiently. This might be an explanation for the anomalously slow shear wave velocity associated with the rupture area of the M_w8.2 2017 Tehuantepec intraslab earthquake.

Our numerical models adjust well several characteristics of the subduction of the Cocos plate in the segment where the Tehuantepec earthquake nucleated. Although the formation of slab segments delimited by damaged zones is observed for a wide range of modeling parameters (S3–S22), only few models fit the present-day slab shape (Figure 4). Additionally, our best-fit model satisfactorily predicts the plate flexure fault spacing in front of the Middle America Trench segment before subduction in southern Mexico (Figure 3). Although our models do not include a yielding criterion under high P-T conditions specific for a ductile regime inside the slab, based on the abovementioned lines of evidence, we hypothesize that the 2017 M_w8.2 Tehuantepec intraslab earthquake was probably facilitated by slab dehydration and ruptured above the hypocenter along a preexistent extensional faulting reactivated in the vicinity of the trench axis. Below the hypocenter, faulting is challenging to be explained. We propose a process by cascading two different mechanisms. Just below the hypocenter, slip might be controlled by transformational faulting that eventually transfers into thermal runaway outside the grain-reduction zone, where temperature exceeds 1,000°C.

Our model can be extended in the future to allow fracture to propagate into the ductile part of the slab but requires a thorough understanding of the failure criterion. As prior weakening is required for thermal runaway, and grain-size reduction is part of our models, the approach of Thielmann (2018) is probably the best candidate. In this model, the viscous part of the rheology is complex and composed of four components, namely, diffusion creep, dislocation creep, low-temperature plasticity, and dislocation-accommodated grain boundary sliding. A possible localized viscosity reduction in the ductile part of the slab depends on the strain rate and temperature and also produces a stress decrease in the affected region. On the other hand, transformational faulting is still not well understood at shallower depths (50–100 km), as in our model where the 2017 Tehuantepec occurred. A recent study by Sindhusuta et al. (2025) has shown that phase transformation of olivine to spinel results in strong strain localization (spinel is mechanically weaker than olivine), depending of the P-T conditions. This model is more applicable for high-pressure conditions specific for the mantle transition zone. However, the phase diagrams of olivine-to-spinel phase transformation (Ross and Navrotsky, 1987; Sindhusuta et al., 2025) shows that this phase transformation is possible even at lower pressure (0.5–1 GPa) and high-temperature conditions (1,200–1,350 K), which are relatively close to our estimates.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

MM: investigation, methodology, software, validation, visualization, writing – original draft, and writing – review and editing. VM: funding acquisition, investigation, methodology, software, supervision, validation, visualization, writing – original draft, and writing – review and editing. SY: formal analysis, funding acquisition, investigation, resources, supervision, validation, and writing – review and editing. EM: formal analysis, investigation, resources, supervision, validation, visualization, and writing – review and editing. NS: formal analysis, investigation, resources, supervision, validation, and writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was partly supported by the MEXT KAKENHI grant (grant number 21H05203), Kobe University Strategic International Collaborative Research Grant (Type B Fostering Joint Research), and FY2021 Invitational Fellowships for Research in Japan (long term). In addition, this work received support from the UNAM-DGAPA-PAPIIT project grant IN114524 and the JST SATREPS Japan grant number JPMJSA2310.

Acknowledgments

Numerical computations benefited from the computational recourses from the National Laboratory for Advanced Scientific Visualization (LAVIS) at Universidad Nacional Autónoma de México (UNAM) (www.lavis.unam.mx) computing facility and Kobe University computing facilities. This work received support from LAVIS software engineers Luis Alberto Aguilar Bautista, Alejandro de León Cuevas, and Alejandro Avalos. Figures were created using ParaView software (version: ParaView 5.4.1, URL link: https://www.paraview.org/download/).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2025.1535039/full#supplementary-material

References

Abe, K. (1972). Mechanics and tectonic implications of the 1966 and 1970 Peru earthquakes. Phys. Earth Planet. Interiors 5, 367–379. doi:10.1016/0031-9201(72)90108-2

CrossRef Full Text | Google Scholar

Abers, G. A. (2005). Seismic low-velocity layer at the top of subducting slabs: observations, predictions, and systematics. Phys. Earth Planet. Interiors 149 (1-2), 7–29. doi:10.1016/j.pepi.2004.10.002

CrossRef Full Text | Google Scholar

Adriano, B., Fujii, Y., Koshimura, S., Mas, E., Ruiz-Angulo, A., and Estrada, M. (2018). Tsunami source inversion using tide gauge and DART tsunami waveforms of the 2017 Mw8.2 Mexico earthquake. Pure Appl. Geophys 175, 35–48. doi:10.1007/s00024-017-1760-2

CrossRef Full Text | Google Scholar

Aguilar-Anaya, D., Mortera-Gutiérrez, C. A., Berndt, C., and Bandy, W. L. (2025). New insights into geomorphological and tectonic processes in the Gulf of Tehuantepec and constraints on tsunami generation. Geomorphology 473, 109612. doi:10.1016/j.geomorph.2025.109612

CrossRef Full Text | Google Scholar

Amante, C., and Eakins, B. W. (2009). ETOPO1 arc-minute global relief model: procedures, data sources and analysis. US Department of Commerce, National oceanic and atmospheric administration, national environmental satellite, data, and information service. Natl. Geophys. Data Cent. Mar. Geol. Geophys. Div. doi:10.1594/PANGAEA.769615

CrossRef Full Text | Google Scholar

Astiz, L., Lay, T., and Kanamori, H. (1988). Large intermediate-depth earthquakes and the subduction process. Phys. Earth Planet. Interiors 53 (Issues 1–2), 80–166. doi:10.1016/0031-9201(88)90138-0

CrossRef Full Text | Google Scholar

Audet, P., Bostock, M., Christensen, N., and Peacock, S. M. (2009). Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing. Nature 457, 76–78. doi:10.1038/nature07650

PubMed Abstract | CrossRef Full Text | Google Scholar

Behr, W. M., Holt, A. F., Becker, T. W., and Faccenna, C. (2022). The effects of plate interface rheology on subduction kinematics and dynamics. Geophys. J. Int. 230 (2), 796–812. doi:10.1093/gji/ggac075

CrossRef Full Text | Google Scholar

Bercovici, D., and Richard, Y. (2012). Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning. Phys. Earth Planet. Inter. 202–203, 27–55. doi:10.1016/j.pepi.2012.05.003

CrossRef Full Text | Google Scholar

Bercovici, D., and Richard, Y. (2014). Plate tectonics, damage and inheritance. Nature 508 (7497), 513–516. doi:10.1038/nature13072

PubMed Abstract | CrossRef Full Text | Google Scholar

Billen, M. I. (2008). Modeling the dynamics of subducting slabs. Annu. Rev. Earth Planet. Sci. 36, 325–356. doi:10.1146/annurev.earth.36.031207.124129

CrossRef Full Text | Google Scholar

Bostock, M., Hyndman, R., Rondenay, S., and Peacock, S. M. (2002). An inverted continental Moho and serpentinization of the forearc mantle. Nature 417, 536–538. doi:10.1038/417536a

PubMed Abstract | CrossRef Full Text | Google Scholar

Calò, M. (2021). Tears, windows, and signature of transform margins on slabs. Images of the Cocos plate fragmentation beneath the Tehuantepec isthmus (Mexico) using Enhanced Seismic Models. Earth Planet. Sci. Lett. 560, 116788. doi:10.1016/j.epsl.2021.116788

CrossRef Full Text | Google Scholar

Čížková, H., and Bina, C. R. (2013). Effects of mantle and subduction-interface rheologies on slab stagnation and trench rollback. Earth Planet. Sci. Lett. 379, 95–103. doi:10.1016/j.epsl.2013.08.011

CrossRef Full Text | Google Scholar

Craig, T. J., Copley, A., and Jackson, J. (2014). A reassessment of outer-rise seismicity and its implications for the mechanics of oceanic lithosphere. Geophys. J. Int. 197 (1), 63–89. doi:10.1093/gji/ggu013

CrossRef Full Text | Google Scholar

Crameri, F., Magni, V., Domeier, M., Shephard, G. E., Chotalia, K., Cooper, G., et al. (2020). A transdisciplinary and community-driven database to unravel subduction zone initiation. Nat. Commun. 11, 3750. doi:10.1038/s41467-020-17522-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellis, S., and Stöckhert, B. (2004). Elevated stresses and creep rates beneath the brittle-ductile transition caused by seismic faulting in the upper crust. J. Geophys. Res. Solid Earth 109, B05407. doi:10.1029/2003JB002744

CrossRef Full Text | Google Scholar

Faccenda, M., Burlini, L., Gerya, T., and Mainprice, D. (2008). Fault-induced seismic anisotropy by hydration in subducting oceanic plates. Nature 455, 1097–1100. doi:10.1038/nature07376

CrossRef Full Text | Google Scholar

Faccenna, C., Heuret, A., Funiciello, F., Lallemand, S., and Becker, W. T. (2007). Predicting trench and plate motion from the dynamics of a strong slab. Earth Planet. Sci. Lett. 257 (Issues 1–2), 29–36. doi:10.1016/j.epsl.2007.02.016

CrossRef Full Text | Google Scholar

Faryad, S. W., Dolejš, D., and Machek, M. (2009). Garnet exsolution in pyroxene from clinopyroxenites in the Moldanubian zone: constraining the early pre-convergence history of ultramafic rocks in the Variscan orogen. J. Metamorph. Geol. 27, 655–671. doi:10.1111/j.1525-1314.2009.00834.x

CrossRef Full Text | Google Scholar

Ferrand, T., Hilairet, N., Incel, S., Deldicque, D., Labrousse, L., Gasc, J., et al. (2017). Dehydration-driven stress transfer triggers intermediate-depth earthquakes. Nat. Commun. 8, 15247. doi:10.1038/ncomms15247

PubMed Abstract | CrossRef Full Text | Google Scholar

Frohlich, C. (2006). Deep earthquakes. Cambridge University Press. doi:10.1017/CBO9781107297562

CrossRef Full Text | Google Scholar

Furusho, M., and Kanagawa, K. (1999). Transformation-induced strain localization in a lherzolite mylonite from the Hidaka metamorphic belt of central Hokkaido, Japan. Tectonophysics 313, 411–432. doi:10.1016/S0040-1951(99)00215-2

CrossRef Full Text | Google Scholar

García, D., Singh, S. K., Herráiz, M., Pacheco, J. F., and Ordaz, M. (2004). Inslab earthquakes of central Mexico: Q, source spectra and stress drop. Bull. Seismol. Soc. Am. 94, 789–802. doi:10.1785/0120030125

CrossRef Full Text | Google Scholar

Gerya, T. V. (2022). Numerical modeling of subduction: state of the art and future directions. Geosphere 18 (2), 503–561. doi:10.1130/GES02416.1

CrossRef Full Text | Google Scholar

Gerya, T. V., Bercovici, D., and Becker, T. W. (2021). Dynamic slab segmentation due to brittle–ductile damage in the outer rise. Nature 599, 245–250. doi:10.1038/s41586-021-03937-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerya, T. V., and Meilick, F. I. (2011). Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts. J. Metamorph. Geol. 29, 7–31. doi:10.1111/j.1525-1314.2010.00904.x

CrossRef Full Text | Google Scholar

Green, H. W., and Houston, H. (1995). The mechanics of deep earthquakes. Annu. Rev. Earth Planet. Sci. 23, 169–213. doi:10.1146/annurev.ea.23.050195.001125

CrossRef Full Text | Google Scholar

Gusman, A. R., Mulia, I. E., and Satake, K. (2018). Optimum Sea surface displacement and fault slip distribution of the 2017 Tehuantepec earthquake (M_w 8.2) in Mexico estimated from tsunami waveforms. Geophys. Res. Lett. 45, 646–653. doi:10.1002/2017GL076070

CrossRef Full Text | Google Scholar

Haberland, C., Rietbrock, A., Lange, D., Bataille, K., and Dahm, T. (2009). Structure of the seismogenic zone of the southcentral Chilean margin revealed by local earthquake traveltime tomography. J. Geophys. Res. 114, B01317. doi:10.1029/2008JB005802

CrossRef Full Text | Google Scholar

Hacker, B. R., Peacock, S. M., Abers, G. A., and Holloway, S. D. (2003). Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions? J. Geophys. Res. Solid Earth 108, 2030. doi:10.1029/2001JB001129

CrossRef Full Text | Google Scholar

Hasegawa, A., and Nakajima, J. (2017). Seismic imaging of slab metamorphism and genesis of intermediate-depth intraslab earthquakes. Prog. Earth Planet. Sci. 4, 12. doi:10.1186/s40645-017-0126-9

CrossRef Full Text | Google Scholar

Hayes, G., Moore, G., Portner, D., Hearne, M., Flamme, H., Furtney, M., et al. (2018). Slab2, a comprehensive subduction zone geometry model. Science 362, 58–61. doi:10.1126/science.aat4723

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirth, G., and Kohlstedt, D. L. (2003). “Rheology of the upper mantle and the mantle wedge: a view from the experimentalists,” in Inside the subduction factory. Editor J. Eiler (American Geophysical Union), 83–105. doi:10.1029/138GM06

CrossRef Full Text | Google Scholar

Hobbs, B. E., Ord, A., and Teyssier, C. (1986). Earthquakes in the ductile regime? pure Appl. Geophys 124, 309–336. doi:10.1007/BF00875730

CrossRef Full Text | Google Scholar

Hummel, N., Buiter, S., and Erdős, Z. (2024). The influence of viscous slab rheology on numerical models of subduction. Solid Earth 15, 567–587. doi:10.5194/se-15-567-2024

CrossRef Full Text | Google Scholar

John, T., Medvedev, S., Rüpke, L. H., Andersen, T. B., Podladchikov, Y. Y., and Austrheim, H. (2009). Generation of intermediate-depth earthquakes by self-localizing thermal runaway. Nat. Geosci. 2, 137–140. doi:10.1038/ngeo419

CrossRef Full Text | Google Scholar

Karato, S. I., and Wu, P. (1993). Rheology of the upper mantle: a synthesis. Science 260 (5109), 771–778. doi:10.1126/science.260.5109.771

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawakatsu, H. (1986). Double seismic zones: kinematics. Jour. Geophys. Res. 91, 4811–4825. doi:10.1029/jb091ib05p04811

CrossRef Full Text | Google Scholar

Kawakatsu, H., and Watada, S. (2007). Seismic evidence for deep-water transportation in the mantle. Science 316, 1468–1471. doi:10.1126/science.1140855

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelemen, P. B., and Hirth, G. (2007). A periodic shear-heating mechanism for intermediate-depth earthquakes in the mantle. Nature 446, 787–790. doi:10.1038/nature05717

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirby, S. H. (1980). Tectonic stresses in the lithosphere: constraints provided by the experimental deformation of rocks. J. Geophys. Res. Solid Earth 85, 6353–6363. doi:10.1029/JB085iB11p06353

CrossRef Full Text | Google Scholar

Kirby, S. H., Engdahl, E. R., and Denlinger, R. (1996). “Intermediate-depth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs,” in Subduction: top to bottom. Editor G. E. Beboutet al. (American Geophysical Union), 195–214. doi:10.1029/GM096p0195

CrossRef Full Text | Google Scholar

Kirby, S. H., and Kronenberg, A. K. (1987). Rheology of the lithosphere: selected topics. Rev. Geophys. 25, 1219–1244. doi:10.1029/RG025i006p01219

CrossRef Full Text | Google Scholar

Kita, S., and Katsumata, K. (2015). Stress drops for intermediate-depth intraslab earthquakes beneath Hokkaido, northern Japan: differences between the subducting oceanic crust and mantle events. Geochem. Geophys. Geosystems 16, 552–562. doi:10.1002/2014GC005603

CrossRef Full Text | Google Scholar

Lee, C., and Kim, Y. (2021). Role of warm subduction in the seismological properties of the forearc mantle: an example from southwest Japan. Sci. Adv. 7, eabf8934. doi:10.1126/sciadv.abf8934

PubMed Abstract | CrossRef Full Text | Google Scholar

Manea, M., and Manea, V. C. (2008). On the origin of El Chichón volcano and subduction of Tehuantepec Ridge: a geodynamical perspective. J. Volcanol. Geotherm. Res. 175 (4), 459–471. doi:10.1016/j.jvolgeores.2008.02.028

CrossRef Full Text | Google Scholar

Manea, M., Manea, V. C., Ferarri, L., Kostoglodov, V., and Bandy, W. L. (2005). Tectonic evolution of the Tehuantepec ridge. Earth Planet. Sci. Lett. 238, 64–77. doi:10.1016/j.epsl.2005.06.060

CrossRef Full Text | Google Scholar

Manea, M., Manea, V. C., and Kostoglodov, V. (2003). Sediment fill in the Middle America Trench inferred from gravity anomalies. Geofísica Int. 42, 603–612. doi:10.22201/igeof.00167169p.2003.42.4.314

CrossRef Full Text | Google Scholar

Masson, D. G. (1991). Fault patterns at outer trench walls. Mar. Geophys. Res. 13, 209–225. doi:10.1007/BF00369150

CrossRef Full Text | Google Scholar

Melgar, D., Ruiz-Angulo, A., Garcia, E. S., Manea, M., Manea, C. V., Xu, X., et al. (2018). Deep embrittlement and complete rupture of the lithosphere during the M_w 8.2 Tehuantepec earthquake. Nat. Geosci. 11, 955–960. doi:10.1038/s41561-018-0229-y

CrossRef Full Text | Google Scholar

Meng, L., Huang, H., Xie, Y., Bao, H., and Dominguez, L. A. (2019). Nucleation and kinematic rupture of the 2017 Mw 8.2 Tehuantepec earthquake. Geophys. Res. Lett. 46, 3745–3754. doi:10.1029/2018GL081074

CrossRef Full Text | Google Scholar

Mulyukova, E., and Bercovici, D. (2019). The generation of plate tectonics from grains to global scales: a brief review. Tectonics 38, 4058–4076. doi:10.1029/2018tc005447

CrossRef Full Text | Google Scholar

Nakajima, J., Hirose, F., and Hasegawa, A. (2009). Seismotectonics beneath the Tokyo metropolitan area, Japan: effect of slab-slab contact and overlap on seismicity. J. Geophys. Res. 114, B08309. doi:10.1029/2008JB006101

CrossRef Full Text | Google Scholar

Nakanishi, M. (2017). Bending-related topographic structures of the incoming oceanic plate at the northwestern pacific ocean trench. J. Geogr. (Chigaku Zasshi) 126, 125–146. doi:10.5026/jgeography.126.125

CrossRef Full Text | Google Scholar

Ogawa, M. (1987). Shear instability in a viscoelastic material as the cause of deep focus earthquakes. J. Geophys. Res. Solid Earth 92, 13801–13810. doi:10.1029/JB092iB13p13801

CrossRef Full Text | Google Scholar

Okuwaki, R., and Yagi, Y. (2017). Rupture process during the M_w 8.1 2017 Chiapas Mexico earthquake: shallow intraplate normal faulting by slab bending. Geophys. Res. Lett. 44 (11), 816–823. doi:10.1002/2017GL075956

CrossRef Full Text | Google Scholar

Ranero, C. R., Phipps Morgan, J., McIntosh, K., and Reichert, C. (2003). Bending-related faulting and mantle serpentinization at the Middle America trench. Nature 425, 367–373. doi:10.1038/nature01961

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranero, C. R., Villaseñor, A., Phipps Morgan, J., and Weinrebe, W. (2005). Relationship between bend-faulting at trenches and intermediate-depth seismicity. Geochem. Geophys. Geosystems 6, Q12002. doi:10.1029/2005GC000997

CrossRef Full Text | Google Scholar

Rebollar, C. J., Quintanar, L., Yamamoto, J., and Uribe, A. (1999). Source process of the Chiapas, Mexico, intermediate-depth earthquake (Mw = 7.2) of 21 october 1995. Bull. Seismol. Soc. Am. 89, 348–358. doi:10.1785/BSSA0890020348

CrossRef Full Text | Google Scholar

Ross, N. L., and Navrotsky, A. (1987). The Mg2GeO4 olivine-spinel phase transition. Phys. Chem. Minerals 14 (5), 473–481. doi:10.1007/bf00628825

CrossRef Full Text | Google Scholar

Sdrolias, M., and Müller, R. D. (2006). Controls on back-arc basin formation. Geochem. Geophys. Geosyst. 7, Q04016. doi:10.1029/2005GC001090

CrossRef Full Text | Google Scholar

Seno, T. (2009). Intraslab seismicity and its generation mechanisms. Zisin J. Seismol. Soc. Jpn. 2nd ser. 61, 357–364. doi:10.4294/zisin.61.357

CrossRef Full Text | Google Scholar

Seno, T., and Yamanaka, Y. (1996). “Double seismic zones, deep compressional trench—outer rise events and superplumes,” in Subduction top to bottom. Editors G. E. Bebout, D. W. Scholl, S. H. Kirby, and J. P. Platt (Washington, DC: AGU), 96, 347–355.

Google Scholar

Shi, F., Wang, Y., Wen, J., Yu, T., Zhu, L., Huang, T., et al. (2022). Metamorphism-facilitated faulting in deforming orthopyroxene: implications for global intermediate-depth seismicity. Proc. Natl. Acad. Sci. 119, e2112386119. doi:10.1073/pnas.2112386119

PubMed Abstract | CrossRef Full Text | Google Scholar

Sindhusuta, S., Chi, S.-W., Foster, C., Officer, T., and Wang, Y. (2025). Numerical investigation into mechanical behavior of metastable olivine during phase transformation: implications for deep-focus earthquakes. J. Geophys. Res. Solid Earth 130, e2024JB030557. doi:10.1029/2024JB030557

CrossRef Full Text | Google Scholar

Singh, S. K., Pérez-Campos, X., Espíndola, V. H., Cruz-Atienza, V. M., and Iglesias, A. (2014). Intraslab earthquake of 16 june 2013 (M_w 5.9), one of the closest such events to Mexico city. Seismol. Res. Lett. 85 (2), 268–277. doi:10.1785/0220130179

CrossRef Full Text | Google Scholar

Singh, S. K., Pérez-Campos, X., Espindola, V. H., Iglesias, A., and Quintanar, L. (2020). An intraslab earthquake at a depth of 100 km in the subducting Cocos plate beneath Nevado de Toluca volcano. Geofísica Int. 59, 5–12. doi:10.22201/igeof.00167169p.2020.59.1.2072

CrossRef Full Text | Google Scholar

Singh, S. K., Ponce, L., and Nishenko, S. P. (1985). The great Jalisco, Mexico, earthquakes of 1932: subduction of the Rivera plate. Bull. Seismol. Soc. Am. 75, 1301–1313. doi:10.1785/bssa0750051301

CrossRef Full Text | Google Scholar

Sippl, C., Dielforder, A., John, T., and Schmalholz, S. M. (2022). Global constraints on intermediate-depth intraslab stresses from slab geometries and mechanisms of double seismic zone earthquakes. Geochem. Geophys. Geosystems 23, e2022GC010498. doi:10.1029/2022GC010498

CrossRef Full Text | Google Scholar

Song, C., and Ge, Z. (2019). 3D model backprojection of the 2017 Mw 8.2 Chiapas earthquake: a two-stage rupture with a barrier-induced velocity increase. Seismol. Res. Lett. 90 (3), 1121–1130. doi:10.1785/0220180268

CrossRef Full Text | Google Scholar

SSN (2017). M8.2–133 km al Suroeste de Pijijiapan, Chis. Mexico City, Mexico: México, Servicio Sismológico Nacional. Available online at: http://www2.ssn.unam.mx:8080/catalogo/ (Accessed November, 2017).

Google Scholar

Suárez, G. (2021). Large earthquakes in the Tehuantepec subduction zone: evidence of a locked plate interface and large-scale deformation of the slab. J. Seismol. 25, 449–460. doi:10.1007/s10950-020-09969-6

CrossRef Full Text | Google Scholar

Suárez, G., Santoyo, M. A., Hjorleifsdottir, V., Iglesias, A., Villafuerte, C., and Cruz-Atienza, V. M. (2019). Large scale lithospheric detachment of the downgoing Cocos plate: the 8 September 2017 earthquake (M 8.2). Earth Planet. Sci. Lett. 509, 9–14. doi:10.1016/j.epsl.2018.12.018

CrossRef Full Text | Google Scholar

Suenaga, N., Yoshioka, S., and Ji, Y. (2021). 3-D thermal regime and dehydration processes around the regions of slow earthquakes along the Ryukyu Trench. Sci. Rep. 11, 11251. doi:10.1038/s41598-021-90199-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, K., Grand, S. P., and Niu, F. (2018). Seismic structure of the upper mantle beneath eastern Asia from full waveform seismic tomography. Geochem. Geophys. Geosystems 19, 2732–2763. doi:10.1029/2018GC007460

CrossRef Full Text | Google Scholar

Thielmann, M. (2018). Grain size assisted thermal runaway as a nucleation mechanism for continental mantle earthquakes: impact of complex rheologies. Tectonophysics 746, 611–623. doi:10.1016/j.tecto.2017.08.038

CrossRef Full Text | Google Scholar

van Keken, P. E., Currie, C., King, S. D., Behn, M. D., Cagnioncle, A., He, J., et al. (2008). A community benchmark for subduction zone modeling. Phys. Earth Planet. Interiors 171 (Issues 1–4), 187–197. doi:10.1016/j.pepi.2008.04.015

CrossRef Full Text | Google Scholar

van Keken, P. E., and Wilson, C. R. (2023). An introductory review of the thermal structure of subduction zones: III—comparison between models and observations. Prog. Earth Planet Sci. 10, 57. doi:10.1186/s40645-023-00589-5

CrossRef Full Text | Google Scholar

Watts, A. B. (2001). Isostasy and flexure of the lithosphere. Cambridge University Press.

Google Scholar

Ye, L., Lay, T., Bai, Y., Cheung, K. F., and Kanamori, H. (2017). The 2017 M_w 8.2 Chiapas, Mexico, earthquake: energetic slab detachment. Geophys. Res. Lett. 44 (11), 824–832. doi:10.1002/2017GL076085

CrossRef Full Text | Google Scholar

Zhang, H., and Brudzinski, M. R. (2019). Evidence for rupture through a double benioff zone during the 2017M_w8.2 Chiapas, Mexico earthquake. Geophys. Res. Lett. 46, 652–660. doi:10.1029/2018GL080009

CrossRef Full Text | Google Scholar

Zhao, D., Hasegawa, A., and Kanamori, H. (1994). Deep structure of Japan subduction zone as derived from local, regional, and teleseismic events. J. Geophys. Res. Solid Earth 99 (B11), 22313–22329. doi:10.1029/94JB01149

CrossRef Full Text | Google Scholar

Zhong, S., and Davies, G. F. (1999). Effects of plate and slab viscosities on the geoid. Earth Planet. Sci. Lett. 170, 487–496. doi:10.1016/S0012-821X(99)00124-7

CrossRef Full Text | Google Scholar