Research ArticleGEOPHYSICS

Semibrittle seismic deformation in high-temperature mantle mylonite shear zone along the Romanche transform fault

See allHide authors and affiliations

Science Advances  09 Apr 2021:
Vol. 7, no. 15, eabf3388
DOI: 10.1126/sciadv.abf3388


Oceanic transform faults, a key element of plate tectonics, represent the first-order discontinuities along mid-ocean ridges, host large earthquakes, and induce extreme thermal gradients in lithosphere. However, the thermal structure along transform faults and its effects on earthquake generation are poorly understood. Here we report the presence of a 10- to 15-kilometer-thick in-depth band of microseismicity in 10 to 34 kilometer depth range associated with a high-temperature (700° to 900°C) mantle below the brittle lithosphere along the Romanche mega transform fault in the equatorial Atlantic Ocean. The occurrence of the shallow 2016 moment magnitude 7.1 supershear rupture earthquake and these deep microearthquakes indicate that although large earthquakes occur in the upper brittle lithosphere, a substantial amount of deformation is accommodated in the semibrittle mylonitic mantle that resides at depths below the 600°C isotherm. We also observe a rapid westward deepening of this band of seismicity indicating a strong lateral heterogeneity.


Mid-ocean ridges (MORs) are segmented by seismogenically active oceanic transform faults (OTFs) (Fig. 1A), first-order discontinuities offsetting ocean spreading centers. As the plates move away from the MOR axis, the oceanic lithosphere cools, subsides, and thickens (1). Consequently, a large-offset OTF separates a cold lithosphere on one side from a hot lithosphere on the other side, producing a complex thermal regime. Although many thermal models have been proposed for OTFs (25), the deep thermal structure along very long-offset OTFs (for example, the Romanche OTF in the equatorial Atlantic) is still not constrained by any geophysical observations.

Fig. 1 Study area.

(A) Location of the study region in the equatorial Atlantic Ocean. Red beach balls indicate three large earthquakes with magnitudes of ≥7.0 (56). The white [black in (B)] lines and numbers indicate the age of the oceanic lithosphere in million years (25). (B) Bathymetric map of the MAR and the eastern Romanche OTF. The white diamonds show the dredged peridotites (17), and the red ones are mylonitic peridotites obtained during the 2019 SMARTIES cruise (37, 45). The triangles indicate the locations of the deployed ocean-bottom seismographs (OBSs). The black triangles indicate stations that generated seismic data; the blue triangles indicate seismographs that were deployed, but did not generate seismic data. Red star and beach ball, the 2016 Mw 7.1 Romanche earthquake. Two red dashed lines mark the transform valley. The white dashed line indicates the seismic refraction profile (28).

The strike-slip motion along OTFs generates seismicity, including large earthquakes (Fig. 1A). Previous studies have indicated that temperature is the primary factor controlling the depth of the brittle to ductile transition along OTFs, approximately defined by the 600°C isotherm (68). As a consequence, the maximum focal depth of seismicity deepens with increasing distance along the transform away from the ridge-transform intersection (RTI) following the 600°C isotherm. However, some thermal models predict shallowing of the 600°C isotherm toward the center of the transform fault (4, 5, 9). Laboratory deformation experiments on olivine indicate that the 600°C isotherm marks the approximate brittle-ductile boundary in the oceanic lithosphere (10, 11). However, recent petrological studies of peridotite mylonites from the Shaka and Prince Edward OTFs, which offset the ultraslow spreading Southwest Indian Ridge, suggest that the brittle-ductile boundary may extend downward to higher temperatures up to 850° to 875°C (12, 13). Peridotite mylonites are usually formed by ductile deformation in a high shear strain regime and at high pressure-temperature conditions (11, 14, 15), giving rise to a semibrittle deformation behavior in mylonite shear zones (12, 13). Their exposures in the ocean may be attributed to the mantle uplift by transpressional motion, for example, at the Romanche and St. Paul OTFs on the Mid-Atlantic Ridge (MAR) (16, 17). Meanwhile, material/physical properties of fault segments (including porosity, mineral alteration, permeability, friction stability, etc.) could also influence seismicity and rupture propagations on OTFs (1821). It is important to note that nearly 85% of plate motion on OTFs is accommodated aseismically (22) or during microearthquakes and earthquake swarms in the mantle (23, 24). Using microseismicity data from the Romanche OTF close to the eastern RTI, here we elucidate thermal states, seismic slips, and rupture propagations on megatransform faults.


