The influence of mantle hydration and flexure on slab seismicity in the southern Central Andes - Communications Earth & Environment

In subduction zones, slab seismicity along the megathrust interface (inter-plate seismicity) and within the oceanic crust and mantle (intra-slab seismicity) represents one of the sources of the greatest seismic hazards on Earth. An understanding of the causative mechanisms behind the spatial distribution of these earthquakes (e.g. extent of the seismogenic zone) allows an improved evaluation of the seismic hazard in these convergent settings. Thermomechanical modelling in various tectonic settings, including rifts, and transform and convergent margins^1,2,3, has shown how the long-term mechanical strength of the lithosphere may control the distribution of seismicity, with the brittle-ductile transition (BDT) generally identified as the downdip limit to the seismogenic zone^4,5. However, oceanic plates in subduction zones generate seismicity at depths greater than 50 km, where lithostatic pressure and temperature conditions should tend to inhibit brittle failure^6,7. In recent years, a growing number of studies has discussed the role of weakening due to fluids released during dehydration processes in the oceanic plate as the cause for such intermediate-depth earthquakes (50–300 km)^8,9,10,11. While the exact physical mechanism is still a matter of debate, the main processes proposed to enable brittle faulting under hydrated conditions include dehydration embrittlement^8,9,10 and dehydration-related stress transfer¹¹. However, geological observations from fossil subduction zones¹² and from rock mechanics experiments¹³ have shown that brittle failure causing intermediate-depth earthquakes can also occur under dry conditions. Thermal runaway slip weakening¹⁴ provides a plausible mechanism to explain their occurrence.

Irrespectively of the exact physical process involved, seismicity is triggered in response to the build-up and distribution of stresses in the slab. These reflect a balance between multiple forces, including buoyancy forces associated with thermal and metamorphic density contrasts, frictional resistance along the subduction interface and slab, and flexural stresses related to changes in the subduction angle^8,9,15,16,17. Quantifying the relative importance of dehydration mechanisms and slab stress in triggering slab seismicity requires a comprehensive integration of available geological, seismological, and geodetic observations.

The convergent margin of the southern Central Andes (29°S–39°S), where the oceanic Nazca Plate subducts under the continental South American Plate, is a subduction environment where a high number of great earthquakes have been recorded (e.g. 2010 M8.8 Maule and 2015 M8.2 Illapel)^18,19. A distinctive feature of this region is the change in the subduction angle of the Nazca Plate between 33°S and 35°S from the Chilean-Pampean flat-slab zone (<5° dip, 27–33°S) in the north to a steeper segment in the south (~30° dip)²⁰ (Fig. 1a). Domains with the highest frequency of slab seismicity²¹ are associated with the oceanic slab subducting sub-horizontally or with a shallow dip. These domains are located either close to the trench or along the flat-slab segment, below the Andean Precordillera fold-and-thrust belt and the thick-skinned Sierras Pampeanas morphotectonic provinces (Fig. 1b).

a Main panel: topography and bathymetry of the southern Central Andes based on the ETOPO1 global relief model⁸⁹. The white rectangle encloses the extent of the area considered in this study (area of Fig. 1b); the black dashed rectangle encloses the area shown in Fig. 1c; the purple rectangle encloses the area shown in Fig. 5. Black-dashed lines in the oceanic domain mark the boundaries between subduction segments (flat slab in the north, steep slab in the south). Fracture zones of the Nazca Plate⁷⁵ are highlighted by solid black lines. Abbreviations of main fracture zones: MFZ = Mocha Fracture Zone; CFZ = Challenger Fracture Zone. White lines show the isobath contours (km below mean sea level) of top of the slab (based on the Slab2 model⁷⁵). The wide transparent white line denotes the position of offshore (full line) and projected onshore (dashed line) tracks of the Juan Fernandez Ridge (JFR)⁹⁰. The red triangles depict locations of active volcanoes⁹¹. The orange lines denote locations of profiles shown in Fig. 2. The boundaries between the main morphotectonic provinces are shown by white dashed lines. Abbreviations of main tectonic provinces: AO Andean orogen, CB Cuyo Basin, ESP Eastern Sierras Pampeanas, EAB extra-Andean basins, FA forearc, NB Neuquén Basin, P Payenia volcanic province, Prc Precordillera, WSP Western Sierras Pampeanas. Left panel: extent of historical seismic ruptures⁹². Seismic gaps (SG) are shown with red lines. b Seismic-event density (logarithmic scale) derived for a circular bin of 75 km radius for slab seismicity²¹. Seismicity marked by red dots. Isobath contours of the top of the slab are drawn in grey. c Degree of locking of the subduction interface⁶⁷.

Full size image

In the southern Central Andes, a recent thermo-mechanical modelling study²² constrained the downdip extent of the brittle-ductile transition to 40–55 km depth (Supplementary Fig. 1), thus implying that the spatial distribution of slab seismicity at depths <50 km is expected from variations in the long-term strength of the lithosphere. However, seismicity also occurs at depths >50–70 km, where ductile rheological conditions are expected. Though numerous studies have contributed to the current understanding of the causative dynamics driving slab seismicity below the brittle-ductile transition in different subduction zones^8,9,11,15, they cannot resolve this apparent inconsistency and a gap in knowledge remains with respect to the southern Central Andes.

In this study, we carry out an integrative analysis of the relationship between the observed seismicity distribution, seismic properties, and the 3D geometrical configuration of the oceanic plate to investigate potential mechanisms influencing the seismic behaviour of the slab. We take advantage of a recent full waveform seismic tomography model of the southern Central Andes²³ to demonstrate the existence of a causative relationship between the state of hydration of the mantle and recorded slab seismicity, both for intra-slab and inter-plate events^21,24,25 (see “Methods” section). We further investigate how in regions of likely slow reaction kinetics and small fluid input to the overriding plate seismic activity can be locally enhanced by flexural stress associated with mapped variations in slab geometry. Overall, although a direct comparison of seismic velocity structure and seismic behaviour is generally well established, the integrative framework presented in our investigation adds thermomechanical and geological constraints to minimise the degree of freedom and non-uniqueness of geodynamical interpretations derived from geophysical modelling.

