Oceanic crust—seismic structure, lithology and the cause of the 2A Event at borehole 504B

SUMMARY

This study focuses on the 3-D velocity structure and thickness of ∼7-Myr-old oceanic crust surrounding borehole 504B, located ∼235 km from the intermediate-spreading Costa Rica Rift (Panama Basin). It investigates how well seismic structure determined by 3-D tomography compares with actual lithology and, consequently, what the origin and cause might be of an amplitude anomaly, the 2A Event, that is observed in multichannel seismic data. Our P-wave model shows an ∼0.3-km-thick sediment layer of velocity between ∼1.6 and 1.9 km s⁻¹ (gradient 1.0 s⁻¹), bound at its base by a velocity step to 4.8 km s⁻¹ at the top of oceanic crustal Layer 2. Layer 2 itself is subdivided into two main units (2A and 2B) by a vertical velocity gradient change at 4.5 km depth, with a gradient of 1.7 s⁻¹ above (4.8–5.8 km s⁻¹) and 0.7 s⁻¹ below (5.8–6.5 km s⁻¹). The base of Layer 2, in turn, is defined by a change in gradient at 5.6 km depth. Below this, Layer 3 has a velocity range of 6.5–7.5 km s⁻¹ and a gradient of ∼0.3 s⁻¹. Corresponding S-wave igneous layer velocities and gradients are: Layer 2A, 2.4–3.1 km s⁻¹ and 1.0 s⁻¹; Layer 2B, 3.1–3.7 km s⁻¹ and 0.5 s⁻¹; Layer 3, 3.7–4.0 km s⁻¹ and 0.1 s⁻¹. The 3-D tomographic models, coupled with gravity modelling, indicate that the crust is ∼6 km thick throughout the region, with a generally flat-lying Moho. Although the P- and S-wave models are smooth, their velocities and gradients are remarkably consistent with the main lithological layering subdivisions logged within 504B. Thus, using the change in velocity gradient as a proxy, Layer 2 is interpreted as ∼1.8 km thick and Layer 3 as ∼3.8 km thick, with little vertical variation throughout the 3-D volume. However, the strike of lateral gradient variation is not Costa Rica Rift-parallel, but instead follows the orientation of the present-day adjacent Ecuador Rift, suggesting a reorientation of the Costa Rica Rift spreading ridge axis. Having determined its consistency with lithological ground-truth, the resulting P-wave model is used as the basis of finite difference calculation of wave propagation to find the origin of the 2A Event. Our modelling shows that no distinct interface, or transition, is required to generate this event. Instead, it is caused by averaging of heterogeneous physical properties by the seismic wave as it propagates through Layer 2 and is scattered. Thus, we conclude that the 2A Event originates and propagates exclusively in the lower part of Layer 2A, above the mean depth to the top of the dykes of Layer 2B. From our synthetic data we conclude that using the 2A Event on seismic reflection profiles as a proxy to determine the Layer 2A/2B boundary's depth will result in an overestimate of up to several hundred metres, the degree of which being dependent on the specific velocity chosen for normal moveout correction prior to stacking.

Composition and structure of the oceanic crust, Controlled source seismology, Crustal imaging, Crustal structure

1. INTRODUCTION

The seismic study of mid-ocean ridges has resulted in a widely accepted, generalized layer-depth model of the oceanic crust (Raitt 1963; Houtz & Ewing 1976), from which ‘type’ average structures have evolved (e.g. White et al. 1992; Grevemeyer et al. 2018b; Christeson et al. 2019). The generalized model comprises three main layers in which Layer 1 is sediment with a velocity ranging between 1.6 and 2.5 km s⁻¹, that is very thin or entirely absent at a ridge axis (zero-age crust), but can be several kilometres thick off-axis (e.g. Shearer & Orcutt 1986). Layer 2 has a velocity between 2.2 and 6.6 km s⁻¹, a high rate of change in velocity with depth (velocity gradient of 1–2 s⁻¹) and variable thickness (e.g. White et al. 1992). Finally, Layer 3 has a velocity between 6.6 and 7.6 km s⁻¹ (e.g. White et al. 1992) or 6.93–7.18 km s⁻¹ (e.g. Grevemeyer et al. 2018b), a low velocity gradient (0.1 s⁻¹) and little thickness variation. A legacy of early models of crustal accretionary processes is that Layer 2 is generally associated with the basaltic extrusive layer of the upper crust (Raitt 1963), while Layer 3 is associated with the intrusive (plutonic) gabbros of the lower crust (White et al. 1992; Carlson & Miller 2004). A further layer, Layer 4, is associated with the upper mantle, with the boundary between Layers 3 and 4 defined as the Moho where a velocity increase to ∼8 km s⁻¹ is generally expected.

As more studies are undertaken with increased resolution and denser sampling, the variability in vertical velocity gradient of Layer 2 has resulted in it being further divided into 2A and 2B, where 2A is ascribed to mainly subhorizontal pillow lavas, massive lava flows and breccias, and 2B to their feeder sheeted dykes (Herron 1982; Christeson et al. 1992; Harding et al. 1993); the boundary between these being defined by a change in velocity gradient (e.g. Grevemeyer et al. 2018b) or a transition zone (e.g. Christeson et al. 2012). Although this simple two subdivision structure (Raitt 1963) is widely accepted for the ridge axis, it is less so off-axis where an alternative three subdivision structure has further developed (Houtz & Ewing 1976) that is alternatively proposed to result from hydrothermal circulation ‘aging’ of the crust as it spreads off-axis. In this model, alteration or the infilling of voids and fractures results in an overall reduction in porosity (e.g. Carlson 2001; Christeson et al. 2007). Here, Layers 2A and 2B are ascribed to less and more altered lavas and 2C to sheeted dykes. This vertical variability in structure is further enhanced by lateral increase in density (and hence velocity) as the crust cools, and as sediment progressively accumulates and seals it from hydrothermal fluid flow as it spreads away from the ridge axis (Christensen 1970; Stein & Stein 1992; Ridley 1995; Alt et al. 1996; Grevemeyer & Weigel 1996; Nedimović et al. 2008; Carlson 2011; Hasterok 2013; Harris et al. 2015; Patten et al. 2016; Wilson et al. 2019; Funnell et al. 2021).

In contrast to the number of seismic surveys that have now been undertaken, there are relatively few in situ sampling studies that enable models based on velocity to be fully understood in terms of lithology. To date, understanding of the lithology has primarily been through the study of ophiolites (e.g. Penrose Conference Participants 1972), dredging of rock samples from the seabed or fault scarps (e.g. Francheteau et al. 1992; Juteau et al. 1995) or the drilling and combined coring and petrophysical logging of a relatively small number of boreholes that have limited subseabed penetration (e.g. Anderson et al. 1982; Wilson et al. 2006). Even then, relating the lithological ground-truth with models and images derived from seismic surveys is challenging given the diversity and complexity of factors that control the seismic imaging response (e.g. Wilkens et al. 1991; Berge et al. 1992; Carlson 2001, 2011, 2014), let alone the contrasting scales of resolution of the various approaches.

1.1 Borehole 504B

Perhaps the best known and most extensively studied borehole yet drilled into the oceanic crust is Deep Sea Drilling Project/Ocean Drilling Program borehole 504B (henceforth simply referred to as 504B), which was drilled into ∼6.9 Ma crust on the southern flank of the Costa Rica Rift spreading centre between 1979 and 1993 (Fig. 1), and is often cited as the standard reference for oceanic crustal structure (e.g. Bratt & Purdy 1984; Carlson & Herrick 1990; Cudrak & Clowes 1993; Grevemeyer et al. 1999; Harris et al. 2015; Christeson et al. 2016). Legs 69 and 70 were first drilled to investigate geothermal processes, penetrating to a depth of 562 m below the overlying sediment (CRRUST 1982). Legs 83, 111, 137, 140 and, in 1993, leg 148 deepened 504B to 1837 m. Becker et al. (1989), Alt et al. (1996) and Gregory (2018) summarize the crustal structure as a 275-m-thick layer of siliceous oozes, chert, limestone and chalk, overlying a 575-m-thick section of pillow lavas, massive lava flows and breccias (hereafter collectively referred to as lavas) which have undergone low-temperature alteration due to seawater circulation. Underneath lies a 209-m-thick transition zone of lavas in which the proportion of dykes increases with depth, and which includes an intensely mineralized zone formed by the interaction between cool, downwelling seawater with hot, upwelling hydrothermal fluids. This, in turn, overlies a 1050-m-thick complex of sheeted diabase dykes which extends to the base of the borehole where gabbros begin to be sampled.

Figure 1.

Bathymetry map, after Funnell et al. (2021), showing the tectonic setting of the Costa Rica Rift within the Panama Basin (GEBCO Compilation Group 2020). Spreading ridges (white lines): Galapagos Spreading Centre (GCS), Ecuador Rift (ER) and Costa Rica Rift (CRR). Transform-fracture zones (black lines): Ecuador Fracture Zone (EFZ), Galapagos Transform (GT), Inca Transform Fault (IT) and Panama Fracture Zone (PFZ). The OSCAR study area is outlined by the blue dashed polygon. South Grid is outlined by the red square (see Fig. 2). Borehole 504B is located by the solid yellow circle. Plate motion directions and rates are shown by the numbered red arrows.

Petrophysical logs show a decrease in porosity from >10 per cent at the top to ∼1 per cent at the bottom of the borehole, a low, but highly variable P-wave velocity zone ∼0.8 km subsediment (4–6 km s⁻¹) and, below that, a less variable but higher velocity region (>6 km s⁻¹) to the base of the borehole. In particular, sonic logs (Cann & Von Herzen 1983; Becker et al. 1988) divide Layer 2 into three layers interpreted as 2A ranging in velocity between 3.0 and 4.0 km s⁻¹, 2B ranging in velocity between 4.5 and 5.5 km s⁻¹ and 2C ranging in velocity between 6.0 and 6.5 km s⁻¹. The final 50 m of the borehole has a velocity more typical of Layer 3 (∼6.8 km s⁻¹; Alt et al. 1993), suggesting that 504B may bottom out in the Layer2/Layer 3 transition zone.

1.2 Velocity versus lithology

Several studies have attempted to correlate velocity–depth structure observed by seismic experiments at or surrounding 504B with the ground-truth it provides. For example, Little & Stephen's (1985) Oblique Seismic Experiment (OSE) showed an increase in velocity with depth to 1.4 km below sediment, with no first-order interfaces (step change in velocity) correlating with either lithological or alteration boundaries. They also observe a ridge-parallel anisotropy in the uppermost 50–100 m subsediment. Vertical seismic profiling (VSP) by Becker et al. (1988) and Swift et al. (1998) show a change in velocity gradient at 875–925 m subsediment, which they interpret as the Layer 2/Layer 3 boundary, while a sonobuoy seismic reflection/refraction study (Hobart et al. 1985) reports a middle crust with a high velocity gradient, a low-velocity zone above the Moho ∼1800-m-thick, an ∼5 km crustal thickness (Collins et al. 1989), and an estimate of Layer 3 velocity of 6.5 km s⁻¹ (Detrick et al. 1994). Finally, Detrick et al. (1998) report the results of a 3-D seismic refraction experiment which revealed a 2-km-deep subsediment, effectively homogeneous velocity structure, with relatively faster and slower velocity regions beneath ridge-parallel basement ridges and basement troughs respectively.

Most recently the OSCAR project (Oceanographic and Seismic Characterisation of heat dissipation and alteration by hydrothermal fluids at an Axial Ridge; Hobbs & Peirce 2015) aimed to investigate the evolution of oceanic crustal structure with age in an attempt to reconcile the 2A/2B versus 2A/2B/2C debate. This study targeted the Costa Rica Rift (CRR) because of the location and proximity of 504B (Fig. 1), lying ∼235 km south along a flowline from the ridge axis, such that it could be used to provide lithological ground-truth to the 2-D and 3-D seismic images acquired. The CRR has been asymmetrically spreading for 11 Myr (Lonsdale & Klitgord 1978), and is currently spreading with a half-spreading rate (HSR) of 30 mm yr⁻¹ for the northern Cocos Plate and 36 mm yr⁻¹ for the southern Nazca Plate (Wilson & Hey 1995). Although the CRR is classified as an intermediate spreading ridge, it has a Mid-Atlantic Ridge (MAR)-type axial valley morphology, and is divided in two by an overlapping spreading centre (OSC) located at 3°20’N, 83°44’W (Fig. 2a). The ∼180-km-long CRR is offset from the adjacent Ecuador Rift (ER) by the Ecuador Fracture Zone (Fig. 1), with these two ridge systems currently having an ∼8–10° difference in their spreading direction from ∼1.2 Ma (Peirce et al. 2023b).

