Constraining the controlling factors of fault rupture is fundamental to understanding the earthquake cycle. Fluids can influence earthquake locations and magnitudes, although the exact pathways of fluids through the lithosphere are not well-known. Ocean transform faults are an ideal laboratory to study the factors controlling fault ruptures and fluid pathways given their relative simplicity. Here, we analyse seismicity recorded by the Passive Imaging of the Lithosphere-Asthenosphere Boundary (PI-LAB) experiment, centred around the Chain Fracture Zone. We find that earthquakes beneath mapped morphological transpressional features occur deeper than the brittle-ductile transition predicted by simple thermal models but elsewhere occur shallower. These features are characterised by multiple parallel fault segments and step overs, high b-values, gaps in large historical earthquakes, and seismic velocity structures consistent with hydrothermal alteration. This suggests that broader fault damage zones preferentially facilitate fluid transport into the lithosphere. Although this cools the mantle, it also reduces the potential for large earthquakes at punctuated locations (barriers). These barriers divide the transform into asperity segments that are shorter, thereby limiting the earthquake magnitudes in these regions.
Article
Broad fault zones enable deep fluid transport and limit earthquake magnitudes
https://doi.org/10.21203/rs.3.rs-1600502/v1
This work is licensed under a CC BY 4.0 License
Journal Publication
published 16 Sep, 2023
Version 1
posted
You are reading this latest preprint version
Ocean transform faults (TF) are ideal places to study the factors controlling seismicity, since they are relatively simple in comparison to their hazardous, continental counterparts. The composition of the oceanic lithosphere is thought to be relatively homogeneous, and the associated fault structures are typically simple i.e., nearly linear and vertical with deformation localised in a narrow, 1–16 km wide, zone9. Plate velocities are also relatively well-known based on magnetic lineations and/or plate motion reconstructions10. Finally, studies of teleseismic earthquake locations11 and earthquake rupture experiments12 show that the depth limit of the region that can potentially slip seismically occurs at the brittle ductile transition, i.e., near the 600°C isotherm, which is well predicted by thermal models.
Despite this relative simplicity, several observational constraints suggest there is still much to learn regarding the factors that control the seismogenic processes. For instance, the amount of moment released seismically on TF systems is smaller than that predicted by plate motions and simple thermal models, ranging from 0 to 83% globally with a mean of 18%7. In addition, oceanic transforms typically do not generate earthquakes much larger than M = 7, i.e., much smaller than the potential size of rupturing the entire transform length5. Finally, earthquakes have been located deeper than the predicted depth of the brittle-ductile transition from thermal models at the Gofar, Blanco, and Romanche TFs13–15.
Fluids likely play an important role in explaining these discrepancies and controlling earthquake rupture location, size, and style. Oceanic transforms are thought to be regions where large amounts of fluids enter the lithosphere4. Specifically, fluids within the fault zone directly affect the pore pressure (i.e., the effective normal stress), potentially advancing the time-to-failure 1,2,16. This process can also reduce the likelihood of large earthquake rupture17,18. Fluids also alter the surrounding rocks, forming serpentinites and other hydrous mineral phases that coat the fault interface and are characterised by weaker rheological and frictional behaviour than their peridotitic and mafic protoliths19. In addition, hydration decreases the strength of olivine, which results in higher strain rates, by up to an order of magnitude, potentially facilitating aseismic creep and/or microseismicity rather than brittle deformation20. Hydrothermal circulation may also cool the lithosphere extending the brittle ductile transition and resulting in deeper seismicity13–15, 21. Variable seawater infiltration and grain size along strike may effectively segment the fault into brittle and ductile sections, as suggested by observations from exhumed peridotites22,23. Geodynamic models suggest that dilatancy mechanisms, whereby fluids become more viscous under pressure may be in effect periodically24,25. Fluids thus play a primary role, however, their exact pathways and relationship with observed earthquakes are not well-known.
Morphology and Seismicity of the Chain Transform Fault
We present the seismicity recorded by the deployment of 39 ocean bottom seismometers (OBS) from March 2016–2017 surrounding the Chain Fracture zone as part of the Passive Imaging of the Lithosphere-Asthenosphere Boundary (PI-LAB) Experiment and the Experiment to Unearth the Rheological Lithosphere-Asthenosphere Boundary (EURO-LAB)8. The seafloor of Chain TF varies in age from 0–20 Myr old across the rupture zone. The active fault zone has a width of up to ~ 15 km, stretching over 300 km in roughly the east-west direction between two spreading segments of the Mid-Atlantic Ridge and slipping at a rate of ~ 30 mm/yr.
The Chain TF contains four lozenge-shaped topographic highs and shows evidence for variability in gravity anomalies and inferred crustal structure along the length of the TF8 (Fig. 1a, b). These features represent crustal blocks exhumed via transpressional stresses, i.e., features typically referred to as positive flower structures (Methods). The eastern section is dominated by the largest flower structure (ELFS; Fig. 1a, Extended Data Fig. 1b), the central section has no flower structures, and the western section is characterised by three small flower structures: western flower structures 1, 2 and 3 (WFS1, WFS2, and WFS3; Fig. 1a, Extended Data Fig. 1a) from west to east.
The Chain TF is also characterized by several fault scarps. There are three relatively long segments: one 30 km long scarp intersecting WFS2 and extending 20 km to the east, one 30 km scarp roughly mid-way between WFS3 and ELFS, and one 50 km long scarp located completely within ELFS (Extended Data Fig. 1). The longer fault segments are separated by much shorter fault scarps typically 2–10 km in length, that often overlap and step over within the fault zone. The overlapping regions occur primarily within and nearby the flower structures, with relatively straight, singular scarps observed in the remaining sections of the TF.
25 MW≥5.6 historical earthquakes along the Chain TF are reported in the relocated catalogue26 which extends back to 1993 (Fig. 1c). Empirical relationships suggest a fault rupture length of ~ 33 km for the largest (MW = 6.6) recorded earthquakes27, in agreement with the scale of two of the longest fault scarps8 that are located on the eastern side of WFS2 and between WFS3 and ELFS (~ 30 km) (Fig. 1). No large earthquakes have been recorded from the longest, 50 km, fault length located within ELFS since 199326. Assuming this rough segmentation scale of 30–50 km, the Chain Transform can be divided into 2 asperity regions (A1, A2) alternating with three barrier regions (B1, B2, B3) where no large (MW≥6.0) earthquakes have been recorded, but are rather characterized by aseismic slip and/or seismicity swarms (Fig. 1, Extended Data Fig. 2). B1 and B3 are located near the ends of the transform. Transform sections near the ridge, e.g., within 20 km26 do not typically host the largest earthquakes given the rapidly changing stress regime28. However, B1 and B3 extend to much further distances from the ridge (> 50 km) and include substantial length scales within regions predicted to be characterised by more uniform stress regimes28. Therefore, seismicity characteristics are not necessarily only a function of proximity to the ridge in these locations. The three barrier regions are co-located with the location of flower structures (WF1, WF3, and ELFS), whereas the asperity regions are co-located with WF2 and the central region of the transform where no flower structures are present. Two of the longest mapped fault scarps also coincide with the asperity regions, occurring on the eastern side of A1 and in A2 (Fig. 1). The amount of historical moment released seismically varies considerably along the fault (Methods, Fig. 1c). Beneath barriers B1, B2, and B3 the moment released seismically is small in comparison to the expected potential (2–6%). Beneath asperities A1 and A2 moment release is much closer to the predicted potential (78% and 67%, respectively).
We detected 626 local events (370 above the magnitude of completeness, MC=2.3) with locatable epicentres31 (Fig. 1; Methods). The number of events per km along strike was highest (2.5 events/km) in segment B3 (Fig. 1b, c). There were 0.54 and 1.36 events/km in B1 and B2, respectively, and 0.98 and 0.82 events/km in A1 and A2, respectively. Depths were sufficiently constrained for 89 events (Methods). These events were distributed from just below seafloor (bsf) to ~ 22 km bsf, mostly at < 15 km depth (82 events). 7 events were located between 19 and 22 km bsf, well below the depth of the predicted brittle-ductile transition based on the 600°C isotherm from simple thermal models (Fig. 2, Extended Data Fig. 3). The deepest events were located directly beneath the WFSs and the ELFS. For 47 events it was possible to determine robust focal mechanisms (Fig. 1a, Methods). Most of those earthquakes (34) occurred with predominantly strike-slip focal mechanisms. However, there were 13 earthquakes with a dominant vertical (normal or reverse) fault component (45° <|rake|<135°), 11 of which are located directly beneath or within 5 km of flower structures (Fig. 1a). Stress inversion results suggest that simultaneous activation of both strike-slip and reverse faults is possible (Methods, Extended Data Fig. 4). Elevated b-values, indicative of greater numbers of smaller earthquakes, are found in regions of shallow bathymetry and low rMBA, i.e., characteristics common to the flower structures (Fig. 1, Fig. 3, Methods; Extended Data Animation 1).
Seismic Velocity and Thermal Structure
We interpret our result considering two seismic velocity models. We present a P-wave refraction model to enhance our interpretation of shallow (< 20 km depth below sea level) depths beneath the ELFS (Fig. 2c, Methods, Extended Data Fig. 5, Extended Data Fig. 6). We also compare our result to thermal structure predicted from the S-wave model derived from local Rayleigh and ambient noise tomography (Methods)29. We preferentially interpret the P-wave model at shallow depths beneath the ELFS, but also considering the S-wave model for deeper velocities and coverage along the entire transform (Methods).
Hydrothermal Alteration and Segmentation along Chain
Several observations beneath barrier B3 and ELFS suggest that another factor besides simple fault length and thermal structure influences this area. Low rMBA at barrier B3 could be consistent with either thickened crust or interaction with fluids in the crust and uppermost mantle because these both decrease the density of the lithosphere as a whole. Hydrothermal alteration of the shallow lithospheric mantle (< 15 km beneath the sea surface) may be suggested, given the anomalously slow sub-crustal P-wave velocities down to at least 15 km (Fig. 2). Alteration caused by fluid infiltration may also explain why the section containing the longest (50 km) fault scarp in the region is not associated with the largest earthquakes, but instead associated with elevated b-values (1.1–1.2) with very low amounts of energy released seismically (4%). The presence of fluids causes a series of mechanical and chemical effects, potentially increasing the pore pressure and reducing the effective normal stresses. This in turn can limit the maximum earthquake magnitude, potentially switching the rheological behaviour from brittle fracture of large earthquakes to aseismic creep and/or earthquakes with very small magnitudes32–35. Fluid migration itself can also result in swarms of very small earthquakes13,36,37. The earthquakes may illuminate fluid pathways, as fluid penetration into the upper mantle can cool the lithosphere and extend the brittle ductile transition to deeper depths21. Therefore, deeper observed seismicity in the OBS data beneath the ELFS may be explained if the cooling of the hydrothermal system has reduced the mantle temperatures depressing the depth of the brittle-ductile transition according to simple thermal models. The S-wave velocity model predicts that the 600°C isotherm associated with the brittle-ductile transition occurs deeper than predicted for a thermal model, supporting the hypothesis that deeper cooling has occurred.
Observations in asperity A2 are consistent with lower degrees of hydrothermal alteration. b-values are lower, and there have been several M ≥ 6.0 historical earthquakes, both consistent with a relatively strong, velocity weakening fault segment, and the lack of small earthquakes caused by fluid migration (Fig. 1c). The moderate rMBA is consistent with a normal crustal thickness and does not require high degrees of alteration. There is a lack of seismicity at > 20 km depth, i.e., the depth of the 600°C from simple thermal models, consistent with a more typical depth to the brittle ductile transition. In addition, the predicted temperatures from the S-wave velocity model are hotter at 25–30 km depth than surrounding regions, suggesting less hydrothermal cooling in this region (Fig. 2b). Resolution testing suggests that the observed variation in the seismic velocity model from which the predicted temperatures are derived along the TF is robust (Methods).
The WFSs observations support the influence of fluids in a similar way to those of ELFS. One exception is that the rMBA in the WFSs region changes from negative to positive from east to west. However, the higher values are likely caused by normal to thinner crustal thickness, which likely dominates over any hydrothermal cooling signature8. The large historical earthquakes beneath WFS2 also distinguish it from the other flower structures. One explanation is that the long fault segment in the region extends 20 km east of the flower structure and may behave more similarly to the long fault in section A2. Another explanation could be that there are temporal variations in rheological behaviour as water is delivered and transported through the system. Episodic fluid delivery and transport and multi-mode rupture models have been previously proposed to explain observed seismicity patterns5,7,38 and have also been predicted by geodynamic modelling25 and inferred from samples from exhumed oceanic TFs23.
Seismicity at depths that exceed the brittle ductile transition predicted by simple thermal models is consistent with previous observations from TF13–15. These studies also find seismicity deeper than predictions from numerical modelling that includes the cooling effects from hydrothermal circulation14,21. Therefore, it was hypothesized that the cooling must be stronger and deeper than allowed in the numerical modelling experiment21. It has also been suggested that ruptures at higher temperatures should be expected in locations of higher strain rates such as transform faults39. Exhumed hydraulically altered high-strain mylonites were also used to argue that brittle and ductile deformation can occur over a broad range of temperatures (300°C – 1000°C) and, therefore, a broad range of depths, owing to variable seawater infiltration and grain sizes along the fault22. Our result does not preclude the possibility that different rheological flow laws22 are in effect in some locations globally. However, the flow laws are not required to explain the seismicity at Chain, which occurs shallower than the 600°C isotherm derived from the S-wave velocity model. In other words, it is more likely that our region has been cooled by hydrothermal circulation. Similar velocity models are not available at most locations where deeper than predicted seismicity has been reported to verify the degree to which deeper and stronger hydrothermal cooling are in effect in comparison to a grain size impacted rheological effect.
Overall, we find that seismicity on the Chain TF is segmented with stronger, locked regions failing seismically (asperities), separated by weaker regions with aseismic slip and/or small earthquake magnitudes (barriers). The weaker regions are located beneath the morphological flower structures, which are also regions characterised by several parallel and overlapping fault strands, arranged as en-echelon structures i.e., much different than other sections characterised by singular, relatively linear faults. The thermal structures inferred from S-wave velocities are also consistent with greater amounts and deeper extents of hydrothermal alteration beneath the barrier regions. Fluids likely cool the plate to greater depths, allowing deep earthquakes to occur, but also weaken the mantle locally, decreasing the likelihood of large earthquakes (M > 7.0), since the ruptures are arrested, preventing the entire transform from failing. Evidence for both deep seismicity from the OBS data and large historical earthquakes beneath anomalous WFS2 indicates that either the larger earthquakes are related to the continuation of the fault strand 20 km east of WFS2, i.e., far from the overlapping fault strands on WFS2 or that barrier-asperity behaviour can vary temporally at any given location. Overall, regions of broad deformation with overlapping fault zones preferentially hydrate the lithosphere, which cools it and deepens the brittle-ductile transition. This process weakens the lithosphere and effectively segments the transform fault and limits the magnitudes of the largest earthquakes.
Data availability
The earthquake catalogue used in the study is available as a pre-print31. The ocean bottom seismic data are archived at the IRIS DMC, as 2016–2017 network XS https://doi.org/10.7914/SN/XS_2016. The relocated data since 199326 are available at https://doi.org/10.5281/zenodo.4646438.
Code availability
Some figures were generated using Generic Mapping Tools v.6.1.1 (www.soest.hawaii.edu/gmt, last accessed March 2022). The code for the Anderson-Darling test can be downloaded from https://git.plgrid.pl/projects/EA/repos/sera-applications/browse/, last accessed March 2022.
Acknowledgements
K. Leptokaropoulos, C. A. Rychert and N. Harmon acknowledge funding from the Natural Environment Research Council (NE/M003507/1) and the European Research Council (GA 638665). J. -M. Kendall was funded by the Natural Environment Research Council (NE/M004643/1). We thank the captain and crew of the R/V Marcus G. Langseth and the RRS Discovery, and also the scientific technicians. The active source seismic data were acquired during the ILPAB-SPARC experiment on board the French R/V Pourquoi Pas funded by the European Research Council Advanced Grant number 339442_TransAtlanticILAB. Refraction data were analysed by Christian Filbrandt.
Author contributions
K.L., C.R. and N.H. developed the main ideas, performed analyses (b-values, statistical analyses velocity/thermal models) and generated figures. K.L. produced the initial manuscript. D.S. analysed the OBS data to determine locations, magnitudes and focal mechanisms. I.G. provided the P-wave refraction model for the eastern part of Chain. All the authors discussed the results and contributed to writing the manuscript.
Competing interests
The authors declare no competing interests.
- Beeler, N. M., Simpson, R. W., Hickman, S. H., & Lockner D. A. Pore fluid pressure, apparent friction and coulomb failure. J. Geophys. Res.105, 25,533-25,542 (2000).
- Cocco, M., & Rice J. R., Pore pressure and poroelasticity effects in Coulomb stress analysis of earthquake interactions. J. Geophys. Res.107, B2, 10.1029/2000JB000138 (2002).
- Mulargia, F., & Bizzarri, A. Anthropogenic triggering of large earthquakes, Sci. Rep., 4: 6100 (2014).
- Hensen C., Duarte J. C., Vannucchi P., Mazzini A., Lever M. A., Terrinha P. et al. Marine Transform Faults and Fracture Zones: A Joint Perspective Integrating Seismicity, Fluid Flow and Life. Front. Earth Sci.7:39 (2019).
- Boettcher, M. S. & Jordan T. H. Earthquake scaling relations for mid-ocean ridge transform faults. J. Geophys. Res.109, B12302 (2004).
- Behn, M. D., Boettcher M. S. & Hirth G. Thermal structure of oceanic transform faults. Geology35, 307-310 (2007).
- Wolfson-Schwehr, M., & Boettcher, M. S. “Chapter 2 – global characteristics of oceanic transform fault structure and seismicity,” in Transform Plate Boundaries and Fracture Zones, ed. J. C. Duarte (Amsterdam: Elsevier), 21–59 (2019).
- Harmon, N., Rychert, C., Agius, M., Tharimena, S., Le Bas, T., Kendall J. M. & Constable, S. Marine geophysical investigation of the Chain Fracture Zone in the equatorial Atlantic from the PI-LAB experiment. J. Geophys. Res. Solid Earth123, 11,016–11,030 (2018).
- Ren, Y., Geersen, J., & Grevemeyer, I. Impact of spreading rate and age-offset on oceanic transform fault morphology. Geophys. Res. Lett. 49, e2021GL096170 (2022).
- Müller, R. D., Sdrolias, M., Gaina, C. & Roest, W. R. Age, spreading rates, and spreading asymmetry of the world's ocean crust. Geochem. Geophys., Geosys., 9, Q04006 (2008).
- Abercrombie R. E. & Ekström G. Earthquake slip on oceanic transform faults. Nature410, 74-77 (2001).
- Boettcher, M. S., Hirth, G. & Evans, B. Olivine friction at the base of oceanic seismogenic zones. J. Geophys. Res. 112, B1, 74 (2007).
- Kuna, V. M., Nabělek J. L. & Braunmiller J. Mode of slip and crust–mantle interaction at oceanic transform faults. Nat. Geosci.12, 138–142 (2019).
- McGuire, J. J., Collins, J. A., Gouedard, P., Roland, E., Lizarralde, D., Boettcher, M. S., Behn, M. D. & van der Hilst, R. D. Variations in earthquake rupture properties along the Gofar transform fault, East Pacific Rise. Nat. Geosci.5, 336–341 (2012).
- Yu, Z., Singh, S. C., Gregory, E. P. M., Maia, M., Wang, Z. & Brunelli, D. Semi-brittle seismic deformation in high-temperature mantle mylonite shear zone along the Romanche transform fault. Sci. Adv.7, eabf3388 (2021).
- Nostro, C., Chiaraluce, L., Cocco, M., Baumont, D. & Scotti, O. Coulomb stress Coulomb stress changes caused by repeated normal faulting earthquakes during the 1997 Umbria-Marche (central Italy) seismic sequence. J. Geophys. Res., 110, B5, (2005).
- 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. Bull. Geol. Soc. Am.70, 115-166 (1959).
- Liu, Y., McGuire, J. J. & Behn, M. D. Aseismic transient slip on the Gofar transform fault. East Pacific Rise, Proc. Natl. Acad. Sci. Unit. States Am.117, 19, 10188-10194 (2020).
- Reinen, L. A., Weeks, J. D. & Tullis T. E. The frictional behavior of serpentinite: Implications for aseismic creep. Geophys. Res. Lett.18, 1921–1924 (1991).
- Hirth, G. & Kohlstedt, D. L. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth Planet. Sci. Lett. 144, 93–108 (1996).
- Roland, E., Behn, M. D. & Hirth G. Thermal‐mechanical behavior of oceanic transform faults: Implications for the spatial distribution of seismicity. Geochem. Geophys. Geosyst.11, Q07001 (2010).
- Kohli, A., Wolfson-Schwehr, M., Prigent, C. & Warren, J. M. Oceanic transform fault seismicity and slip mode influenced by seawater infiltration. Nat. Geosci., 14, 606-611 (2021).
- Cox, S., Fagereng, A., & MacLeod, C. J. Shear zone development in serpentinized mantle: Implications for the strength of oceanic transform faults. J. Geophys. Res. Solid Earth126, e2020JB020763 (2021).
- Géli, L., J.-M. Piau, R. Dziak, V. Maury, D. Fitzenz, Q. Coutellier & P. Henry (2014), Seismic precursors linked to super-critical fluids at oceanic transform faults, Nat. Geosci., 7, 757-761.
- Liu, Y., McGuire, J. J. & Behn, M. D. Frictional behavior of oceanic transform faults and its influence on earthquake characteristics. J. Geophys. Res.117, B04315 (2012).
- Shi, P., Meng W. & Pockalny, R. The ubiquitous creeping segments on oceanic transform faults, Geology50, 199-204 (2021).
- Wells, D.L. & Coppersmith, K.J. New empirical relationships among magnitude, rupture length, rupture width, rupture area and surface displacement. Bull. seism. Soc. Am.84, 974–1002 (1994).
- Parmentier, E. M. & Morgan, J. P. Spreading rate dependence of three‐dimensional structure in oceanic spreading centres. Nature 348, 325–328 (1990).
- Saikia, U., Rychert, C., Harmon, N. & Kendall, J. M. Upper mantle anisotropic shear velocity structure at the equatorial Mid-Atlantic ridge constrained by Rayleigh wave group velocity analysis from the PI-LAB experiment. Geochem. Geophys. Geosyst.22, e2020GC009495, (2021).
- Stixrude, L. & Lithgow-Bertelloni, C. Thermodynamics of mantle minerals – I. Physical properties, Geophys. J. Int. 162, 610-632 (2005).
- Schlaphorst, D., Rychert, C., Harmon, N., Hicks, S., Bogiatzis, P., Kendall, J. M. & Abercrombie, R. Local seismicity around the Chain Transform Fault at the Mid-Atlantic Ridge from OBS observations. preprint submitted in Geophys. J. Int. doi: 10.1002/essoar.10511147.1.
- Froment, B., McGuire, J. J., van der Hilst, R. D., Gouedard, P., Roland, E. C., Zhang, H. & Collins J. A. Imaging along-strike variations in mechanical properties of the Gofar transform fault, East Pacific Rise. J. Geophys. Res. Solid Earth119, 7175–7194 (2014).
- Roland, E., Lizarralde, D., McGuire, J. J. & Collins J. A. Seismic velocity constraints on the material properties that control earthquake behavior at the Quebrada-Discovery-Gofar transform faults, East Pacific Rise. J. Geophys. Res.: Solid Earth117, B11102 (2012).
- Cox, S., Ikari, M. J., MacLeod, C. J., & Fagereng, A. Frictional characteristics of oceanic transform faults: Progressive deformation and alteration controls seismic style. Geophys. Res. Lett.48, e2021GL096292 (2021).
- Renard, F. Reaction-induced fracturing: when chemistry breaks rocks. J. Geophys. Res. Solid Earth.126, e2020JB021451 (2021).
- El-Isa, Z.H. & Eaton D. W. Spatiotemporal variations in the b-value of earthquake magnitude-frequency distributions: classification and causes. Tectonophysics615-616, 1-11 (2014).
- Tilmann, F., Grevemeyer, I., Flueh, E.R., Dahm, T. & Goßler, J. Seismicity in the outer rise offshore southern Chile: indication of fluid effects in crust and mantle, Earth Planet Sci. Lett.269, 41-55, 2008.
- Aderhold, K. & Abercrombie, R. E. The 2015 Mw 7.1 earthquake on the Charlie-Gibbs transform fault: Repeating earthquakes and multimodal slip on a slow oceanic transform. Geophys. Res. Lett.43, 6119–6128 (2016).
- Molnar, P. The Brittle‐Plastic Transition, Earthquakes, Temperatures, and Strain Rates. J. Geophys. Res. Solid Earth. 125, e2019JB019335 (2020)
Origin of the lozenge-shaped topographic highs. One possibility is that these features represent intra-transform volcanism40. However, the observed low backscatter intensity of these structures8 does not support recent volcanism. A more likely scenario is that these features represent crustal blocks exhumed via transpressional stresses underlain by altered mantle, i.e., features typically referred to as positive flower structures. Similar morphological features, are frequently observed and fault plane locations are sometimes available to substantiate a transpressional origin. This is supported by the highly tectonised mafic and ultra-mafic compositions of rocks from features with similar morphology and gravity signatures at the nearby St. Paul TF that are also thought to have a transpressional origin41.
Seismic catalogue and focal mechanisms. The analysis of the broadband ocean bottom seismic deployment revealed 972 events, 812 of which are located within the area surrounded by the OBS network31, which belongs to the Chain TF and the adjacent ridge spreading centres (Fig. 1b). The location of the events is performed by NonLinLoc software42 using a 1-D velocity model of the Chain TF region based on CRUST1.043 and earthquake arrival times31. In order to focus on the TF seismicity, we discard the events that unequivocally occurred at the ridge segments and at the inside corners after considering their location and focal mechanisms (north-south striking, normal faulting events). This selection leaves 626 events along the Chain transform valley, within a local magnitude (ML) range between 1.1 and 5.4. We further remove events below the completeness threshold MC=2.3 (see below), resulting to a dataset of 370 events along the transform valley, which is used for our analysis.
The epicentral coordinates of these events are sufficiently constrained with median lateral uncertainties of 2.5 km. Vertical uncertainties are larger (median ~18 km) and several events are above 20 km. For this reason, the hypocentre depths as well as focal mechanisms of well-recorded events are re-evaluated using the Grond software44, which carries out a Bayesian-bootstrap time-domain deviatoric moment tensor inversion. We calculated 119 focal mechanism solutions. For the stress inversion we select the 47 focal mechanisms classified as having good-fits (Fig. 1). We include 42 additional events (89 events total) with sufficiently well-determined depth (mean vertical error of 6km).
We supplement our seismicity data recorded by the OBS with historical earthquakes from the Global Centroid Moment Tensor catalogue45. Due to the relatively large epicentral uncertainties, we only consider the strong events (MW≥5.6) after 1993 relocated by Shi et al.26.
OBS seismicity with respect to bathymetry and rMBA. The bathymetry and residual mantle Bouguer anomaly (rMBA) data acquisition and analysis are described in Harmon et al.8. Here, we investigate whether there is preference for seismic events in the local OBS data occurring at sites with particular rMBA/ bathymetry values. In doing so, we compare the rMBA/bathymetry values at the epicentral location of events with the background rMBA/bathymetry data. Background data are defined as all the rMBA/bathymetry values calculated within a 15km wide zone, centred on the transform valley axis and bounded in the along-strike direction by our seismicity data (roughly between the ridge spreading centres). In such way, we consider only the values corresponding to the active deformation zone along the transform valley, where the vast majority of seismicity in the local OBS data occurs, discarding the essentially aseismic areas elsewhere. If the earthquake epicentres show no preference for particular parameter values, then the red and blue histograms in Fig. 3d and Fig. 3e should be similar. We see, however, that there is a clear preference of events occurring at negative rMBA values and an event deficit for positive rMBA (Fig. 3d). There is also a preference of events occurring at very shallow ocean depths (~<3700m, Fig. 3e), i.e., directly beneath the ELFS.
Magnitude distribution. We investigate the seismicity magnitude distribution dependence on bathymetry and rMBA, after assigning a bathymetry and rMBA value to each earthquake by nearest neighbour interpolation. Then the events are sorted by bathymetry and rMBA values in an ascending order and the b-values from the corresponding event magnitudes are calculated, together with standard error. Calculations are performed for overlapping 75-event windows, sliding by 1 event (different windows of 50 to 100 events were also tested, as shown in Extended Data Animation 1).
The exponentiality of magnitude distribution is investigated by the Anderson-Darling Test46,47 and further established by goodness of fit test48,49, applying the Maximum Likelihood Estimation (MLE) of b-value50. Both techniques suggest that the magnitude distribution can be sufficiently modelled by an exponential distribution for ML≥2.3. A total of 370 events with ML≥2.3 comprise the complete data for the study area (circles, Fig. 1b), which roughly corresponds to an average rate of ~1 event/day. However, the MLE returns considerably unstable b-values, even well above Mc, fluctuating between 0.76 and 0.89 for 2.3≤ML≤3.0 (Extended Data Table 1, Extended Data Fig. 7). This is incompatible with the unique b-value predicted by the GR law and may lead to artifacts and misinterpretation of the derived outcome. For this reason, we apply alternative methods which are more suitable for small datasets and are less sensitive to MC selection51,52. The comparison of the results from the aforementioned methods, assuming different MC are shown in Extended Data Fig. 7. We choose to estimate b-values by applying the repeated medians53 (RM), a non-parametric regression technique, which provides more robust results for small data sets, since it is highly resistant to observational uncertainties and outliers51. According to the RM approach, we consider n magnitude intervals between two points, i and j, having Ni and Nj events, respectively, and calculate n-1 slopes (b-values) between these points as:
With Mj≠Mi. For each point, i, we calculate the median slopes, thus n median values. Then, the bRM is estimated as the median of the n median values:
bRM=-median(median(bij))
The standard errors, SE, of the corresponding bRM are derived by bootstrap resampling process as described in Amorèse et al.51. The significance of the b-value difference for two datasets, e.g. A and B, can be estimated by a bootstrap t-test, with the t statistic defined as:
A t=1.96 corresponds to a significance at 0.05 level. This approach returns b=0.83±0.09, with this value being insensitive to the applied magnitude cut-off (b=0.82-0.84 for 2.0≤M≤3.0) as shown in Extended Data Fig. 7 and Extended Data Table 1. We therefore apply the RM technique for all b-value calculations in this study, since it provides the most consistent results, regardless the MC choice.
The highest b-values are highest beneath the ELFS (1.1-1.2), they are lowest (<0.8) in the central parts of the fault (~14.7°E – 13.8°E) and intermediate (0.9 – 1.1) in the western area (Extended Data Animation 1). Shallower water depths generally correspond to higher b-values in comparison to the deeper bathymetry areas, where b-values gradually decline with increasing depth at ~0.15 units/km (Fig. 3a). There are higher b-values (b=1.04±0.13) in regions of low rMBA (<-3mgal) than regions of high rMBA (<-3mgal) where b=0.72±0.11 (Fig. 3b, c) with the difference between these two b-values being statistically significant at 0.05 level, indicated by the calculated t statistic equal to 1.96. Events also occur preferentially in areas with negative rMBA (Fig. 3d). These areas are mainly located in the regions of shallowest bathymetry (<3700 m), i.e., beneath the ELFS (Fig. 1; Fig. 3e).
Seismic Vs Predicted Moment. We compare the seismic moment release, M0, in each asperity (A1, A2) and barrier (B1, B2, B3) using the formula54 , where M0 is given in N∙m. We express the results in terms of percentage (a) of seismically released moment (SM0s) over the potential moment release assuming full seismic coupling (SM0p), i.e., a= SM0s/SM0p. In doing so, we estimate SM0p by summing the seismic moment of all MW≥5.6 events occurred since 1993 in each patch (barrier or asperity). We then estimate M0p as5:
M0p=GLwsDT
where, G is the shear modulus [40-60 GPa], L is the patch length, w, is the patch width from sea floor to the 600 °C isotherm, s is the slip rate [28-33mm/yr] and DT is the time period duration since 1993, equal to 28 years. Our results range from a=0.25 to a=0.46 with an average value of a=0.32 (for G=50GPa and s=32mm/yr). This range agrees well with previous findings (a=0.33) for the entire Chain Transform7. However, we find the proportion of seismically released moment varies along the length of the transform from 0.02-0.06 in barrier regions to 0.67-0.78 in asperity regions (Fig. 1c; Fig. Extended Data Fig. 2).
Stress Inversion. The applied stress inversion algorithm55,56 is used to invert focal mechanisms for a set of earthquakes into tectonic stress, under the assumptions that a uniform regional stress field applies and that the earthquakes occur on faults with varying orientations, but they do not interact with each other. Given that there are only 2 ML>5.0 events, with the stronger earthquake being ML=5.4, the effects of co-seismic stress changes can be neglected57,58 in our study. Under these assumptions, we invert focal mechanisms to estimate the orientation of the principal stress axes (and consequently, the principal focal mechanism) and the shape ratio, R, which determines the relative amplitude of the principal stresses. A typical problem in stress inversion from focal mechanisms concerns the selection of the actual fault plane from the two nodal planes, a choice which may severely bias the results57 mainly the shape ratio determination. To avoid this effect, an algorithm which inverts jointly for stress and fault orientations is considered59. The fault orientation determination is based on the evaluation of fault instability index, I:
where τ and σ are the normalized shear and normal tractions, respectively and m, the friction coefficient. The fault plane is selected as the one having the highest I value, after applying an iterative process. The iteration process performs a linearized inversion55, each time introducing a different noise realization, which randomly rotate the given focal mechanisms. We applied 1000 such realizations in each case tested. The final stress tensor is derived as the mean of the results from all realizations. The uncertainty of the principal stress axes orientations is demonstrated by the scatter of the derived values (Extended Data Fig. 4a), whereas the shape ratio uncertainty is determined as the 95% confidence interval. The instability of principal stress axes shows the impact of the uncertainties from random perturbation of the solution, however it is not clear whether this instability arises exclusively from the spatial stress variations or uncertainties in the focal mechanism determination.
We use the available focal mechanisms along the Chain transform valley to determine the stress field. The vast majority of the reported GCMT solutions45 suggest strike-slip faulting along Chain TF. However, our 1-year seismicity analysis of the OBS data reveals that a proportion of smaller magnitude events (3.0<ML<5.0) demonstrate a considerable vertical component, plausibly related with the positive flower structures (Fig. 1a). We derived 47 best constrained focal mechanisms that are considered for the stress inversion. The principal focal mechanism has strike=275°, dip=62° and rake=165° (Extended Data Fig. 4b). We also derived a shape ratio, R=0.86, with 95% confidence interval [0.69-0.89]. Dividing the region into smaller areas for a spatial analysis would significantly lower the significance of the results due to the decreasing sample size, therefore, we are not able to quantify the spatial variation of the stress along Chain.
The high R value also suggests that the intermediate, s2 and minimum, s3 principal stresses, both deviate considerably from the vertical direction (plunge is 60° for s2 and 30° for s2) and they have similar amplitudes (Extended Data Fig. 4). Therefore, they are virtually indistinguishable and might switch with each other. With s1 being clearly sub-horizontal (2° plunge), this physically means that under such stress field the simultaneous activation of both strike-slip and reverse faults is plausible60. This suggests a non-negligible reverse movement, related to the flower structures, especially within the ELFS, indicating the existence of active transpressional features along Chain.
S-wave model from local Rayleigh and ambient noise tomography. We compare our result to the S-wave velocity structure derived from local Rayleigh wave and ambient noise tomography29. Group velocity measurements were inverted from 16-40 s period to determine 2-D velocity maps for each period of interest. Then the group velocities at each pixel across all of the maps were inverted for S-wave velocity structure. Additional details can be found in the original work. The lateral smoothing or correlation length in the group velocity map is 100 km, and checkerboard tests from the study suggest that anomalies at this length scale are well resolved by the method, particularly in the region around the Chain Transform Fault as it is within the centre of the array, with many crossing raypaths with a wide range of azimuths. The formal resolution for the inversion for S-wave velocity with depth is 0.1 for 1 km thick layers in the upper 10 km beneath the seafloor and 20 km at greater depths. This suggests that the S-wave velocity inversion can uniquely resolve the average velocity over the upper 10 km of the model beneath the seafloor, the average velocity from 10-30 km beneath the seafloor and the average value from 30-50 km beneath the seafloor. The lateral and vertical resolution suggests that the variation in lithospheric thickness/velocity we observe across the transform fault at ~100 km scale (WFS1-WF3, A2 vs ELFS) and over 50 km depth is robust. The broader thermal structure of the lithosphere in the region is well within the resolving power of the model presented here. The low P-wave velocities just beneath the Moho of the ELFS from the active source model do not appear in the S-wave velocity model derived from the surface waves, although this is expected given the different resolution of the waveforms. The scale of the crustal and shallow mantle variations observed in our P-wave refraction model is~25-30 km laterally along strike of the transform, 10-20 km across the transform and occurring over ~ 5-10 km in depth. This is below the resolution of the tomographic model used here as described above. In other words, the velocity anomalies in the refraction model would be averaged out in the S-wave velocity model with the relatively faster lithosphere on either side of the fault zone.
The 600°C isotherm derived from the S-wave velocity derived model is closer to 30 km depth across the entire transform, deeper and different in shape than that of simple thermal models (Fig 2). The scale of this variation, 100 km laterally, is well resolved as described above. Isotherms with shapes that deviate from those of simple thermal models can be caused by upwelling in the centre of the transform, as predicted by geodynamic modelling with a visco-elastic-plastic rheology6. While such a rheology cannot be precluded here, we also require an additional factor, given that our isotherms are also deeper than predicted by those models.
S-wave model translated into thermal model. We examine two methods for translating seismic S-wave velocity to temperature based on the assumption of a peridotitic mantle. In the main text we use the empirical fit to predictions for a pyrolite mantle30:
where Vs is S-wave velocity in km/s, P is pressure in GPa and T is temperature in Kelvin. We also investigate the effects of different composition assumptions on the conversion from S-wave velocity to temperature. We use the Abers and Hacker61 MATLAB toolbox to calculate the predicted S-wave velocity for three end member compositions. We use a depleted upper mantle (harzburgite), an undepleted mantle (lherzolite) and, for reference, an idealized unmelted mantle composition (pyrolite). The specific mineral modes we use are as follows: harzburgite (fo 72.48% volume, fa 7.52 % , en 18.24%, fs 1.76%), lherzolite (fo 45.41%, fa 5.07%, en 22.48%, fs 2.50% di 20.03, sp 4.50%) and pyrolite (alm 2.10%, gr 1.10%, py 10.80%, fo 54.50%, fa 6.70%, en 14.70%, fs 1.50%, di 7.50%, hed 1.2%). We do not account for phase changes of the aluminious phases in the calculations although the spinel – garnet transition is encompassed somewhat between the lherzolite and pyrolite models. We apply a frequency dependent attenuation correction to the S-wave velocity output from the codes using both the approach and the 1-D attenuation structure62. The velocities were calculated for a range of temperatures at a given pressure corresponding to the depth in the model and the temperature corresponding to each velocity in the model was determined by interpolation.
The overall effect of choosing a different starting composition shallows the isotherms predicted from the S-wave velocity model relative to the Stixrude and Bertelloni29 parameterization, particularly the 600°C isotherms, although the effect is not substantial. With the Abers and Hacker61 parameterization the 700°C isotherm lies at nearly the same depth as the 600°C isotherm from the Stixrude and Bertelloni30 parameterization. The difference between the compositions from the Abers and Hacker61 calculations is small, typically the isotherms are within a few km of each other. If the compositions and calculations from Hacker are used all the well-located events are shallower than the 700°C isotherm.
Seismic tomography of active source airgun profiling. We used data collected at the easternmost section of the Chain Fracture Zone as part of the ILAB-SPARC experiment63 (Extended Data Fig. 5, Extended Data Fig. 6). Seismic data were acquired in October of 2018 aboard the French research vessel Pourquoi Pas?. A 140 km long seismic profile was acquired at the eastern Ridge-Transform Intersection running from the pressure ridge into the adjacent fracture zone (Extended Data Fig. 5). Five OBS spaced every 12 km recorded airguns shots fired with the Pourquoi Pas?; two sub-arrays containing eight G-guns each provided a total volume of 82 litres and were towed at a depth of 10 m.
Travel times of first arrival P-waves and secondary wide-angle-refection interpreted to be crust/mantle boundary reflection (PmP) arrivals are hand-picked. In general, picking uncertainties were 20-30 ms for short-offset P-waves (Pg) and reach 60 ms for far-offset P-waves (Pn) and secondary PmP arrivals.
Seismic refraction travel time data were used to derive 2-D velocity models using a seismic tomography approach64. The method employs a hybrid ray‐tracing scheme combining the graph method with further refinements utilizing ray bending with the conjugate gradients method. Smoothing and damping constraints regularize the iterative inversion. Details of the procedure are given elsewhere23. Picking errors and starting velocity models may control inversion results. We therefore chose a nonlinear Monte Carlo-type error analysis to derive model uncertainties. The approach consists of randomly perturbing the velocity values of an initial average 1-D model to create a set of 100 2-D reference models, providing a well constrained seismic velocity model (Extended Data Fig. 6).
In the centre of the flower structure the <7.0 km/s contour, corresponding to the Moho, is uplifted. On the eastern end of the flower structure, in the mantle beneath the Moho velocities may be as low as 7.0 – 7.5 km/s, i.e., 4 – 10 % slower in comparison to velocities >7.8 km/s on the western side of the ELFS.
References
40. Grevemeyer, I., Rüpke, L.H., Morgan, J.P., Iyer, K, and Devey, C.W. Extensional tectonics and two-stage crustal accretion at oceanic transform faults. Nature591, 402–407, (2021).
41. Maia, M., Sichel, S., et al. (2016) Extreme mantle uplift and exhumation along a transpressive transform fault. Nature Geosci9, 619–623.
42. Lomax, A., Virieux, J., Volant, P., & Berge-Thierry, C. Probabilistic earthquake location in 3D and layered models. Advances in Seismic Event Location, 18, 101– 134 (2000).
43. Laske, G., Masters, G., Ma, Z. & Pasyanos M. Update on CRUST1.0—a 1-degree global model of Earth’s crust. Geophys. Res. Abstr. 15, 2658 (2013).
44. Heimann, S., Isken, M., Kühn, D., Sudhaus, H., Steinberg, A., Vasyura-Bathke, H., et al. Grond - A probabilistic earthquake source inversion framework. V. 1.0. Potsdam: GFZ Data Services (2018).
45. Ekström, G., Nettles, M. & Dziewoński, A. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes. Phys. Earth. Planet. Int.200–201, 1–9 (2012).
46. Marsaglia, G. & Marsaglia J. Evaluating the Anderson-Darling distribution. J. Stat. Soft.9, 1-5 (2004).
47. Leptokaropoulos, K. Magnitude distribution complexity and variation at The Geysers geothermal field. Geophys. J. Int. 222, 893-906 (2020).
48. Wiemer, S. & Wyss M. Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the Western United States, and Japan. Bull. Seismol. Soc. Am.90 (4), 859–869 (2000).
49. Leptokaropoulos K., Karakostas, V., Papadimitriou, E., Adamaki, A., Tan, O. & İnan S. A Homogeneous Earthquake Catalogue for Western Turkey and Magnitude of Completeness Determination. Bull. Seismol. Soc. Am.103(5), 2739-2751 (2013).
50. Aki, K. Maximum likelihood estimate of b in the formula logN = a - bM and its confidence limits. Bull. Earthq. Res. Inst. Tokyo Univ. 43, 237–239 (1965).
51. Amorèse D, Grasso J.-R. & Rydelek P.A. On varying b values with depth: results from computer-intensive tests for Southern California. Geophys. J. Int.180, 347–360 (2010).
52. van der Elst, N. J. B-positive: A robust estimator of aftershock magnitude distribution in transiently incomplete catalogs. J. Geophys. Res.: Solid Earth126, e2020JB021027 (2021).
53. Siegel, A. F. Robust regression using repeated medians. Biometrika69, 242–244, (1982).
54. Hanks, T. C. & Kanamori H. A moment magnitude scale. J. Geophys. Res.84(B5), 2348–50 (1979).
55. Michael, A. J. Determination of stress from slip data: faults and folds. J. Geophys. Res.89, 11,517-11,526 (1984).
56. Vavryčuk, V. Iterative joint inversion for stress for stress and fault orientations from focal mechanisms. Geophys. J. Int.199, 69-77 (2014).
57. Martínez-Garzón, P., Ben-Zion, Y., Abolfathian, N., Kwiatek, G. & Bohnhoff, M. A refined methodology for stress inversions of earthquake focal mechanisms. J. Geophys. Res. Solid Earth121, 8666–8687 (2016).
58. Fojtíková, L. & Vavryčuk, V. Tectonic stress regime in the 2003-2004 and 2012-2015 earthquake swarms in the Ubaye Valley, French Alps. Pure Appl. Geophys.175, 1997-2008 (2018).
59. Vavryčuk, V. Moment tensor decompositions revisited. Pure Appl. Geophys.19, 231-252 (2015).
60. Hallo, M., Opršal, I., Asano, K. & Gallovič, F. Seismotectonics of the 2018 northern Osaka M6.1 earthquake and its aftershocks: joint movements on strike‑slip and reverse faults in inland Japan. Earth Planets Space71:34 (2019).
61. Abers, G. A., & Hacker, B. R. A MATLAB toolbox and Excel workbook for calculating the densities, seismic wave speeds, and major element composition of minerals and rocks at pressure and temperature. Geochem. Geophys. Geosyst.,17, 616–624 (2016).
62. Saikia, U., Rychert, C., Harmon, N. & Kendall, J.M. Seismic attenuation at the equatorial Mid-Atlantic Ridge constrained by local Rayleigh wave analysis from the PI-LAB experiment. Geochem. Geophys. Geosyst., 22, e2021GC010085. (2021).
63. Marjanović, M., Singh, S. C., Gregory, E. P. M., Grevemeyer, I., Growe, K., Wang, Z., et al. Seismic crustal structure and morphotectonic features associated with the Chain Fracture Zone and their role in the evolution of the equatorial Atlantic region. J. Geophys. Res. Solid Earth, 125, e2020JB020275 (2020).
64. Korenaga, J., Holbrook, W. S., Kent, G. M., Kelemen, P. B., Detrick, R. S., Larsen, H. C., et al. Crustal structure of the southeast Greenland margin from joint refraction and reflection seismic tomography. J. Geophys. Res.: Solid Earth, 105, 21591-21614 (2000).
There is NO Competing Interest.
- Leptokaropoulosetal.ExtendedData.docx
Extended Data Table 1,Extended Data Animation 1,Extended Data Figures 1-7,