To infer the driving physical mechanisms responsible for brittle failure requires knowledge of parameters such as the state of hydration of the mantle and slab geometry. Seismic tomography, due to the sensitivity of P- and S-wave velocities to lithology, fluid/melt-content and temperature, can be used as a proxy to investigate the state of hydration of the mantle^26,27,28 (Supplementary Fig. 2). Specifically, low Vp/Vs (1.65–1.8) characterise dry peridotite²⁹, whereas high Vp/Vs (>2) indicate a hydrated mantle with free fluids^30,31. Intermediate Vp/Vs values between 1.8 and 1.9 likely correspond to a serpentinized mantle^29,32,33,34. In the oceanic mantle, serpentinization may result from water penetration along normal faults cross-cutting the slab;³⁵ conversely, in the mantle wedge, serpentinite is formed from fluids migrating upwards from the subducting slab³⁶. Ascent of fluids from the subducting slab through the mantle of the overriding plate may drive partial melting below the volcanic arc, likewise leading to an increase in Vp/Vs (>1.9)³⁷. Finally, very low Vp/Vs (<1.65) may be indicative of a hydrated and anisotropic mantle, the foliation of antigorite or the presence of a network of preferentially oriented cracks^30,38,39.

Here, we analyse the Vp/Vs distribution from the full-waveform inversion tomography model²³. In the following, we discuss our results relying on (i) three cross sections perpendicular to the trench, considered representative of the different subduction segments (Fig. 2a–c); (ii) a N-S cross section along the forearc (Fig. 2d); and (iii) a map along a constant depth of 90 km, which corresponds to the mantle of the overriding plate immediately above the slab in the flat-subduction segment (Fig. 3). Although the tomographic inversion solves for Vs radial anisotropy (i.e., Vsv and Vsh) to avoid introducing artifacts into the isotropic model, it is not robustly constrained⁴⁰ and we thus refrain from interpreting it.

a–d Vp/Vs ratio along cross-sections A–D (see Figs. 1 or 3 for profile locations). Areas with poor resolution are masked with semi-transparent white shading. Dots show the location of hypocentres of slab seismicity²¹ within 5 km distance from each profile. Pale pink dots denote seismicity due to brittle stick-slip assisted by dehydration, magenta dots mark dehydration-related seismicity, while black dots show events attributed to flexure-assisted seismicity. The interfaces of the lithospheric model units are marked by black lines, except for the Moho discontinuity, which is shown by red lines. On top of each profile, topography of the upper plate is displayed with a vertical exaggeration of 10:1. Right panels in a–c and bottom panel in d show the histogram of slab seismicity and the average Vp/Vs within the slab as a function of depth and latitude, respectively.

Full size image

Additionally shown are profile locations of Fig. 2 (black lines), depth contours to the top of the slab at 30-km intervals (white thin lines, except 90 km-depth contour in yellow), the dominant orientation of the T-axes of focal mechanisms^24,25 projected horizontally (white thick bars), the projected onshore (thick grey line) tracks of the Juan Fernandez Ridge (JFR)⁹⁰, the active magmatic arc (red triangles), and the coastline (blue line). For abbreviations of morphotectonic provinces, see Fig. 1.

Full size image

The W-E cross sections depict an across-strike variability in Vp/Vs indicative of progressive eastward dehydration of the oceanic mantle (Fig. 2). The southernmost part of the flat-slab segment (32°S–33°S) shows a very low ratio (<1.65) at depths shallower than 50 km, which can arise in an anisotropic mantle, where anisotropy could be due to the foliation of antigorite or the presence of an oriented network of fluid-filled cracks^30,38,39 (S1,O1; Fig. 2a). Between 50 and 100 km, the ratio within the slab and the mantle of the overriding plate increases from intermediate (1.65–1.8; S2, O2; Fig. 2a) to high values (1.8-1.9; S3, O3; Fig. 2a), indicating mantle serpentinization. At ~80 km, the Vp/Vs ratio reaches values >1.9 (S4, O4; Fig. 2a), suggesting hydrated conditions with a free fluid phase. Deeper than 100 km, the ratio decreases from values around 1.8 (S5; Fig. 2a) in the central part of the slab to intermediate and low values (1.65–1.8; S6, O5, O6; Fig. 2a) within the top part of the slab and within the overlying mantle. While the Vp/Vs values in the central part of the slab likely indicate a serpentinized mantle, the O6 anomaly seems to indicate either dry or anisotropic hydrated conditions. A dry slab that underwent complete eclogitization and deserpentization in its topmost domain would imply high densities; this, in turn, is difficult to reconcile with the buoyancy needed to maintain the sub-horizontal subduction angle of the slab in these domains. Our preferred interpretation is that the slab is still hydrated, but the decrease in Vp/Vs may be partly due to wave-propagation artefacts in anisotropic minerals. Indeed, antigorite is stable at upper-mantle pressures and temperatures <650 °C³⁶, which is consistent with the conditions modelled in this region of the flat slab segment (~600 °C)^41,42. However, Anderson et al.⁴³ reported moderate anisotropy (~4%) for this part of the slab, which is inconsistent with the very strong seismic anisotropic signature (>30%)³² of foliated antigorite. Instead, these authors attributed the medium-scale seismic anisotropy to sheared olivine. In this context, the combined observation of low to intermediate Vp/Vs ratios and moderate anisotropy in the flat slab could indicate hydrated peridotite with isotropic antigorite and anisotropic olivine⁴⁴.

For the mantle of the continental plate above the flat slab, there is little information on seismic anisotropy to draw more extensive conclusions on its hydration state. We can speculate that dehydration reactions within the slab being kinetically delayed due to low temperatures would imply less fluid release detectable by seismic methods in the overlying lithospheric mantle, which could explain the observed low Vp/Vs in our study. An additional factor contributing to the smaller fluid content within the mantle of the overriding plate could be the pressure gradient generated by the slab re-steepening, which would cause the fluids released by the dehydration reactions in the central part of the slab to be diverted from their expected vertical pathways further to the east⁴⁵. The increase in Vp/Vs (~1.7) from low values at 68.5°W to high (~1.8) values at 67.5°W (Fig. 2a) lends support to this hypothesis.

Unlike the flat-slab segment, the transitional segment to the steeper subduction angle displays a more uniform Vp/Vs distribution. Below 50 km, an increase from low and intermediate (1.65–1.9; S7, S8, O7, O8) to high (>1.9; S9, O9) values marks the transition from serpentinized³² rocks to a mantle with a free fluid phase³⁰ (Fig. 2b). At depths >100–150, Vp/Vs decreases to values <1.8 (S10, S11, O10, O11; Fig. 2b), which indicates a slab that has already been dehydrated²⁹.