Figure 2.

South Grid and its setting within the OSCAR project. (a) Swath bathymetry with the Costa Rica Rift spreading axis represented by the white-red dashed line. The red box shows the footprint of the South Grid. CRR, Costa Rica Rift; NG, North Grid (Robinson et al. 2020); SAP, Synthetic Aperture Profiles (Wilson et al. 2019; Funnell et al. 2021) and SG, South Grid. (b) Layout of the South Grid seismic acquisition with profile and OBS names annotated. The map footprint shows the extent of the 3-D inversion models (45 km × 45 km). In both parts 504B is located by the solid yellow circle, black lines represent seismic profiles and white triangles OBS locations.

The geophysical objectives of the OSCAR study (Figs 1 and 2) were to:

1. Resolve the velocity–depth structure of zero-age crust via a 3-D seismic refraction tomography study centred at the ridge axis; ‘North Grid’ (NG; Zhang et al. 2017; Robinson et al. 2020).

2. Resolve the velocity–depth structure of ∼7-Myr-old crust off-axis via a 3-D seismic refraction tomography study sited at borehole 504B; ‘South Grid’ (SG; as part of which Gregory (2018) focuses solely on the lithology and anisotropy of Layer 2).

3. Resolve the velocity–depth structure of 0 to ∼7-Myr-old crust via three coincident flowline Synthetic Aperture multichannel seismic (MCS) reflection and 2-D long-offset wide-angle (WA) refraction profiles (SAP_A, SAP_B and SAP_C; Wilson et al. 2019; Funnell et al. 2021).

4. Resolve the structure of the Ecuador Fracture Zone via a network of 13 MCS reflection profiles (EX; Tedd 2021; Peirce et al. 2023b).

5. Investigate the origin of the 2A Event interpreted in MCS reflection images and often used to determine the depth and geometry of the Layer 2A/2B boundary (e.g. Harding et al. 1993; Christeson et al. 1996, 2007, 2010; Newman et al. 2011; Arnulf et al. 2012; Audhkhasi & Singh 2019; Estep et al. 2019).

These geophysical objectives were complemented by physical oceanography and heat flow measurements and modelling, to determine the influence of crustal heat on the bottom boundary layers of the ocean (Kolandaivelu et al. 2017; Barnes et al. 2018; Banyte et al. 2018a,b; Lowell et al. 2020).

1.3 This study

This study addresses objectives (2) and (5) above and focuses on the 3-D velocity structure and thickness of the ∼7-Myr-old crust of the South Grid. In addition, it also aims to investigate how well, or otherwise, the seismic structure determined by 3-D tomography compares with the lithological ground-truth provided by 504B and, thus, provide a better understanding and framework within which to interpret seismic studies of the oceanic crust in the absence of any lithological control. The resulting 3-D model of Layer 2 structure is further used to inform modelling of the origin of the 2A Event and demonstrate how MCS surveys, and the processing of their data, should be used to constrain Layer 2 structure and thickness. The geological setting is described in detail in each of the OSCAR papers cited above and will, therefore, not be repeated here where we focus on the specific part of the project centred on 504B.

2. DATA

Swath bathymetry and gravity data were acquired throughout the South Grid during RRS James Cook cruise JC114, and the details of their acquisition can be found in Hobbs & Peirce (2015). The seismic component of OSCAR aimed to image the crust surrounding 504B using a 3-D tomographic approach, with an array of 25 ocean-bottom seismographs (OBSs) deployed in a 5 km × 5 km grid centred on 1°14.017′N, 83°46.664′W (OBS SG_13; Fig. 2). This South Grid was connected to a matching North Grid (Robinson et al. 2020), to enable comparison with the zero-age crust currently forming at the CRR, via flowline profiles (e.g. SAP_B; Wilson et al. 2019; Funnell et al. 2021) designed to determine how crustal structure has changed with age in between. The 3-D tomography seismic velocity–depth model of the South Grid can, thus, be lithologically ground-truthed by 504B, and this seismic reference propagated to all other seismic images (reflection profiles or refraction velocity–depth models) via their intersections.

2.1 MCS and WA data acquisition

Seismic reflection and refraction data were acquired contemporaneously along five E–W (SG_A to SG_E) and five N–S (SG_F to SG_J) profiles through the OBS array (Fig. 2b), using two seismic source configurations (GI and Bolt) fired in a flip-flop fashion. Two further profiles (SG_M and SG_N) where shot to the north and east of the South Grid respectively, to enhance deeper crust and uppermost mantle sampling. The shot cycle was repeated every 30 s, with a higher dominant frequency GI airgun array fired first for MCS imaging to determine the thickness and internal layering structure of the sediment cover. A lower dominant frequency Bolt airgun array was fired 6 s later, for WA crustal imaging and to record the 2A Event, to allow the wave train from the GI array to dissipate sufficiently so as not to obscure the Bolt array first arrivals. The 1620 in³ (26.55 l) seven Bolt airgun array was towed at 8 m depth and the 420 in³ (6.88 l) two GI airgun array at 5 m depth. Strong southerly surface currents resulted in profiles shot from north-to-south having a shot interval of 75 m and those south-to-north 95 m, with east–west profiles averaging 70 m. For both the South Grid and SAP_B profiles the shots were recorded by a 4.5-km-long Sercel Sentinel streamer towed by the Cook, which had 360 channels spaced at 12.5 m. Example MCS 4.5 km streamer record sections, for both Bolt and GI shots, are shown in Figs 3(a)–(d), and the general processing scheme is outlined in the Supporting Information.

Figure 3.

Multichannel seismic images along coincident flowlines from the Costa Rica Rift, traversing the South Grid (profiles SG_I and SAP_B). (a) Higher frequency GI airgun source profile for SG_I (4.5 km streamer), primarily used to image the sediment thickness and layering. (b) Lower frequency Bolt airgun source profile for SG_I (4.5 km streamer), primarily used to image the sediment-basement interface depth and geometry. (c) Lower frequency Bolt airgun source profile for SAP_B (4.5 km streamer), equivalent to (b), except shot under different sea state and ambient current conditions. (d) As (c) with processing aimed at Moho reflection imaging, if present. This event should arrive at a two-way traveltime within the dashed black outline but is not observed. (e) Synthetic aperture supergathers (8.7 km streamer simulation), shot with both Bolt and G airguns, processed to enhance imaging of the 2A Event, which is observed consistently along flowline through the South Grid footprint.

2.2 MCS imaging of the 2A Event

Profile SAP_B was also shot using the two-ship synthetic aperture approach (Buhl et al. 1982) to generate the sufficiently long subgather offsets required to image the 2A Event. The second vessel, FS Sonne, sailed 8.7 km behind the Cook, which towed the 4.5 km streamer, to provide a 0.3 km recording overlap between shots fired from the Cook and shots fired from the Sonne (Hobbs & Peirce 2015). The Cook fired both a GI array (at time 0 s) and a Bolt array (at +6 s), while the Sonne fired a G array (at +30 s), with the entire shot firing cycle repeated every 60 s. The 1120 in³ (18.35 l) six Bolt airgun array was towed at 8 m depth, the 420 in³ (6.88 l) two GI airgun array at 5 m depth, and the 4280 in³ (70.14 l) nine G airgun array at 10 m depth. The shot interval ranged between 130 and 160 m due to variable surface currents. The MCS synthetic 8.7 km streamer record section imaging the 2A Event for the part of SAP_B running through the South Grid is shown in Fig. 3(e), and the general synthetic aperture processing scheme is outlined in the Supporting Information.

2.3 OBS traveltime picking

Each OBS recorded the output of a three-component geophone package and a hydrophone using a sampling rate of 500 Hz. Traveltimes of P-wave first arrivals were picked using hydrophone data (Figs 4, 5, S1 and S2) as it demonstrated the highest signal-to-noise ratio. Traveltimes of S-wave arrivals were picked primarily from hydrophone data, whilst comparing with geophone data to avoid phase misidentification. We interpret the picked S-waves as having first travelled as a P-wave through the seawater and sediment, and then as an S-wave for the remainder of their travel path, converting between P and S at the sediment–basement interface. Approximately 144 000 P-wave and 141 800 S-wave arrival traveltimes were picked and each assigned an uncertainty to account for picking and instrument location errors; P-waves 40 ms and S-waves 60 ms.

Figure 4.

Example wide-angle seismic data from OBS SG_12 (red triangle). (a) P-wave traveltime pick locations (white), highlighting those from profile SG_C [red (d)] and SG_I [blue (f)]. Borehole 504B is located by the solid yellow circle. (b) As (a) for S-wave traveltime pick locations. (c) Hydrophone record section for ridge-parallel profile SG_C, reduced at 6.0 km s⁻¹ and bandpass filtered between 2–4 and 40–60 Hz. (d) As (c), showing the P- and S-wave traveltime picks (red) and base of crust (Moho) reflections (PmP) calculated using the 7.5 km s⁻¹ contour from the 3-D inversion model (orange; Fig. 7) as a proxy for the Moho, and a Moho depth and geometry derived from gravity modelling (purple; Fig. 9). (e) As (c) for flowline profile SG_I. (f) As (d) for flowline profile SG_I. The pick density through the South Grid is shown in Fig. S3. The intermittent nature of the PmP arrival calculated using the 7.5 km s⁻¹ contour proxy is caused by its undulating shape (cf. Figs 9a and b).

Figure 5.

Example wide-angle seismic data from OBS SG_05 (red triangle). (a) P-wave traveltime pick locations (white), highlighting those from profile SG_A [red (d)] and SG_F [blue (f)]. Borehole 504B is located by the solid yellow circle. (b) As (a) for S-wave traveltime pick locations. (c) Hydrophone record section for ridge-parallel profile SG_A, reduced at 6.0 km s⁻¹ and bandpass filtered between 2–4 and 40–60 Hz. (d) As (c), showing the P- and S-wave traveltime picks (red) and base of crust (Moho) reflections (PmP) calculated using the 7.5 km s⁻¹ contour from the 3-D inversion model (orange; Fig. 7) as a proxy for the Moho, and a Moho depth and geometry derived from gravity modelling (purple; Fig. 9). (e) As (c) for flowline profile SG_F. (f) As (d) for flowline profile SG_F. The pick density through the South Grid is shown in Fig. S3. The intermittent nature of the PmP arrival calculated using the 7.5 km s⁻¹ contour proxy is caused by its undulating shape (cf. Figs 9a and b).

Within the South Grid the maximum shot-receiver offset is less than ∼40 km. When set in the context of a crustal 2-D transect (e.g. offsets up to and beyond 100 km), this would be regarded as near or short offset and it would be expected that it should be possible to pick traveltimes to effectively the same degree of accuracy within that offset range. The OSCAR data is of particularly good quality, in that arrivals are clear and of an amplitude well above the background noise at all offsets within the 3-D grid (cf. Figs 4 and 5 with Figs S1 and S2). Thus, a variation in pick uncertainty with offset was not deemed necessary, and an effectively consistent density of P-wave and S-wave picks was achieved (Fig. S3), with no apparent bias in the north–south or east–west directions.

3. TOMOGRAPHIC INVERSION