Microearthquake experiment and data

Our study region lies along the eastern part of the Romanche OTF in the equatorial Atlantic Ocean, which offsets the slow-spreading MAR by ~880 km, with a maximum age contrast of ~45 million years (Ma) old (Fig. 1) (25). In this area, the transform valley seafloor is 10- to 20-km wide, up to 7300-m deep and the northern bounding fault wall rises up to 2000 m below the sea level (Fig. 1B). Thermal modeling indicates that the horizontal thermal gradient across the Romanche OTF could exceed >30°C/km at 20 km depth (26) beneath the RTI. Such a strong across-strike thermal gradient would result in a marked reduction in crustal thickness near the RTI and thickening of the brittle lithosphere beneath the OTF (17). Previous studies suggested that there are a large number of peridotites south of the OTF and the absence of basalts from dredges in this area (Fig. 1B) (17, 26), indicating a thin crust here.

In July and August 2019, we deployed 19 four-component short-period ocean-bottom seismographs (OBSs) in the study region, with an average instrument spacing of ~30 km (Fig. 1B). The average continuous recording period for 17 OBSs is ~21 days but ~29 days for OBS12. Nine OBSs are located near the Romanche transform valley, covering a distance of ~100 km along its strike (Fig. 1B); five OBSs were deployed around the RTI and MAR (Fig. 1B); three OBSs were located on the northern suspended transform valley and transverse ridge (Fig. 1B) (17).

Both P- and S-wave arrivals picked from the seismometer records were used to locate earthquakes using a probabilistic, nonlinear earthquake location algorithm (27), and station corrections were calculated and iteratively updated to suppress the impact of the inhomogeneous three-dimensional (3D) structures (figs. S1 and S2) (see Materials and Methods). Travel times were calculated using a tested 1D velocity model (figs. S3 and S4 and table S1) constrained by a seismic refraction profile across the Romanche OTF (Fig. 1B) (28). We carried out extensive resolution tests by changing the velocity models (figs. S3 and S4) and initial depths (fig. S5), which show that microseismicity depths are robust (see Materials and Methods and figs. S3 to S5). In this study, we removed S-wave delays caused by low-velocity unconsolidated sediments (fig. S6) (29). The S-wave delays range from ~0.1 s on OBS03 to ~1.3 s on OBS06 (table S2) in the valley, and the hypocentral parameters were remarkably improved after removing these delays (figs. S7 and S8). We determined 209 earthquakes using at least six P- and S-wave arrivals, with depth uncertainties of ~1.8 km (fig. S4F) and local magnitudes of ML = 0.5 to 3.5 (figs. S9 and S10). In total, 155 of these events were also relocated using the double-difference location method (Fig. 2) (30), with a station gap of <270°. Both the average horizontal and depth uncertainties of these well-located events are ~1.6 km (Fig. 2A). Focal mechanism solutions for three earthquake swarms were obtained (Fig. 3), which show normal faulting (fig. S11) instead of strike-slip faulting expected along the transform fault (6). The focal mechanism solutions were not sensitive to changes in the source depth (fig. S12) and the absolute velocity values (table S3), suggesting that the focal mechanism solutions are robust.

Fig. 2 Seismicity along the Romanche transform fault.