The steep-slab segment depicts a region of intermediate Vp/Vs values (1.8–1.9; S12, O12; Fig. 2c) within the slab domain proper and the overriding plate mantle extending down to 60 km, which again indicates a serpentinized³² mantle. Between 75 and 120 km, we notice a local Vp/Vs low (1.7–1.8) within the slab (S13) and the mantle of the overriding plate (O13), suggesting dry conditions²⁹ (Fig. 2c). However, to interpret a fully dehydrated mantle at these depths is difficult, especially given the observed Vp/Vs increase (>1.8) below 120 km (S14, O14, O15; Fig. 2c) and the presence of the active volcanic arc immediately above in the upper plate. Our preferred alternative scenario to reconcile the observed intermediate Vp/Vs with the overlying volcanic activity is the presence of hydrated but anisotropic mantle.

The southern Central Andes also display large along-strike variability in Vp/Vs (Fig. 2d–e). For instance, at 90 km depth, the orogenic domain displays intermediate (1.8–1.9) to high (>1.9) values only in its central part (between 31°S–35°S, S8, S12, O8, O12; Fig. 2d), while in the northern and southern regions Vp/Vs decreases to values between 1.65-1.8 (S1, S15, O1, O15; Fig. 2d). We can explain this enhanced hydration in the central orogenic domain with its proximity to the Juan Fernandez Ridge, a seamount chain subducting at 33°S (Fig. 1a). This ridge is associated with pervasive faulting⁴⁶, which enables fluids to reach the and hydrate the oceanic mantle., These fluids are released through dehydration reactions after subduction and ultimately enhance hydration of the overriding mantle⁴⁷. South of 33°S, the active magmatic arc is located near the western edge of an extended region of high (>1.8–1.9) Vp/Vs values in the mantle wedge at 90 km depth (Fig. 3). This suggests that the observed increase in Vp/Vs can indeed be associated with the presence of melts due to the reopening of the mantle wedge and the re-steepening of the Nazca Plate^42,48,49. Finally, Vp/Vs is high (1.84–1.9) beneath the Precordillera and the western Sierras Pampeanas (29°S–31°S, Fig. 2e), which spatially correlates with a zone of low Vs (Supplementary Fig. 3) interpreted by Gao et al. ²³ as a region of slab thinning and tearing along the inland projection of the Juan Fernandez Ridge.

Intra-slab seismicity, fluids and flexure

We now compare the Vp/Vs distribution with the distribution of inter-plate and intra-slab seismicity. First, for the subduction segment in the transition zone (33°S–35°S), slab seismicity at depths shallower than 120 km is limited to the slab domain with high Vp/Vs (>1.8, S8–S10, Fig. 2b). Secondly, seismicity also occurs below a region of hydrated continental mantle with a high velocity ratio (O8–O10, Fig. 2b, d; O12, Fig. 2d). This points to a spatial correlation between slab seismicity and mantle hydration in this region. Seismicity ceases quite abruptly where the slab undergoes (nearly) complete dehydration (as indicated by a lower Vp/Vs) and is not able to release any more fluids to the mantle of the overriding continental plate (e.g. S11, O11; Fig. 2b; S15, O15, Fig. 2d). Likewise, the observed change from medium (1.7–1.8) to high (>1.9) Vp/Vs in the overriding mantle spatially correlates with the location of the active volcanic arc (Figs. 2d and 3) in the transition zone and steep-slab segments, at 69°W–70°W, which confirms the presence of fluid-related melts. In short, fluids responsible for the generation of partial melts, which feed the volcanic arc south of 33°S, are released west of 70°W and the downdip limit of intra-slab seismicity is actively controlled by the availability of fluids in the slab mantle.

The spatial relationship between the regions of mantle hydration and the occurrence of slab seismicity showcases how fluids play a primary role in the promotion of intermediate depth seismicity. The release of fluids increases pore pressure and enhances frictional weakening of the slab, thereby triggering intra-slab seismicity even in regions where brittle failure should be inhibited given local conditions of lithostatic pressure and temperature. In this context, metamorphic reactions of constitutive minerals in the oceanic crust (e.g. amphibole and epidote breakdown during eclogitization^9,50) and in the mantle (e.g. serpentine breakdown^10,28) can lower the rock frictional strength, thereby offering an explanation for the observed intra-slab seismicity at intermediate depths (50–300 km). Fluids from these reactions can migrate trenchward, enabling intraplate earthquakes at shallower depths (~40 km) through a variety of mechanisms^51,52: (i) pore-pressure increase (dehydration embrittlement)^8,9,10,53,54; (ii) grain-size reduction¹³; and (iii) dehydration-related stress transfer from the hydrous phases to the dry host rock¹¹. Since these reactions all affect the state of hydration of the slab and overlying lithospheric mantle, it should be expected that hydrated regions are seismically more active than dry regions^55,56,57,58, as observed in our analysis. Similar relationships between seismic velocities and the distribution of seismicity have been found in other subduction zones, such as Cascadia^55,56, Hirukangi⁵⁷ and the north-western Pacific margin⁵⁸.