Both P- and S-wave traveltimes were inverted using FAST (First-Arrival Seismic Tomography; Zelt & Barton 1998), which produces a smooth velocity model without any discrete interfaces below the bathymetry, which was held static during inversion. Above the bathymetry interface, the seawater layer remained fixed throughout inversion. All inversion models were parametrized on a 0.1 km × 0.1 km × 0.1 km uniform (cubic) forward grid (to be consistent with Wilson et al. 2019; Robinson et al. 2020; Funnell et al. 2021), within a 45 km × 45 km model footprint, extending to a depth of 15 km below sea surface (bss). The 0,0 km (in x,y) origin of the model space is located at 1°02.628′N, 83°58.433′W, and model distance was assigned to increase northwards and eastwards. The seabed (bathymetry interface) was constructed by sampling the high-resolution swath bathymetry data at the forward model node spacing.

3.1 Initial models

Flowline SAP_B traverses 504B and runs through the South Grid along profile SG_I (Fig. 2). Wilson et al. (2019) report the result of higher-resolution P-wave 2-D tomographic modelling of Layer 2 along SAP_B (Fig. 6a), while Funnell et al. (2021) report the result of both P- and S-wave 2-D tomographic modelling of the entire crust and uppermost mantle (Figs 6c and d). Consequently, these models were sampled at the location of 504B (27.1, 20.4 km in model space) to create a 1-D velocity–depth profile which was subsequently added to the seabed depth at each model node to, in effect, create seabed-following, starting point 3-D initial models for both P- and S-wave inversions (Figs 6g and h). The rationale for this approach was based on:

1. The linearization assumption embedded in FAST (Zelt & Barton 1998) that requires that ‘small perturbations to the starting model are determined’.

2. That the initial model should be sufficiently different from a model that fits the traveltime picks within error, such that the inversion iterates away from the starting point to a point of minimum misfit.

3. That the initial model should be free from any imposed modeller bias or preconceived idea of structure (velocity with depth or laterally).

Figure 6.

Initial model. (a) Wilson et al.’s (2019) P-wave velocity–depth model for Layer 2 along SAP_B, based on both MCS gather and OBS upper crustal first arrival traveltime picks. Blue dashed lines show the flowline extent of the South Grid. Red dashed line locates 504B. (b) P- and S-wave 1-D velocity–depth profiles sampled through Funnell et al.’s (2021) models at 504B (black dashed) plotted relative to the seabed (red). These models (c and d) are based on OBS traveltime picks for all observed first arrivals. The solid black line shows the equivalent profiles sampled through the initial inversion models of this study. (c) Funnell et al.’s (2021) P-wave velocity–depth model for SAP_B. (d) Funnell et al.’s (2021) S-wave velocity–depth model for SAP_B. (e) Funnell et al.’s (2021) P-wave velocity–depth model for the flowline extent of the South Grid. (f) As (e) for the S-wave model. (g) Initial P-wave velocity–depth model of this study. (h) As (g) for the initial S-wave model.

Based on these P-wave and S-wave starting points and assumptions, a suite of initial models was created that varied in velocity by up to ±10 per cent. These initial models each formed the basis of a series of inversions whose outcomes produced very similar best-fitting models. Consequently, an end-member initial model (SAP_B plus 10 per cent) of this suite was selected for both the P-wave (Fig. 7b) and S-wave (Fig. 8b) initial models.

Figure 7.

Slices through the 3-D P-wave inversion model. (a) Bathymetry map showing the location of 504B (solid yellow circle). Dashed lines in (a) and all horizontal (depth) slices indicate the locations of the vertical slices through the velocity model. (b) Vertical slice through the initial model at y = 21 km; ridge-parallel through 504B. (c) Depth slice through the inversion model at z = 4.5 km bss. (d, f, h and j) Vertical slices through the inversion model at the annotated model distances. (e, g and i) Depth slices through the inversion model at the annotated model depths bss. In all inversion model parts, the slices are masked by the inversion model ray coverage. The 7.5 km s⁻¹ contour is used as the proxy for the crust-mantle transition and as the basis of PmP traveltime predictions shown in Figs 4 and 5. Depth slices are illuminated by the seabed topography. Model parameters, the number of traveltime picks included and the resulting χ² and RMS misfit are indicated.

Figure 8.

Slices through the 3-D S-wave inversion model. (a) Bathymetry map showing the location of 504B (solid yellow circle). Dashed lines in (a) and all horizontal (depth) slices indicate the locations of the vertical slices through the velocity model. (b) Vertical slice through the initial model at y = 21 km; ridge-parallel through 504B. (c) Depth slice through the inversion model at z = 4.5 km bss. (d, f, h and j) Vertical slices through the inversion model at the annotated model distances. (e, g and i) Depth slices through the inversion model at the annotated model depths bss. In all inversion model parts, the slices are masked by the inversion model ray coverage. Depth slices are illuminated by the seabed topography. Model parameters, the number of traveltime picks included, and the resulting χ² and RMS misfit are indicated.

3.2 Inversion

All modelling parameters were systematically tested for both P- and S-wave modelling, to find a combination that resulted in minimal modelling artefacts, good lateral and vertical resolution, and a dense and even ray coverage within the model that reflected the sampling provided by the traveltime pick distribution (Fig. S3). The parameters consequently selected for inversion of both the P- and S-wave traveltime pick sets are summarized in Table 1. A single-pass inversion process was applied comprising a series of five iterations, with each iteration testing five values of the trade-off parameter (λ), which controls the minimization of the data misfit compared to the minimum structure required to achieve the fit. Both P- and S-wave inversions converged to a χ² of ∼1 within two iterations (Fig. S4). The root-mean-square (RMS) misfit distributions of the initial and final models, for both P- and S-wave traveltime pick sets, are shown in Fig. S5, which also shows that the inversion reduces the misfit for the majority of picks to within the pick error, and that the misfit associated with both the P- and S-wave final models has little-to-no dependence on shot-receiver offset. Vertical and horizontal slices through the best-fitting P-wave model are shown in Fig. 7 and this model has an overall RMS misfit of 39.7 ms, while the best-fitting S-wave model is shown in Fig. 8 which has an RMS misfit of 60.0 ms.

Table 1.

Inversion parameters, traveltime picks and goodness of fit.

Model	Forward model node spacing (km)	Inversion grid spacing (km)	Perturbation weighting factor, α	Initial trade-off parameter, λ0	Reduction factor for trade-off parameter, λ	Smoothness factor, Sz	Edge constraint, Sedge
Both	0.1 × 0.1 × 0.1	0.5 × 0.5 × 0.5	0.95	100	1.414	0.125	10
	Traveltime picks (total)	Pick uncertainty (ms)	TRMS (ms)	χ2
P-wave	144 041	40	39.7	1.0
S-wave	141 778	60	60.0	1.0

3.3 Resolution

To determine how well the observed velocity anomalies are resolved, checkerboard testing of both P- and S-wave models was undertaken using the approach of Zelt & Barton (1998). A ±5 per cent velocity perturbation was added to each inversion model (Figs S6 and S7), taking the form of a columnar pattern with depth (Zelt 1998). This pattern was selected based on Zelt's (1998) assertion that, for 3-D tomographic inversion, horizontal resolution is primarily controlled by the ray coverage through the model, while vertical resolution is controlled by both the vertical model parametrization (forward cell size) and the inversion cell size. On this basis, we define the vertical resolution of our models as the larger of these two; 0.5 km, the inversion vertical cell size.

The basis of checkerboard testing was a set of synthetic traveltimes calculated by ray tracing through each perturbed model, to which Gaussian noise was added which was scaled to the pick uncertainty. These synthetic P- and S-wave pick sets were then inverted using the same parametrization (Table 1) used to achieve the corresponding best-fitting inversion model, with the latter model used as the checkerboard inversion starting point. The resulting P- and S-wave checkerboard inversions both have a χ² of ∼1.

The resolvability of particular length scales within a model is quantified by the semblance (Figs S6 and S7; Zelt 1998). We adopted Zelt's (1998) semblance threshold of 0.7 to indicate where each model is well resolved. Three checkerboard pattern sizes in the x- and y-directions were tested (2, 3 and 5 km), with half-cell shifts also applied to test for anomaly location and bias (Zelt 1998; Zelt & Barton 1998). Consequently, we tested nine unique checkerboard patterns for each pattern size, which were averaged to create the semblance volume and quantify model recoverability.

Checkerboard testing suggests that both P- and S-wave inversion models can resolve 2 km × 2 km sized, ±5 per cent velocity anomalies to a depth of ∼5 km bss, 3 km × 3 km to a depth of ∼6 km bss and 5 km × 5 km to a depth of ∼7 km bss along the length of the model in both x and y dimensions, and 5 km × 5 km to a depth of ∼10 km bss between 10 and 35 km model offset associated with the deepest travelling ray coverage (Figs S8 and S9). The latter represents ∼6.5 km depth below seabed and, based on standard models of crustal thickness (e.g. White et al. 1992; Grevemeyer et al. 2018b), suggests that the depth and geometry of the Moho should be resolved within a 25 km × 25 km sized region beneath the centre of the South Grid.

3.4 Moho and crustal thickness

To determine the depth and geometry of the Moho throughout the South Grid, and whether it is a distinct interface or transition zone at the seismic wavelength scale of this study, a Moho proxy velocity contour was defined as inversion models resulting from FAST (Zelt & Barton 1998) are smooth and interface-free. Here we chose the 7.5 km s⁻¹ P-wave contour because it matches, in depth, the crustal thickness along SAP_B at 504B determined by Funnell et al. (2021), based on 2-D forward ray-trace modelling of the traveltimes of longer offset uppermost mantle first arrivals (Pn) and the later arriving Moho event, which we assume to be a reflection (PmP) for the purposes of validating our velocity and gravity models.

As a check on the veracity of this proxy, we cut two north-south and two east–west slices through the 3-D P-wave model (along profiles SG_A and SG_C and SG_F and SG_I respectively; Fig. 2b) and from these, created 2-D forward models using the 4.5, 5.5, 6.5 and 7.5 km s⁻¹ velocity contours; depth markers of changes in gradient within the 3-D P-wave model that might be used to define the subdivisions of Layer 2 (lavas, lava-to-dyke transition and dykes from 504B logs), the Layer 2-Layer 3 transition and the base of crust, respectively. Using rayinvr (Zelt & Smith 1992), PmP arrivals were predictively traced through these models and compared to the observed data record sections (Figs 4 and 5). In general, a reasonable qualitative fit was achieved.

As the deeper crustal and uppermost mantle velocity–depth structure is only constrained within a 25 km × 25 km region beneath the centre of the South Grid, the 2-D forward models were then converted into density models to, primarily, determine crustal thickness outside of this region. Typical densities for seawater (1030 kg m⁻³), Layer 2 (2500, 2600 and 2700 kg m⁻³) and Layer 3 (2900 kg m⁻³) were assigned based on Carlson & Raskin (1984), and 3300 kg m⁻³ assigned to the uppermost mantle to be consistent with Funnell et al. (2021). Calculation of the free-air anomaly (FAA), when compared with the ship-derived observed (Hobbs & Peirce 2015) and the global data set (v32) of Sandwell et al. (2014), shows that the apparent thinning that could be interpreted from the shallowing of the 7.5 km s⁻¹ P-wave contour proxy Moho in the ∼5-km-wide region around the periphery of the South Grid is, in fact, an artefact of the ray coverage (cf. Fig. 9 with Figs S8 and S9). In other words, a lack of deep enough travelling arrivals are included in the inversion to provide constraint on crustal thickness outside of the central region that is regarded as well constrained. Instead, modelling of the FAA suggests an effectively constant thickness crust, with an effectively flat-lying Moho at a depth of between 9.4 and 10.2 km bss. Such a flat-lying Moho was then built into the ray-trace models and the PmP predictively recalculated. A good fit between calculated and observed PmP was then achieved (Figs 4 and 5).

Figure 9.

Modelling of the gravity free-air anomaly (FAA). (a) Initial density model constructed along profile SG_C using the velocity contours extracted from the 3-D inversion model. See text for details. Location of 504B is shown by the red triangle. Orange dashed line is the 7.5 km s⁻¹ contour initially used as the base of crust proxy. This model results in the anomaly shown by the orange line in (c) and (d). (b) Preferred density model for profile SG_C. Black dashed lines show the minimum and maximum likely depth of the Moho determined from the 3-D inversion model where constrained around 20 km offset. The average depth was used as the basis of PmP traveltime predictions shown in Figs 4 and 5. (c) Comparison with the ship acquired (black) and Sandwell et al. (2014) v32 (green) FAAs, with grey shading showing the modelling error bound. (d) Difference between observed and calculated. (e, f, g and h) Same for profile SG_I. RMS modelling misfits for both preferred models are annotated.