(A) Epicenters of microearthquakes on the transform fault, with a major-semi axis of <10 km, are shown by magenta circles (black ones have station gaps of >270°). The trace of active transform fault is shown by white dashed lines. Blue beach balls show three calculated focal mechanism solutions shown in Fig. 3, and the red one shows the 2016 Mw 7.1 Romanche earthquake (33). The average horizontal uncertainty (~1.6 km) is marked by the black plus sign. (B) Depth of the earthquakes along-strike line A′A below the seafloor. Zero corresponds to the RTI location. The longitude is shown on the top for reference. Isotherms are calculated using a plate cooling model (53, 54) and averaging both sides of the transforms (see Materials and Methods). The solid green line indicates the 800°C isotherm extracted from the thermal model by Ligi et al. (3). Solid color rectangles show eight earthquake swarms in Fig. 3. The magenta line denotes the base of the seismogenic zone. Three shaded areas with numbers 1 to 3 indicate three distinct zones discussed in the text. The gray histograms indicate the numbers of earthquakes along the profile. The red arrows indicate the rupture direction and extent of the 2016 Mw 7.1 earthquake (33). (C) The seafloor depths along nine transects across the OTF [red lines in (A)] and their distances in the north (positive) and south (negative) of the transform trace line in (A). (D) Earthquake depths across the OTF. The red dashed lines indicate the inferred fault planes (~80°). The beach balls indicate focal mechanism solutions in the side-on projection view.

Fig. 3 Earthquake temporal distribution, earthquake swarms, and focal mechanisms.

(A) Focal depth distribution of 209 located earthquakes as a function of dates from 20 July to 16 August 2019. The red line denotes the cumulative event number. Magnitude scales are shown on the upper right. Green histograms on the left show the depth distribution. Eight colored rectangles with numbers denote earthquake swarms shown in Fig. 2. They repeatedly occurred in Zones 1 and 2, with shifting focal depths (<5 km) and dates. (B and C) Waveforms of two stations (OBS03 and OBS04) on the vertical components are aligned on the P-arrivals for eight earthquake swarms. The numbers of events in each earthquake swarm are shown by brown histograms between (B) and (C). The blue and red arrows denote the P and S arrivals, respectively. (D to F) Focal mechanism solutions for three earthquake swarms (nos. 3, 5, and 7). P-wave polarities are shown in black (top) and white dots (bottom).


In this study, most microearthquakes occurred at depths of ~10 to 34 km below the seafloor (bsf), and only a few shallow events with depths <10 km are observed near the RTI (Fig. 2). All events are located beneath the southern portion of the transform valley, within a 10-km-wide belt toward the warmer lithosphere (Fig. 2A); no events are observed beneath the central or northern portions of the valley (Fig. 2A). Their focal depths rapidly increase westward from ~10 km at the RTI to 34 ± 2 km at 70 ± 5 km distance along the transform fault (Fig. 2B). To the east (20 to 60 km), these events are focused beneath a bathymetric high in the transform valley (Figs. 1B and 2A) in a 10 to 20 km depth range, whereas in the west (70 to 90 km) beneath a deep valley at an 18 to 34 km depth range, with an abrupt transition in between. Peridotites have been sampled from the seafloor in a ~50-km-wide zone south of the Romanche OTF close to the RTI (Fig. 1B) (31, 32), indicating that the bathymetric high here might be a peridotite hill formed near the ridge axis. Microearthquakes align along two steep (~80°) south-dipping planes extending down to 34 km (Fig. 2D), which are consistent with the dip inferred for the 1994 [moment magnitude (Mw = 7.0)] and 2016 (Mw = 7.1) earthquakes (6, 33), but extend much deeper than their hypocentral depths (~17 km bsf) (33), indicating that the deformation extends to much greater depths. However, three focal mechanism solutions show normal faulting (Figs. 2 and 3, D to F), implying extensional stress here. Their dips vary from 40° to 50°, which is common for normal earthquakes but much lower than that for the 2016 Mw 7.1 strike-slip earthquake (Fig. 2D) (33). These earthquake swarms are all located beneath the bathymetric highs (Fig. 2A), indicating that they might be linked with the uplift. The presence of a wide and deep transform valley supports the existence of normal faulting.


On the basis of the seismicity pattern, the study area can be divided into three zones, Zones 1, 2, and 3 (Fig. 2). These results highlight four key observations: (i) The absence of seismicity in the uppermost 5 to 18 km from east to west; (ii) the seismicity lies in a band, 11 to 18 km bsf in Zone 1 and 18 to 34 km in Zone 3, defining the seismogenic zone; (iii) the seismogenic zone in Zone 1 rapidly deepens more than 7 km in a narrow (25 km) transition Zone 2 to 18 to 34 km depth in Zone 3 (Fig. 2B); (iv) about 60% of the seismicity and most of the swarms occurred in the transition Zone 2.