In the flat-slab segment (29°S–33°S), the relationship between slab seismicity and hydration of the mantle is not straightforward. For the flat-slab region between 31°S and 32°S, seismicity shallower than 80 km occurs where the slab is hydrated (S1–S3, Fig. 2a). However, this is not the case for the seismicity at larger depths (80–100 km), east of the mantle wedge corner, beneath the Precordillera and the Sierras Pampeanas (S6, Fig. 2a). Here, the upper cluster of seismicity close to the top surface of the slab is localised below a domain in the continental mantle and the subduction interface that is characterised by low Vp/Vs (<1.8, O5, O6) and, therefore, by low fluid content (Figs 2a, 4 and 5). Unlike the southern part of the flat-slab segment, the overriding plate mantle, and the subduction interface in the northern part (29°S–31°S) are characterised by high Vp/Vs (1.8); these areas are thus likely hydrated, while the underlying subducting plate is aseismic (Figs. 3–5). These two observations (i.e., seismic and partially anhydrous overriding plate mantle in the southern flat segment versus aseismic and hydrated overriding plate mantle in the northern flat segment) make it difficult to link the observed intra-slab seismicity in the flat-slab segment to dehydration reactions^59,60. When examining the earthquake focal-mechanism solutions in the region^24,25 (Fig. 5a), we observe a lack of compressional mechanisms, which also rules out the effect of plate convergence as the main trigger for these earthquakes. Seismicity in these domains could instead be related to enhanced bending and unbending stresses in the presence of a non-fluid related weakening mechanism where the slab resumes its steep subduction angle east of 65°W, beneath the Sierras Pampeanas and, to the south, beneath the Cuyo Basin. A possible mechanism is thermal shear instability, where a positive feedback between strain and temperature reduces friction and favours brittle rupture^13,14. This hypothesis agrees with the findings of previous studies^15,16 suggesting that, in regions of intermediate to low pore-fluid pressure, flexural stresses can be more relevant than fluids in triggering seismicity. Variations in the geometry of the subducting slab may have triggered seismicity within its bending and unbending parts, if the effective strain rates associated with slab flexure exceed those associated with slab pull^15,16,61,62. Additional evidence for such a control by slab flexure on the occurrence of seismicity comes from the observed spatial correlation between the flexural regime and the orientation of the earthquake moment tensors^15,16. For example, slab segments with positive curvature (i.e. concave down) gradients are usually characterised by earthquake moment tensors with downdip tensional axes in the upper part of the slab, whereas those with negative curvature gradients display seismicity associated with downdip compressional axes (Supplementary Fig. 4)^4,16. Indeed, we found a rather striking spatial correlation between seismicity and slab geometry in the southern Central Andes based on the orientation of the P- and T-axes of focal mechanisms^24,25 (Fig. 5b). The focal mechanism solutions with a down-dip orientation of the P-axis spatially correlate spatially with the portions of the slab characterised by a negative curvature gradient due to unbending in the transition from flat to steeper subduction angles (Fig. 5b), and equivalently for down-dip T-axes indicating bending. This supports the notion that inter-plate seismicity in these regions is triggered by internal deformation of the slab, as required by the flat-slab geometry.

a 3D diagram of Vp/Vs distribution reflecting the state of hydration within the mantle and the subduction interface, and slab seismicity colour-coded by hypocentral depth. The extent of labelled dry or hydrated slab only pertains to the cross-sections in Fig. 2a, b and strongly varies along the strike of the subduction zone. The dashed white line shows the brittle ductile transition (BDT)²². Magenta-coloured contours show the areas with a high degree of locking⁶⁷. The 3D lithospheric layers used for this graph are based on the structural model of Rodriguez Piceda et al.⁹³. b Vp/Vs distribution in the region of the flat-slab to steep-dip subduction transition overlain with focal mechanism solutions of inter-plate and intra-slab earthquakes^24,25, colour-coded by depth and with white and black dots indicating P and T axes, respectively. The inter-plate events are marked with a magenta circular outline. The isobath contours of the top of the slab⁷⁵ are shown in white lines.

Full size image

Curvature gradient of the top of the slab in the transition between flat and steep-slab subduction (Supplementary Note 2), superposed with a focal mechanism solutions of inter-plate and intra-slab earthquakes^24,25 colour-coded by depth (dots show the orientation of the P- (black) and T-(white) axes); and b P-(red) and T-(blue) axes projected horizontally. Only the sub-horizontal axes (plunge <45°) are shown. The inter-plate events in a are bordered in magenta. The isobath contours of the top of the slab are shown in grey lines. The geometry of the top of the slab was taken from Rodriguez Piceda et al.⁹³, based on Hayes et al.⁷⁵.

Full size image

In the steep-slab segment (35°S–39°S), seismicity is generally limited to domains of hydrated mantle with high Vp/Vs (S12, S14; Fig. 2c). However, similarly to events below the Precordillera and Sierras Pampeanas, the largest seismic cluster in the steep-slab segment correlates at 35°S with a local low in Vp/Vs at 80–120 km (S13, Fig. 2c). Unlike the flat-slab cluster, the steep-slab segment lacks unbending portions below 20 km, making it difficult to relate these events to flexural stresses. An alternative explanation to reconcile the low Vp/Vs ratio with the seismic activity is to consider that the mantle here is hydrated but anisotropic. To date there is no direct evidence of anisotropy in this part of the slab, hence further work targeting the anisotropic signature of the steep-slab segment is needed to fully explain this feature.

Interplate seismicity, fluids and locking degree

A wide range of inter-plate seismic phenomena have been attributed to the effects of high pore-fluid pressure, including deepening of the updip limit of background seismicity and trenchward propagation of megathrust earthquakes⁶³. Typically, an increase in pore pressure is linked to release of fluids, which can lead to enhanced frictional weakening of the subduction interface. In this case, fluid sources encompass: (i) pore water in the sediments released by compaction and (ii) bound water within hydrous minerals of the sediments, oceanic crust, and mantle released by dehydration reactions⁶³. Dewatering due to sediment compaction is restricted to the shallow part of the fault interface (3–7 km depth)⁶⁴, while fluids from dehydration reactions become predominant below this depth. For the Chilean margin between 37°S and 42°S, previous studies have identified a correlation between the state of hydration of the mantle and pore-fluid content, the spatial extent of locked zones along the subduction interface⁶⁵, and the recurrence time of great earthquakes⁶⁶. They found that domains of high Vp/Vs, interpreted as being dominated by high pore-fluid pressure, correlate with areas of a low degree of locking⁶⁵, which are able to slip predominantly aseismically. In contrast, interfaces with low Vp/Vs and thus inferred low pore-fluid pressure, display a high degree of locking⁶⁵ and tend to cause seismic slip less frequently. Our findings support these observations and expand their spatial coverage. For the domains of the southern Central Andes between 29°S and 39°S we have been able to detect a similar correlation between low Vp/Vs and locked areas⁶⁷ along the plate interface (Fig. 4a and Supplementary Fig. 5). This suggests that variations in pore-fluid pressure should be considered as the main factor controlling pre-seismic locking along the plate interface in the southern Central Andes. In agreement with previous studies⁶⁵, we can correlate the presence of N–S variations in pore-fluid pressures with the distribution of fracture zones in the Nazca plate (Fig. 1a). Due to this structural control, slab regions that underwent a higher degree of fracturing are prone to contain more fluids and are characterised by a higher degree of hydrologic connectivity. This in turn translates into higher Vp/Vs values and a lower degree of plate locking compared to regions where fracturing is less intense.