4. MODEL DESCRIPTION

We follow the standard nomenclature for the layering of the oceanic crust to describe our models. Both the P- and S-wave models reveal a generally consistent crustal thickness and velocity–depth structure throughout the South Grid. Four vertical slices, with one intersecting 504B in the north–south and one in the east–west direction, and four horizontal (constant depth) slices nominally through Layers 2 and 3, demonstrate the main features revealed by the P-wave (Fig. 7) and the S-wave (Fig. 8) tomographic modelling, while corresponding gravity modelling (Fig. 9) indicates that the crust is ∼6 km thick throughout the South Grid region, with a generally flat-lying (within the modelling depth error) Moho. Predictive forward ray tracing calculates arrivals that match clear and distinct PmP arrivals observed on each OBS record section (Figs 4 and 5), regardless of profile orientation or location within the South Grid. Thus, within the seismic wavelength at this depth, based on the 3-D tomographic modelling the Moho would be regarded as a distinct interface, rather than a gradient zone.

4.1 Velocity and structure—1-D

To describe the general velocity layering of each model, we take a 1-D velocity–depth profile through both P- and S-wave models (Fig. 10) at the location of 504B (27.1, 20.4 km in x, y model space), to enable direct comparison not only with the 2-D velocity–depth models of Wilson et al. (2019) and Funnell et al. (2021), but also the sonic logs from 504B (Anderson et al. 1982; Guerin et al. 2008) and 1256D (Wilson et al. 2006; Guerin et al. 2008), located in ∼15 Ma crust formed at the East Pacific Rise during a period of ultrafast spreading, both of which are directly relatable to lithology via downhole sampling. Borehole 1256D sampled ∼0.250 km of sediment, ∼0.750 of basaltic lava flows, a 0.06 km transition zone and ∼0.350 km of sheeted dykes, before first sampling gabbro at ∼1.16 km subsediment (Expedition 309 and 312 Scientists 2006). We also compare our P-wave model with the standard oceanic crust models of White et al. (1992) and Grevemeyer et al. (2018b).

Figure 10.

Comparison of velocity–depth structure derived from traveltime inversions, boreholes and standard envelopes. (a) P- and S-wave 1-D velocity–depth sample at 504B, comparing the results from this study (in 3-D; black) with those of Wilson et al. (2019) (in 2-D; green) and Funnell et al. (2021) (in 2-D; blue). The Layer 2 structure is derived from the lithological logs of 504B, while the base of crust is determined from the gravity modelling (Fig. 9). The vertical grey dashed line marks 7.5 km s⁻¹, used as the proxy for base of crust in the 3-D inversion model (Fig. 7). (b) Velocity–depth profiles from this study compared with the sonic logs (Guerin et al. 2008) of 504B (purple) and 1256D (orange), and the standard envelopes of White et al. (1992) for 0–7 Ma oceanic crust (pink shading) and Grevemeyer et al. (2018b) for Mid-Atlantic Ridge (lilac) and Pacific (turquoise) crust. (c and d) Enlargement of Layer 2 showing the better fit to the borehole-derived sonic logs, for both P-wave (c) and S-wave (d) structure, achieved by the South Grid 3-D crustal modelling. The node locations of the 3-D model are shown by the red dots. See text for discussion. Note that both the P-wave and S-wave vertical velocity gradient change (intra-Layer 2) occurs at the base of the ∼200-m-thick (deep) transition zone identified in 504B. This would be known as the Layer 2A/2B boundary using the traditional ridge-axis-based naming convention, and the Layer 2B/2C boundary off-axis following the convention of, for example, Carlson (2011, 2018).

A limitation of any tomographic model is the constraint imposed by defining the vertical and lateral variation in velocity by the forward cell size (0.1 km in this case), with velocity defined at each node within the regular 3-D mesh, and not in between nodes. In contrast, the sonic logs from 504B and 1256D were reprocessed by Guerin et al. (2008) to a vertical sampling interval of 0.015 km and the lithological logging has a similar resolution. Our P-wave model shows an ∼0.3-km-thick sediment layer of velocity ranging between ∼1.6 and 1.9 km s⁻¹ (gradient 1.0 s⁻¹), bound at its base by a one model node thick velocity step (interpreted as a first-order interface) to 4.8 km s⁻¹ at the top of Layer 2. Layer 2 itself is subdivided into two main units by a vertical velocity gradient change at 4.5 km depth bss (interpreted as a second-order interface), with a gradient of 1.7 s⁻¹ above (4.8–5.8 km s⁻¹) and 0.7 s⁻¹ below (5.8–6.5 km s⁻¹). The base of Layer 2 is defined by a change in vertical velocity gradient at 5.6 km depth bss. Below this, Layer 3 has a velocity range of 6.5–7.5 km s⁻¹ and a gradient of ∼0.3 s⁻¹. On this basis, Layer 2 is interpreted as ∼1.8 km thick and Layer 3 as ∼3.8 km thick. Corresponding S-wave igneous layer velocities and gradients are 2.4–3.1 km s⁻¹ and 1.0 s⁻¹ for Layer 2A, 3.1–3.7 km s⁻¹ and 0.5 s⁻¹ for Layer 2B, and 3.7–4.0 km s⁻¹ and 0.1 s⁻¹ for Layer 3.

Table 2 and Fig. 10(a) show a comparison with the P- and S-wave 2-D velocity models of Funnell et al. (2021) along SAP_B. The latter models are fixed at the sediment-basement interface based on a sediment layer thickness measured from the coincident MCS record section and a velocity derived from its stacking velocity and sonic logs run in 504B, while the models derived from this study are fixed above the seabed. All of the models are parametrized with a forward cell size of 0.1 km horizontally and vertically. We adopted this approach for two reasons:

1. To determine if an inversion with a forward cell size similar to the sediment thickness could resolve its dimensions and velocity.

2. What the effect and consequence are of fixing a model multiple nodes-deep beneath an inversion ray source.

Table 2.

P-wave and S-wave layer velocities and gradients for the South Grid 3-D tomographic volume compared with those for 2-D inversion of SAP_B (Funnell et al. 2021), both sampled at 504B.

This study						Funnell et al. (2021)
	P-wave velocity and gradient		S-wave velocity and gradient		Thickness	P-wave velocity and gradient		S-wave velocity and gradient
	(km s−1)	(s−1)	(km s−1)	(s−1)	(km)	(km s−1)	(s−1)	(km s−1)	(s−1)
Sediment	1.6–1.9	1.0	1.6	-	0.3	-	-	-	-
Layer 2–upper	4.8–5.8	1.7	2.4–3.1	1.0	0.7	5.0–5.5	0.7	2.7–3.1	0.6
Layer 2–lower	5.8–6.5	0.7	3.1–3.7	0.5	1.1	5.5–6.5	0.9	3.1–3.7	0.5
Layer 3	6.5–7.5	0.3	3.7–4.0	0.1	3.8	6.5–7.8	0.3	3.7–4.4	0.2
Crustal thickness					5.9

Although all achieving a χ² fit of ∼1, there are differences in outcome, most notably for the SAP_B models that show a higher-to-lower velocity trade-off between the upper and lower parts of Layer 2 (Figs 10a, c and d). However, the best fit to the velocity and gradient of the sonic logs is achieved by solely fixing the seabed and allowing the model freedom to vary in velocity immediately beneath the ray source (immediately under the seabed beneath an OBS position in this case), with the vertical offset between observed and calculated trends dictated by where the fixed surface lies in relation to the depth of a node. Not only is the velocity of the sediment layer recovered, so is its thickness to a resolution of the forward node spacing in depth.

The 1-D velocity–depth models of this study (Fig. 10b) show a remarkably good fit to the 504B and 1256D sonic logs of oceanic Layer 2 structure for both P- and S-wave velocity, and also recover a Layer 3 structure that is consistent with standard model envelopes (e.g. White et al. 1992; Grevemeyer et al. 2018b). The primary uncertainty, however, lies in the base of crust determination. Although a change in gradient could be interpreted from the 3-D tomographic models to mark the deepest the Moho could be, we conclude that inversion alone is insufficient to make a reliable determination and, instead, forward modelling of second arriving PmP phases is required. Such ray-trace modelling (Figs 4 and 5) shows that the Moho lies at a depth of ∼10 km bss beneath the central region of the South Grid, and further gravity modelling suggests that it is generally flat-lying at this depth throughout the entire 45 km × 45 km study region.

4.2 Velocity and structure—3-D

The correspondence in depth of the change in vertical velocity gradient with the main lithological layering subdivisions logged within 504B (Fig. 10), allows calculation of the velocity gradient throughout the 3-D tomographic volume to be used as a proxy to investigate the lateral variation in structure and thickness of these primary seismic layers within the 3-D volume. Inherently crust formed at spreading ridges would also be expected to have a ridge-parallel fabric, for example due to tectonism associated with spreading. Therefore, any change in spreading style [magma-dominated (magma rich) to tectonism-dominated (magma poor)] or orientation may also be manifest in crustal structure.

As might be expected, variation in horizontal gradient (both P-wave and S-wave) is only observed in the flowline direction (Figs 11c and e), mirroring the seabed bathymetry highs and lows (Fig. 11a) which, in turn, correspond to basement highs and lows beneath which coincident MCS record sections image (Figs 3a–c). However, laterally, the strike direction of the gradient variation is not CRR-parallel, but instead lies at an azimuth of 8–10° to that, following the current spreading direction of the ER (Fig. 1). Peirce et al. (2023b) investigate the structure and dynamics of the Ecuador Fracture Zone (EFZ; Fig. 1), a system comprising three transform-fractures with two short ridge segments captured in between, which separates the CRR and ER. In this system, the most recent eastern transform offset to form ∼1.5 Ma, is orientated CRR-perpendicular, 8–10° anticlockwise to the ER-trend of remainder of the EFZ. We interpret the gradient trend as reflecting the CRR spreading direction when the crust of the South Grid was formed ∼7 Ma which was, thus, 8–10° clockwise from that observed at the present-day.

Figure 11.

Slices through the 3-D P- and S-wave vertical and lateral velocity gradient models. The lateral gradient was calculated in the north-south flowline direction. (a) Bathymetry map showing the location of 504B (solid yellow circle). Black dashed lines in (a) and all horizontal (depth) slices indicate the locations of the vertical slices through the models (profiles SG_C and SG_I). The ridge-strike directions of the Costa Rica Rift (CRR) and Ecuador Rift (ER) are annotated (arrows). Red dashed line shows the 7 Ma magmatic-to-tectonic-dominated spreading transition of Wilson et al. (2019), based on magnetic anomaly modelling. (b and f) Vertical slices through the P-wave vertical gradient model along profiles (b) SG_C (CRR-parallel) and (f) SG_I (CRR-flowline). 1-D velocity–depth sample (white solid line) and the lithological layering at 504B (horizontal white dash) are annotated, together with its location (vertical white dash). (c) Depth slice through the lateral P-wave horizontal gradient model at 5.5 km bss, (mid-Layer 2B). (d and g) As (b) and (d) for the S-wave model. (e) As (c) for the S-wave model. Note that both the P- and S-wave lateral gradient models show a fabric orientation parallel to the current spreading direction of the ER and not the present-day CRR at which the crust formed, suggesting a plate reorganization has occurred in the ∼6.9 Myr since this crust formed at the CRR spreading axis (Peirce et al. 2023b). In all model parts, the slices are masked by the inversion model ray coverage.