Microseismicity and thermal structure

Our results indicate that the microseismicity occurred within three distinct zones: (i) Eastern Zone 1 between 0 and 30 km distance range, where the seismicity lies at ~11 to 18 km depth (bsf) roughly bounded by the 400° and 700°C isotherms, (ii) Transition Zone 2 between 40 to 65 km distance, constrained by the 500° and 800°C isotherms with the upper bound of seismicity rapidly deepening from 12 km to 18 km, and (iii) Western Zone 3 between 65 and 90 km distance with seismicity at 18 to 34 km depths limited by the isotherms of 520° and 900°C (Fig. 4). The temperature limit for earthquakes, defined by the maximum depth of earthquakes from global seismic observations, remains a matter of debate, ranging between 600° and 800°C (2, 6, 29, 34, 35). However, our results indicate that the maximum depths of microseismicity along the OTF are not consistent with a fixed temperature isotherm but vary remarkably along the OTF within 90 km of the RTI.

Fig. 4 Schematic diagram showing the thermal structure along the eastern Romanche transform fault.

The frontal plane is dipping ~80° southward, corresponding to the fault plane of the 2016 Mw 7.1 earthquake (red star) and these are shown in Fig. 2D. The maximum hypocentral depth, the base of the semibrittle zone, is bounded by the 700° to 900°C isotherms (magenta line). The upper boundary of the seismogenic zone is shown in a blue dashed line. The rupture propagations of the 2016 Mw 7.1 earthquake (33) are shown in gray with black arrows. Mantle mylonite shear zone is shown in brown (see the main text).

A similar observation has been noted for the Gofar (18), Discovery (36), and Blanco (24) OTFs between rupture patches and rupture barriers. There could be three possibilities to explain such a variation. (i) The classical thermal models used here (Fig. 2B) do not portray the true thermal structure of the lithosphere. More sophisticated numerical models (25, 9) indicate a somewhat deeper 600°C isotherm away from the RTI. When the viscoplastic rheology and brittle weakening of the lithosphere are incorporated, the calculated 600°C isotherm would reach 22 km (9), which is still not deep enough for this deep microseismicity. On the other hand, some thermal models (4, 5, 9) show the shallowing of the 600°C isotherm toward the center of OTFs. Figure 2B clearly shows that even when complicated parameters are incorporated in the thermal models (3), these deep microearthquakes still occur at depths below the 800°C isotherm. However, it is possible that our short study region, which is just ~100 km long close to the vicinity of the RTI, may not be able to present the shallowing microseismicity corresponding to the elevated temperatures (4, 5, 9). (ii) Because of the previous thermal control on the maximum depth of earthquakes is based on the global catalog of large earthquakes (6, 7), it is plausible that the depth limit for large seismic ruptures is physically different from that for microseismicity and the associated transient slip. (iii) Furthermore, it is well known that the microseismicity is strongly influenced by local dynamic processes, e.g., material properties of fault segment (18), dilatancy strengthening (21), hydrothermal circulation (4), and brittle-fracturing within the mylonite zone (13), and these processes could introduce strong heterogeneities along the OTFs.

The absence of shallow microseismicity

In Zones 1 and 2, the absence of seismicity in the uppermost 5 to 16 km (Fig. 2B) indicates that the lithosphere is either completely locked or slipping aseismically at these depth ranges. The absence of seismicity here is located beneath the two bathymetric highs (~200 to 300 m in height) on the seafloor (Figs. 1, 2C, and 4). Peridotites were sampled here during the SMARTIES cruise (37) and dredges (17, 26) south of the transform valley (Fig. 1B), which suggest that the upper parts of Zones 1 and 2 might consist of serpentinized peridotites. However, a recent seismic study just west of Zone 2 indicates the presence of a 6-km-thick normal oceanic crust (Fig. 1B) (28), suggesting the existence of magmatic crust. It is difficult to say the lateral extent of this magmatic crust, but the few events observed in Zone 1 at ~6 km depth might be at the crust-mantle boundary (Fig. 4). Altered gabbro under supercritical water conditions would change from velocity-weakening to velocity-strengthening behavior (38) and may not be able to host microearthquakes. The mantle underneath is likely to be serpentinized peridotite, which results in velocity-strengthening behavior that inhibits earthquake nucleation (19), possibly linked to hydrothermal circulation (Fig. 4).

