Broad fault zones enable deep fluid transport and limit earthquake magnitudes

Domain: earth_science
Published: 2023-09-16
Editorial outcome: minor_revision

Abstract
Constraining the controlling factors of fault rupture is fundamentally important. Fluids influence earthquake locations and magnitudes, although the exact pathways through the lithosphere are not well-known. Ocean transform faults are ideal for studying fault and fluid pathways given their relative simplicity. We analyse seismicity recorded by the Passive Imaging of the Lithosphere-Asthenosphere Boundary (PI-LAB) experiment, centred around the Chain Fracture Zone. We find earthquakes beneath 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, higher proportions of smaller events, gaps in large historical earthquakes, and seismic velocity structures consistent with hydrothermal alteration. Therefore, broader fault damage zones preferentially facilitate fluid transport. This cools the mantle and reduces the potential for large earthquakes at localized barriers that divide the transform into shorter asperity regions, limiting earthquake magnitudes on the transform.

Introduction
The role of fluids is fundamental to understanding earthquake rupture location, size and style1–4. Oceanic transform faults are ideal places to study the factors controlling seismicity, since they are relatively simple in comparison to their hazardous, continental counterparts5–7. 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 localized in a narrow, 1–15 km wide, zone8. Plate velocities are also relatively well-known based on magnetic lineations and/or plate motion reconstructions9. Finally, studies of teleseismic earthquake locations10 and earthquake rupture experiments11 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 seismogenic processes. For instance, the amount of moment released seismically on transform fault systems is smaller than that predicted by plate motions and simple thermal models, ranging from 0 to 83% 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. For instance, at the Chain Transform Fault the magnitude of rupturing the entire transform is predicted to be M > 8.012 based on empirical relationships. Finally, earthquakes have been located deeper than the predicted depth of the brittle-ductile transition from thermal models at the Gofar, Blanco, and Romanche Transform Faults13–15.

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-failure1,2,16. This process can 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 characterized by weaker rheological and frictional behavior 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 micro-seismicity rather than large-scale ruptures20. Hydrothermal circulation may also cool the lithosphere extending the brittle ductile transition and resulting in deeper seismicity13,14,21, and the degree to which this occurs may vary along the fault22. Geodynamic models suggest that dilatancy mechanisms, whereby fluids become more viscous under pressure may be in effect periodically23,24. Fluids thus play a primary role. However, their exact pathways and relationship with observed earthquakes are not well-known.

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)25–30. The seafloor of Chain Transform Fault 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 is slipping at a rate of ~30 mm/yr.

Results
Historical seismicity
25 M ≥ 5.6 historical earthquakes along the Chain Transform Fault are reported in the relocated catalog31 which extends back to 1993 (Fig. 1c). Empirical relationships12 suggest a fault rupture length of ~33 km for the largest (M = 6.6) recorded earthquakes, in agreement with the scale of two of the longest fault scarps25 that are located on the eastern side of WFS2 and between WFS3 and ELFS (~30 km) (Figs. 1 and 2). No large earthquakes have been recorded from the longest, 50 km, fault length located within ELFS since 1993 (green lines, Figs. 1 and 2)31. 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 (M ≥ 6.0) earthquakes have been recorded but are rather characterized by aseismic slip and/or seismicity swarms (Fig. 1, Supplementary information Fig. 1). B1 and B3 are located near the ends of the transform. Transform sections near the ridge, e.g., within 20 km31 do not typically host the largest earthquakes given the rapidly changing stress regime32. However, B1 and B3 extend to much further distances from the ridge (>50 km) and include substantial length scales within regions predicted to be characterized by more uniform stress regimes32. 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 (WFS1, WFS3, and ELFS), whereas the asperity regions are co-located with WFS2 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, d). 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).

Local seismicity recorded by the OBS
We detected 626 local events (370 above the magnitude of completeness, M = 2.3) with locatable epicenters33,34 (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. 3, Supplementary Fig. 2). 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 component (45° < |rake| < 135°), 2 of which were normal and 11 of which were reverse, and 11 of these are located directly beneath or within 5 km of flower structures (Fig. 1a, Methods). Stress inversion results suggest that simultaneous activation of both strike-slip and reverse faults is possible (Methods, Supplementary Fig. 3). Elevated b-values are found in regions of shallow bathymetry and low rMBA, i.e., characteristics common to the flower structures (Figs. 1 and 4, Methods; Supplementary Movie 1). The length scale of heterogeneity in seismicity is 10s of km, which is much greater than the scale of bathymetric resolution (100 m), allowing us to make meaningful comparisons.

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. 3a, Methods, Supplementary Fig. 4, Supplementary Fig. 5). We also compare our result to thermal structure predicted from the S-wave model derived from local Rayleigh and ambient noise tomography (Methods, Fig. 3c)35. We preferentially interpret the P-wave model at shallow depths beneath the ELFS, but also consider the S-wave model for deeper velocities and coverage along the entire transform (Methods).

Seismic velocities can be influenced by a variety of factors such as composition including degree of alteration/hydration, pore fluid, or temperature. Typically, the bulk composition of depleted mantle is expected to be relatively homogeneous, but even large variations in depletion would yield very minor amounts of seismic velocity anomalies36. Hydrothermal alteration and serpentinization are expected to reduce seismic velocities. However, the S-wave velocities we observe in the lithosphere are relatively high, 4.55 km/s on average, and inconsistent with these factors37. Also, porosity effects on S-wave speeds are not likely to be important at the mantle depths where surface waves have sensitivity (>~20 km depth)38. Therefore, we proceed interpreting the S-wave velocities in terms of temperature. Observed slow P-wave velocities could be explained by thickened crust, enhanced porosity, and/or hydrothermal alteration, which we will discuss in greater detail in the next section.

Methods
Origin of the lozenge-shaped topographic highs
One possibility is that these features represent intra-transform volcanism44. However, the observed low backscatter intensity of these structures25 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. This is also supported by the 11 reverse fault mechanisms recorded by the OBS along the transform (Fig. 1). There were only 2 events characterized by a dominantly normal focal mechanism, potentially explained by flank collapses. 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 transform fault that are also thought to have a transpressional origin42.

Seismic catalog and focal mechanisms
The analysis of the broadband ocean bottom seismic deployment has revealed 972 events, 812 of which are located within the area surrounded by the OBS network33, which belongs to the Chain transform fault and the adjacent ridge spreading centers (Fig. 1b). The location of the events is performed by NonLinLoc software45 using a 1-D velocity model of the Chain transform fault region based on CRUST1.046 and earthquake arrival times47. In order to focus on the transform fault 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, with a local magnitude (M) range between 1.1 and 5.4. We further remove events below the completeness threshold M = 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 software48, 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 6 km).

We supplement our seismicity data recorded by the OBS with historical earthquakes from the Global Centroid Moment Tensor catalog49. Due to the relatively large epicentral uncertainties, we only consider the strong events (M ≥ 5.6) after 1993 relocated by Shi et al.31.

OBS seismicity with respect to bathymetry and rMBA
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 epicenter location of events with the background rMBA/bathymetry data. Background data are defined as all the rMBA/bathymetry values calculated within a 15 km wide zone, centered on the transform valley axis and bounded in the along-strike direction by our seismicity data (roughly between the ridge spreading centers). 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 epicenters show no preference for particular parameter values, then the red and blue histograms in Fig. 4d and Fig. 4e 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. 4d). There is also a preference of events occurring at very shallow ocean depths (~<3700 m, Fig. 4e), i.e., directly beneath the ELFS. We perform two statistical tests, to quantify the difference between the two distributions. The Wilcoxon rank sum test is used to test the null hypothesis that the data in the two sets (i.e., epicenters and background data) are drawn from distributions with equal medians. In addition, the two sample Kolmogorov-Smirnov test was applied to test the hypothesis that the data in the two sets (i.e., epicenters and background data) are drawn from the same continuous distributions. Both Wilcoxon rank sum test and Kolmogorov-Smirnov test reject the corresponding null hypotheses: at p < 10−6 for rMBA and at p < 10−4 for bathymetry.

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 neighbor 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 Supplementary Movie 1).

The exponentiality of magnitude distribution is investigated by the Anderson-Darling Test50,51 and further established by goodness of fit test52,53, applying the Maximum Likelihood Estimation (MLE) of b-value54. Both techniques suggest that the magnitude distribution can be sufficiently modeled by an exponential distribution for M ≥ 2.3. A total of 370 events with M ≥ 2.3 comprise the complete data for the study area (circles, Figs. 1b and 2), which roughly corresponds to an average rate of ~1 event/day. However, the MLE returns considerably unstable b-values, even well above M, fluctuating between 0.76 and 0.89 for 2.3 ≤ M ≤ 3.0 (Supplementary information Table 1, Supplementary 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 M selection55,56. The comparison of the results from the aforementioned methods, assuming different M are shown in Supplementary Fig. 7.

We choose to estimate b-values by applying the repeated medians57 (RM), a non-parametric regression technique, which provides more robust results for small datasets, since it is highly resistant to observational uncertainties and outliers55. According to the RM approach, we consider n magnitude intervals between two points, i and j, having N and N events, respectively, and calculate n-1 slopes (b-values) between these points as:

b ij = log 10 (N j) - log 10 (N i) / M j - M i (1)

With M j ≠ M i. 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(b ij))

The standard errors, SE, of the corresponding bRM are derived by bootstrap resampling process as described in Amorèse et al.55. 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:

t = |bRM(A) - bRM(B)| / sqrt(SE2 A + SE2 B) (2)

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 Supplementary Fig. 7 and Supplementary 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 M choice.

The 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 (Supplementary Movie 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. 4a). The apparent cyclic relationship between the b-values and bathymetry is an artefact related to the low number of events in each calculation (100). b-values calculated using a small number of data are more sensitive to the larger magnitude events, which coincide with the x-axis locations where jumps in b-value are apparent (Fig. 4a). There are higher b-values (b = 1.04 ± 0.13) in regions of low rMBA (< −3 mgal) than regions of high rMBA (> −3 mgal) where b = 0.72 ± 0.11 (Fig. 4b, 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. 4d). These areas are mainly located in the regions of shallowest bathymetry (< 3700 m), i.e., beneath the ELFS (Fig. 1; Fig. 4e).

Seismic vs predicted moment
We compare the seismic moment release, M 0, in each asperity (A1, A2) and barrier (B1, B2, B3) using the formula58 M 0 = 10 1.5 M W + 9.105, where M 0 is given in N∙m. We express the results in terms of percentage (α) of seismically released moment (ΣM 0s) over the potential moment release assuming full seismic coupling (ΣM 0p), i.e., α = ΣM 0s / ΣM 0p. In doing so, we estimate ΣM 0p by summing the seismic moment of all M ≥ 5.6 events occurred since 1993 in each patch (barrier or asperity). We then estimate M 0p as5:

M 0p = G L w s ΔT (3)

where, G is the shear modulus [40–60 GPa], L is the patch length, w, is the patch width from seafloor to the 600°C isotherm, s is the slip rate [28–33 mm/yr] and ΔT is the time period duration since 1993, equal to 28 years. Our results range from α = 0.25 to α = 0.46 with an average value of α = 0.32 (for G = 50 GPa and s = 32 mm/yr). This range agrees well with previous findings (α = 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 to 0.06 in barrier regions to 0.67 to 0.78 in asperity regions (Fig. 1c, d; Supplementary Fig. 1).

Stress inversion
The applied stress inversion algorithm59,60 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 M > 5.0 events, with the stronger earthquake being M = 5.4, the effects of co-seismic stress changes can be neglected61,62 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 results63 mainly the shape ratio determination. To avoid this effect, an algorithm which inverts jointly for stress and fault orientations is considered64. The fault orientation determination is based on the evaluation of fault instability index, I:

I = τ + μ(σ - 1) / sqrt(μ + 1 + μ2) (4)

where τ and σ are the normalized shear and normal tractions, respectively and μ, 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 inversion59, 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 (Supplementary Fig. 3a), 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 solutions49 suggest strike-slip faulting along Chain Transform Fault. However, our 1-year seismicity analysis of the OBS data reveals that a proportion of smaller magnitude events (3.0 < M < 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° (Supplementary Fig. 3b). 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, σ 2 and minimum, σ 3 principal stresses, both deviate considerably from the vertical direction (plunge is 60° for σ 2 and 30° for σ 3) and they have similar amplitudes (Supplementary Fig. 3a, c). Therefore, they are virtually indistinguishable and might switch with each other. With σ 1 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 plausible64. 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 tomography35. Group velocity measurements were inverted from 16 to 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 center of the array, with many crossing ray paths 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 to 30 km beneath the seafloor and the average value from 30 to 50 km beneath the seafloor. Therefore, we do not interpret S-wave velocity model variations at any particular depth, rather in an average sense, over some depth range. The lateral and vertical resolution suggests that the variation in lithospheric thickness/velocity we observe across the transform fault at ~100 km scale (WFS1-WFS3, 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. 3). 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 center of the transform, as predicted by geodynamic modeling 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 mantle47:

V s = 4.77 + 0.038 P - 0.00038 (T - 300) (5)

where V s 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 Hacker65 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 structure66. 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 Lithgow-Bertelloni47 parameterization, particularly the 600°C isotherms, although the effect is not substantial. With the Abers and Hacker65 parameterization the 700°C isotherm lies at nearly the same depth as the 600°C isotherm from the Stixrude and Lithgow-Bertelloni47 parameterization. The difference between the compositions from the Abers and Hacker65 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 experiment67 (Supplementary Data 1; Supplementary Figs. 4–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 (Supplementary Fig. 4). 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 approach68. 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 elsewhere40. 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 (Supplementary Fig. 5).

Our P-wave refraction study illuminates anomalous crust and mantle structure beneath ELFS. There is a clear Moho reflector visible across the region, which is relatively flat at 10 km depth, especially beneath ELFS (Supplementary Fig. 6, Fig. 4a). Within the crust of ELFS there is evidence of a high velocity core, while the crustal structure away from ELFS has relatively flat velocity contours (Fig. 4a). In the uppermost mantle beneath the Moho reflector we observe a broad region of slower than expected mantle velocities (<7.5 km/s), i.e., 4–10% slower in comparison to velocities >7.8 km/s on the western side of the ELFS. The slow velocity anomaly extends to at least 15 km depth and is centered just east of ELFS.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper describes a study that investigates the characteristics of seismicity that occurs in the\nChain Fracture Zone, an oceanic transform fault. In particular, it looks at how the magnitude\ndistribution changes spatially with characteristics such as morphological features identified in\nbathymetry and gravity data. It finds that the Chain can be divided into 2 asperity and 3 barrier\nzones. The barrier zones are characterized by smaller moment release than expected and higher\nb-values (i.e., relatively more smaller earthquakes than larger ones). They also tend to occur\nwhere the flower structures and therefore broader fault zones are located, which suggests that\nthese features facilitate fluid transport deeper into the lithosphere and weakening it, thus reducing\nthe potential for large earthquakes to occur.\nI very much enjoyed reading this paper; I thought it was well-written and well-reasoned, and easy\nto follow. I appreciated the thoroughness with which the b-values were investigated, especially\nconsidering multiple methods. It is certainly on appropriate topics for this journal, combining\nearthquake statistics with seafloor geophysics observations and modeling, and I think a broad\nmultidisciplinary range of readers will appreciate it. What controls earthquake rupture size is of\ncourse a topic of great interest with implications for seismic hazard, and this study could be a\nvaluable contribution to that topic, by examining it in the relatively simpler oceanic transform fault\nenvironment. I have just a few minor comments and suggestions listed below.\np. 1, 2nd paragraph: What magnitude would the rupture of the entire transform fault length\npotentially result in? (to compare it numerically with M7)\np. 3, Fig. 1 caption: The caption for panel (b) states that black dots are in low rMBA and white dots\nare in high rMBA, but the figure itself shows the reverse. I think the caption might have just\naccidentally switched them.\np 4, Fig. 2: I think it could help to label the isotherms that are shown in panels (b) and (c). I’m not\nexactly sure which one is the 600C isotherm, but I suppose I could guess from the location of the\nred dots. Also, the caption states that “the events located deeper than the 600C isotherm from\nsimple thermal models are highlighted with a red dot.” Is this the thermal model shown in panels\n(b) and (c)? Finally, the caption for (b) states that “the dashed red isolines indicate the thermal\nstructure…” but I think these are actually shown in blue, correct?\np. 11, section on OBS seismicity (and Fig. 3d and e): The red and blue distributions in these\npanels do look qualitatively different. Have you done any statistical testing to show that they are\nquantitatively significantly different?\nExtended Data Fig. 7: It might be informative to indicate in this figure where the nominal Mc that\nyou use in this study lies (Mc=2.3), perhaps with a vertical dashed line.\np. 12: I believe that the numerator for the equation for bij should actually be “log10(Nj)-log10(Ni)”", "reviewer_2": "The authors present a seismological study of the Chain transform fault in the mid-Atlantic ocean\nridge. Combining several different datasets together, the authors argue that morphological\nfeatures of the fault zone, primarily flower structures linked to step overs in the fault system,\nfacilitate the circulation of fluids into the fault zone, pushing the brittle-ductile transition zone\ndeeper and acting as barriers to major seismic ruptures. The manuscript is well structured and\neasy to read (commendable given the amount of different datasets involved!) and the figures\ninformative. While the interpretation is not necessarily that surprising, the authors’ argument is\nwell justified from the unique combination of data types exploited here. The manuscript is thus of\nsure interest to the community.\nThe technical approach and analysis seem sound to me, as the methods are for the most part\nstandard techniques. One main comment would be on the proposed interpretation of the data that\nfails to address the question of fault heterogeneity. The prominence of microseismicity at the\napparent barriers to seismic slip would suggest to me that the distribution of heterogeneity length\nscales is significantly different between “asperities” (where large earthquakes have occurred) and\n“barriers”. The text makes reference to “mechanical and chemical effects…increasing pore\npressure” (L183) caused by the presence of fluids in the fault zone, but this vague explanation is\nlikely only part of the story. A reduction of effective normal stress would indeed favor seismic\nrupture, but how the rupture plays out after nucleation is controlled by heterogeneity (whether\nstructural or in the stress state). Given that we know that fluids (water) can alter lithologies, it\nwould be logical to me that it would also affect the distribution of fault heterogeneity, potentially\nleading to favoring a distributed pattern of seismic release (microseismicity) rather than large\nfault-wide ruptures.\nBeyond this point of discussion that I think should be addressed, I have a few other minor\ncomments to make.\nL23 - punctuated locations (barriers) -> localized barriers\nL52 - isn’t microseismicity brittle rupture?\nFig 1 - vertical alignment of axes labels is quite strange; should likely just rotate 90 degrees as is\ncommon\nFig 1 - It would be quite helpful to show a panel with cumulative seismic moment as a function of\nalong-strike position. This would clearly highlight where potential asperities and barriers are. This\nwould also complement well the seismic moment percentages\nL131 - red -> blue (?)\nFig 3a - why are there cyclic variations in the b-values as a function of bathymetric depth?\nFig 3b and c - Why not use the same color scheme for these two panels? It would make it clearer\nthat the annotations in b are the fits to the distributions in c", "reviewer_3": "This article outlines a series of observations that the authors propose support an evolving theory\nthat hydrothermal processes lead to earthquake segmentation of oceanic transform faults. Much of\nthe reasoning in the article leans heavily on recent work done at other fault systems, particularly\nthat there are distinct rupture patches on transform faults adjacent to rupture barriers, and that\nthis leads to the overall low seismic coupling seen at oceanic transofrms. New observations from\nthe Chain transform are presented, which is a valuable comparison to add to a growing number of\nregional studies using marine seismic and OBS data.\nThe authors suggest that the barrier regions can be associated on the Chain transform fault with\ntranspressional flower structures, interpreted in (undefined resolution) multibeam bathymetry\ndata. They also reference evidence for hydrothermal processes in these zones from a few\nobservations including 1. Several (8?) locatable deep earthquakes that occur within the barrier\nzones, 2. low RMBA gravity anomalies within the barrier zones, possibly associated with\nalteration?, and 3. seismic velocity models and estimated thermal structures from these models.\nAlthough I appreciate the new set of observations that the authors present in this paper, and I\nbelieve much of the logic they present is good for arguing for the link between earthquake rupture\nproperties and fluid-rock interactions, I do not feel that the data set shows the correlations they\noutline clearly as it is presented. My primary concerns stem from a few specific things: 1. I am\nconcerned that the multibeam bathymetry images are over interpreted, and that ‘flower structures’\nidentified by the authors may be course-resolution expressions of a more complex structural\nconfiguration that is accommodating variation in\nhydrogeology or some other source of material\nheterogeneity, 2. Inferences made regarding along-strike changes in thermal structure based on\nseismic velocities seem to be strangely presented and interpreted. Seismic velocity models shown\nin the paper do not seem to show striking low velocity zones where the authors claim there are\naltered zones. (Or I may just be having a difficult time making the needed correlation b/w seismic\nvelocity and rupture zones) 3. Geophysical trends in gravity from RMBA are presented as evidence\nof hydrothermal processes, but it’s not clear to me that they necessarily support this\ninterpretation. Forward modeling of how low RMBA could relate to (for example)\nserpentinization or porosity (not mentioned) as well as thickened crust (how does this\npromote hydrothermal cooling?) would be helpful to strengthen this interpretation and clarify how\nlow RMBA could link to material variability that results from the proposed hydrothermal processes.\nBelow I have outlined some specific sections of the manuscript that I found to be confusing or in\nneed of reconsideration:\n- Line 70: “The Chain TF contains four lozenge-shaped topographic highs and shows evidence for\nvariability in gravity anomalies and inferred crustal structure along the length of the TF8 (Fig 1a,\nb). These features represent crustal blocks exhumed via transpressional stresses, i.e., features\ntypically referred to as positive flower structures (Methods). “\nI am not certain, as the it is currently presented that the presence of these flower structure is very\ncompelling. Figure 1a shows bathymetry (displayed at an unlabeled ? resolution) with the\ntopographic highs interpreted. My skepticism comes from the seeming likelihood that these\nfeatures are being over-interpreted based on the apparently low-resolution dataset and the\nhorizontal earthquake uncertainties as well. What is the resolution of the seafloor morphology that\nis being presented in the bathymetry? Are we sure that topographical highs and earthquake\npatters are really correlating with the uncertainties of both datasets? Could there be finer scale\nfeatures (i.e. fault step-overs, multiple strands, as the authors point out, but also extensional\nbends or wide damage zones?) that are equally responsible for the observed heterogeneity in\nearthquake behavior along strike? It’s not clear how robust the morphology interpretations are,\nand it seems like they do not perfectly correlate with either the RMBA or the earthquake locations\nof more-frequent smaller events.\n- Line 192: “The S-wave velocity model predicts that the 600°C isotherm associated with the\nbrittle-ductile transition occurs deeper than predicted for a thermal model, supporting the\nhypothesis that deeper cooling has occurred. “\nIt seems like the authors may need to spend more time convincing us that the interpretation of S-\nwave velocity is a function of temperature, and not other physical properties. The S-wave velocity\nstructure does not predict that it is cooler in the barrier zone relative to the asperitiy zone, just\nthat the temperatures are cooler than predicted by a (very simple!) half space cooling model. More\nadvanced numerical models predict a cooler thermal structure both near the transform fault ends\n(flatter isotherms) and if hydrothermal cooling is present… this is not a surprise anywhere along a\ntransform fault, but doesn’t seem to correlate with earthquake barrier zones shown here, except\npossibly below 30 km, and well below the 900° C isotherm based on S-wave velocities (where we\nwould not expect brittle failure?)\n- Line 201: “In addition, the predicted temperatures from the S- wave velocity model are hotter at\n25 – 30 km depth than surrounding regions, “\nSimilar to comment above: This doesn’t’ seem to be true based on either Figure 2b or 2c? I don’t\nsee how the S-wave velocity model shows an anomaly in this range, or that there is a significant\nthermal anomaly (predicted from the S-wave velocity model) that is consistent with this. I\nwouldn’t expect that S-wave anomalies would necessarily be caused by changes in temperature as\nopposed to lithologic variability or variation in porosity.\n- Line 177: “Hydrothermal alteration of the shallow lithospheric mantle (<15 km beneath the sea\nsurface) may be suggested, given the anomalously slow sub-crustal P-wave velocities down to at\nleast 15 km”\nI am finding it difficult to interpret Figure 2a showing the P-wave velocity. Even as someone who\nmodels crustal velocity structure, it seems like this P-wave model is only a little bit slower than\ntypical crust at ~8 km below the surface (sub-crustal depths?) – with no reference, it’s difficult to\nbe sure of what the velocities mean in here. Also where is this model exactly relative to the Chain\nTF? Does it extend past the ridge – to the east of the fault shown in b and c panels? Although I\nthink it is useful to have a velocity model if the authors are arguing for the impact of hydrothermal\nalteration, in this form, I am not sure I gain much in the way of evidence for hydrothermal\nprocesses from this figure.\n- Line 180: “Alteration caused by fluid infiltration may also explain why the section containing the\nlongest (50 km) fault scarp in the region is not associated with the largest earthquakes, but\ninstead associated with elevated b-values (1.1 – 1.2) with very low amounts of energy released\nseismically (4 %).“\nThis is quite confusing. Are the authors saying here that hydrothermal alteration is prevalent\neverywhere along the fault – in the asperity and barrier zones? Where specifically are they\nreferring to when they reference the 50 km long fault scarp (is this within the AX or BX zones?).\n- Re Methods - S-wave model translated into thermal model:I am not an expert in using S-wave\nvelocities to estimate thermal structure, but it does seem to me that there could be physical\nfactors other than just temperature that are influencing S-wave velocity at relatively shallow\ndepths, such as the 5-20 km below the seafloor depths that are shown in the paper overlapping\nwith seismicity."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThis paper describes a study that investigates the characteristics of seismicity that occurs in the Chain\nFracture Zone, an oceanic transform fault. In particular, it looks at how the magnitude distribution\nchanges spatially with characteristics such as morphological features identified in bathymetry and\ngravity data. It finds that the Chain can be divided into 2 asperity and 3 barrier zones. The barrier zones\nare characterized by smaller moment release than expected and higher b-values (i.e., relatively more\nsmaller earthquakes than larger ones). They also tend to occur where the flower structures and therefore\nbroader fault zones are located, which suggests that these features facilitate fluid transport deeper into\nthe lithosphere and weakening it, thus reducing the potential for large earthquakes to occur.\nI very much enjoyed reading this paper; I thought it was well-written and well-reasoned, and easy\nto follow. I appreciated the thoroughness with which the b-values were investigated, especially\nconsidering multiple methods. It is certainly on appropriate topics for this journal, combining earthquake statistics\nwith seafloor geophysics observations and modeling, and I think a broad multidisciplinary range of\nreaders will appreciate it. What controls earthquake rupture size is of course a topic of great interest with\nimplications for seismic hazard, and this study could be a valuable contribution to that topic, by\nexamining it in the relatively simpler oceanic transform fault environment. I have just a few minor\ncomments and suggestions listed below.\np. 1, 2nd paragraph: What magnitude would the rupture of the entire transform fault length\npotentially result in? (to compare it numerically with M7)\nWe added the following sentence:\n“For instance, at the Chain TF the magnitude of rupturing the entire transform is predicted to be M > 8.013 based\non empirical relationships ”\nWith 13: Wells and Coppersmith, 1994\np. 3, Fig. 1 caption: The caption for panel (b) states that black dots are in low rMBA and white dots are in\nhigh rMBA, but the figure itself shows the reverse. I think the caption might have just accidentally\nswitched them.\nYes, fixed.\np 4, Fig. 2: I think it could help to label the isotherms that are shown in panels (b) and (c). I’m not exactly\nsure which one is the 600C isotherm, but I suppose I could guess from the location of the red dots.\nWe increased the size of the contour labels and framed them within blue boxes to make them more visible.\nAlso, the caption states that “the events located deeper than the 600C isotherm from simple thermal\nmodels are highlighted with a red dot.” Is this the thermal model shown in panels (b) and (c)?\nWe added the parenthetical reference in the following sentence in the figure caption:\n“The events located deeper than the 600°C isotherm from simple thermal models (dashed blue lines in panels b\nand c) are highlighted with a red dot.”\nFinally, the caption for (b) states that “the dashed red isolines indicate the thermal structure…” but I\nthink these are actually shown in blue, correct?\nWe changed ‘red dashed’ to ‘blue dashed.’\np. 11, section on OBS seismicity (and Fig. 3d and e): The red and blue distributions in these\npanels do look qualitatively different. Have you done any statistical testing to show that they are\nquantitatively significantly different?\nWe added the following to the text of the methods section:\n“We perform two statistical tests, to quantify the difference between the two distributions. The Wilcoxon rank sum\ntest is used to test the null hypothesis that the data in the two sets (i.e., epicentres and background data) are\ndrawn from distributions with equal medians. In addition, the two sample Kolmogorov-Smirnov test was applied to\ntest the hypothesis that the data in the two sets (i.e., epicentres and background data) are drawn from the same\ncontinuous distributions. Both Wilcoxon rank sum test and Kolmogorov-Smirnov test reject the corresponding null\nhypotheses: at p<10-6 for rMBA and at p<10-4 for bathymetry.”\nExtended Data Fig. 7: It might be informative to indicate in this figure where the nominal Mc that you use\nin this study lies (Mc=2.3), perhaps with a vertical dashed line.\nDone (now Supplementary material Fig. 8).\np. 12: I believe that the numerator for the equation for bij should actually be “log10(Nj)-log10(Ni)”\nCorrected.\nReviewer #2 (Remarks to the Author):\nThe authors present a seismological study of the Chain transform fault in the mid-Atlantic ocean ridge.\nCombining several different datasets together, the authors argue that morphological features of the fault\nzone, primarily flower structures linked to step overs in the fault system, facilitate the circulation of fluids\ninto the fault zone, pushing the brittle-ductile transition zone deeper and acting as barriers to major\nseismic ruptures. The manuscript is well structured and easy to read (commendable given the amount of\ndifferent datasets involved!) and the figures informative. While the interpretation is not necessarily that\nsurprising, the authors’ argument is well justified from the unique combination of data types exploited\nhere. The manuscript is thus of sure interest to the community.\nThe technical approach and analysis seem sound to me, as the methods are for the most part\nstandard techniques. One main comment would be on the proposed interpretation of the data that\nfails to address the question of fault heterogeneity. The prominence of microseismicity at the\napparent barriers to seismic slip would suggest to me that the distribution of heterogeneity length\nscales is significantly different between “asperities” (where large earthquakes have occurred) and\n“barriers”. The text makes reference to “mechanical and chemical effects…increasing pore\npressure” (L183) caused by the presence of fluids in the fault zone, but this vague explanation is\nlikely only part of the story. A reduction of effective normal stress would indeed favor seismic\nrupture, but how the rupture plays out after nucleation is controlled by heterogeneity (whether\nstructural or in the stress state). Given that we know that fluids (water) can alter lithologies, it\nwould be logical to\nme that it would also affect the distribution of fault heterogeneity, potentially leading to favoring a\ndistributed pattern of seismic release (microseismicity) rather than large fault-wide ruptures.\nWe revised the final paragraph to add more details: “The Chain TF provides a view into how segmentation on\noceanic TF can develop, which may in turn explain why oceanic TF do not generate larger earthquakes for their\nfault length. Multiple overlapping faults located on the flower structures are associated with broader damage zones,\nwhich provide a wider pathway for fluids into the crust and uppermost mantle. The presence of fluids causes a\nseries of mechanical and chemical effects, potentially increasing the pore pressure, reducing the effective normal\nstresses, lubricating the fault, and/or creating weaker regions of altered material. The hydrothermal circulation that\noccurs at shallow depths also cools the underlying mantle deepening the brittle ductile transition. However, the\ndamaged/altered crust and upper mantle are weaker, permitting only low magnitude earthquakes and/or aseismic\naccommodation of tectonic strain accumulation in these barrier regions beneath the flower structures. The barrier\nregions separate stronger, locked regions (asperities) beneath singular, linear fault strands with narrower damage\nzones, where hydrothermal circulation is lower, the brittle ductile transition and observed earthquake depths are\nshallower. The faults in the asperity regions are stronger, permitting larger events. In other words, broader fault\nzones can effectively segment the fault.”\nBeyond this point of discussion that I think should be addressed, I have a few other minor comments to\nmake.\nL23 - punctuated locations (barriers) -> localized barriers\nFixed\nL52 - isn’t microseismicity brittle rupture?\nChanged to ‘rather than large scale ruptures’\nFig 1 - vertical alignment of axes labels is quite strange; should likely just rotate 90 degrees as is\ncommon\nWe added the following explanation to the figure caption: We use an oblique Mercator projection so that the\ntransform appears horizontal, which is the reason for angled ticks marking longitude and latitude.\nFig 1 - It would be quite helpful to show a panel with cumulative seismic moment as a function of along-\nstrike position. This would clearly highlight where potential asperities and barriers are. This\nwould also complement well the seismic moment percentages\nDone. We added an extra panel (d) in figure 1, showing the along-strike cumulative seismic moment of the\nhistorical events.\nL131 - red -> blue (?)\nFixed\nFig 3a - why are there cyclic variations in the b-values as a function of bathymetric depth?\nThe b-value drops are caused by the stronger events, which makes sense since the windows have only 100 events,\nmaking the calculations sensitive to the existence of strong events. What really matters is the general decreasing\ntrend of the b-value. We added the following to the text:\n“The apparent cyclic relationship between the b-values and bathymetry is an artefact related to the low number of\nevents in each calculation (100). b-values calculated using a small number of data are more sensitive to the\nlarger magnitude events, which coincide with the x-axis locations where jumps in b-value are apparent (Fig. 3a).”\nFig 3b and c - Why not use the same color scheme for these two panels? It would make it clearer that the\nannotations in b are the fits to the distributions in c\nFixed\nReviewer #3 (Remarks to the Author):\nThis article outlines a series of observations that the authors propose support an evolving theory that\nhydrothermal processes lead to earthquake segmentation of oceanic transform faults. Much of the\nreasoning in the article leans heavily on recent work done at other fault systems, particularly that there\nare distinct rupture patches on transform faults adjacent to rupture barriers, and that this leads to the\noverall low seismic coupling seen at oceanic transofrms. New observations from the Chain transform are\npresented, which is a valuable comparison to add to a growing number of regional studies using marine\nseismic and OBS data.\nWe agree that our observations are complimentary to some previous work. However, our reasoning and\nconclusions do not rely on it, and we have altered the text to better clarify this. Much previous work relies on 1 or\n2 observational constraints, e.g., seismicity or seismicity and gravity, etc. We are showing new independent\nvelocity models that provide independent support for hydrothermal circulation and cooling beneath regions with\nwide fault deformation zone. This kind of association and process are not part of previous work that we are aware\nof. Our work also provides a possible explanation for the cause of the barrier regions, in other words that broader\nfault damage zones preferentially facilitate fluid transport into the lithosphere.\nWe reworded the following paragraph to increase clarity:\n“A few other studies have reported TF seismicity at deeper depths than the 600oC isotherm14-16. Geodynamic\nmodels predict a deepened 600oC isotherm if there is hydrothermal alteration, and this is true particularly near the\nends of the transform22. However, seismicity has also been reported at deeper depths than predictions including\nhydrothermal alteration15,22. Therefore, it has been suggested that this can be explained by higher strain rates at\ntransform faults38. Alternatively, different rheological flow laws may be in effect, and brittle and ductile deformation\nmay occur over a broad range of temperatures (300oC – 1000oC) and, therefore, a broad range of depths,\nowing to variable seawater infiltration and grain sizes along the fault23. Our result does not preclude the possibility that\neither or both of these two scenarios are in effect in some locations globally. However, what is different about our\nwork is that our S-wave velocity model offers an independent constraint that is not often available. It shows that\ndifferent flow laws and/or high strain rates are not required at Chain, since seismicity occurs shallower than the\n600(cid:31)C isotherm derived from the S-wave velocity model. In other words, our S-wave wave velocity model\ndemonstrates that hydrothermal circulation has cooled our region more efficiently and to deeper depths beneath\nthe flower structures with multiple fault strands and enhanced damage and the likelihood of large earthquakes,\nand this model is model is likely in effect elsewhere.”\nWe also added the following sentence:\n“Although previous studies have suggested that transforms may be segmented, and it has been hypothesized\nthat variable water infiltration into the lithosphere plays a role, our study provides new constraints on the Earth\nproperties that can result in such phenomena.”\nThe authors suggest that the barrier regions can be associated on the Chain transform fault with\ntranspressional flower structures, interpreted in (undefined resolution) multibeam bathymetry data. They\nalso reference evidence for hydrothermal processes in these zones from a few observations including 1.\nSeveral (8?) locatable deep earthquakes that occur within the barrier zones, 2. low RMBA gravity\nanomalies within the barrier zones, possibly associated with alteration?, and 3. seismic velocity models\nand estimated thermal structures from these models. Although I appreciate the new set of observations\nthat the authors present in this paper, and I believe much of the logic they present is good for arguing for\nthe link between earthquake rupture properties and fluid-rock interactions, I do not feel that the data set\nshows the correlations they outline clearly as it is presented. My primary concerns stem from a few\nspecific things: 1. I am concerned that the multibeam bathymetry\nimages are over interpreted, and that ‘flower structures’ identified by the authors may be course-\nresolution expressions of a more complex structural configuration that is accommodating variation in\nhydrogeology or some other source of material heterogeneity, 2. Inferences made regarding along-strike\nchanges in thermal structure based on seismic velocities seem to be strangely presented and\ninterpreted. Seismic velocity models shown in the paper do not seem to show striking low velocity zones\nwhere the authors claim there are altered zones. (Or I may just be having a difficult time making the\nneeded correlation b/w seismic velocity and rupture zones) 3. Geophysical trends in gravity from RMBA\nare presented as evidence of hydrothermal processes, but it’s not clear to me that they necessarily\nsupport this interpretation. Forward modeling of how low RMBA could relate to (for example)\nserpentinization or porosity (not mentioned) as well as thickened crust (how does this\npromote hydrothermal cooling?) would be helpful to strengthen this interpretation and clarify how\nlow RMBA could link to material variability that results from the proposed hydrothermal processes.\nWe thank the reviewer for this useful query. We have responded to the content of this paragraph in the more\ndetailed review sections below.\nBelow I have outlined some specific sections of the manuscript that I found to be confusing or in need of\nreconsideration:\n- Line 70: “The Chain TF contains four lozenge-shaped topographic highs and shows evidence for\nvariability in gravity anomalies and inferred crustal structure along the length of the TF8 (Fig 1a, b).\nThese features represent crustal blocks exhumed via transpressional stresses, i.e., features typically\nreferred to as positive flower structures (Methods). “\nI am not certain, as the it is currently presented that the presence of these flower structure is very\ncompelling. Figure 1a shows bathymetry (displayed at an unlabeled ? resolution) with the topographic\nhighs interpreted. My skepticism comes from the seeming likelihood that these features are being over-\ninterpreted based on the apparently low-resolution dataset and the horizontal earthquake uncertainties\nas well. What is the resolution of the seafloor morphology that is being presented in the bathymetry? Are\nwe sure that topographical highs and earthquake patters are really correlating with the uncertainties of\nboth datasets?\nWe added the horizontal resolution of the bathymetry (100 m) and vertical resolution of the bathymetry (~10s of\nm) of meters to the figure caption. These uncertainties are an order of magnitude smaller than the size of the\nstructures interpreted in the bathymetry. For instance, the width of the transform is ~15 km.\nThe argument for transpression is presented in the Supplementary material section. We added some more\ninformation, commenting on the earthquake focal mechanisms (itallics, below) in more detail:\n“Origin of the lozenge-shaped topographic highs. One possibility is that these features represent intra-\ntransform volcanism39. However, the observed low backscatter intensity of these structures8 does not support\nrecent volcanism. A more likely scenario is that these features represent crustal blocks exhumed via\ntranspressional stresses underlain by altered mantle, i.e., features typically referred to as positive flower structures.\nThis is also supported by the 11 reverse fault mechanisms recorded by the OBS along the transform (Fig. 1). There\nwere only 2 events characterized by a dominantly normal focal mechanism, potentially explained by flank collapses.\nSimilar morphological features are frequently observed and fault plane locations are sometimes available to\nsubstantiate a transpressional origin. This is supported by the highly tectonised mafic and ultra-mafic compositions\nof rocks from features with similar morphology and gravity signatures at the nearby St. Paul TF that are also thought\nto have a transpressional origin40.”\nCould there be finer scale features (i.e. fault step-overs, multiple strands, as the authors point out, but\nalso extensional bends or wide damage zones?) that are equally responsible for the observed\nheterogeneity in earthquake behavior along strike?\nYes, we agree that there are multiple fault strands visible in the bathymetry (Supplementary Material Fig. 1). They\nsupport our hypothesis that the flower structures are wide damage zones that contribute to the heterogeneity.\nWe added the following to the text, “There could be additional fault structures beneath the resolution of our\nbathymetry, but there is no reason to believe that they would be different or more important than the large-scale\nstructures on which we have based our arguments.”\nWe added the following sentence, “The length scale of heterogeneity in seismicity is 10s of kms, which is much\ngreater than the scale of bathymetric resolution (100 m), allowing us to make meaningful comparisons.”\nWe are interpreting the features at the scale that we can resolve. We see no evidence for transtension in the\nfocal mechanisms of the earthquakes along strike.\nIt’s not clear how robust the morphology interpretations are, and it seems like they do not perfectly\ncorrelate with either the rMBA or the earthquake locations of more-frequent smaller events.\nPlease see other responses that describe the resolution of the bathymetry.\n- Line 192: “The S-wave velocity model predicts that the 600°C isotherm associated with the brittle-\nductile transition occurs deeper than predicted for a thermal model, supporting the hypothesis that\ndeeper cooling has occurred. “\nIt seems like the authors may need to spend more time convincing us that the interpretation of S-wave\nvelocity is a function of temperature, and not other physical properties.\nWe added the following to the text:\n“Seismic velocities can be influenced by a variety of factors such as composition including degree of\nalteration/hydration, pore fluid, or temperature. Typically, the bulk composition of depleted mantle is expected to\nbe relatively homogeneous, but even large variations in depletion would yield very minor amounts of seismic\nvelocity anomalies33. Hydrothermal alteration and serpentinization are expected to reduce seismic velocities.\nHowever, the S-velocities we observe in the lithosphere are relatively high, 4.55 km/s on average, and inconsistent\nwith these factors34. Also, porosity effects on S-wave speeds are not likely to be important at the mantle depths\nwhere surface waves have sensitivity35. Therefore, we proceed interpreting the S-wave velocities in terms of\ntemperature. Observed slow P-wave velocities could be explained by thickened crust, enhanced porosity, and/or\nhydrothermal alteration, which we will discuss in greater detail in the next section.”\nThe S-wave velocity structure does not predict that it is cooler in the barrier zone relative to the asperitiy\nzone, just that the temperatures are cooler than predicted by a (very simple!) half space cooling model.\nMore advanced numerical models predict a cooler thermal structure both near the transform fault ends\n(flatter isotherms) and if hydrothermal cooling is present… this is not a surprise anywhere along a\ntransform fault, but doesn’t seem to correlate with earthquake barrier zones shown here, except\npossibly below 30 km, and well below the 900° C isotherm based on S-wave velocities (where we\nwould not expect brittle failure?)\nWe revised the following paragraph to improve clarity:\n“A few other studies have reported TF seismicity at deeper depths than the 600°C isotherm14-16. Geodynamic\nmodels predict a deepened 600°C isotherm if there is hydrothermal alteration, and this is true particularly near the\nends of the transform22. However, seismicity has also been reported at deeper depths than predictions including\nhydrothermal alteration15,22. Therefore, it has been suggested that this can be explained by higher strain rates at\ntransform faults38. Alternatively, different rheological flow laws may be in effect, and brittle and ductile deformation\nmay occur over a broad range of temperatures (300°C – 1000°C) and, therefore, a broad range of depths, owing\nto variable seawater infiltration and grain sizes along the fault23. Our result does not preclude the possibility that\neither or both of these two scenarios are in effect in some locations globally. However, what is different about our\nwork is that our S-wave velocity model offers an independent constraint that is not often available. It shows that\ndifferent flow laws and/or high strain rates are not required at Chain, since seismicity occurs shallower than the\n600°C isotherm derived from the S-wave velocity model. In other words, our S-wave wave velocity model\ndemonstrates that hydrothermal circulation has cooled our region more efficiently and to deeper depths beneath\nthe flower structures with multiple fault strands and enhanced damage, and this model is model is likely in effect\nelsewhere.”\nWe added the italicized sentence to the following paragraph in the methods section:\n“The formal resolution for the inversion for S-wave velocity with depth is 0.1 for 1 km thick layers in the upper 10\nkm beneath the seafloor and 20 km at greater depths. This suggests that the S-wave velocity inversion can\nuniquely resolve the average velocity over the upper 10 km of the model beneath the seafloor, the average\nvelocity from 10-30 km beneath the seafloor and the average value from 30-50 km beneath the seafloor.\nTherefore, we do not interpret shear velocity model variations at any particular depth, rather in an average sense,\nover some depth range.”\nGiven surface wave sensitivity it doesn’t matter whether the 600 or 900 isotherm dips lower. We added the\nfollowing in the discussion:\n“Indeed, surface wave-derived shear velocities indicate higher velocities on average over the shallow mantle\n(down to 40 km depth) beneath the flower structures (Fig. 2). The predicted thermal contour depths (e.g., 600 °C)\nand shapes (e.g., 900 °C) are much different than the predictions for a simple conductive cooling thermal model,\naccounting for the depth sensitivity of surface wave inversion (Methods).”\n- Line 201: “In addition, the predicted temperatures from the S- wave velocity model are hotter at 25 – 30\nkm depth than surrounding regions, “\nSimilar to comment above: This doesn’t’ seem to be true based on either Figure 2b or 2c? I don’t see how\nthe S-wave velocity model shows an anomaly in this range, or that there is a significant thermal anomaly\n(predicted from the S-wave velocity model) that is consistent with this. I wouldn’t expect that S-wave\nanomalies would necessarily be caused by changes in temperature as opposed to lithologic variability or\nvariation in porosity.\nThis is similar the previous response. We also added the italicized phrase below to clarify:\n“In addition, the predicted temperatures from the S-wave velocity model averaged over depth, given the\nsensitivity of surface waves (Methods), are hotter than surrounding regions which could be consistent with less\nhydrothermal cooling in this region (Fig. 2b).”\n- Line 177: “Hydrothermal alteration of the shallow lithospheric mantle (<15 km beneath the sea surface)\nmay be suggested, given the anomalously slow sub-crustal P-wave velocities down to at least 15 km”\nI am finding it difficult to interpret Figure 2a showing the P-wave velocity. Even as someone who models\ncrustal velocity structure, it seems like this P-wave model is only a little bit slower than typical crust at ~8\nkm below the surface (sub-crustal depths?) – with no reference, it’s difficult to be sure of what the\nvelocities mean in here. Also where is this model exactly relative to the Chain TF? Does it extend past the\nridge – to the east of the fault shown in b and c panels? Although I think it is useful to have a velocity\nmodel if the authors are arguing for the impact of hydrothermal alteration, in this form, I am not sure I\ngain much in the way of evidence for hydrothermal processes from this figure.\nWe have tempered our interpretation of the active source study and provided additional information to support\nour interpretation. We also added Supplementary material Figure 7, which shows the Moho reflections, which are\nrelatively simple and clear, illuminating a flat and simple structure. For this reason, we interpret the velocity\nvariations beneath the Moho as related to altered mantle material.\nWe added the following text to the methods: “Our P-wave refraction study illuminates anomalous crust\nand mantle structure beneath ELFS. There is a clear Moho reflector visible across the region, which is relatively\nflat at 10 km depth, especially beneath ELFS (Supplementary Material Fig 7, Fig 3a). Within the crust of ELFS\nthere is evidence of a high velocity core, while the crustal structure away from ELFS has relatively flat velocity\ncontours (Fig. 3 a). In the upper most mantle beneath the Moho reflector we observe a broad region of slower than\nexpected mantle velocities (< 7.5 km/s), i.e., 4 – 10 % slower in comparison to velocities >7.8 km/s on the western\nside of the ELFS. The slow velocity anomaly extends to at least 15 km depth and is centered just east of ELFS.”\nWe added the following text to the main text: “Observed slow P-wave could be explained by thickened\ncrust, enhanced porosity, and/or composition including hydrothermal alteration, which we will discuss in greater\ndetail in the next section.”\nand\n“Clear Moho reflections are visible across the region from the active source experiment, and these arrivals\nsupport relatively constant crustal thickness at 11 km beneath the sea surface (Supplementary material, Fig. ED7);\nhowever, anomalously slow P-wave velocities in the subcrustal lithosphere are present down to at least 15 km\ndepth beneath the sea surface, which are consistent with fluids and damaged/porous altered low density mantle\nmaterial.”\nSupplementary material Figure 5 shows the location of the active source study.\n- Line 180: “Alteration caused by fluid infiltration may also explain why the section containing the longest\n(50 km) fault scarp in the region is not associated with the largest earthquakes, but instead associated\nwith elevated b-values (1.1 – 1.2) with very low amounts of energy released seismically (4 %).“\nThis is quite confusing. Are the authors saying here that hydrothermal alteration is prevalent\neverywhere along the fault – in the asperity and barrier zones? Where specifically are they referring to when\nthey reference the 50 km long fault scarp (is this within the AX or BX zones?).\nWe amended the section to clarify which fault we are referring to:\n“Fluids may also reduce the frictional properties of the TF beneath ELFS and explain why the longest (50 km)\nfault scarp in the region (Fig. 1, green line on ELFS, Supplementary Material Fig. 1) is not associated with the\nlargest earthquakes, but instead associated with elevated b-values (1.1 – 1.2) with very low fractions of\nseismically released moment (4 %).”\n- Re Methods - S-wave model translated into thermal model:I am not an expert in using S-wave velocities\nto estimate thermal structure, but it does seem to me that there could be physical factors other than just\ntemperature that are influencing S-wave velocity at relatively shallow depths, such as the 5-20 km below\nthe seafloor depths that are shown in the paper overlapping with seismicity.\nPlease see previous response regarding the possible effects of Earth properties on S-wave velocities."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am satisfied with how my comments from the previous version of the manuscript have been\naddressed. In particular, I appreciate the additional statistical testing that was done, and the\nupdated labeling of the figures. I believe this manuscript is now acceptable for publication. I only\nnote one very minor thing that should be corrected: on l. 388, “<-3mgal” should be “>-3mgal)”.", "reviewer_2": "Thank you for the detailed response to both my comments and those of the other reviewers. My\nconcerns and questions have been addressed by the authors, and I now recommend this\nmanuscript for acceptance.\nMy last two remarks are:\n- rereading this manuscript for the first time since the original submission, there are quite a few\nacronyms that are difficult to quickly process. Any reduction of acronyms (such as the simple TF\nfor transform fault) would be greatly appreciated and would improve the readability\n- a tracked changes version of the manuscript would have made the understand what was done in\nrevision much easier...!", "reviewer_4": "This paper looks into the factors that controls the rupture length of large earthquakes on oceanic\ntransform faults. They focused on the Chain Fracture Zone in the mid Atlantic, where the PI-Lab\nproject collected OBS data. They presented correlation between several geophysical and geological\nobservations. They observed that morphological transpressional features is correlated with rupture\nbarriers and proposed that fluids going down through the flower structures and altered the fault.\nThe paper is well written. The data is mostly supporting the conclusions. I suggest publication with\nminor revisions. Below are my comments.\nLine 28: Why use TFs? OTFs is usually used for oceanic transform faults.\nLine 68: The four topographic highs are difficult to find in Figure 1.\nLine 198-200: This seems out of place. Or there is not enough information.\nLine 265: \"model is\" is duplicated.\nLine 543: The most critical part of this paper is the high resolution bathymetry data. It was not\nmentioned in the Data availability. I am not familiar with all the rules but shouldn't it be provided\nfor people to verify the geological structures in Figure 1a.\nFigure 1: The seismic segment and the flower structure do not aligned well. Part of the WFS2 is\nwithin A1. Isn't this contradict to the conclusion?\nFigure 2 is very hard to follow. (a) is showing a different section as (b) and (c). It should be\npointed out in the caption in the beginning and also provide an explanation. The P-wave model\nseems to be based on active source survey and only covers the east segment. This information is\nin the Method section but should be mentioned here. This will clarify things and help readers to\nunderstand the figure.\nAn additional supplemental figure in map view might be helpful. What is the horizontal location\nuncertainty of the OBS data? In a map view figure, do all the events follow the main fault trace?\nOr do any occur on the side faults in the flower structure?\nFaults do not operate in this simple way. I suggest to make it clear that the main conclusion is\nbased on Chain only and might not work for all. There are many oceanic transform faults that\ntransit from seismic to barrier without any surface structure change. Of course, some does such as\nGarrett."}, "author_response": "Response to Reviewers NCOMMS-22-23238B\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nI am satisfied with how my comments from the previous version of the manuscript have been\naddressed. In particular, I appreciate the additional statistical testing that was done, and the\nupdated labelling of the figures. I believe this manuscript is now acceptable for publication. I\nonly note one very minor thing that should be corrected: on l. 388, “<-3mgal” should be “>-\n3mgal)”.\nWe thank the reviewer for comments and corrections. The typo was corrected.\nReviewer #2 (Remarks to the Author):\nThank you for the detailed response to both my comments and those of the other reviewers.\nMy concerns and questions have been addressed by the authors, and I now recommend this\nmanuscript for acceptance.\nMy last two remarks are:\n- rereading this manuscript for the first time since the original submission, there are quite a\nfew acronyms that are difficult to quickly process. Any reduction of acronyms (such as the\nsimple TF for transform fault) would be greatly appreciated and would improve the\nreadability\n- a tracked changes version of the manuscript would have made the understand what was\ndone in revision much easier...!\nWe thank the reviewer for all comments and suggestions. We replaced TF with transform fault\nthroughout the manuscript. We decided however to keep the abbreviations for ELFS (East Large\nFlower Structure), WFS (western Flower Structures) and rMBA (residual Mantle Bouguer Anomaly),\nwhich are all shown in figure 1 and are explained in the text as well as in the figure captions.\nReviewer #4 (Remarks to the Author):\nThis paper looks into the factors that controls the rupture length of large earthquakes on\noceanic transform faults. They focused on the Chain Fracture Zone in the mid Atlantic, where\nthe PI-Lab project collected OBS data. They presented correlation between several\ngeophysical and geological observations. They observed that morphological transpressional\nfeatures is correlated with rupture barriers and proposed that fluids going down through the\nflower structures and altered the fault.\nThe paper is well written. The data is mostly supporting the conclusions. I suggest\npublication with minor revisions. Below are my comments.\nLine 28: Why use TFs? OTFs is usually used for oceanic transform faults.\nWe removed all abbreviation of TF from the manuscript.\nLine 68: The four topographic highs are difficult to find in Figure 1.\nIt was very hard to find an appropriate colour scale for the topographic highs to be visible, and at the\nsame time keep the rest of the regional features visible. For that reason, we highlight these areas\nwithin the yellow curves in Fig. 1. In addition, we moved previous Supplementary Figure 1 into the\nmain text (new Fig. 2), which provides a larger, zoomed in version in which the topography can be\nseen.\nLine 198-200: This seems out of place. Or there is not enough information.\nWe thank the reviewer for pointing this out. We removed this sentence from the manuscript.\nLine 265: \"model is\" is duplicated.\nDuplicate removed.\nLine 543: The most critical part of this paper is the high resolution bathymetry data. It was\nnot mentioned in the Data availability. I am not familiar with all the rules but shouldn't it be\nprovided for people to verify the geological structures in Figure 1a.\nWe added all resources in the Data Availability section.\nFigure 1: The seismic segment and the flower structure do not aligned well. Part of the WFS2\nis within A1. Isn't this contradict to the conclusion?\nThis is discussed in lines 258 – 271 of the revised (annotated) manuscript:\n“The WFSs observations share some similarities to those of ELFS, but also some differences.\nSimilarities include high topography, broad deformation zones, and relatively high S-wave velocities\nand likely cooler temperatures on average which could be consistent with larger degrees of\nhydrothermal circulation in the region. The rMBA in the WFSs region changes from negative to\npositive from east to west. The higher values are likely caused by normal to thinner crustal thickness,\nwhich likely dominates the gravity signature25. The large historical earthquakes beneath WFS2 also\ndistinguish it from the other flower structures. One explanation may be that the long fault segment in\nthe region extends 20 km east of the flower structure and may behave more similarly to the long fault\nin section A2. Another explanation could be that there are temporal variations in rheological behaviour\nas water is delivered and transported through the system. Episodic fluid delivery and transport and\nmulti-mode rupture models have been previously proposed to explain observed seismicity\npatterns5,7,40 and have also been predicted by geodynamic modelling24 and inferred from samples\nfrom exhumed oceanic transform faults41.”\nFigure 2 is very hard to follow. (a) is showing a different section as (b) and (c). It should be\npointed out in the caption in the beginning and also provide an explanation. The P-wave\nmodel seems to be based on active source survey and only covers the east segment. This\ninformation is in the Method section but should be mentioned here. This will clarify things\nand help readers to understand the figure.\nWe added reviewer’s recommendations in the caption of our new figure 2 (previously Supplementary\nFigure 1).\nAn additional supplemental figure in map view might be helpful. What is the horizontal\nlocation uncertainty of the OBS data? In a map view figure, do all the events follow the main\nfault trace? Or do any occur on the side faults in the flower structure?\nThis is now shown in the new Fig. 1 (previously Supplementary Figure 1), where we added the\nlocation of the OBS seismicity with the median and 95% quantile of the horizontal uncertainty.\nFaults do not operate in this simple way. I suggest to make it clear that the main conclusion is\nbased on Chain only and might not work for all. There are many oceanic transform faults that\ntransit from seismic to barrier without any surface structure change. Of course, some does\nsuch as Garrett.\nWe added the following phase (Lines 287-289 of the revised manuscript) to be more specific: “This\nalso likely occurs elsewhere, particularly in transpressional environments, for example, St. Paul and\nOwen Transform Faults43,44”. (Maia et al., 2016; Janin et al., 2023). We also changed some other\ninstances of “transforms” to “the transform”."}]

Decision letter