Our findings therefore provide new insights into the general characteristics of regional seismic hazard of the southern Central Andes. At the latitudes of the flat-slab segment, large megathrust earthquakes exhibit a recurrence interval of ~60–90 year⁶⁸, based on the frequency of the most recent historical events with M ≥ 8 (AD 1822, 1906, and 1985 in the Valparaiso sub-segment and 1880, 1943 and 2015 in the Illapel sub-segment). The segment north of 29°S, which last experienced a major rupture in 1922 represents a seismic gap, the so-called Atacama gap (Fig. 1a). Geodetic measurements furthermore indicate that the near-trench subduction interface at these latitudes is locked⁶⁹ (Figs. 1c and 4a). This raises the question as to how a seismogenic rupture will be produced along these portions of the subduction interface and if stresses will be released either in the form of either several large but not great earthquakes (M7-8) or during one great earthquake (M > 8.5). Regions that are subjected to high fluid pressures are generally weaker and more unstable under the stress build-up than faults under dry conditions^63,70. As such, weak faults are more prone to seismic rupture and favour stress release by several earthquakes of moderate magnitude⁶⁶. Conversely, regions with a drier interface can accumulate elastic energy for longer periods of time, which can then be released by less frequent, but very high-magnitude earthquakes. The locked region, the subduction interface displays low Vp/Vs, thereby indicating dry conditions (Fig. 4a, b). It thus follows that this domain is likely to be able to generate a great earthquake in the future, a conclusion that calls for increased geodetic monitoring of the seafloor to better assess the updip limit of locking and thus provide more accurate estimates of the size of a future earthquake.

Seismic catalogue

We derived the seismicity distribution for the study area between 1964 and 2018 from the reviewed bulletin of the International Seismological Centre^71,72,73,74, from which we removed all events with fixed depth or a depth error >5 km. We then filtered the events related to the oceanic slab (‘slab seismicity’), according to their proximity to the subduction interface, which is represented by the top of the slab from the Slab2 subduction zone geometry model⁷⁵. All events below the top of the slab, with a tolerance given by the uncertainty of the Slab2 model were extracted, a total of 9125 slab earthquakes, which comprises both plate interface events and intraslab events. To classify the events into ‘intra-slab’ and ‘inter-plate’ seismicity, we relied on the focal mechanism solutions from the GCMT catalogue^24,25. An earthquake was considered as inter-plate seismicity if it fulfilled four criteria⁷⁶: (i) it occurs in the interval defined by the depth of the top slab ± the uncertainty of the Slab2 model⁷⁵; (ii) the slip vector orientation represents thrusting; (iii) the strike of one of the nodal planes is equal to the strike of the top slab surface from Slab2⁷⁵ (±15°); and (iv) the dip angle of one of the nodal planes differs <15° from the dip angle of the top of the slab surface⁷⁵. The remaining earthquakes were considered as intra-slab seismicity. There were 873 events with GCMT mechanisms of which 114 events were classified as intra-slab events and 759 as inter-plate events.

P- and S- wave seismic tomography

We derived the Vp/Vs model from the isotropic P and S velocity models presented in Gao et al.²³. Here, we only provide a brief summary of the methodology followed to construct the tomography model; more technical details about the inversion workflow can be found in the original publication²³.

The seismic velocity model is the result of a multi-scale full waveform inversion⁶ based on the adjoint methodology^77,78, starting from the 3D seismic velocity model SP12RTS⁷⁷. Solutions of the visco-elastic wave equation in a radially anisotropic Earth model are obtained from Salvus^77,79.

Long-wavelength surface topography from the EGM2008 Geoid⁸⁰ and Earth2014 global topography model⁸¹ with Earth ellipticity according to WGS84 and the Moho depth from the global model CRUST1.0⁸² were implemented by deforming the mesh grid vertically. The surface topography and Moho depth have been filtered with a maximum angular order of 128, equivalent to a spatial resolution of 155 km. In order to constrain the deep structure of the upper mantle, we initialised our inversion from long-period surface wave data at 60–120 s and progressed in seven stages to a final period range of 12–120 s. The progressive extension to shorter periods mitigates the risk of falling into local minima. The model updates are driven by the Limited-Memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS)⁸³. We employ the Time-Frequency Phase Shift and Cross-Correlation-Coefficient misfits as misfit functions during stages 1–5 and 6–7, respectively, with the assistance of the Large-scale Seismic Inversion Framework 2.0^84,85.

1. Anikiev, D., Cacace, M., Bott, J., Gomez Dacal, M. L. & Scheck-Wenderoth, M. Influence of lithosphere rheology on seismicity in an intracontinental rift: the case of the Rhine Graben. Front. Earth Sci. 8, 492 (2020).

   Article  Google Scholar 

2. Gholamrezaie, E. et al. Lithospheric strength variations and seismotectonic segmentation below the Sea of Marmara. Tectonophysics 815, 228999 (2021).

   Article  Google Scholar 

3. Ibarra, F., Prezzi, C. B., Bott, J., Scheck‐Wenderoth, M. & Strecker, M. R. Distribution of temperature and strength in the central andean lithosphere and its relationship to seismicity and active deformation. J. Geophys. Res. Solid Earth 126, e2020JB021231 (2021).

4. Sibson, R. H. Fault zone models, heat flow, and the depth distribution of earthquakes in the continental crust of the United States. Bull. Seismol. Soc. Am. 72, 151–163 (1982).

   Google Scholar 

5. Scholz, C. H. The brittle-plastic transition and the depth of seismic faulting. Geol. Rundsch 77, 319–328 (1988).

   Article  Google Scholar 

6. Brace, W. F. & Kohlstedt, D. L. Limits on lithospheric stress imposed by laboratory experiments. J. Geophys. Res. Solid Earth 85, 6248–6252 (1980).

   Article  Google Scholar 

7. Tullis, J. & Yund, R. International Geophysics. Chapter 4. (eds. Evans, B. & Wong, T.) vol. 51 89–117 (Academic Press, 1992).

8. Kirby, S. Interslab earthquakes and phase changes in subducting lithosphere. Rev. Geophys. 33, 287–297 (1995).

   Article  Google Scholar 

9. Hacker, B. R., Peacock, S. M., Abers, G. A. & Holloway, S. D. Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions? J. Geophys. Res. Solid Earth 108, 1–26 (2003).

10. Peacock, S. M. Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle? Geology 29, 299–302 (2001).

    Article  CAS  Google Scholar 