In Zone 3, the absence of seismicity down to 18 km depth coincides with the occurrence of the 2016 Mw 7.1 earthquake that initiated at 17 km depth bsf (33), propagated eastward up to 10 km depth at the western boundary of Zone 2, and then reversed direction and propagated westward at a supershear speed (Fig. 4). The presence of gabbroic crust as revealed by seismic refraction study might have enabled the supershear rupture during the 2016 earthquake. The absence of seismicity could be explained by the release of stress during the 2016 earthquake, and the region is now in the interseismic stage of the seismic cycle following the 2016 Mw 7.1 earthquake.

Model of deep microearthquakes

The eastern Zone 1, close to the RTI, has a subhorizontal band of seismicity in an 11 to 18 km depth range (Figs. 2 and 4) that extends to 700°C isotherm. This depth extend and the flat base of the seismogenic zone could be influenced by two different processes: (i) heat diffusion from the young (0 to 4 Ma) extremely hot lithosphere to the 45-Ma-old thick and cold lithosphere (Figs. 1 and 4) (26, 31, 39) or (ii) efficient hydrothermal circulation that will steepen the isotherm (40) downward near the ridge axis over a short distance, and consequently, the brittle lithosphere will thicken and flatten (40) in Zone 1 (Figs. 2 and 4). A recent seismic study from the equatorial Atlantic Ocean suggests that the lithosphere younger than 4 Ma is influenced by hydrothermal circulation (41), and therefore, we interpret that the band of seismicity between 11 and 18 km depth in Zone 1, which lies in the 400° to 700°C temperature range, is being affected by active hydrothermal circulation and subsequent serpentinization and mylonitization (Fig. 4).

In Zone 2, the steeply dipping seismicity band contains the largest number of earthquakes and swarms (Fig. 2B). The earthquake swarms in Zone 2 have similar epicenter locations but different depths (Figs. 2A and 3A), indicating that they may be produced by brittle failure through the creep of the aseismic serpentinized mantle activating frictionally velocity-weakening asperities as observed on the Blanco OTF (22, 24). However, in this study, these microearthquake swarms mostly occurred around the 800°C isotherm instead of the speculated 600°C isotherm (24), indicating higher temperature and deeper deformation. Zone 2 might have acted as a rupture barrier for the 2016 Mw 7.1 earthquake propagation (33). Both the potential extension indicated by the normal earthquake focal mechanisms (Figs. 2 and 3) and the wide transform valley (Fig. 2B) would support the dilatancy strengthening in Zone 2, resulting in a barrier for rupture propagation (21). Similar barriers have also been observed along the Discovery and Gofar OTFs (18, 19, 36) and have been interpreted as evidence for the presence of an active damage zone. We suggest that the earthquake swarms in Zone 2 are induced by active faulting, resulting in a large amount of fractures, similar to damage zones observed on continental strike-slip faults (42, 43), allowing the percolation of enhanced hydrothermal fluids to a greater depth (Fig. 4) (4, 12, 13). The seawater interacting with mantle peridotite would foster the mantle serpentinization and the active faulting would extend it deeper.

