Upper crustal Vp/Vs ratios along the northern East Pacific Rise derived from downward-continued streamer data

SUMMARY

Multi-channel seismic (MCS) imaging has been extensively used to investigate fast-spreading East Pacific Rise (EPR) crustal compressional wave velocity (Vp) structure and tectono-magmatic behaviors. However, its upper oceanic crust’s shear wave velocity (Vs) profile has remained a rarity. We first confirm that additional offset ranges can be derived for traveltime picking from downward-continued MCS data in fast-spreading tectonic settings for both early-arrival P and S waves. We then inverse independent 2-D Vp and Vs structures along a ∼80-km-long along-axis stretch of the northern EPR. The resulting Vp/Vs ratio exhibit that the upper crust comprises pillow lavas, transition zone and sheeted dykes. The average thickness of pillow lavas is ∼125 m, with Vp increasing from ∼2.8 to 3.2 km s⁻¹ and Vs from ∼1.2 to 1.5 km s⁻¹. The lava unit with a transition zone has high Vp/Vs ratios (∼2.1 ± 0.2), indicating that fracturing and alteration are variable but pervasive. The average thickness of the transition zone is ∼400 m, with Vp increasing from ∼3.2 to 5.3 km s⁻¹ and Vs from ∼1.5 to 2.8 km s⁻¹. The pillow lavas and the transition zone constitute the layer 2A with an average thickness of ∼525 m. The boundary of layer 2A/2B can be defined using a Vp/Vs ratio contour of 1.9. The layer 2B exhibits lower Vp gradients (∼1.51 s⁻¹), Vs gradients (∼1.30 s⁻¹) and Vp/Vs ratios (∼1.8–1.9) compared to the layer 2A (∼4.65 s⁻¹, ∼2.98 s⁻¹ and ∼2.1 ± 0.2). Porosity variation and crack morphology are critical in controlling the seismic velocities of layer 2A. The strong lateral heterogeneity of the Vp/Vs ratios indicates hydrothermal signatures in the upper crust. The high Vp/Vs anomalies indicate fluid pathways into and out of the oceanic upper crust. This study demonstrates that the Vp/Vs ratio can be obtained from seismic tomography of downward-continued streamer data and used as a reference to investigate the crustal structure and hydrothermal activities along fast-spreading ridges.

1 INTRODUCTION

The East Pacific Rise (EPR) 9°–10° N is one of the most dynamic and magmatic mid-ocean ridges (MORs) with long-lived hydrothermal activities. Along fast-spreading MORs, the accretion of igneous oceanic crust is predominantly governed by magmatic processes, often resulting in a so-called Penrose model, which comprises a layered structure of pillow lavas (layer 2A), sheeted dykes (layer 2B) and layered/massive gabbros (layer 3). At the top of the igneous section, layer 2A typically exhibits relatively low compressional wave velocity (Vp ∼2.0 km s⁻¹) due to associated cracks and fissures, increasing to ∼4.5 km s⁻¹ as porosity decreases with plate divergence (Houtz 1976; Carlson 2011). Upper lavas are subhorizontal, undeformed or weakly fractured, whereas deeper regions are more fractured and altered (Karson et al. 2023). The transition zone is located between pillow lavas and sheeted dykes (Becker et al. 1989; Karson et al. 2002a, b, 2023), though this layer is not widely observed in seismic reflection data (Vera & Diebold 1994; Christeson et al. 2007). Layer 2B, characterized by higher Vp (∼4.5–6.0 km s⁻¹) values, overlies mid-crustal magma lenses at the ridge axis (Vera et al. 1990; Harding et al. 1993; Christeson et al. 1994; Estep et al. 2019). However, previous studies have primarily correlated the 2A/2B boundary with either the extrusive-to-sheeted dyke transition (e.g. Herron 1982; Toomey et al. 1990; Christeson et al. 1992, 1994, 1996; Harding et al. 1993; Vera & Diebold 1994; Hooft et al. 1996, 1997; Carbotte et al. 1997; Hussenoeder et al. 2002) or an alteration front (e.g. McClain et al. 1985; Becker et al. 1989; Burnett et al. 1989; Carlson & Herrick 1990; Fisher et al. 1990; Jacobson 1992; Wilcock et al. 1992; Christeson et al. 2007). Therefore, the depth, lateral distribution and seismic properties of the 2A/2B boundary need to be better constrained.

The shear wave velocity (Vs) and Vp/Vs ratios are highly sensitive to the physical state of rocks (Domenico 1984; Christensen & Mooney 1995; Christensen 1996), crack distribution (Collier & Singh 1998) and melt and fluid activity (Taylor & Singh 2002; Kim et al. 2019). In the young oceanic crust, particularly the upper oceanic crust, seismic velocities are primarily influenced by porosity structure and crack morphology (Carlson 2010, 2014a), whereas the impact of mineral alteration is comparatively minor. Measurement of rock samples from Ocean Drilling Program (ODP) Hole 504B indicated that alteration reduces Vp by only ∼0.2 km s⁻¹, whereas cracks decrease the average in situ Vp by ∼1.1 km s⁻¹ (Carlson 2014b). The correlation between Vp/Vs ratios and lithology is well established through both laboratory measurements (e.g. Domenico 1984; Christensen & Mooney 1995; Christensen 1996; Wang et al. 2009; Juneja & Endait 2017) and case studies (e.g. Christensen 1972; Christensen et al. 1974, 1980; Grevemeyer et al. 2018). These findings provide a solid foundation for petrographic interpretation based on seismic velocities, enabling Vp/Vs ratios to investigate the crustal structure and hydrothermal systems along MORs.

Velocity models derived from multi-channel seismic (MCS) data sets acquired over the past decades predominantly comprise Vp models obtained through traveltime tomography (e.g. Canales et al. 2008; Xu et al. 2009; Arnulf et al. 2011), 2-D and 3-D acoustic/elastic full waveform inversions (e.g. Arnulf et al. 2012, 2014a; Harding et al. 2016; Marjanović et al. 2017, 2019b; Xu et al. 2020; Zhang et al. 2022). However, obtaining and analysing converted shear waves with streamer-based seismic data is more challenging compared to wide-angle ocean bottom seismometer (OBS) data (e.g. Zhao et al. 2010; Grevemeyer et al. 2018; Li et al. 2021) due to the limitations of the marine environment, acquisition techniques, and severe factors (e.g. sediment thickness, vertical velocity gradient) in generating and receiving converted shear waves. Consequently, the characteristics of Vs and Vp/Vs structures remain weak in MCS studies.

In this paper, we first use synthetic modelling to demonstrate that our MCS data, after downward continuation (DC), should contain S waves amenable to traveltime tomography. We then integrate the DC method with traveltime tomography to obtain Vp/Vs ratios for an along-axis MCS profile at the northern segment of the East Pacific Rise (EPR). The Vs result reveals distinct features in layers 2A and 2B, suggesting a novel approach to delineate them. Correlations between modelled Vp/Vs variations, venting locations, ridge segmentation (White et al. 2006), numerical modelling (Marjanović et al. 2019a) and microseismicity (Waldhauser & Tolstoy 2011) help constrain the spatial distribution of fluid pathways, offering a new perspective on studying hydrothermal signatures along MORs.

2 BACKGROUND OF THE STUDY AREA

The EPR is an archetype for fast-spreading MORs with a spreading rate of ∼110 mm a⁻¹ (Carbotte & Macdonald 1992). Between 8° and 10°N, the EPR forms the boundary between the Pacific and Cocos Plates at 8° and 10°N (Fig. 1). The Ridge 2 K Program has selected it as an integrated study site because of the presence of all orders of ridge-axis discontinuities and a number of geological features, events and geophysical characteristics, such as continuous Clipperton transform fault (TF) and segmented Siqueiros TF (Gregg et al. 2009; Hebert & Montési 2011), 9°03′N overlapping spreading centre (OSC; Bazin et al. 2003; Arnulf et al. 2014b), the occurrence of repeated eruptions (Haymon et al. 1993; Tolstoy et al. 2006; Soule et al. 2007), low-velocity zones within the crust and upper mantle (Dunn et al. 2000; Zha et al. 2014), on-axis multilevel melt lens (Marjanović et al. 2014), off-axis magma bodies (Canales et al. 2012; Han et al. 2014; Aghaei et al. 2017), skewness of mantle supply (Toomey et al. 2007) and widely distributed hydrothermal activities (Haymon et al. 1991, 1993; Baker et al. 2016).

Before this study, most active-source seismic velocity modelling studies at the EPR focused on investigations of Vp structure using MCS (Marjanović et al. 2017, 2019b; Vithana et al. 2019) and OBS data (e.g. Orcutt et al. 1976; Grevemeyer et al. 1998; Canales et al. 2003; Hammond & Toomey 2003). Since the pioneering work of Vera et al. (1990), only a few studies have constrained Vs structure in this segment. Vera et al. (1990) used the reflectivity method to determine crust-scale Vs profiles at ∼9°N using traveltime calculation, waveform modelling, and the forward modelling of expanded spread profiles (ESP) data. The Vp/Vs structure near 9°30′N was constrained through amplitude modelling of seismic refraction data by placing both the source and receiver close to the seafloor (Christeson et al. 1994). Christeson et al. (1997) used 1-D waveform inversion of refracted arrivals from OBS data to model Poisson's ratio structure along the 9°–10°N stretch of the ridge. Seafloor compliance techniques were also utilized to provide constraints on Vp/Vs ratios at 9°48′N (Crawford et al. 1999), 9°33′N (Hulme et al. 2003) and 9°–10°N (Zha et al. 2014). However, in the above studies, the Vp/Vs ratios were either restricted to one specific site or obtained as average values throughout the entire crust. This study offers a 2-D heterogeneous upper crustal Vp/Vs ratios structure along the EPR at ∼9°–10°N.

3 SYNTHETIC TESTS OF DOWNWARD CONTINUATION WITH MCS DATA

3.1 Synthetic data generation