11. Ferrand, T. P. et al. Dehydration-driven stress transfer triggers intermediate-depth earthquakes. Nat. Commun. 8, 15247 (2017).

    Article  CAS  Google Scholar 

12. Scambelluri, M. et al. Fossil intermediate-depth earthquakes in subducting slabs linked to differential stress release. Nat. Geosci. 10, 960–966 (2017).

    Article  CAS  Google Scholar 

13. John, T. et al. Generation of intermediate-depth earthquakes by self-localizing thermal runaway. Nat. Geosci. 2, 137–140 (2009).

    Article  CAS  Google Scholar 

14. Kelemen, P. B. & Hirth, G. A periodic shear-heating mechanism for intermediate-depth earthquakes in the mantle. Nature 446, 787–790 (2007).

    Article  CAS  Google Scholar 

15. Isacks, B. & Molnar, P. Distribution of stresses in the descending lithosphere from a global survey of focal-mechanism solutions of mantle earthquakes. Rev. Geophys. 9, 103 (1971).

    Article  Google Scholar 

16. Sandiford, D., Moresi, L. M., Sandiford, M., Farrington, R. & Yang, T. The fingerprints of flexure in slab seismicity. Tectonics 39, e2019TC005894 (2020).

17. Scholz, C. H. Earthquakes and friction laws. Nature 391, 37–42 (1998).

    Article  CAS  Google Scholar 

18. Costa, C. et al. Hazardous faults of South America; compilation and overview. J. South Am. Earth Sci. 104, 102837 (2020).

    Article  Google Scholar 

19. Anderson, M., Alvarado, P., Zandt, G. & Beck, S. Geometry and brittle deformation of the subducting Nazca Plate, Central Chile and Argentina. Geophys. J. Int. 171, 419–434 (2007).

    Article  Google Scholar 

20. Barazangi, M. & Isacks, B. L. Spatial distribution of earthquakes and subduction of the Nazca plate beneath South America. Geology 4, 686–692 (1976).

    Article  Google Scholar 

21. International Seismological Centre. On-line Bulletin. https://doi.org/10.31905/D808B830 (2021).

22. Rodriguez Piceda, C., Scheck-Wenderoth, M., Cacace, M., Bott, J. & Strecker, M. R. Long-term lithospheric strength and upper-plate seismicity in the Southern Central Andes, 29°–39°S. Geochem. Geophys. Geosyst. 23, 22 (2022).

    Article  Google Scholar 

23. Gao, Y. et al. Impact of the juan fernandez ridge on the pampean flat subduction inferred from full waveform inversion. Geophys. Res. Lett. 48, e2021GL095509 (2021).

    Article  Google Scholar 

24. Ekström, G., Nettles, M. & Dziewoński, A. M. The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. Phys. Earth Planet. Inter. 200–201, 1–9 (2012).

    Article  Google Scholar 

25. Dziewonski, A. M., Chou, T.-A. & Woodhouse, J. H. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. Solid Earth 86, 2825–2852 (1981).

    Article  Google Scholar 

26. O’Connell, R. J. & Budiansky, B. Seismic velocities in dry and saturated cracked solids. J. Geophys. Res. 79, 5412–5426 (1974).

    Article  Google Scholar 

27. Eberhart‐Phillips, D., Han, D.-H. & Zoback, M. D. Empirical relationships among seismic velocity, effective pressure, porosity, and clay content in sandstone. Geophysics 54, 82–89 (1989).

    Article  Google Scholar 

28. Hacker, B. R., Abers, G. A. & Peacock, S. M. Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H 2O contents: subduction zone mineralogy and physical properties. J. Geophys. Res. 108, 1–20 (2003).

29. Christensen, N. I. Serpentinites, peridotites, and seismology. Int. Geol. Rev. 46, 795–816 (2004).

    Article  Google Scholar 

30. Bloch, W. et al. Watching dehydration: seismic indication for transient fluid pathways in the oceanic mantle of the subducting nazca slab. Geochem. Geophys. Geosyst. 19, 3189–3207 (2018).

    Article  Google Scholar 

31. Takei, Y. Effect of pore geometry on VP/VS: from equilibrium geometry to crack. J. Geophys. Res. Solid Earth 107, 6–12 (2002).

    Article  Google Scholar 

32. Bezacier, L., Reynard, B., Bass, J. D., Sanchez-Valle, C. & Van de Moortèle, B. Elasticity of antigorite, seismic detection of serpentinites, and anisotropy in subduction zones. Earth Planet. Sci. Lett. 289, 198–208 (2010).

    Article  CAS  Google Scholar 

33. Christensen, N. I. Elasticity of ultrabasic rocks. J. Geophys. Res. 71, 5921–5931 (1966).

    Article  Google Scholar 

34. Christensen, N. I. The abundance of serpentinites in the oceanic crust. J. Geol. 80, 709–719 (1972).

    Article  CAS  Google Scholar 

35. Rüpke, L. H., Morgan, J. P., Hort, M. & Connolly, J. A. D. Serpentine and the subduction zone water cycle. Earth Planet. Sci. Lett. 223, 17–34 (2004).

    Article  Google Scholar 

36. Ulmer, P. & Trommsdorff, V. Serpentine stability to mantle depths and subduction-related magmatism. Science 268, 858–861 (1995).

    Article  CAS  Google Scholar 

37. Hirschmann, M. M. Mantle solidus: experimental constraints and the effects of peridotite composition. Geochem. Geophys. Geosyst. 1, 1–26 (2000).

38. Healy, D., Reddy, S. M., Timms, N. E., Gray, E. M. & Brovarone, A. V. Trench-parallel fast axes of seismic anisotropy due to fluid-filled cracks in subducting slabs. Earth Planet. Sci. Lett. 283, 75–86 (2009).

    Article  CAS  Google Scholar 

39. Reynard, B., Nakajima, J. & Kawakatsu, H. Earthquakes and plastic deformation of anhydrous slab mantle in double Wadati-Benioff zones. Geophys. Res. Lett. 37, 1–6 (2010).

40. Huang, Z., Tilmann, F., Comte, D. & Zhao, D. P wave azimuthal anisotropic tomography in Northern Chile: insight into deformation in the subduction zone. J. Geophys. Res. Solid Earth 124, 742–765 (2019).

    Article  Google Scholar 