Similarly, in Zone 3, the deepest band of seismicity lies underneath the 2016 earthquake rupture zone at 18 to 34 km depth (Fig. 2B), bounded by the 520° and 900°C isotherms (Fig. 4), far exceeding the 600°C isotherm as the suggested boundary for the brittle-ductile transition (6, 10). The downward propagation of fracturing into the high-temperature (HT) mantle mylonitic shear zone below the 2016 earthquake rupture zone (12, 13) could explain the band of microseismicity (Fig. 4). Large amounts of peridotite mylonites have been collected on the OTF and ridge systems (11). Recent studies found that semibrittle deformation did occur at HT (> 850° to 875°C) mylonitic shear bands at the Shaka and Prince Edward OTFs (13), leading to the formation of fractures overprinting the ductile deformation. It is possible that the mantle at 18 to 34 km depth range may not only deform ductilely but have brittle fractures in the HT peridotite mylonites. This means that the brittle-ductile transition can extend into mantle mylonite at a HT of ~900°C due to high strain rates within the shear bands (13) rather than the currently assumed 600°C (11). The presence of mylonitic peridotites (Fig. 1B) (37) supports this possibility. Therefore, we propose that Zone 3 has an 18-km-thick seismically coupled lithosphere that supports the nucleation and rupture of large earthquakes, underlain by large amounts of mylonites with brittle fractures that contribute to deep microearthquakes (Fig. 4), creating a semibrittle zone.

This model can also support the deep microseismicity above the 800°C isotherm in Zone 2, which is influenced by both the downward extend of serpentinization of peridotites and HT fracturing in the semiductile mylonites (Fig. 4). Our observations of earthquakes in the deep mantle associated with the mylonite shear zone may be prevalent on other OTFs. Deep microseismicity has also been reported for the Blanco (24), Chain (44), Discovery (36), and Gofar OTFs (18) and has been interpreted to be caused by enhanced cooling, serpentinization, or hydrothermal circulations (36) but may be caused by the semibrittle deformation within the mantle mylonitic shear zone as suggested here.

We also found that the microearthquakes are all located on the younger side of the transform valley on a southward-dipping fault plane (Fig. 2, A and D), indicating that earthquake slip and rupture propagation along the long-offset OTF are more influenced by the young (hot) plate rather than the old (cold) plate (Fig. 4). It is possible that the colder part of the OTF in the north is either highly fractured or hydrothermal alterations have filled the fractures, leading to velocity strengthening, and microseismicity occurs at the transition. Together, these findings suggest that earthquake slip and rupture processes on an OTF can be affected by the combined influence of temperature, fluid circulation, and lithology of the lithosphere.

It is important to note that our relatively short time period dataset, compared to the time scale of a seismic cycle, can only provide a snapshot of the microseismicity here. Over a longer period, the long-term strain release dynamics might be significantly different. However, the link between the band of deep microseismicity, mantle shear deformation, and HT mylonitization might extend all along the OTF, which could only be addressed by future longer-period transform-scale studies.


Data acquisition and arrival detection

During the SMARTIES cruise (45) in July and August 2019, an array of 19 OBSs was deployed at the intersection of the eastern part of the Romanche Transform Fault and the MAR in the equatorial Atlantic Ocean. Eighteen OBSs were successfully recovered, one of which recorded no data because the logger tube was flooded (Fig. 1B). The OBS positions on the seafloor were determined by a triangulation exercise (46) during the deployment, and the linear clock drifts were applied and removed after the OBS recovery.

Initial P- and S-wave arrivals were detected automatically in the continuous time series using a short-term-average/long-term-average trigger algorithm in the SEISAN package (47), and arrival times were refined with a kurtosis-based picking tool written in MATLAB (48). The arrival times of events were then checked manually. An event was considered to be an earthquake when five or more stations were triggered coincidently. The estimated picking errors were assigned at 0.05 to 0.1 s and 0.1 to 0.2 s for P- and S-wave arrivals, respectively.

Earthquake location

For the initial location, five travel time arrivals, including one S-arrival, are required. We used the nonlinear oct-tree search algorithm: NonLinLoc program (27). The maximum likelihood location is chosen as the preferred hypocenter for each event. NonLinLoc estimates a 3D error ellipsoid (68% confidence) from the posterior density function scatter samples (27). Iterative calculations were used to find the best solution when the average RMS misfit yields a minimum, which was ~0.08 s after nine iterations in this study (fig. S1). At each iteration, the average residuals for P- and S-phase at each station were calculated, and then they were used as the station corrections for the later run of the location algorithm (27). The cumulative delay times were given by the sum of the average residuals calculated and the input station corrections at each iteration (fig. S2).