In recent years, the use of long and ultra-long (>12 km) streamers in scientific research has led to continuous improvements in velocity modelling (e.g. Qin & Singh 2018; Audhkhasi & Singh 2019; Kardell et al. 2021). Consistent with actual data acquisition, we used a 6-km-long streamer containing 468 receivers spaced at 12.5 m, with a 200 m distance between the first channel and the source. The shot interval was 37.5 m, totalling 368 shots. We utilized the general oceanic crustal architecture, that is the Penrose model, along the fast-spreading ridges (Fig. 2a) for data simulation. The structural model spanned 20 km in width and 8 km in depth, with 12.5 m grid spacing in both directions. The seafloor depth was fixed at 2 km, and the Vp of seawater was 1.5 km s⁻¹ (Fig. 2b). The Vs below the seafloor was converted from Vp using variable Vp/Vs ratios. The Vp/Vs ratio of layer 2A is 2.3256 (Vera et al. 1990), with a thickness of 0.5 km, and the ratio for the remainder of the crust is 1.7321 (Fig. 2c). Ccrust density distributions were calculated using Gardner's equation (Gardner et al. 1974; Fig. 2d). We used DENISE, an open-source 2-D time-domain isotropic viscoelastic finite difference modelling and full waveform inversion code for P/SV waves, to generate synthetic data (Köhn et al. 2012). All shots and receivers were placed at 10 m depth for simulation.

3.2 Downward continuation of the synthetic data

The simulated synthetic shot gathers (Figs 3a and c) reveal that the crustal P-wave refractions only appear at far offsets (>3.8 km). We used the Kirchhoff integral formulation of the wave equation (Berryhill 1979) to downward continue all shots and receivers directly to the seafloor (Fig. 3b). This method not only brings P-wave refractions as first arrivals in a wider offset range (e.g. Harding et al. 2007; Arnulf et al. 2012; Marjanović et al. 2017; Jimenez-Tejero et al. 2022) but also makes it possible to identify S-wave arrivals. To test this conjecture, we performed ray tracing with Vp and Vs models (Fig. 2) based on the after-DC geometry configuration, respectively. The synthetic arrival times of P and S refractions in the after-DC geometry accurately match the downward-continued synthetic surface data (Fig. 3d). Like P waves, the traveltime of S waves could be picked within a certain offset range (∼1–5 km), despite the interference of other phases and the boundary effect of the DC process (Zhang et al. 2022). This synthetic test confirmed that the MCS data after DC can be used for shear wave tomography in fast-spreading tectonic settings. Additional synthetic tests with simple one-gradient models are presented in the supporting information (Text S1, Figs S1 and S2).

4 DATA PROCESSING AND TRAVELTIME TOMOGRAPHY

The 3-D MCS experiment (Cruise MGL0812) was conducted in 2008 to image the magmatic-hydrothermal system at the EPR (Mutter et al. 2009). Four 6-km-long digital hydrophone streamers were towed behind the R/V Marcus G. Langseth to record air-gun shots. Each streamer contained 468 channels with a spacing of 12.5 m. The shot interval was 37.5 m. The sources and receivers were placed at 7.5 m below the sea surface. In this study, we extracted data from a ridge-parallel profile axis2r1 (Fig. 1) recorded along streamer 2 to investigate the Vp/Vs ratio structure along the ridge axis. The southern and northern ends mark the northernmost limit of the 9°03′N OSC, and the 2005–2006 seafloor eruption lava flows, respectively. All shots and receivers were downward continued to the seafloor with variable bathymetry for traveltime picking (Fig. 4).

Traveltime tomography of Vp and Vs structures was performed independently for the traveltime picks obtained from the DC data. Unlike P waves, the S waves always interfere with other phases (Figs 3d and 4), especially at the near offsets, so the S waves phase picking is the biggest challenge throughout this study (Fig. S3). A total of 152 406 P wave and 115 902 S-wave traveltimes were semi-automatically picked with manual correction for every fourth shot gathers (∼150 m). Here we only briefly report the Vs inversion process. All S-wave traveltimes were picked at 1.0–4.5 km offsets (Figs 4, S4 and S5). To avoid the bias from the Vp structure, we followed the initial model selecting a method of Canales et al. (2008). We chose a one-gradient model as our starting model, with a seafloor velocity of 1.25 km s⁻¹ and a vertical gradient of 2.50 s⁻¹ (Fig. S6). The open-source traveltime tomography code FAST (Zelt & Barton 1998) was used to obtain Vs model. The traveltime calculation and ray tracing were solved in regular grids with 25 m, and tomography was conducted using a regularized, damped least squares inversion on a 50 m grid. The final model was chosen when the weighted misfit function (chi-square) was closest to 1 (Figs S7 and S8). Further details of traveltime picking, estimation of traveltime picking uncertainty, the Vp inversion process, the uncertainty of the Vs model, and checkerboard resolution tests can be found in the supporting information (Text S2, Figs S9–S13).

5 RESULTS

Fig. 5 displays the final inverted Vp, Vs and Vp/Vs ratios models masked by ray path coverage. The final Vp model, derived from first-arrival traveltime tomography, generally aligns with the model obtained through wave-equation tomography (Marjanović et al. 2017). Up to approximately 200 m below seafloor (mbsf), the final Vs model appears more homogenous than the Vp model, which exhibits low-velocity sections (e.g. at model distances of ∼5–13, 43–46 and 65–77 km) and notable along-axis variations. At depths exceeding approximately 500 mbsf, the Vp and Vs models reveal alternating high and low-velocity structures. The Vp/Vs ratios demonstrate significant horizontal and vertical variations. Checkerboard resolution tests of Vp and Vs models suggest that features as small as 1.0 km × 0.5 km are well resolved at a depth less than 1.0 km, and 0.5 km × 0.25 km are resolved at a depth less than 0.5 km (Figs S12 and S13).

The along-axis upper-crustal averaged Vp structure consists of three primary velocity gradient layers: ∼2.91 s⁻¹ (∼125 m thick), ∼5.25 s⁻¹ (∼400 m thick) and ∼1.51 s⁻¹, from top to bottom. These layers are interpreted as the pillow lavas, transition zone, and sheeted dykes, respectively (Fig. 6a). The Vp values are ∼2.8 and ∼3.2 km s⁻¹ at the top and bottom of the pillow lava layer, whose thickness is consistent with numerous studies (e.g. Han et al. 2014; Marjanović et al. 2019b). The averaged Vs model also exhibits three distinct velocity gradient layers of 2.18, 3.25 and 1.30 s⁻¹, respectively (Fig. 6b), with values of ∼1.2 and ∼1.5 km s⁻¹ at the top and bottom of the pillow lava layer. The transition zone in both averaged Vp and Vs models display the most significant velocity gradients, reflecting extrusive-to-sheeted dyke transition and/or possible comprehensive fracturing/alteration in this layer. In the Vp/Vs ratio model, pillow lavas, transition zone, and sheeted dykes exhibit entirely distinct characteristics (Fig. 6c). The Vp/Vs ratios of layer 2A vary within a higher and broader range (∼1.9–2.2) and display an increasing trend except for the very shallow part, while the layer 2B concentrates within a narrow range (∼1.8–1.9), and shows a decreasing trend. An increase occurs between the pillow lava and the top of the transition zone, which then decreases gradually. Moreover, the most prominent features of the Vp/Vs ratios are the strong heterogeneities along the ridge axis (Fig. 5c). This phenomenon may result from magma intrusion, fracturing and fluid alteration.

6 DISCUSSION

6.1 Comparison with previous results

We present here, for the first time, 2-D upper crustal Vs and Vp/Vs ratios structures along the northern EPR, derived from tomographic inversion of streamer data (Fig. 5). The application of the DC technique enables better resolution of shallow subseafloor structures. Although the final models exhibit similar characteristics as previous studies (Figs 6a and b), all published velocity models display a wide range of values.

Discrepancies among earlier studies can be attributed to the varying data sets and limitations of modelling methods. In addition to the pioneering work of Vera et al. (1990) and Christeson et al. (1997) utilized shear waves recorded by OBS to model the regional Vs structure of the 9°–10°N segment. Still, the shallow structure (depth <0.5 km) could not be well resolved. Hulme et al. (2003) jointly inverted seismic and compliance data at 9°33′N, discovering that Vs from seismic data were higher than from seafloor compliance data due to anelasticity and anisotropy. Zha et al. (2014) inferred Vs profiles and melt distributions at 9°–10°N from 3-D seafloor compliance modelling. Generally, in the upper section (depth <1.0 km), our Vs model exhibits similar features to those of previous studies, including velocity drops (Fig. 6b). The deeper part (depth >1.0 km) displays higher values than those reported in earlier studies, which may result from the influence of the starting model and limited data coverage in depth.

6.2 Upper crustal structure demarcation from Vp/Vs ratios

Ocean drilling is the most direct and effective method for investigating the upper crustal structure. However, due to significant technical and programmatic challenges associated with scientific drilling, only a few holes, such as ODP Holes 504B (e.g. Becker et al. 1989), 735B (e.g. Bach et al. 2001), 1256D (e.g. Wilson et al. 2006) and IODP (International Ocean Discovery Program) Hole U1309D (e.g. Beard et al. 2009), have been successfully cored deeper than 1 km into the oceanic basement. Two of these holes (504B and 1256D), located at fast-spreading ridges, detected a transition zone between layers of pillow lavas and sheeted dykes (Becker et al. 1989; Wilson et al. 2006), suggesting that the transition zone may be widespread.

Before this study, the most prevalent method of identifying layer 2A involved picking the pseudo-reflection phase generated by the strong velocity gradient at its base (e.g. Han et al. 2014; Marjanović et al. 2018). However, these seismic depth images were often migrated using a simple 1-D Vp profile (i.e. ESP profiles; Vera et al. 1990) or velocity models with a limited shallow resolution, leading to overstatement or underestimation of layer 2A thicknesses. Previous studies on layer 2A thickness at EPR 9°–10°N reported 0.1–0.3 km at the ridge axis and 0.3–0.6 km at the ridge flank (Vera et al. 1990; Harding et al. 1993; Christeson et al. 1994; Aghaei et al. 2017; Marjanović et al. 2018). These works roughly divided the upper crust into pillow lavas and sheeted dykes without discerning the transition zone. Considering only the pillow lavas as layer 2A, our thickness result (∼125 m) is consistent with the aforementioned seismic studies. However, suppose we regard the pillow lavas combined with the transition zone as layer 2A. In that case, our result indicates a thickness of ∼525 m for layer 2A, which is thinner than Hole 504B (∼780 m) in the intermediate-spreading crust (Becker et al. 1989).

The nature of the layer 2A/2B boundary has been debated. Previous studies have interpreted it as a strong velocity gradient marking the transition from lavas to dykes (petrological boundary; e.g. Harding et al. 1993; Vera & Diebold 1994; Smallwood & White 1998), although this is not consistently observed (Christeson et al. 2007). At EPR 14°S, simultaneous 1-D inversion of Vp and Vs enabled assessment of crustal porosity (Collier & Singh 1998), which decreases from >30 per cent in layer 2A to 6–7 per cent in the layer 2A/2B transition and 5 per cent within layer 2B. This change appears consistent with a transition from lavas to dykes and the closure under pressure or infilling by sediment or hydrothermal minerals of numerous gaps and thin cracks within a lava pile (e.g. Carlson 2011). Christeson et al. (2007) argued that the bottom of seismic layer 2A corresponds to an alteration boundary, which can be either within the pillow lavas or near the top of the sheeted dykes. In our Vp and Vs models (Figs 6a and b), the upper crust consists of three distinct velocity gradient layers, which we interpret as pillow lavas, transition zone, and sheeted dykes, respectively. The pillow lavas, along with the transition zone, likely constitute the entire layer 2A, suggesting that the seismic layer 2A identified from the pseudo-reflection phase (e.g. Marjanović et al. 2018) cannot be used reliably to distinguish lavas and dykes at the ridge axis.

The average Vp/Vs ratio for layer 2A is estimated to be 2.1 ± 0.2, as illustrated in Fig. 6(c). This value is slightly higher than the laboratory measurements for metabasalts, which typically do not exceed 1.9 (Christensen 1996; Wang et al. 2009; Juneja & Endait 2017). Furthermore, it surpasses the ratios observed in oceanic rock samples (Fig. 6d; Christensen et al. 1974, 1980). Nonetheless, the estimated ratio aligns with the findings of Vera et al. (1990) at ∼2.36. However, it remains lower than the elevated ratios ranging from 3.6 to 5.0 identified by Christeson et al. (1997) and Collier & Singh (1998). The Vp/Vs ratios measured on Greenschist facies basalts collected from oceanic drilling were higher than those of laboratory samples containing the same minerals (Christensen 1996). High Vp/Vs ratios are often attributed to considerable porosity, hydration, and an extensive population of cracks (Contreras-Reyes et al. 2008; Moscoso & Grevemeyer 2015). Brantut & David (2019) used differential effective medium theory to confirm that the Vp/Vs ratios increase with higher fluid-saturated porosity. Additionally, permeability plays a significant role in Vp/Vs evolution. Numerous velocity-permeability relationships (e.g. Fabricius et al. 2007) demonstrate a positive correlation between the Vp/Vs ratios and permeability, with high Vp/Vs ratios resulting from elevated pore fluid pressure and crack anisotropy (Wang et al. 2012). Consequently, we hypothesize that low rock consolidation and high fluidity primarily account for the higher Vp/Vs ratios within layer 2A.

The Vp/Vs ratios of layer 2B are smaller than 1.9 (Fig. 6c) and exhibit characteristics similar to those found in previous studies (Fig. 6d) (Christensen et al. 1974, 1980), corresponding to unaltered or minimally altered basalt (Vp/Vs = 1.5–1.9; Christensen 1996; Grevemeyer et al. 2018). Although both lavas and dykes consist of basaltic material, dykes should have higher seismic velocities and lower Vp/Vs ratios due to their reduced fracture density and volume of void spaces. Our findings establish three distinct linear relationships (layers 2A, 2B, 2A + 2B) between Vp and Vs (Fig. 6d) and confirm the Vp/Vs ratio as a valuable tool for delineating oceanic crustal structure.

6.3 Porosity estimation of the layer 2A

Porosity is widely recognized as a crucial factor influencing both seismic velocity and the permeability in the upper crust and is a major determinant of fluid activities (e.g. Rosenberg et al. 1993; Alt et al. 1996, 2010; Coogan 2008; Carlson 2014a). Numerous studies have attempted to derive crustal porosity through various methods, such as laboratory measurements (e.g. Yu et al. 2016), borehole logs (e.g. Carlson 2010), effective medium theory (e.g. Taylor & Singh 2002) and relationships with Vp (Carlson 2014a). Carlson (2010) used effective medium models on borehole log data to demonstrate that variations in crack morphology and porosity primarily govern seismic velocities in the upper crust, whereas seismic velocities in the lower crust are entirely controlled by rock composition. However, some studies also highlighted the challenge of accounting for the extremely low velocities and high porosities in layer 2A (e.g. Whitmarsh 1978; Purdy 1987; Wilkens et al. 1991). In this study, we use various methods to estimate the porosity of the zero-age layer 2A at the EPR using Vp and Vp/Vs ratios, respectively.

Based on the Vp, we use effective medium theory with self-consistent approximation (SCA; Budiansky 1965; Hill 1965; Wu 1966) and differential effective medium (DEM; McLaughlin 1977; Cleary et al. 1980; Norris 1985; Marjanović et al. 2019a) to derive porosity-Vp relationships (Fig. 7a). The SCA determines the elastic moduli of the two-phase medium, consisting of background basalt and water inclusions in the pores, by introducing a single inclusion anomaly and approximating the interaction between the inclusion anomalies that replace the background matrix. The parameters for basalts are Vp = 6.0 km s⁻¹, Vs = 3.5 km s⁻¹ and ρ = 2.7 g cm⁻³; the parameters for water are Vp = 1.5 km s⁻¹, Vs = 0 and ρ = 1.0 g cm⁻³. We utilize the sphere-like pores in our SCA modelling. The DEM calculates macroscopic elastic moduli of the two-phase medium by incrementally adding the presence of inclusion anomalies to the background matrix. The elastic model parameters are the same as for the SCA. In contrast to the results from the SCA tests, the DEM yields relatively reasonable values for all aspect ratios except for the sphere-like pores (curves b, d and e in Fig. 7a). In addition to effective medium theory, we also examine the porosity–Vp relationship derived for layer 2A by Carlson (2014a).

Based on the Vp/Vs ratios or Poisson's ratios, we adopt the method of Brantut & David (2019), utilizing the DEM scheme to calculate the elastic properties of solids containing voids (i.e. dry pores) and Gassmann's equations (Gassmann 1951) to compute the effect of a fluid filling the voids (Fig. 7b). Overall, models based on effective medium theory demonstrate that both pore geometry and fluid compressibility exert a strong influence on the variations in Vp/Vs ratios with increasing fluid content (e.g. O'Connell & Budiansky 1974; Zimmerman 1994; Berryman et al. 2002). Low aspect ratios (α < 1) correspond to oblate spheroids (crack-like shapes), high aspect ratios (α > 1) correspond to prolate spheroids (needle-like shapes) and α = 1 corresponds to spheres. Low fluid compressibility ratios (ζ ≤ 10⁻²) indicate highly compressible fluids. From all the test results, we find that both α and ζ close to 0.1 provide a reasonable solution (Fig. 7b), which corresponds to a crack-like pores model with a high content of fluids and aligns well with the numerical modelling result of Marjanović et al. (2019a).

Both porosity estimations derived from individual Vp and Vp/Vs ratios exhibit a substantial range of variations within layer 2A (0.45–0.02 and 0.43–0.05), with higher maximum values than borehole data (Holes 504B and 1256D; Carlson 2010) due to low velocities (Fig. 7a). The simulation results also suggest that the layer 2A is more highly hydrated than previously assumed. Wilkens et al. (1991) modelled the rapid Vp increase in layer 2A by invoking the sealing of only thin cracks. Kim et al. (2019) used an OBS data set to present a joint inversion for Vp and Vs structures in the upper crust of the Endeavour segment of the Juan de Fuca Ridge. They attributed the high Vp/Vs ratios in layer 2B to a high density of thin cracks. Unfortunately, due to data limitations, we are unable to discuss the porosity variations within layer 2B. Porosity and crack morphology are crucial in controlling the seismic velocities in layer 2A, and high porous fluid content accounts for low Vp, Vs and high Vp/Vs ratios in this region.

6.4 Hydrothermal signatures from lateral heterogeneity of Vp/Vs ratios

The most striking feature of the Vp/Vs ratios model is the channel-like anomalies (Figs 5c and 8b), for which several potential explanations can be attributed. The effect of temperature is likely to be minimal at shallow depths (<1.0 km). Christensen (1996) demonstrated that rocks’ Vp/Vs ratios do not change significantly from 0 to 700 °C at 600 MPa, which was subsequently confirmed by Kern et al. (2001). Rock composition could influence Vp/Vs ratios; however, compositional variations are insufficient to account for the substantial lateral variations observed in the oceanic upper crust (e.g. Spudich & Orcutt 1980; Carlson 2014b). Barclay & Wilcock (2004) concluded that when earthquake activities shift to shallower depths, Vp and Vs increase, but Vp/Vs ratio decreases. Although large earthquake events are relatively rare in this region, microseismic events are present in the 9°50′N area. However, their border distribution is not fully constrained due to limitations in data coverage (Waldhauser & Tolstoy 2011). Consequently, lateral Vp/Vs ratio variations within the upper crust are primarily related to porosity/fracture characteristics and hydrothermal activities.

Ridge segmentation reflects differing degrees of mantle upwelling and melts delivery to the overlying crust (Macdonald et al. 1991). Third-order ridge-axis discontinuities are small-scale, with ∼30–100 km lengths and ∼0.5–1.0 km offsets. They resemble smaller OSCs and only exist for hundreds to thousands of years. Fourth-order discontinuities are characterized by small offsets (∼0.05–0.5 km) in the axial summit trough or axial pillow ridges and have a short lifespan (Macdonald et al. 1988; Smith et al. 2001; White et al. 2002). Previous studies have confirmed that the fourth-order discontinuities could represent hydrothermal recharge zones (e.g. Tolstoy et al. 2008). Marjanović et al. (2017) interpreted pipe-like low Vp anomalies along this profile as evidence of hydrothermal pathways, basing their conclusions on the spatial distribution of seafloor venting (Baker et al. 2016) and ridge segmentation patterns (White et al. 2006). Combining these results with our Vp/Vs ratios profile (Fig. 8b), we find that our Vp/Vs structure correlates well with their interpretations but offers more details (e.g. ∼9°19′N, ∼9°30′N, ∼9°37′N, ∼9°50′N). Additionally, we note that hydrothermal signatures identified from the high-precision Vp models (Marjanović et al. 2017; Vithana et al. 2019) are primarily located in the transition zone. The evolution of the upper crust is influenced mainly by mineral precipitation due to hydrothermal circulation, with increased and reduced velocities alternately marking the locations of hydrothermal discharge and recharge zones, respectively (Newman et al. 2011). Therefore, we propose that hydrothermal activities are primarily responsible for the broad Vp/Vs variations (∼2.1 ± 0.2) within the transition zone.

The ∼9°50′N area was the site of two documented volcanic eruptions in 1991–1992 (Haymon et al. 1993; Rubin et al. 1994) and 2005–2006 (Tolstoy et al. 2006). The origin of the pipe-like cluster of microseismicity between 9°49′ and 9°50′N was interpreted as a thermal contraction of the solid when in contact with cold fluid, occurring within the nearest fourth-order discontinuity (Fig. 8b; Waldhauser & Tolstoy 2011). Furthermore, numerous numerical simulations have confirmed that a systematic hydrothermal cell could develop at the ridge axis (e.g. Coumou et al. 2008, 2009; Lowell et al. 2013), and one of them, derived from the Vp model of Marjanović et al. (2019a), is also visible in our results. Cold seawater penetrates the oceanic crust through ridge discontinuities (cracks/fractures), causing high Vp/Vs anomalies. The seawater is heated above the underlying magma bodies and returns to the seafloor through pores, causing high Vp/Vs anomalies. Thus, the pipe-like Vp/Vs ratios structure could be interpreted as hydrothermal pathways. However, due to the limited streamer length, data quality and presence of the AMLs, we could not resolve deeper Vp/Vs anomalies.

7 CONCLUSIONS

1. Based on synthetic data modelling and analysis of 6-km MCS data, we demonstrate that a wider range of the compressional and shear waves can be observed and picked in downward-continued MCS data at fast-spreading ridges. We present a 2-D Vp/Vs ratio structure derived from independent P- and S-wave traveltime tomography along an ∼80-km-long stretch of the EPR 9°–10°N.

2. The average velocity structure we derive consists of three distinct velocity gradient layers, interpreted as pillow lavas, transition zone and sheeted dykes, respectively. The average thickness of pillow lava is ∼125 m, and its base coincides with the Vp contour of ∼3.2 km s⁻¹ and Vs contour of ∼1.5 km s⁻¹.

3. The approximately 400-m-thick transition zone exhibits the most pronounced velocity gradient and strongest heterogeneities in Vp/Vs ratios, indicating low rock consolidation and extensive fracturing/alteration within this layer. In conjunction with the transition zone, the pillow lavas likely constitute layer 2A, which has an average thickness of ∼525 m.

4. The Vp/Vs ratios of layers 2A (∼2.1 ± 0.2) and 2B (∼1.8–1.9) exhibit distinct characteristics. The 2A/2B boundary can be approximately identified by the Vp/Vs contour of 1.9 with Vp of ∼5.3 km s⁻¹ and Vs of ∼2.8 km s⁻¹. The pseudo-reflection of layer 2A cannot be reliably used to differentiate lavas and dykes at the ridge axis.

5. Porosity estimations derived from individual Vp and Vp/Vs ratios exhibit a substantial range of variations within layer 2A, spanning from 0.02 to 0.45, which signifies pronounced heterogeneities. Seismic velocities within layer 2A are predominantly governed by porosity and crack morphology.

6. Integrating our findings with a previous Vp interpretation (Marjanović et al. 2017), numerical simulation results (Marjanović et al. 2019a), ridge-axis discontinuity locations (White et al. 2006) and microseismicity distribution (Waldhauser & Tolstoy 2011), we interpret the lateral Vp/Vs anomalies within the upper crust as signatures of the hydrothermal system. Both fluid pathways into and out of the oceanic crust display elevated Vp/Vs anomalies. Our results suggest that the Vp/Vs ratios can be a robust reference parameter for investigating the oceanic crustal structure and hydrothermal system.

SUPPORTING INFORMATION

Figure S1. 1-D structure models for synthetic tests. (a) 1-D Vp profiles with the seafloor Vp are set to 2.5 km s⁻¹. (b) 1-D Vs profiles, where the Vs values below the seafloor are derived from Vp using a fixed Vp/Vs ratio of 1.7321. (c) 1-D density profiles, where the density distributions below the seafloor are calculated using Gardner's equation (Gardner et al. 1974) based on Vp profiles.

Figure S2. Synthetic downward continuation (DC) testing results for one-gradient velocity models. (a) Vertical Vp gradient = 1.0 s⁻¹. (b) Vertical Vp gradient = 1.5 s⁻¹. (c) Vertical Vp gradient = 2.0 s⁻¹. (d) Vertical Vp gradient = 2.5 s⁻¹. (e) Vertical Vp gradient = 3.0 s⁻¹. The green and red lines represent the synthetic traveltimes of P and S waves, respectively. Solid lines indicate pickable data, whereas dotted lines denote data that cannot be picked. TWTT, two-way traveltime.

Figure S3. Observed data from line axis2r1 at the East Pacific Rise (EPR). (a) Shot# 2019 at 11.95 km. (b) Shot# 3519 at 68.20 km.

Figure S4. Observed shot gathers used for traveltime picking of P and S refractions from line axis2r1 at the EPR. (a) Shot# 2087 at 14.50 km. (b) Shot# 2487 at 29.50 km. (c) Shot# 2887 at 44.50 km. (d) Shot# 3287 at 59.50 km. Shot locations are indicated by blue triangles in Fig. 1. The offset range for traveltime picking is based on waveform consistency within each shot gather. Traveltime uncertainties of P and S waves are 10 and 20 ms, respectively. The red line represents observed (picked) traveltimes. Green and red lines indicate ray-traced traveltimes for the initial and final models, respectively.

Figure S5. P- and S-wave traveltime picking for profile axis2r1. (a, b) Observed traveltimes of P and S waves are displayed as a function of shot number (#) and channel. 490 shots with 152 406 traces are picked for P waves (a) and 115 902 for S waves (b). A vertical column indicates the time picking along the streamer of a single shot. Due to data quality or phase interference, white regions represent the absence of traveltime picks for specific shot-receiver pairs. TWTT, two-way traveltime.

Figure S6. Selection of initial models for P and S waves. (a, b) The contour of the root-mean-squared traveltime residuals calculated for P- and S-wave models with various seafloor velocities and vertical velocity gradients. Red stars represent the models selected as the initial models for traveltime tomography.

Figure S7. Traveltime residuals predicted from initial models and traveltime tomography models. (a) P-wave tomography. (b) S-wave tomography. N represents the number of picked traveltimes, RMS denotes root-mean-squared traveltime residuals, and χ² indicates the weighted misfit function.

Figure S8. Traveltime residual histograms of initial and final models. (a) P-wave tomography. (b) S-wave tomography.

Figure S9. Velocity perturbations with respect to the average 1-D structure. (a) Vp perturbation. The one-gradient model with a seafloor velocity of 3.25 km s⁻¹ and vertical gradient of 3.00 s⁻¹ was chosen as the initial Vp model. (b) The Vs perturbation. The one-gradient model with a seafloor velocity of 1.25 km s⁻¹ and vertical gradient of 2.50 s⁻¹ was chosen as the initial Vs model.

Figure S10. Vertical velocity gradient distributions. (a) Vp vertical velocity gradient. (b) Vs vertical velocity gradient.

Figure S11. Uncertainty analysis of Vs model. (a) Initial Vs models with varying Vp/Vs ratios for traveltime tomography. (b) The standard deviation of inverted Vs models.

Figure S12. Checkerboard resolution tests for P-wave traveltime tomography. The true anomaly pattern (left-hand column) is created using a sine function with a maximum amplitude of ±10 per cent. Recovered anomalies are displayed in the right-hand column. Checkerboard sizes include 2.0 km × 1.0 km, 1.5 km × 0.5 km, 1.0 km × 0.5 km and 0.5 km × 0.25 km.

Figure S13. Checkerboard resolution tests for S-wave traveltime tomography. The true anomaly pattern (left-hand column) is created using a sine function with a maximum amplitude of ±10 per cent. Recovered anomalies are shown in the right-hand column. Checkerboard sizes include 2.0 km × 1.0 km, 1.5 km × 0.5 km, 1.0 km × 0.5 km and 0.5 km × 0.25 km.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGMENTS

We thank Editor Jenny Collier, Associate Editor Louise Alexander, Reviewer Helene Carton and an anonymous referee for their valuable input in improving the original manuscript. We also extend our gratitude to the crew and science parties of the MGL0812 experiment. Synthetic seismic data simulation was performed using the DENISE code (https://github.com/daniel-koehn/DENISE-Black-Edition). The downward continuation code, developed by Alistair Harding, is available upon request from the corresponding author. Traveltime tomography was conducted using the FAST (First Arrival Seismic Tomography program) code (https://terra.rice.edu/department/faculty/zelt/fast.html). The microseismicity locations at EPR 9°50′N were obtained from Marine Geoscience Data System (MGDS; https://www.marine-geo.org/tools/datasets/24032). All figures in this paper were plotted using GMT6 (https://www.generic-mapping-tools.org), NumPy (https://numpy.org) and Matplotlib (https://matplotlib.org). This research was supported by the National Natural Science Foundation of China (91858207, 41676044) and the Guangdong Basic and Applied Basic Research Foundation (2021B1515020023).

Author contributions: HD: investigation, methodology, software, data curation, original draft writing and revising. WX: review and editing. MX: conceptualization, review and editing, funding acquisition.

CONFLICT OF INTEREST

Declare no competing interests.

DATA AVAILABILITY

Supplementary text and figures supporting the conclusions are provided in Supporting Information. The MCS data of MGL0812 can be accessed through Marine Geoscience Data System (MGDS; https://www.marine-geo.org/tools/search/entry.php?id=MGL0812). The downward-continued MCS data, P- and S-wave traveltime picking files, and velocity models generated in this study are available on Zenodo at https://doi.org/10.5281/zenodo.7309649.

References