41. Rodriguez Piceda, C. et al. Controls of the lithospheric thermal field of an ocean-continent subduction zone: the southern Central Andes. Earth ArXiv https://doi.org/10.31223/X5B05D (2021).

42. Marot, M. et al. Flat versus normal subduction zones: a comparison based on 3-D regional traveltime tomography and petrological modelling of central Chile and western Argentina (29°–35°S). Geophys. J. Int. 199, 1633–1654 (2014).

    Article  Google Scholar 

43. Anderson, M. L., Zandt, G., Triep, E., Fouch, M. & Beck, S. Anisotropy and mantle flow in the Chile-Argentina subduction zone from shear wave splitting analysis. Geophys. Res. Lett. 31, 1–4 (2004).

44. Boudier, F., Baronnet, A. & Mainprice, D. Serpentine mineral replacements of natural olivine and their seismic implications: oceanic lizardite versus subduction-related antigorite. J. Petrol. 51, 495–512 (2010).

    Article  CAS  Google Scholar 

45. Faccenda, M., Gerya, T. V., Mancktelow, N. S. & Moresi, L. Fluid flow during slab unbending and dehydration: implications for intermediate-depth seismicity, slab weakening and deep water recycling. Geochem. Geophys. Geosyst. 13, 1–23 (2012).

46. Kopp, H., Flueh, E. R., Papenberg, C. & Klaeschen, D. Seismic investigations of the O’Higgins Seamount Group and Juan Fernández Ridge: Aseismic ridge emplacement and lithosphere hydration. Tectonics 23, 1–21 (2004).

47. Grevemeyer, I., Ranero, C. R. & Ivandic, M. Structure of oceanic crust and serpentinization at subduction trenches. Geosphere 14, 395–418 (2018).

    Article  Google Scholar 

48. Gutscher, M., Spakman, W., Bijwaard, H. & Engdahl, E. R. Geodynamics of flat subduction: seismicity and tomographic constraints from the Andean margin. Tectonics 19, 814–833 (2000).

    Article  Google Scholar 

49. Kay, S. M. & Mpodozis, C. Magmatism as a probe to the Neogene shallowing of the Nazca plate beneath the modern Chilean ¯at-slab. J. South Am. Earth Sci. 19, 39–57 (2002).

50. Incel, S. et al. Laboratory earthquakes triggered during eclogitization of lawsonite-bearing blueschist. Earth Planet. Sci. Lett. 459, 320–331 (2017).

    Article  CAS  Google Scholar 

51. Wang, H., Huismans, R. S. & Rondenay, S. Water migration in the subduction mantle wedge: a two-phase flow approach. J. Geophys. Res. Solid Earth 124, 9208–9225 (2019).

    Article  Google Scholar 

52. Halpaap, F. et al. Earthquakes track subduction fluids from slab source to mantle wedge sink. Sci. Adv. https://doi.org/10.1126/sciadv.aav7369 (2019).

53. Jung, H., Green, H. W. II & Dobrzhinetskaya, L. F. Intermediate-depth earthquake faulting by dehydration embrittlement with negative volume change. Nature 428, 545–549 (2004).

    Article  CAS  Google Scholar 

54. Raleigh, C. B. & Paterson, M. S. Experimental deformation of serpentinite and its tectonic implications. J. Geophys. Res. 70, 3965–3985 (1965).

    Article  Google Scholar 

55. Guo, H., McGuire, J. J. & Zhang, H. Correlation of porosity variations and rheological transitions on the southern Cascadia megathrust. Nat. Geosci. 14, 341–348 (2021).

    Article  CAS  Google Scholar 

56. Gao, H. Three-dimensional variations of the slab geometry correlate with earthquake distributions at the Cascadia subduction system. Nat. Commun. 9, 1204 (2018).

    Article  Google Scholar 

57. Mochizuki, K., Henrys, S., Haijima, D., Warren-Smith, E. & Fry, B. Seismicity and velocity structure in the vicinity of repeating slow slip earthquakes, northern Hikurangi subduction zone, New Zealand. Earth Planet. Sci. Lett. 563, 116887 (2021).

    Article  CAS  Google Scholar 

58. Fujie, G. et al. Controlling factor of incoming plate hydration at the north-western Pacific margin. Nat. Commun. 9, 3844 (2018).

    Article  Google Scholar 

59. Wagner, L. S., Caddick, M. J., Kumar, A., Beck, S. L. & Long, M. D. Effects of oceanic crustal thickness on intermediate depth seismicity. Front. Earth Sci. 8, 244 (2020).

    Article  Google Scholar 

60. Porter, R. et al. Shear wave velocities in the Pampean flat-slab region from Rayleigh wave tomography: implications for slab and upper mantle hydration. J. Geophys. Res. 117, 1–21 (2012).

61. Engdahl, E. R. & Scholz, C. H. A double Benioff Zone beneath the central Aleutians: An unbending of the lithosphere. Geophys. Res. Lett. 4, 473–476 (1977).

    Article  Google Scholar 

62. Wang, K. Unbending combined with dehydration embrittlement as a cause for double and triple seismic zones. Geophys. Res. Lett. 29, 36–1 (2002).

    Article  Google Scholar 

63. Saffer, D. M. & Tobin, H. J. Hydrogeology and mechanics of subduction zone forearcs: fluid flow and pore pressure. Annu. Rev. Earth Planet. Sci. 39, 157–186 (2011).

    Article  CAS  Google Scholar 

64. Bekins, B. A. & Dreiss, S. J. A simplified analysis of parameters controlling dewatering in accretionary prisms. Earth Planet. Sci. Lett. 109, 275–287 (1992).

    Article  Google Scholar 

65. Moreno, M. et al. Locking of the Chile subduction zone controlled by fluid pressure before the 2010 earthquake. Nat. Geosci 7, 292–296 (2014).

    Article  CAS  Google Scholar 

66. Moreno, M. et al. Chilean megathrust earthquake recurrence linked to frictional contrast at depth. Nat.Geosci. 11, 285–290 (2018).

    Article  CAS  Google Scholar 

67. Métois, M., Vigny, C. & Socquet, A. Interseismic coupling, megathrust earthquakes and seismic swarms along the chilean subduction zone (38°–18°S). Pure Appl. Geophys. 173, 1431–1449 (2016).

    Article  Google Scholar 

68. Dura, T. et al. Coastal evidence for Holocene subduction-zone earthquakes and tsunamis in central Chile. Quat. Sci. Rev. 113, 93–111 (2015).

    Article  Google Scholar 