Reference velocity model and focal depth resolution

Our earthquake locations were made using a 1D P-wave velocity model, which was constructed through a large number of tests, using various crustal thicknesses as well as velocities in the crust and mantle (fig. S3). The crustal thicknesses ranged from 3 to 7 km. The velocity of the top 1 km and the depth (50 km) of the lower end of the velocity profile were kept fixed (fig. S3). The P-wave velocity structure in these models was derived from an active-source wide-angle seismic refraction profile (Fig. 1B) (28). We constructed a subset of 85 earthquakes that occurred on at least 10 stations, with RMS of <0.4 s, which subsequently was repeatedly located using various velocity models. We then selected the best one (table S1) exhibiting the best possible combination of a large number of located events, a low average RMS residual, a small average semimajor axis, and a large number of phases used in localization, which was then used to locate the hypocenters by the NonLinLoc (27), hypoDD (30) programs, as well as to calculate mechanism solutions. The selected velocity model is continuous, without the Moho discontinuity. The statistic table S1 shows that the continuous models seem to result in a better result but not obvious. Also, fewer PmP phases were identified beneath the transform valley from the seismic refraction profile (28), all probably indicating a weak Moho interface here.

For the full dataset, the selected velocity model also leads to more well-located events, fewer depth uncertainties, and lower RMS residuals than the fastest and slowest ones (fig. S4, C to E). However, when using the end-member models, absolute focal depths of most events changed between 20 and 30 km, which is significant (fig. S4, C to E). Notice that the fastest and slowest models are not very reasonable due to the too-thin crust (3 km) or too-slow velocity at the base of the lower crust (5.7 km/s) (fig. S4A). The shallow events show larger vertical uncertainties than the deep events. To validate the robustness of the focal depths of these deep earthquakes, we constructed a subset of 90 events at depths of 20 to 25 km (fig. S5). We forced their depths to be fixed at shallower depths (5 and 10 km), and then the results show that the RMS residuals were significantly increased (fig. S5A), which, in turn, demonstrates that the synthetic shallow depths are not reasonable. In addition, the S-P times of these earthquakes with vertically propagating rays can testify the reliability of our earthquake depth. Unfortunately, in our dataset, we lack earthquakes occurring directly beneath any of the stations (Fig. 2). But waveforms of OBS04 show small S-P time differences for earthquake swarms (nos. 1, 4, 5, 7, and 8) with similar focal depths (20 to 25 km) and epicenters (Figs. 2 and 3), which share similar ray paths. All the evidence indicates the reliability of our earthquake depth determination. Even using very different 1D velocity models, these deep earthquakes between 70 and 90 km distance were well located, with small uncertainties (fig. S4, C to E). Therefore, we believe that the deep events are not an artifact of the network geometry. All 209 earthquakes have RMSs less than 0.3 s with the mean semimajor axis of ~3.9 km (fig. S4).

S-wave delays

The unconsolidated sediments at MORs can cause large delays in S-arrivals (29, 35), which usually affect the focal depths in the earthquake location process. Here, we assume that the direct S-wave delays caused by sediments share the same values with that between the converted S-wave from the basement sediment interface and first P-wave arrival times, where the P-arrivals are strong on the vertical component and the corresponding S-arrivals on the horizontal component (fig. S6). The S-wave delays range from 130 ms (OBS03) to 1280 ms (OBS06) (table S2). In comparison with the initial location results, after removing the S-wave delays from the original S-onsets, Wadati diagrams show a better fit (fig. S7), more phases are used (2609 versus 1621) (fig. S8), and lower residuals are obtained (fig. S8). The S-wave velocity can be approximated from the Vp/Vs ratio by the Wadati diagram (Vp/Vs ~1.7) (fig. S7).

Double-difference hypocenter relocation

To improve the relative location accuracy, we relocated the NonLinLoc hypocenters using the hypoDD double-difference earthquake location algorithm (30) for the 177 microearthquakes detected on more than six OBSs, with RMS residuals of <0.25 s and azimuthal gaps of <270°. The differential travel times for both P- and S-phases were calculated from the catalog. Station corrections calculated by the NonLinLoc program were successfully applied. A minimum of eight catalog links per event pair was required to form a continuous cluster, with a solution obtained using a method of LSQR (49). Five iterations were carried out, with a maximum event separation of 4 km.


Earthquake magnitudes were calculated using the local magnitude scale ML (50): ML = lgA + 1.11 lg(D) + 0.00189 D − 2.09. The maximum amplitude A is measured on a seismogram simulating the original Wood-Anderson seismogram using the SEISAN program (47). D is the hypocentral distance in kilometers. Figure S9 shows the most frequent magnitudes are ~1 to 2. We calculated the residual magnitudes, which were given by subtracting the average magnitude of all stations from the absolute magnitudes at each station (fig. S10). Magnitude results at stations OBS12, OBS16, and OBS19 were removed on the basis of the observed bias in individual magnitude distribution diagrams (fig. S10).

Focal mechanism solutions

We determined P-phase first-motion polarities from unfiltered earthquake waveform data on the vertical component, which were subsequently used to calculate the focal mechanism solutions using the HASH software package (51). We believe that those eight small earthquake swarms show similar focal mechanisms, so we determined focal mechanism solutions for each earthquake swarm. Using a selection criterion based on P-wave polarities of >8, an azimuthal gap of <120°, an RMS fault plane uncertainty of <40°, average misfit <20%, station distribution ratio of >0.4, and mechanism probability >70%, we obtained well-constrained focal mechanisms for three earthquake swarms (Fig. 3, D to F, and table S3). But fig. S11 shows that the solutions are not very clustered, which may be caused by the limitation of our observation network. To validate the influence of the depth uncertainty, we forced their focal depths to be at shallower or deeper depths (fig. S12). The angular difference will change larger at shallow depths than that at deeper depths, which also indicates that the results for our deep earthquakes are robust. To validate the influence of the 1D velocity models, we calculated mechanism solutions using different 1D velocity models in fig. S4A. Table S3 shows that the solutions did not change much, which agrees with the fact that the focal mechanism calculations are not sensitive to the absolute velocity values (51).

Estimates of isotherms

The temperature contours in Fig. 2B were calculated using a plate cooling model (5254), with a lithospheric thickness of 106 km (8). We assumed the mantle potential temperature is 1350°C (54). The half spreading rate of the lithosphere in this region is 16 mm year−1. We calculated the temperature of ridge segments north and south of the Romanche transform fault using the plate cooling model, and then the temperature within the fault zone was estimated by averaging the temperature on both sides. A 1D profile along the middle of the fault valley was extracted and used as the temperature in the transform fault, as shown in Fig. 2B.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to the officers, crew, and scientific technicians of the 2019 SMARTIES cruise for their hard work. We thank P. Moimeux and his crew for the assistance during the cruise. We are very grateful to W. C. Crawford for helping with automatic phase picks and OBS data preprocessing. The discussions with C. Prigent on petrology and M. Cuffaro on thermal modeling were extremely useful. We would like to thank S. Hicks and two anonymous reviewers for constructive reviews, which improved the paper significantly. This is an Intitut de Physique du Globe de Paris contribution number 4206. Most of the figures are made using GMT software (55). Funding: The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC Advance Grant agreement no. 339442_TransAtlanticILAB. This work was also supported by the ISblue project, Interdisciplinary graduate school for the blue planet (ANR-17-EURE-0015), the Italian Programma di Rilevante Interesse Nazionale (PRIN_2017KY5ZX8), and cofunded by a grant from the French government under the program “Investissements d’Avenir.” The ship time for the SMARTIES cruise was funded through the TGIR French Oceanographic Fleet. Author contributions: Z.Y. processed the microseismicity data and wrote the paper. S.C.S. developed the project, supervised the data processing and interpretation, and wrote the paper. M.M., D.B., and E.P.M.G. designed the SMARTIES project. M.M., E.P.M.G., Z.W., and D.B. participated in the data collection. Z.W. carried out the thermal modeling. All authors discussed results and contributed to the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. The earthquake catalog obtained by this study and 1D velocity models can be found in Zenodo:

Stay Connected to Science Advances

Navigate This Article