An alternative way to investigate lateral variation is to calculate the difference between each of the P- and S-wave inversion models and the corresponding P- and S-wave 1-D velocity–depth profiles (Fig. 10) sampled at the location of 504B (27.1, 20.4 km in x, y model space). These 1-D references were draped beneath the seabed interface to create the 3-D reference models which were then subtracted from the corresponding P- and S-wave inversion models (Figs 7 and 8). The resulting difference models are shown in Figs S10 and S11. For the parts of the inversion models that are considered well constrained, what little variation that is revealed (±0.3 km s⁻¹) follows, unsurprisingly, the same 8–10° anticlockwise ER-trend.

The sediment, Layer 2 and Layer 3 boundaries are evident in both the P- and S-wave vertical velocity gradients, and suggest that there are no significant variations in thickness or velocity structure in any layer either ridge-strike (Figs 11f and g) or ridge-parallel (Figs 11b and d). Consequently, the conclusion to be drawn from both the gradient and difference investigations is that the oceanic crust throughout the South Grid is rather uniform in layer thickness and velocity structure, both ridge-parallel and ridge-perpendicular, making the location of 504B a good choice as a representative example of oceanic crustal lithology.

4.3 Crustal accretion and ridge structure inheritance

The crustal structure imaged at the South Grid represents a snapshot of crustal formation occurring at the CRR ∼7 Myr ago. This crust has been aged (porosity and permeability infilling) by hydrothermal circulation as it has spread off-axis and has been covered by a ∼275-m-thick layer of sediment that should, eventually, seal fluid flow pathways (e.g. Wilson et al. 2019). However, MCS seismic imaging (Fig. 3 and Peirce et al. 2023b) shows that, even in crust now located 250–300 km off-axis, active faulting is still taking place, with basement block-bounding faults observed to propagate through the sediment layer to the present-day seabed. Consequently, the legacy of structure embedded in the crust at its formation still plays an active role in determining its properties and behaviour millions of years later.

At the present-day, an active overlapping spreading centre (OSC) forms a non-transform offset at the CRR that divides the ridge axis in two (Fig. 12). Robinson et al.’s (2020) modelling of the North Grid shows that the OSC is associated with a P-wave velocity low in Layer 2 that is 10-km-long along-axis and 5-km-wide across-axis, which they interpret to reflect a region of higher bulk porosity and permeability due to the fracturing and faulting associated OSC opening and propagation. Robinson et al. (2020) further postulate that these fracture pathways promote hydrothermal circulation that rapidly cools the ridge axis and any accumulated magma. Thus, the velocity and thickness characteristics consequent of magma rich or magma poor crustal formation are set at zero-age and, as Wilson et al. (2019) and Funnell et al. (2021) show, are preserved and carried in the crust as it spreads off-axis. However, neither P-wave (Figs 7 and S10) or S-wave (Figs 8 and S11) velocity or difference models show (within their resolution) any relic of the OSC which leads to the conclusion that either the OSC is younger than 6–7 Myr or that it has propagated from a more eastwards location outside the South Grid footprint over that time frame.

Figure 12.

Comparison between the on-axis North Grid (top) and off-axis South Grid (bottom) locations relative to present-day tectonic features along the CRR ridge axis. The CRR is shown by the red dashed line, together with the overlapping spreading centre (OSC) between two of its constituent ridge segments. Wilson et al.’s (2019) crustal ages, based on magnetic anomaly modelling, are marked by black dashed lines. From analysis of the spreading rate Wilson et al. (2019) conclude that, between 6 and 7 Ma, crust was forming at the CRR under magma rich conditions.

Wilson et al.’s (2019) magnetic modelling along SAP_A and SAP_C identifies phases of magma rich and magma poor spreading over 7 Myr of spreading at the CRR and predicts that the crust of the South Grid formed during a magma rich phase (Fig. 12). The Vp/Vs ratio can be used as a proxy to distinguish between crust formed under magma rich and magma poor conditions, with Grevemeyer et al.(2018a) proposing a Vp/Vs ratio of <1.9 to indicate magma poor and >1.9 magma rich crustal formation while, alternatively, Peirce et al. (2019, 2020, 2023a) propose 1.85. Here we set the magma poor/magma rich discriminator at 1.85, and note that this corresponds to a Poisson's ratio (σ) of 0.29. The Vp/Vs ratio for the South Grid (Fig. S12) is consistent with Wilson et al.’s (2019) predictions in the flowline direction and also suggests, given the lack of ridge-parallel variability (within the resolution), that magma supply to the ridge-axis was also consistent along-axis. Consequently, it is reasonable to conclude that magma supply and spreading processes were most likely consistent over the ∼1 Myr sampled by the South Grid survey.

5. SMOOTH VERSUS HETEROGENEOUS MODELS

The 3-D tomography models, particularly the P-wave, show a change in vertical gradient velocity associated with the transition from lavas to feeder dykes, and such a change in velocity gradient has been used to determine not only Layer 2A thickness, but also how Layer 2 as a whole evolves as newly formed oceanic crust ages as it spreads off-axis (e.g. Audhkhasi & Singh 2019; Estep et al. 2019). At 504B, the transition between lavas and dykes is not interpreted as a distinct interface lithologically (Carlson 2011); instead it is viewed as a transition zone within which the percentage of dykes increases with depth which has, in turn, brought about the alternative, off-axis, terminology for the layering as 2A (lavas), 2B (transition zone) and 2C (dykes) (Houtz 1976; Houtz & Ewing 1976).

The lava-dyke transition is associated with an amplitude anomaly within individual MCS gathers, which is variously called the 2A Reflection, 2A Retrograde or 2A Caustic (Harding et al. 1993; Christeson et al. 2007; Arnulf et al. 2012 respectively). To avoid pre-empting this paper's interpretation of the cause of this anomaly, we will simply refer to it as the 2A Event. By optimising the processing of MCS data, this amplitude anomaly can be ‘imaged’ on the resulting MCS record section, where it has been used to infer the location of the boundary between the lavas and dykes and/or an alteration front, with authors citing evidence for both interpretations (e.g. Harding et al. 1993; Christeson et al. 2010; Newman et al. 2011; Audhkhasi & Singh 2019; Estep et al. 2019). Although it is generally agreed that the 2A Event is associated with the velocity change at the 2A/2B transition, the exact mechanism of genesis remains undetermined. The 2A Event is consistently observed across the north-south extent of the South Grid along the SAP MCS profiles (Fig. 3e).

The ambiguity in cause is further fuelled as the 2A Event's location (two-way traveltime and lateral geometry) on an MCS record section is also a function of the velocity model used for the normal moveout correction as part of stacking (Harding et al. 1993; Christeson et al. 1996). In oceanic crust typically devoid of any internal reflectivity with which to undertake velocity analysis prior to stacking, velocity models from coincident refraction data, if they exist, are often used. If not, stacking has a more qualitative basis of improving visual appearance (Figs 3a–c). MCS data along profile SAP_B were acquired using a synthetic aperture approach (Buhl et al. 1982) to produce receiver gather offsets up to 8700 m, with the specific purpose of imaging the 2A Event, which modelling based on 504B predicted would appear within individual gathers at ∼5500 m shot-receiver offset. Profile SAP_B is coincident with profile SG_I (maximum MCS shot-receiver offset 4500 m) where it runs through the South Grid (Fig. 2).

Here, the optimum stacking velocity for the 2A Event (2200 m s⁻¹ in the case of the SAP_B—Fig. S13), that gives the strongest response on an MCS profile, is not what would be predicted based on a refraction-derived velocity model due to the non-hyperbolic behaviour of this arrival (Harding et al. 1993). More recently, full-waveform inversion (FWI), which simultaneously matches the traveltime, amplitude and phase of the arrivals on a seismogram (Christeson et al. 1996, 2012; Arnulf et al. 2012), has been used with both 1-D and 2-D smooth velocity models to show a correlation between the 2A Event and the base of a rapid increase in velocity.

Such smooth models are also at odds with the observed velocity variation with depth in well logs (e.g. Anderson et al. 1982; Adamson 1985; Gilbert & Salisbury 2011; Chadwick et al. 2013; Ingale & Singh 2021) and geological observations from ophiolites (Kidd 1977; Pallister 1981) that show significant velocity heterogeneity, especially so in the upper lava layer (Fig. 10). The mismatch can be visualized using a scattering regimes diagram (Fig. S14; Wu & Aki 1988). Ray tracing imposes a simple model that can be described by the physics of geometric optics (simple refraction and reflection) whereas, in reality, the geological model suggests that the oceanic crust lies within the weak-to-strong scattering regimes. This is particularly true for the Layer 2A lavas which, due to their episodic emplacement, have a complex velocity structure depending on the lithology and internal structure of the flows (Planke 1994; Single & Jerram 2004; Nelson et al. 2009). Thus, a propagating seismic wave will be scattered by this heterogeneity as it is small in size relative to the dominant signal wavelength and the total length of the propagation path. On the scattering regime diagram such a propagation path lies outside of the region where geometric optics is valid (Fig. S14), and upon which all ray tracing-based modelling is founded. This leads to two primary questions:

1. How applicable are the smooth velocity models derived from WA data forward modelling and inversion for stacking MCS data acquired in oceanic crustal settings?

2. What is the actual origin and cause of the arrival being stacked into MCS record sections as the 2A Event?

To answer the first question, we start by using a 1-D sample at 504B in the 3-D tomography model (Figs 2, 7 and 10) to create a laterally consistent 2-D model (henceforth the smooth model). We then create a laterally and vertically variable model (henceforth heterogeneous model) by perturbing the smooth model to simulate the effects of scattering caused by the lavas (Layer 2A), the feeder dykes beneath (Layer 2B), and an irregular boundary between them. We then compute the seismic response of both models using the finite-difference wave propagation modelling programme SOFI2D (Saenger & Bohlen 2004) to create synthetic receiver gathers to compare with those acquired along profiles SG_I and SAP_B and identify arrivals. Finally, we invert the synthetic first arrival traveltimes for both smooth and heterogeneous models in the same manner as the real WA data, and compare the resulting inversion models with a slice along profile SG_I through the South Grid 3-D tomography model (henceforth real 3-D model) that is based on observed arrivals and an equivalent inversion in 2-D of the real traveltime picks along profile SG_I only (henceforth real 2-D model). To answer the second question, we analyse snapshots in time of wave front propagation through both the smooth and heterogeneous models, from which we identify wave propagation modes and where the 2A Event forms.

5.1 Model construction

To investigate the effect of small-scale heterogeneity, first we build a simple structural framework for the six primary layers [seawater, sediment (Layer 1), lavas (Layer 2A), dykes (Layer 2B), gabbro (Layer 3) and mantle] using a slice along profile SG_I (Fig. 2) through the P-wave 3-D velocity model of this study (Fig. 7), which is consistent with the P-wave 2-D inversion models of Wilson et al. (2019) and Funnell et al. (2021) along SAP_B (Fig. 6). The corresponding S-wave velocity (V_S) and density (ρ) models were derived using the Brocher (2005) relationships:

$$\begin{eqnarray} {{{V}}}_{{S}}{{ = 1000 \times }}\left( {{{0}}{{.7858 - 1}}{{.2344}}{{{V}}}_{{P}}{{\ +\ 0}}{{.7949V}}_{{P}}^{{2}}{{\ -\ 0}}{{.1238V}}_{{P}}^{{3}}{{\ +\ 0}}{{.0064V}}_{{P}}^{{4}}} \right) \end{eqnarray}$$

$$\begin{eqnarray} {{\rho}} = {{ 1000\ \times }}\left( {{{1}}{{.6612}}{{{V}}}_{{P}}{{\ - \ 0}}{{.4721V}}_{{P}}^{{2}}{{ \ + \ 0}}{{.0671V}}_{{P}}^{{3}}} {{{\ - \ 0}}{{.0043V}}_{{P}}^{{4}}{{\ + \ 0}}{{.000106V}}_{{P}}^{{5}}} \right), \end{eqnarray}$$

where the P-wave velocity (V_P) is defined in km s⁻¹. The models are discretized with a node spacing of 2 m both horizontally and vertically. This spacing is chosen as a compromise between being sufficient to capture the heterogeneity, avoiding numerical dispersion in the finite difference simulation, and limiting the computational resources necessary to calculate the seismic response using SOFI2D (Saenger & Bohlen 2004).

For the smooth model (Fig. 13a), the P-wave velocity is extrapolated from a 1-D sample, at the location of 504B, through the 3-D tomographic inversion of the observed traveltime picks from the South Grid, and the corresponding S-wave velocity and density are computed using the relationships above. The S-wave velocity of seawater is set to zero and its density set to 1000 kg m⁻³. The seawater and sediment layers are assigned P-wave velocities of 1.5 and 1.665 km s⁻¹, respectively, and the interface (the seabed) between the two layers is set to be planar at 3.5 km depth below sea surface (bss); the average seabed depth around 504B. The Layer 3/mantle boundary is modelled as a gradient zone at a depth of 9 km to 10 km bss, based on the South Grid modelling of this study. Besides the steps in the velocity and density at the seabed and base of sediment, the crustal velocity structure has a monotonic increase of physical properties with depth with the layers defined by changes in gradient or P-wave velocity.

Figure 13.

Synthetic models used to determine the origin of the 2A Event. See text for a description of their construction. (a) Smooth P-wave velocity model based on a 1-D sample of the 3-D P-wave inversion model at 504B. (b) Comparison between 1-D velocity–depth profiles from the 3-D P- (black solid) and S-wave models (black dashed), 504B (purple), 1256D (orange), the smooth model [(a); blue] and the heterogeneous model [(c); red]. (c) Heterogeneous P-wave velocity model, with black dashed boxes showing model selections shown in (d), (f) and (h) as annotated. (e) Heterogeneous S-wave velocity model. (g) Heterogeneous density model.

The heterogenous model (Figs 13c, e and g) is built in a top down fashion by first setting the seawater and sediment layers to have the same homogeneous characteristics as the smooth model. Layer 2A (Fig. 13d) is then built as a series of successive units where each is randomly assigned a lithology, thickness and surface topography (Table S1) based on Nelson et al. (2015), who used 3-D mapping and well logs to synthesise the internal structure of flood basalts along the northwest margin of Europe, and multibeam surveys of ocean floor lava flows (Chadwick et al. 2013, 2016). The percentage of each type of lithology in the layer is adjusted, together with the velocity, until a vertical velocity variation similar to that in 504B logs (Anderson et al. 1982; Alt et al. 1996) is achieved. Note that there is no direct lithological correlation between model and samples taken from 504B as there was only 20 per cent core recovery, so most of the actual lithology is unknown.

Layer 2B (Fig. 13f), the dykes, is constructed by varying the depth to the top, the width and velocity of each vertical dyke (Table S1). The depth of the dyke top has a log-normal distribution that is scaled such that the shallowest dyke tops reach the top of Layer 2A, and the mean depth is set to a pre-assigned value of 4500 m which is 1000 m below seafloor based on the velocity inversions from this study and Detrick et al. (1998). The dyke width has a Gaussian distribution with a mean of 10 m and a standard deviation of 4 m. The model node spacing precludes dykes narrower than 2 m, even though this is wider than those mapped by Pallister (1981) in the Samail ophiolite in Oman and from studies of dykes on Iceland (Gudmundsson 1984; Galindo & Gudmundsson 2012), where the mean is less than 1 m. However, individual intrusions narrower than 2 m penetrating into Layer 2A would cool rapidly and so be more likely to be fractured, increasing their porosity and lowering their seismic velocity (Wilkens et al. 1991), so we have assumed that they would have similar properties to that of the heterogeneity of Layer 2A. The assigned dyke velocity has a mean of 6100 m s⁻¹, which is perturbed by the depth of its top such that dykes that intrude to a shallower subsurface depth have a lower velocity (proxy for lower density and higher porosity; Galindo & Gudmundsson 2012). A vertical velocity gradient of 0.3 s⁻¹ is imposed (Carlson 2001) to match the trend in the sonic log from 504B, to which a further 3 per cent Gaussian perturbation was added to mirror heterogeneity in dyke composition (Becker et al. 1989). Although the resultant representation of Layer 2 is still a simplification of field observations, we believe it captures the important geophysical characteristics of a distributed and variable velocity structure, both laterally and vertically, that averages to produce the velocity gradient within Layers 2A and 2B, and the change in gradient at the boundary between them, observed in both the South Grid and SAP_B velocity–depth models (and hence the smooth model; Fig. 13b).

For Layer 3 (Fig. 13f), we adopt a multiple sill-like unit construction for modelling simplicity with a weak velocity perturbation (Table S1). The transition between dykes (Layer 2B) and gabbro (Layer 3) is not well constrained and not fully sampled by 504B, although isolated gabbroic intrusions are observed in the vicinity of the base of the borehole (Wilson et al. 2006). We have, therefore, assumed that the Layer 2/Layer 3 transition occurs over 1000 m depth between 5500 and 6500 m, with the probability of a sill-like gabbro intrusion linearly increasing from zero at the top to one at the bottom of this interval.

Based on studies of ophiolites and models inferred from seismic reflection images, the transition between Layer 3 and the lithospheric mantle (the Moho) may not be an abrupt interface, but may instead occur over a finite thickness (e.g. Karson et al. 1984; Brocher et al. 1985; Jousselin & Nicolas 2000; Nedimović et al. 2005; Ildefonse et al. 2010). As processing of OSCAR MCS reflection data did not reveal any reflectivity above the background noise at the anticipated Moho two-way traveltime (Fig 3d), and the 3-D South Grid model reveals no unambiguous step in velocity or significant change in velocity gradient at Moho depth (Figs 10 and 11), we conclude that this is the case around 504B. Thus, like the Layer 2/Layer 3 transition, the Layer 3/mantle boundary is constructed with an increasing density of sill-like units (Fig. 13h) within a 1000 m transition zone between 9000 and 10 000 m depth, each with a velocity randomly assigned from a normal distribution with a mean of 8.0 km s⁻¹ and standard variation of 0.05 km s⁻¹ which results in a velocity gradient average of 0.75 s⁻¹, similar to the hybrid of models presented by Ingale & Singh (2021). Once below the transition zone, the mean velocity was fixed at 8.0 km s⁻¹.

An example from the heterogeneous model set is shown in Fig. 13. In all but the seawater and sediment layers there are multiple complex interfaces, some with significant impedance contrasts that will produce a strong scattering response. Overall, though, its layer velocities are consistent with those compiled by Christeson et al. (2019), the OSCAR WA models for SAP_B (Fig. 6) and the South Grid (Fig. 7), as well as the sonic logs from 504B (Fig. 13b). However, it is also significantly heterogeneous when compared to the smooth model equivalent (Fig. 13a) input as the starting point or resulting from forward ray tracing or tomographic inversion.

Thus, a suite of models was generated based on the range of properties described above and the resulting simulated wavefield for each was compared to the surface and seabed observations made during OSCAR around the site of 504B. The model included here (Fig. 13) is one of a subset which result in the best match to the 2A Event observed at 5–7 km offset in the MCS data (Fig. 14), that occurs behind the 2A and 2B refracted arrivals.

Figure 14.

Synthetic MCS gathers calculated for the smooth and heterogeneous models (Fig. 13). (a) Observed 8.7 km streamer supergather from profile SAP_B at the location of 504B and, adjacent to the right, overlain by arrivals calculated using the heterogeneous model. See legend for identification. Black line is the 2200 m s⁻¹ normal moveout curve that fits the 2A Event (solid red line). (b) As (a) for the observed 4.5 km streamer gather from profile SG_I. (c) Synthetic gather calculated based on the heterogeneous model, with phases annotated. The vertical black dashed line marks the 8.7 km maximum observed supergather offset (a). (d) Heterogeneous model synthetic gather annotated to show where the Layer 2A refracted phase emerges as a first arrival at ∼5.0–5.5 km offset, and the incoherent scattered energy throughout the gather caused by the subwavelength-scale heterogeneity. (e) Smooth model synthetic gather annotated to show the arrival of top-bottom sediment layer pegleg multiples. Note the lack of scattered energy and 2A Event. All gathers are plotted with a reduction velocity of 6.0 km s⁻¹.

5.2 Seismic wavefield simulation

Seismic waveform propagation through both the smooth and heterogeneous models was undertaken with SOFI2D which uses a finite difference approach (Saenger & Bohlen 2004). Each simulation took 2 hr to run on a 120-node high performance computer cluster. The run parameters are summarized in Table S2. The seabed depth of 3500 m meant that free-surface (sea surface) multiples would not interfere with the primary wavefield and were, thus, not included in the simulation. Edge effects from the sides and base of the model were suppressed by using perfectly matched layers 50 nodes wide (Komatitsch & Tromp 2003). Modelling was set to be elastic so that attenuation of the higher frequencies was solely due to scattering (Maresh et al. 2006; Shaw et al. 2008).

The acoustic source is a zero phase Ricker wavelet with a centre frequency of 10 Hz. The grid node spacing of 2 m is more than sufficient to avoid numerical dispersion for the slowest S-wave velocity in the igneous crust (0.9 km s⁻¹). However, it is marginal for the sediment layer although visible dispersion artefacts are not observed in the simulation because the most significant P-to-S-wave conversion occurs at the sediment/basement boundary, such that S-wave energy is largely trapped in the igneous part of the crust and propagation in the sediment layer is dominated by P-waves, which are not dispersive with a 2 m grid node spacing.

To visualize seismic wave front propagation and aid analysis of the effect of scattering and interaction with the various interfaces, snapshots of the simulation were taken every 50 ms, the pressure wavefield being the most useful, as this is what is recorded by hydrophones and is most easily compared with the field recordings. Snapshots at 3.2 s for both the smooth and homogeneous models are shown in Fig. 15. Finally, the succession of snapshots was merged into a movie sequence for each model (Figs SA1 and SA2).

Figure 15.

Synthetic wave propagation at 3.2 s after shot. (a) Smooth model. (b) Model masked by phases as annotated. (c) Wave propagation showing phase amplitude. (d, e and f) Same as (a), (b) and (c) for the heterogeneous model. Note that the 2A Event is created by the geological heterogeneity in the crust, and originates and propagates solely in the lower part of Layer 2A. Time-lapse views of seismic wave propagation through the smooth and rough models are shown in Figs SA1 and SA2. Red star marks the transition between the 2A Event and the refracted wave in Layer 2B.

5.3 Synthetic versus observed gathers

To simulate MCS gathers both the shot and receivers were located at the surface of both models and the direct water-wave, and refracted and reflected arrivals calculated for all layers in each model (Figs 14d and e). The similarity between the synthetic gathers calculated from both the smooth and heterogeneous models and observed gathers (Fig. 14) is remarkable, and demonstrates that both lithology-based models can produce seismic data that is geophysically plausible. However, the smooth model (Fig. 14e) does not result in an arrival from Layer 2A, but does create a series of pegleg multiples generated within the sediment layer (repeated reflection between the top and bottom of this layer), It also demonstrates how the scatter generated by the small-scale heterogeneity adds what might be regarded as background noise to the heterogeneous model-derived synthetic gather (Fig. 14d).

The arrivals calculated from the heterogeneous model show a remarkable consistency with the observed, both for the ∼8.7 km gathers of profile SAP_B (Fig. 14a) and the ∼4.5 km gathers of profile SG_I (Fig. 14b). Due to the water depth, the refraction from Layer 2A is only just recorded by the longest offset traces in the profile SG_I gather and, apart from the seabed reflection (Top Layer 1), no other arrivals are recorded by this length of streamer. However, the profile SAP_B longer streamer records the Layer 2B refraction first arrival and a second arrival that arrives just after it. This arrival is what is generally stacked as the 2A Event.

5.4 Traveltime picks and inversion

Having established the viability of model construction, a second set of gathers was calculated in which the receiver was located at 3.5 km depth (the seabed) and the shots located at the surface to simulate an OBS record section. The traveltimes of the first arrivals were picked (Fig. 16), and mirrored about each OBS position along profile SG_I. A comparison between the various traveltime pick sets is shown in Fig 16, where those for OBS SG_12 along profile SG_I (Fig. 2) are used as the real data reference, as this OBS is located closest to 504B and this profile runs through 504B in a north–south direction (flowline). Examining the P-wave arrivals the following observations can be made:

1. Both sets of synthetic picks match the observed (Figs 16a and b).

2. The smooth model generates arrivals that are generally earlier (faster) than the heterogeneous model equivalent (Fig. 16c; red line).

3. Both the heterogeneous and smooth models generate picks that generally arrive after the observed P-wave picks (Fig. 16c; orange and blue lines).

4. Almost all of the mismatches are within the observed (real) arrival pick uncertainty (±40 ms) and are, thus, equivalent.

Figure 16.

Comparison between observed and synthetic traveltime picks. (a) Smooth model synthetic OBS gather annotated with observed P- and S-wave picks [blue; cf. (e)] and calculated picks (red). (b) Heterogeneous model synthetic OBS gather annotated with observed P- and S-wave picks [blue; cf. (e)] and calculated picks (red). (c) Differences between pick sets relative to the observed P-wave traveltime pick uncertainty of 40 ms. Note that, generally, the heterogeneous model generates picks that arrive later than the smooth model equivalent and that the observed picks arrive earlier than either synthetic picks. However, the synthetic picks generally arrive within the observed pick uncertainty and are, thus, equivalent within the error. (d) Observed record section for OBS SG_12, located closest to 504B, and shots from profile SG_I (Fig. 2) for comparison. (e) As (d) with observed traveltime picks annotated (blue).

The later S-wave arrivals show a more noticeable difference between the synthetic smooth and heterogeneous models (Figs 16a and b). The modelled S-wave arrival for the heterogeneous model is consistent with the actual traveltime picks from OBS SG_12, whereas the equivalent from the smooth model is consistently early and this mismatch increases with offset even though the same P-to-S-wave velocity relationship was used for both models. The reason for this mismatch is the stronger interaction of the shorter-wavelength S-waves with the complex velocity structure in the heterogeneous model that is absent in the corresponding smooth model, which increases with the length of the wave's travel path through the model. Both the smooth and heterogeneous model pick sets, together with the observed pick sets for profile SG_I (OBSs SG_02, SG_09, SG_12, SG-19 & SG-22; Fig. 2) were then tomographically inverted in 2-D using FAST (Zelt & Smith 1992; Zelt & Barton 1998) with the initial model (Figs 7b and 8b) and inversion parameters (Table 1) set equivalent to those used for the 3-D tomographic inversion of the South Grid.

The primary difference between the synthetic and observed picks sets is the maximum shot-receiver offset. For the synthetic picks this is ∼37.5 km; sufficient in theory to constrain the base of crust as the pick sets contain upper mantle (Pn) arrivals as first arrivals. For the observed picks the maximum offset along profile SG_I is ∼29 km, but only for SG_22 (Fig. 2) which is the only OBS to record Pn arrivals given the fixed profile length and the various OBS positions along profile. However, comparison between inversion models (Figs 17a, d and g) shows what role these longer offset arrivals play in constraining the crust below ∼6 km bss. Even then, it is also clear that both the smooth and heterogeneous model inversions only constrain the crust deeper than ∼6 km bss between ∼10 and 35 km model offset, an outcome not dissimilar to that of the real 3-D model, while the real 2-D model has little-to-no constraint below ∼6 km bss. This is shown by calculating the difference between the real 3-D and real 2-D models (Fig. 17h). Perhaps the most striking outcome is demonstrated by calculation of the difference between the smooth and heterogeneous inversion models and that obtained for the real 2-D inversion (Figs 17c and f). In effect, there is no difference between real and heterogeneous models (Fig. 17f) shallower than ∼5 km bss (Layer 2A), in contrast to the difference between the real and smooth models (Fig. 17c). This result implies that although the real 2-D model would be regarded as smooth, or smoothed, relative to the scale of the heterogeneity, it nevertheless contains information on that heterogeneity via averaging of it at the seismic wavelength scale.

Figure 17.

Comparison between 2-D inversion models based on the synthetic and observed (real) traveltime picks. (a) Smooth model picks. (b) Difference between the smooth model inversion and a slice through the 3-D South Grid (real) model. (c) Difference between the smooth model inversion and a 2-D inversion model (g) using only the observed (real) picks from OBSs located along profile SG_I (OBSs SG_02, SG_09, SG_12, SG_19, SG_22; Fig. 2). (d, e and f) As (a), (b) and (c) for the heterogeneous model. (h) Difference between inversions based on observed traveltime picks undertaken in 2-D (profile SG_I OBS picks only) and 3-D (using all South Grid OBS picks with the 3-D model sliced along profile SG_I). For parts (a), (d) and (g), the 5.0, 6.0, 7.0 and 7.5 km s⁻¹ velocity contours are annotated as solid black lines. In all parts, the 5.0, 6.0, 7.0 and 7.5 km s⁻¹ velocity contours from the 3-D model slice along profile SG_I are annotated as dashed black lines. A good correlation is observed between contour sets, except for the 7.0 and 7.5 km s⁻¹ contours for the observed (real) 2-D and 3-D models, where the 2-D contours lie ∼1 km shallower, indicating that the maximum shot-receiver offset inline of each profile is insufficient to constrain the lower crust and below. The constraint at these depths achieved by the 3-D inversion, results from picks obtained largely from profiles SG_M and SG_N (Fig. 2).

5.5 Wave propagation

Comparison of snapshots taken at 3.2 s after the shot instant (Fig. 15) shows that the subhorizontal heterogeneous structure of Layer 2A in the heterogeneous model creates a leaky waveguide that channels seismic energy across the model that scatters back into the overlying water layer. This mode of propagation can result in en-echelon arrivals (Sanford et al. 2018) if the source frequency is high enough to distinguish alternating zones of a higher or lower average velocity. Leaky waveguides have been modelled to demonstrate rapid attenuation of first arrivals in gas hydrates (Zanoth et al. 2007), whereby energy is scattered out of the leading wave front into lower velocity zones where propagation is slower. This observation was first reported in young ocean crust by Lewis & Jung (1989). At 3.2 s a high-amplitude arrival has formed in the lower part of Layer 2A (above the mean depth to the top of the vertical dykes; Figs 15e and f). In the depth range of 4.2–4.5 km the wave front is delayed due to the presence of a lower velocity channel in Layer 2A, even though there is an increase in the average velocity due to the increasing number of higher velocity dykes (Fig. 13b), and its amplitude is locally increased (Fig. 15f). With increasing elapsed time (Fig. SA2), this amplitude anomaly propagates upwards until it arrives at the receiver array between 5 and 8 km offset (Fig. 14c), where it forms a secondary arrival and appears as a fold in the wave front.

Although wave propagation modelling through both the smooth and heterogeneous models predicts the observed traveltimes (Fig. 16) within their uncertainty, there are some notable differences in modes of propagation. The first is the nature of the wave propagation in Layer 2A. In the heterogeneous model the wave front propagates horizontally, guided by the horizontal layering, with sections advanced in time where there is a concentration of higher velocity layers and retarded where there is a concentration of lower velocity layers. Here, the average velocity of the model has a low gradient of 0.7 s⁻¹ that increases the velocity from 4.5 to 4.95 km s⁻¹ over a thickness of 640 m. Below this, the heterogeneous model has a 200-m-thick layer of higher gradient of 5.25 s⁻¹ that increases the velocity from 4.95 to 6.0 km s⁻¹ due to the increasing number of dykes. This generates the caustic in the modelled traveltimes, here termed the 2A Event (Fig. 14c). Then the arrivals (Layer 2B refraction in Fig. 14c) from the higher velocity dyke layer, Layer 2B, overtake the 2A Event and form the first arrival out to offsets of 22 km. The second notable difference is the significant amplitude anomaly second arrival between 20 and 28 km offset. This is a second folded wave front created by the velocity gradient in the crust–mantle transition zone. Traditionally, in WA seismic forward ray tracing, this arrival is modelled as a critical angle PmP (Moho) reflection caused by an abrupt step change in velocity (first-order interface).

The primary observation is that the 2A Event is generated and propagates exclusively in the lower section of Layer 2A. In our models its cause is a combination of strong scattering encountered by the wave front where the increasing number of higher velocity dykes intrude the lower velocity lavas which delays the peak energy, coupled with energy scattered upwards from the sharp velocity change generated by the dyke tops. Thus, equating the 2A Event to the top of the transition from lavas to dykes to enable determination of the depth of this transition, is challenged by the velocity (here 2200 m s⁻¹) that is used to image (stack) it into an MCS profile (e.g. Fig 3e). In Fig. 3(e), the 2A Event arrives ∼100 ms later than the heterogeneous model predicts, placing the Layer 2A/2B boundary ∼300 m below its true position. This result is consistent with those from full-waveform inversion (e.g. Christeson et al. 2012) and tomographic inversion (e.g. Christeson et al. 2010) which propose that the 2A Event is located at or below the bottom of the high velocity gradient zone.

In lithological terms, this would place the source of the 2A Event just above the abrupt porosity change observed in 504B (Carlson 2011). However, although we conclude that the 2A Event is created by the geological heterogeneity of the crust, correlating its location on an MCS profile or velocity model with this geological ‘boundary’ is problematic. Our modelling suggests that the 2A Event mapped in MCS profiles should be treated as a maximum depth for the transition from lower to higher velocity caused by the change in lithology and/or degree of alteration. There is a caveat to this conclusion; the model presented here does not include the effects of local dip and faulting of the layers which will complicate the issue further.

The heterogeneous model shown here (Fig. 13) is one of several tens of similar models that all show equivalent behaviour of the wavefield, with a guided horizontally travelling wave in the horizontally layered Layer 2A, and with the 2A Event generated at the base of this layer. They also show that its amplitude is sensitive to the velocity structure and is most prominent if there are lower velocity channels in the horizontal layering in the lower part of Layer 2A, over the depth range where there is an increasing number of vertically intruding dykes.

5.6 Velocity model for MCS stacking

Stacking MCS data in the oceanic crustal setting is challenging, especially so for younger crust where there is effectively no sediment cover which results in significant seabed scattering. The igneous crust normally lacks first-order (reflection generating) interfaces to enable velocity analysis to underpin normal moveout correction, and the velocity of any sediment cover (even though it may be internally layered and thus reflection generating) is only a few 100 m s⁻¹ greater than that of seawater. Consequently, at abyssal ocean depths, stacking and migrating at 1500 m s⁻¹ generally does as good a job of imaging the seabed and the sediment-basement interface and minimizing diffraction hyperbola (Figs 3a–c), as the more complicated approaches, such as using a stacking velocity that varies with two-way traveltime or applying dip moveout for example. To test this observation, we reprocessed the MCS data from profiles SG_I and SAP_B, for both GI and Bolt source arrays (Figs 18a–h), using a 1-D velocity–depth profile extracted from both the smooth and heterogeneous models (Fig. 18i). Comparison of Figs 3 and 18 shows a marginal gain in the clarity of the sedimentary layering. The most notable improvement is a reduction in the ringing following the sediment-basement reflection and a marked improvement in the resolution of the lateral variation in the geometry of this interface. The same outcome is achieved regardless of smooth or heterogeneous model source of the 1-D stacking velocity function.

Figure 18.

(a)–(h) Reprocessing of the MCS record sections along profiles SG_I and SAP_B (Fig. 2) using velocity functions (i) for normal moveout correction derived from the smooth and heterogeneous synthetic models. The outcome of standard processing using a constant velocity of 1500 m s⁻¹ is shown in Fig. 3. Note that a marginal gain in the clarity of the sedimentary layering, a notable improvement in the ringing following the sediment-basement reflection, and a marked improvement in the resolution of the lateral variation in the geometry of this interface is achieved regardless of a smooth or heterogeneous model source of the 1-D stacking velocity function.

6. DISCUSSION

This study focuses on the 3-D velocity structure and thickness of the ∼7-Myr-old crust of the South Grid, aiming to investigate how well the seismic structure determined by 3-D tomography compares with the lithological ground-truth provided by 504B and, thus, provide a better understanding with which to interpret seismic studies of the oceanic crust. We further use our 3-D P-wave model to determine the origin of the 2A Event and demonstrate how multichannel seismic surveys can be used to constrain upper crustal structure.

6.1 Structure of the oceanic crust—seismic versus lithological

Our tomographically determined P-wave velocity model has an ∼0.3-km-thick sediment layer of between ∼1.6 and 1.9 km s⁻¹ (gradient 1.0 s⁻¹), bound at its base by a velocity step to 4.8 km s⁻¹ at the top of Layer 2. Layer 2 itself is subdivided into two main units by a vertical velocity gradient change at 4.5 km depth bss, with a gradient of 1.7 s⁻¹ above (4.8–5.8 km s⁻¹) and 0.7 s⁻¹ below (5.8–6.5 km s⁻¹). The base of Layer 2, in turn, is defined by a change in gradient at 5.6 km depth. Below this, Layer 3 has a velocity range of 6.5–7.5 km s⁻¹ and a gradient of ∼0.3 s⁻¹. The S-wave model has corresponding layer velocities and gradients for Layer 2A of 2.4–3.1 km s⁻¹ and 1.0 s⁻¹, Layer 2B of 3.1–3.7 km s⁻¹ and 0.5 s⁻¹ and Layer 3 of 3.7–4.0 km s⁻¹ and 0.1 s⁻¹. Together, the 3-D tomographic models, when coupled with gravity modelling, indicate that the crust is ∼6 km thick throughout the region, with a generally flat-lying Moho. Thus, overall, Layer 2 is found to be 1.8 km thick and Layer 3 to have a thickness of 3.8 km.

Although the P- and S-wave models are smooth when compared to the lithology, and parametrized at a vertical node spacing (forward cell size) of 0.1 km, their velocity–depth structure and vertical gradients are remarkably consistent with the main lithological layering subdivisions logged within 504B (Fig. 10) and these layers show little vertical variation throughout the entire 3-D volume (Fig. 11). Consequently, 504B can be regarded as representative of magmatic intermediate spreading generated crust, at least for a 45 km × 45 km square area around it.

Traditionally, the seismic layering model of the oceanic crust first developed at spreading ridge axes (Raitt 1963) has also been applied to aged crust off-axis that has cooled, undergone hydrothermal alteration and infilling of the permeability formed between lava flows and along faults and fractures. Houtz (1976) and Houtz & Ewing (1976) were the first to further extend the 2A/2B subdivision of Layer 2 to include a Layer 2C off-axis, which Salisbury et al. (1985) supported by the interpretation of borehole sonic logs. To date, this 2A/2B versus 2A/2B/2C disagreement has not been resolved by more recent studies (e.g. Christeson et al. 1992, 2007), even by using more sophisticated modelling and data processing approaches on higher quality and denser seismic data sets.

This poses an interpretation challenge between on- and off-axis studies in terms of lithological identification, how the crust is altered as it spreads away from the ridge axis, and how changes between magmatic and tectonic-dominated spreading manifest in seismic data and their resulting models. OSCAR's South Grid is the highest resolution and most densely sampled seismic survey yet undertaken at any borehole in the oceanic crust. What the modelling of both its WA and MCS data show is that there is no need to invoke a third subdivision in Layer 2, but instead that the Layer 2A/2B boundary can be viewed as a transition zone whose thickness is less than the seismic wavelength at its depth, and within which the density of incursion of dykes into the overlying lavas increases (Figs 13d and 15d). There is also no need to invoke a gradient change that corresponds to the top of this transition zone as the changes in physical properties due to the increase in dyke density are simply averaged over a seismic wavelength until, at its base, effectively 100 per cent dyke concentration occurs. From a seismic modelling perspective, this places the Layer 2A/2B boundary gradient change at the base of the lava-dyke transition zone identified in 504B (e.g. Alt et al. 1996). In our models of the oceanic crust of the South Grid that surrounds 504B, Layer 2A is therefore 0.7 km thick and includes the transition zone, and Layer 2B is 1.1 km thick. Neither the top or the base of the transition zone result in an observable MCS reflection event (Figs 14 and 15) and, consequently, there is no direct seismic reflection marker of its depth.

The boundary between Layer 2 and Layer 3 is commonly associated seismically with a change to a low velocity gradient (∼0.1 s⁻¹; Carlson 2001) or higher absolute velocity (6.6–7.6 km s⁻¹, White et al. 1992; 6.93–7.18 km s⁻¹, Grevemeyer et al. 2018b). Various estimates have been placed on the depth and characteristics of the Layer 2/3 boundary at 504B. Swift et al. (1998) interpreted an intraborehole vertical seismic profile to place the boundary at 4.7 km below seabed (bsb), which correlates it with the top of the sheeted dykes traditionally interpreted as the Layer 2A/2B boundary. Detrick et al. (1994) used a sonobuoy study (Collins et al. 1989) to estimate the Layer 3 velocity at 504B to be 6.5 km s⁻¹ and place the Layer 2/3 boundary at 1.5 km bsb, lying within the dykes drilled in 504B. Salisbury et al. (1996) further argue that the Layer 2/3 boundary occurs within the sampled sheeted dyke section near the base of the borehole at ∼2.1 km bsb, and that it is primarily a metamorphic front within the dykes. Interpretation of both the WA and synthetic wave propagation modelling undertaken in this study shows that the Layer 2/3 transition occurs when, over a seismic wavelength, the horizontal fabric inherent of sill-like formation of the lower crust starts to be encountered (Figs 13f and 14). On this basis, 504B samples the lowermost reaches of Layer 2, and is unlikely to have sampled true Layer 3.

In summary, at least for the magma-dominated accretion occurring ∼7 Ma at the intermediate-spreading Costa Rica Rift (Wilson et al. 2019), the 3-D inversion models can be interpreted to locate the boundaries between layers to an accuracy of the vertical model node spacing. They can also be interpreted to identify the Layer 2A/2B boundary as the base of the transition from 100 per cent lavas to 100 per cent dykes, and the Layer 2/3 boundary as the depth at which sill-like fabric starts to become prevalent. The wave propagation modelling further demonstrates that the seismic response of the oceanic crust is most likely a consequence of its fabric orientation (vertical for Layer 2; horizontal for Layer 3), and that it is the contrast between the two inherent orientations that dictates where the seismic gradient changes occur that are associated with the main layering definitions.

6.2 Cause and origin of the 2A Event

Our modelling of wave propagation through smooth and heterogeneous synthetic crustal models shows that no distinct interface, or transition, is required to generate the 2A Event. Instead, it is caused by averaging of heterogeneous physical properties by the seismic wave as it propagates through Layer 2A and is scattered. Thus, we conclude that the 2A Event originates and propagates exclusively in the lower part of Layer 2A, entirely above the mean depth to the top of the dykes of Layer 2B (Fig. 15). Our modelling, that matched the 2A Event observed at 5–7 km offset in the MCS data (Fig. 14), shows that it occurs behind (arrives later) the Layer 2A and Layer 2B refracted arrivals. For the model shown in Fig. 15, the Layer 2A refraction results from a leaky guided wave trapped in the upper 250 m of the Layer 2A lava sequence, and the Layer 2B refraction is produced by the leading coherent wave front in the weakly scattering dyke sequence below 4500 m, equivalent to the base of the 504B Layer 2A/2B transition zone.

However, the key characteristic of the heterogeneous model is the presence of a region of lower velocity in the deepest part of Layer 2A that coincides with the increasing density of dykes (504B's Layer 2A/2B transition zone). The juxtaposition of the strongly scattering interfaces across the direction of wave propagation and the higher velocity layer below created by the dykes, together with the focusing created by the lower velocity region, creates a delayed high amplitude event in the depth range between 4200 m and 4500 m (i.e. above the mean depth of the top of the dykes), which then leaks upward to be recorded by both the seabed (OBS) and sea surface (MCS) receivers. Regardless, modelling does show that its mapping on MCS profiles is a reasonable proxy (i.e. equivalent to the resolution of the 3-D tomographic models) for the maximum depth of the Layer 2A/2B boundary.

Furthermore, a transition zone between Layer 2B and Layer 3 that occurs between 5500 and 6500 m depth, correctly emulates the higher amplitude of the refracted arrival at around 17 km offset and, similarly, a gradient transition at the Moho generates the high amplitude second arrival observed in the WA data (compare Fig. 14 with Figs 4 and 5), which is most likely to be a folded wave front or caustic.

6.3 Value of smooth versus rough velocity models

Although smooth at the scale of the underlying lithological heterogeneity, the contrast in wave propagation through the smooth and heterogeneous models (cf. Figs 14d and e and cf. Figs 15b and e) not only shows that the characteristics of the observed data are influenced by it, even though averaged, it also shows that the outcome of 3-D tomographic inversion also retains the fingerprint of the inherent heterogeneity, albeit also averaged at the model node spacing (cf. Layer 2A velocity–depth structure in Fig. 17c with Fig. 17f). So, although the exact subwavelength heterogeneity cannot be extracted by such modelling, it can be determined to be present; that is, that the crust being modelled is not one-dimensionally homogeneous within layers. A comparison between Figs 17(g) and (h) shows that, even for the same traveltime picks for the same shots, recorded by the same instruments, data acquisition and subsequent inversion modelling is better done in 3-D in the presence of structural or lithological variation-imposed heterogeneity.

7 CONCLUSIONS

We draw the following conclusions:

1. 504B can be regarded as a representative model of magma-dominated, intermediate-spreading generated crust.

2. The Costa Rica Rift had an orientation equivalent to the present-day adjacent Ecuador Rift, when the South Grid crust was formed ∼7 Ma.

3. The velocity–depth structure resulting from 3-D tomographic inversion modelling shows that the crust conforms to the traditional 2A/2B/3 nomenclature derived from ridge-axis studies, despite having undergone cooling and hydrothermal permeability-infilling as it has spread off-axis.

4. The Layer 2A refracted arrival results from a leaky guided wave trapped in the upper part of the Layer 2A.

5. The Layer 2B refracted arrival is produced entirely below the weakly scattering Layer 2A/Layer 2B transition zone.

6. The 2A Event originates at the base of Layer 2A due to the increasing density of dyke incursion (504B's Layer 2A/Layer 2B transition zone).

7. If used as a proxy for the Layer 2A/Layer 2B boundary, the velocity used for optimal stacking of the 2A Event results in its depth being over predicted by several 100 m.

8. Gradient transition at the Moho generates a high amplitude fold in the wave front or caustic.

9. Although smooth at the scale of the underlying lithological heterogeneity, 3-D tomographic inversion models do retain information about that inherent heterogeneity, albeit averaged at the node spacing.

10. Data acquisition and subsequent inversion modelling is better done in 3-D in the presence of structural or lithological variation-imposed heterogeneity.

ACKNOWLEDGEMENTS

This research project was funded by the Natural Environmental Research Council (NERC) grants NE/I027010/1, NE/J022551/1 and NE/J021741/1. We would like to thank all those involved in the planning and acquisition of data during research cruises JC113, JC114 and SO238, including the officers, engineers and crew of the RRS James Cook and the FS Sonne, the scientific party, and all seagoing technicians and engineers. The NERC Ocean-Bottom Instrumentation Facility (Minshull et al. 2005) also provided technical support at sea. Seismic data were processed and manipulated for plotting using Seismic Unix (Cohen & Stockwell 2010). All figures were prepared using the Generic Mapping Tools (GMT) package (Wessel et al. 2013). The final accepted version of this paper is available through Durham Research Online (dro.dur.ac.uk). For the purposes of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence to any author accepted manuscript arising from this work. Finally, we thank Tim Minshull and an anonymous reviewer for their helpful and supportive comments on this paper.

DATA AVAILABILITY

Data from cruises JC113, JC114 and SO238 are archived at the NERC's British Oceanographic Data Centre and are available from the following links:

JC113–www.bodc.ac.uk/resources/inventories/cruise_inventory/report/15035/,

JC114–www.bodc.ac.uk/resources/inventories/cruise_inventory/report/15036/, and

SO238–www.bodc.ac.uk/resources/inventories/cruise_inventory/report/15045/.