69. Moreno, M., Rosenau, M. & Oncken, O. 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone. Nature 467, 198–202 (2010).

    Article  CAS  Google Scholar 

70. Hubbert, M. K. & Rubey, W. W. Role of fluid pressure in mechanics of overthrust faulting: I. Mechanics of fluid-filled porous solids and its application to overthrust faulting. GSA Bull. 70, 115–166 (1959).

    Article  Google Scholar 

71. International Seismological Centre. http://www.isc.ac.uk/.

72. National Earthquake Information Center. https://www.usgs.gov/programs/earthquake-hazards/national-earthquake-information-center-neic (2021).

73. Instituto Nacional de Prevención Sísmica. https://www.inpres.gob.ar/desktop/ (2021).

74. Centro Sismológico Nacional. https://www.sismologia.cl/ (2021).

75. Hayes, G. P. et al. Slab2, a comprehensive subduction zone geometry model. Science 362, 58 (2018).

    Article  CAS  Google Scholar 

76. Tichelaar, B. W. & Ruff, L. J. Depth of seismic coupling along subduction zones. J. Geophys. Res. Solid Earth 98, 2017–2037 (1993).

    Article  Google Scholar 

77. Fichtner, A., Kennett, B. L. N., Igel, H. & Bunge, H.-P. Full waveform tomography for radially anisotropic structure: New insights into present and past states of the Australasian upper mantle. Earth Planet. Sci. Lett. 290, 270–280 (2010).

    Article  CAS  Google Scholar 

78. Tape, C., Liu, Q., Maggi, A. & Tromp, J. Seismic tomography of the southern California crust based on spectral-element and adjoint methods. Geophys. J. Int. 180, 433–462 (2010).

    Article  Google Scholar 

79. Afanasiev, M. et al. Modular and flexible spectral-element waveform modelling in two and three dimensions. Geophys. J. Int. 216, 1675–1692 (2019).

    Article  Google Scholar 

80. Pavlis, N. K., Holmes, S. A., Kenyon, S. C. & Factor, J. K. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). J. Geophys. Res. Solid Earth 117, 1–38 (2012).

81. Hirt, C. & Rexer, M. Earth2014: 1 arc-min shape, topography, bedrock and ice-sheet models – available as gridded data and degree-10,800 spherical harmonics. Int. J. Appl. Earth Obs. Geoinf. 39, 103–112 (2015).

    Google Scholar 

82. Laske, G., Masters, G., Ma, Z. & Pasyanos, M. Geophysical Research Abstracts. Vol. 15. p. 2658 (EGU General Assembly Vienna, 2013).

83. Liu, D. C. & Nocedal, J. On the limited memory BFGS method for large scale optimization. Math. Program. 45, 503–528 (1989).

    Article  Google Scholar 

84. Krischer, L. et al. ObsPy: a bridge for seismology into the scientific Python ecosystem. Comput. Sci. Discov. 8, 014003 (2015).

    Article  Google Scholar 

85. Thrastarson, S., van Herwaarden, D.-P., Krischer, L. & Fichtner, A. Large-scale seismic inversion framework, an updated version. Earth ArXiv. https://doi.org/10.31223/X5NC84 (2021).

86. QGIS Development Team. QGIS Geographic Information System (Open Source Geospatial Foundation, 2009).

87. Wessel, P. et al. The generic mapping tools version 6. Geochem. Geophys. Geosyst. 20, 5556–5564 (2019).

    Article  Google Scholar 

88. Henderson, A. ParaView Guide, A Parallel Visualization Application (2007).

89. Amante, C. & Eakins, B. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. https://doi.org/10.7289/V5C8276M (2009).

90. Yáñez, G. A., Ranero, C. R., von Huene, R. & Díaz, J. Magnetic anomaly interpretation across the southern central Andes (32–34 S): The role of the Juan Fernández Ridge in the late Tertiary evolution of the margin. J. Geophys. Res. Solid Earth 106, 6325–6345 (2001).

    Article  Google Scholar 

91. Smithsonian Institution. Global GIS: Volcanoes Of The World; Volcano Basic Data (American Geological Institute, 2003).

92. Drápela, J., Calisto, I. & Moreno, M. Locking-derived tsunami scenarios for the most recent megathrust earthquakes in Chile: implications for tsunami hazard assessment. Nat. Hazards 107, 35–52 (2021).

    Article  Google Scholar 

93. Rodriguez Piceda, C. et al. Lithospheric density structure of the southern Central Andes constrained by 3D data-integrative gravity modelling. Int. J. Earth Sci. https://doi.org/10.1007/s00531-020-01962-1 (2020).

Download references

Authors and Affiliations

1. HelmholtzZentrum GFZ – German Research Centre for Geosciences, Potsdam, Germany

   Constanza Rodriguez Piceda, Ya-Jian Gao, Mauro Cacace, Magdalena Scheck-Wenderoth, Judith Bott & Frederik Tilmann

2. University of Potsdam, Institute of Geosciences, Potsdam, Germany

   Constanza Rodriguez Piceda & Manfred Strecker

3. University of Plymouth, School of Geography, Earth and Environmental Sciences, Plymouth, UK

   Constanza Rodriguez Piceda

4. Freie Universität Berlin, Berlin, Germany

   Ya-Jian Gao & Frederik Tilmann

5. Karlsruhe Institute of Technology, Geophysical Institute (GPI), Karlsruhe, Germany

   Ya-Jian Gao

6. RWTH Aachen University, Aachen, Germany

   Magdalena Scheck-Wenderoth

Contributions

C.R.P. conceived and designed the main conceptual ideas of this work, processed the seismic catalogue, provided the 3D lithospheric model, performed the analysis and wrote the main parts of the manuscript. Y.G. conducted the seismic tomography imaging, and contributed to the analysis of the results and the writing of the manuscript. M.C. helped devising the main conceptual ideas and proof outline of this work, and contributed to the analysis of the results and the writing of the manuscript. J.B. and F.T. contributed to the analysis of the results and the writing of the manuscript. M.S.W. and M.S. supervised the study and contributed to the interpretation of results and the writing of the paper.

Corresponding author

Correspondence to Constanza Rodriguez Piceda.

Competing interests

The authors declare no competing interests.

Peer review information

Communications Earth & Environment thanks Diana Comte, Wasja Bloch and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Luca Dal Zilio and Joe Aslin. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions