Does intermediate spreading-rate oceanic crust result from episodic transition between magmatic and magma-dominated, faulting-enhanced spreading?—The Costa Rica Rift example

SUMMARY

Ocean-bottom seismograph and multichannel streamer wide-angle seismic data are jointly analysed and compared with reflection images, bathymetric maps and potential field data, to reveal the detailed structure of layer 2 of the oceanic crust formed at the intermediate spreading Costa Rica Rift (CRR). Separate modelling of each wide-angle data set independently reveals a gradual increase in P-wave velocity with distance (hence crustal age) from the ridge axis, with a model derived from their joint inversion, in turn, displaying a pattern of shorter-wavelength structural complexity in addition to a background flow-line trend. Normalizing against a ridge-located reference velocity–depth model reveals that, off-axis, velocity perturbations are correlated with trends in basement roughness and uplift; regions of rougher and uplifted basement correlate with slower layer 2 velocity, <0.5 km s⁻¹ faster than at the ridge axis, and thinner sediment cover, while smoother basement and locations where sediment cover forms a continuous seal over the oceanic basement, are mirrored by regions of relatively higher velocity, 1.0–1.4 km s⁻¹ faster than at the CRR. These velocity variations are interpreted to reflect periodic changes in the degree of magma supply to the ridge axis.

Using a combination of global and shipboard magnetic data, we derive a spreading history model for the CRR which shows that, for the past 5 Ma, spreading has been asymmetric. Comparing the seismic model structure with variations in full spreading rate over this period, reveals a correlation between periods of slower spreading and slower layer 2 velocity, basement roughness and uplift, and faster spreading, higher velocity and smoother basement structure. Zones of slower velocity also correlate with lows in the residual mantle Bouguer anomaly, interpreted as most likely reflecting corresponding regions of lower density in the lower crust or upper lithospheric mantle. Using ODP borehole 504B as ground-truth, we show that periods of faster spreading are associated with phases of magmatic accretion, interspersed by phases of increased asymmetric tectonic extension that likely facilitates fluid flow to the deeper crust and results in metamorphic alteration, manifest as the modelled density anomalies.

Overall, our study shows that the mode of CRR crustal formation is sensitive to relatively small changes in full spreading rate within the range of 50–72 mm yr⁻¹, that tips the balance between magmatic and magma-dominated crustal formation and/or tectonic stretching, as characterized by significant variation in the fabric and physical properties of layer 2. We further hypothesize that this inherited structure has a direct influence on the subsequent evolution of the crust through secondary alteration. We conclude that descriptive phrases like ‘ocean crust formed at an intermediate-spreading rate’ should no longer be used to describe an actual crustal formation process or resulting crustal structure as, over the full range of intermediate spreading rates, a fine tipping-point dictates an episodic transition between primarily magmatic accretion and magma-dominated crustal formation coupled with enhanced faulting, with their asymmetry recorded in either ridge flank.

Composition and structure of the oceanic crust, Controlled-source seismology, Mid-ocean ridge processes

1 INTRODUCTION

In the simplest model of oceanic crustal formation at mid-ocean ridges, growth occurs consistently over timescales of several million years, even though at most ridges there is evidence for temporal fluctuations in spreading rate (Müller et al.2008). For ridges spreading at the slowest rates, these are interpreted as having the poorest magma supply and, consequently, undergo significant tectonic extension, expressed in terms of a ratio, M, of magmatic crustal accretion to amagmatic or tectonic extension, where M = 0.5 acts as the discriminator between these types of spreading (e.g. Buck et al.2005). Thus, the structure and morphology of the oceanic crust along a flow-line away from a ridge axis can be used as a proxy for determining changes in the mode of formation (e.g. Kappel & Ryan 1986; Phipps Morgan & Chen 1993). If the ridge axis fluctuates between magmatic and reduced magma supply conditions on timescales of a few million years, as may be the case at intermediate spreading rates, studying this ridge type provides an opportunity to understand the tipping-points between the different modes of crustal formation, as well as the effects of variable spreading rate on processes that age the crust as it spreads off-axis.

Based on the M-ratio classifications of Buck et al. (2005), we first define a set of terms that we use to describe the different spreading modes that may occur at spreading ridges, together with the transitions between them. These comprise:

• magmatic, where M ≈ 1, which we define as meaning that almost all of the plate separation is accommodated by eruption of lavas and intrusion of dykes at the ridge axis;

• magma-dominated, where 0.5 < M < 1, which we define as meaning eruption of lavas and dyke intrusion are still the principal processes accommodating spreading, but do not account for all of the plate separation;

• faulting-enhanced, which lies at the low-M end of the magma-dominated range, and where lower M values correspond to an increasing component of spreading being accommodated by faulting or tectonic extension, and;

• tectonic-dominated, where M < 0.5, and accretion of new magmatic material is no longer the principal process in plate separation, which, instead, primarily occurs through faulting/tectonic extension. This range of M-ratio is restricted to mainly slow and ultra-slow spreading ridges and is, therefore, not applicable to the intermediate spreading Costa Rica Rift studied here, however we include it for completeness.

A series of models have been proposed to account for the different spreading styles commonly observed at slow- and ultraslow mid-ocean ridges, with half-spreading rates (HSR) <25 mm yr⁻¹, or full-spreading rates (FSR) <50 mm yr⁻¹. Under increasingly tectonic-dominated conditions (M < 0.5), crustal growth is predominantly asymmetric (e.g. Cannat 1993; Ranero & Reston 1999; Cannat et al.2006; MacLeod et al.2009). Depending on the degree and focus of any magma supply, one flank of the ridge axis may be accreted with a ‘normal’ layer of volcanic rocks fed by a system of dykes, while the other flank either develops into a series of volcanic fault blocks (Tucholke et al.2008; Olive et al.2010), while, at the lowest values of M, a low-angle detachment fault initiates to produce oceanic core complex-style structures (e.g. Searle et al.2003; Smith et al.2006; Dannowski et al.2010), exhuming sections of lower crustal gabbro and ultimately ultramafic mantle rocks. Both of these possibilities result in a rougher basement surface, and a potential for alteration of peridotite assemblages to serpentinite by fluid-rock interaction.

In contrast, at fast spreading ridges (70 mm yr⁻¹ FSR) crust is accreted under magmatic conditions (M ≈ 1) and a continuous magma supply is frequently observed (e.g. Detrick 1987; Kent et al.1990; Sinton & Detrick 1992; Hooft et al.1997). A relatively smooth top basement surface results, and the subsurface is composed of basaltic lava flows and dykes of oceanic crustal layer 2, accreted relatively symmetrically about the ridge axis (Sinton & Detrick 1992). The seismic velocity–depth structure of layer 2 is well documented at fast spreading ridges (e.g. Christeson et al.1992; Harding et al.1993; Vera & Diebold 1994; Grevemeyer et al.1998), and is commonly subdivided into two seismically distinct sublayers, defined at the ridge axis as layers 2A, which comprises extrusive, high porosity pillow basalts and 2B, comprising sheeted dykes (e.g. Herron 1982; Christeson et al.1992; Harding et al.1993). As the crust ages, there is a trend towards an overall increase in seismic velocity, related to cooling and hydrothermal alteration processes (Houtz & Ewing 1976; Christensen 1979; Carlson 1998), which may alter the appearance and/or manifestation of the layer 2A/2B boundary, largely as a result of the progressive closure or infilling of the bulk porosity (Christensen 1978; Vera et al.1990; Detrick et al.1994; Christeson et al.2007).

Between these two end-members, crust formed at intermediate spreading ridges (50–70 mm yr⁻¹ FSR) might be expected to demonstrate some form of average of the two end-member structures, although the significant heterogeneity of crust formed at the slower spreading rates suggests that this may not be the case.

Studies of intermediate spreading ridge axes and flanks are relatively limited. Wide-angle seismic imaging has been used to characterize the off-axis velocity structure of the intermediate spreading (56 mm yr⁻¹ FSR) Juan de Fuca Ridge (JdFR) flanks, over several million years of spreading. Close to the Blanco Transform, layer 2A has a thickness of 0.49–0.54 km and velocity of 3.0–3.2 km s⁻¹, and is located above a 0.23–0.28-km-thick high-velocity gradient transition zone into layer 2B (Christeson et al.2012). Uppermost layer 2B velocities at the JdFR vary depending on location, with values of 4.7–4.9 km s⁻¹ (Christeson et al.2012) and 5.2 ± 0.3 km s⁻¹ (Newman et al.2011) obtained by different studies. In both cases, however, these values are elevated relative to those of 4.4 ± 0.3 km s⁻¹ at the ridge axis (Newman et al.2011), with the increase occurring within ∼0.5–1.0 Myr of crustal formation. This increase occurs at a rate of up to four times faster than would be expected purely by cooling alone (Christensen 1979), suggesting that layer 2B velocity evolution is controlled by active, rather than passive, hydrothermal circulation. Overall, the observed upper layer 2B velocity in 1.3–6.0 Ma magmatically generated crust at the JdFR (Horning et al.2016) is not significantly different from that in 0.5–8.3 Ma crust formed at the superfast spreading southern East Pacific Rise (EPR; Grevemeyer et al.1998), indicating that the magmatic crustal accretion processes are fundamentally similar in their action and result.

Variable along-ridge axial morphology has been documented at the intermediate spreading: JdFR (Canales et al.2005), which displays an overall bathymetric high with a shallow central axial rift, similar to the fast spreading EPR (Phipps Morgan & Chen 1993); the Southeast Indian Ridge (SEIR; ∼75 mm yr⁻¹ FSR), which shows an overall increase in axial depth from 2.3 to 4.8 km between 88°E and 118°E (Cochran et al.1997; Ma & Cochran 1997); and the Galapagos Spreading Centre (GSC; Detrick et al.2002; Blacic et al.2004), which changes from a developed axial rift valley in the west, in the region around Hess Deep Rift and the Galapagos Triple Junction, to an EPR-type axial high in the east, where it approaches the Galapagos hot spot. In turn, each of these ridge systems also show evidence for shallow axial melt lenses along their axes (Detrick et al. 1987; Blacic et al.2004; Baran et al.2005; Canales et al.2006; van Ark et al.2007). Inter-relationships between axial morphology, evidence for active magma supply, and the thickness of layer 2A suggest that these features are controlled by distinct modes of crustal accretion (Phipps Morgan & Chen 1993). EPR-like axial highs tend to be associated with the presence of shallow melt lenses and a thinner on-axis layer 2A, whereas rifted axial highs correlate with deeper melt lenses and thicker on-axis layer 2A. Segments with well-developed axial valleys do not exhibit melt lenses but have a thicker layer 2A (e.g. Small 1994; Baran et al.2005; Canales et al.2005; van Ark et al.2007).

Whilst layer 2A seismic velocity increases away from the ridge (Houtz & Ewing 1976; Christensen 1979; Carlson 1998), systematic studies of correlations between this and layer 2A thickness, layer 2A/2B transition thickness, or layer 2B seismic velocity with distance from the ridge axis do not exist. Consequently, our current understanding of how the oceanic crust evolves over time (e.g. how and when the porosity closes, and how layer 2 structure changes as a result) assumes that the growth style at the ridge axis does not vary through time. In this paper we test this assumption, by examining the shallow structure of young oceanic crust that has been formed at the Costa Rica Rift, which has a variable intermediate spreading rate (Wilson & Hey 1995). We use the combined proxy of basement roughness and seismic velocity characteristics of the crust to determine magma supply variation at the ridge axis at the time of formation, and use this proxy to investigate how the characteristics generated at formation may affect subsequent crustal evolution, or ageing.

2 TECTONIC SETTING

The Panama Basin is a small ocean basin located in the Eastern Pacific between 0–7°N and 78–90°W, proximal to the subduction zones of Central and South America (Fig. 1a). The oceanic crust in this basin is generated along the eastern limb of the Galapagos Spreading Centre (GSC), which separates the Cocos and Nazca Plates. The Cocos and Carnegie Ridges bound the basin to the west and the south, respectively. These are shallow aseismic ridges of thickened crust that formed as the plates moved over the Galapagos hotspot (e.g. Sallarès et al.2003, 2005). The GSC is divided into a number of ridge segments, offset by transform faults. The easternmost spreading ridge segment is known as the Costa Rica Rift (CRR), which is bounded to the west by the Ecuador Fracture Zone (EFZ) and to the east by the Panama Fracture Zone (PFZ).

Figure 1.

Tectonic setting of the Costa Rica Rift and ODP borehole 504B in the Panama Basin. (a) GEBCO bathymetry (GEBCO 2008) of the Panama Basin, with solid red lines showing west–east trending ridge segments that mark the boundary between the Cocos and Nazca Plates, offset by north–south trending fracture zones (dark red dashed lines). Motion of the Cocos and Nazca Plates (NRR-MORVEL56; Argus et al.2011) is shown by labelled arrows, with rates in mm yr⁻¹. GSC, Galapagos Spreading Centre; ER, Ecuador Ridge; CRR, Costa Rica Rift; EFZ, Ecuador Fracture Zone; PFZ, Panama Fracture Zone. Black dashed box shows the location of the OSCAR seismic experiment shown in part (b). ODP borehole 504B is indicated with a green dot, SAP profiles by black lines, and SAP_B by the orange line. (b) Swath bathymetry-derived seabed topography illuminated from the north, highlighting the axis-parallel fault planes. Inverted black triangles mark OBS locations, with instrument numbers labelled. OBSs with example data record sections shown in Fig. 3 are coloured red. SAP profile locations are labelled. ODP borehole 504B is indicated with a green circle. Labelled arrows at the ridge axis indicate present-day half-spreading rates (mm yr⁻¹).

The CRR is ∼140 km long and has been spreading for 11 Myr (Lonsdale & Klitgord 1978). Based on shipboard magnetic anomaly data, Wilson & Hey (1995) characterized the ridge as an asymmetric intermediate spreading ridge, with current HSRs of 34 mm yr⁻¹ to the south and 28 mm yr⁻¹ to the north (Fig. 1b), and also showed that the spreading rate has fluctuated significantly over the past 8 Myr, encompassing the full range of HSRs defined for intermediate spreading ridges, of 25–35 mm yr⁻¹ (Dick et al.2003). There is evidence for an active magma supply to the ridge axis currently, in the form of an ∼5-km-long axial melt lens reflection (Buck et al.1997; Lowell et al.2016; Zhang et al.2017), although its small size indicates that it would require regular replenishment to prevent ‘freezing’ (Liu & Lowell 2009; Lowell et al.in review).

The oceanic crust in the Panama Basin has been sampled and logged at Ocean Drilling Program (ODP) borehole 504B, located ∼233 km south of the CRR (Becker et al.1989, and references therein), where the crust has been dated at 6.9 Ma (Wilson et al.2003). The 275 m sediment succession (Cann et al.1983) implies an average sedimentation rate of ∼40 m Myr⁻¹, which is double that for the Pacific Ocean over the same time period (Davies et al.1977). The lithological structure and hydrothermal alteration of the igneous crust at 504B have been characterized in detail throughout an ∼1830 m section, sampling the volcanic lavas (2A) and sheeted dykes (2B) (Alt et al.1996). This geological reference point can be used to ground-truth geophysical interpretations, and a number of previous seismic surveys detect an ∼600-m-thick layer 2A with P-wave velocity varying between ∼4.5 and 5.2 km s⁻¹, and an ∼1200-m-thick layer 2B with P-wave velocity between ∼5.8 and 6.5 km s⁻¹ (Detrick et al.1998) with, more recently, equivalent values obtained from coupled 2-D forward modelling and 3D inversion by Gregory et al. (2017).

In this study, we use newly acquired geophysical data along a flow-line from the CRR to ODP borehole 504B to investigate changes in the morphology and velocity structure of layer 2 in young oceanic crust. From these observations, we infer the formation history and subsequent ageing of the crust formed at the CRR over the past 8 Myr.

3 DATA ACQUISITION

The geophysical component of the OSCAR (Oceanographic and Seismic Characterization of heat dissipation and alteration by hydrothermal fluids at an Axial Ridge) project was carried out during a pair of research expeditions in January–March 2015 (Hobbs & Peirce 2015; Morales Maqueda 2015), onboard the RRS James Cook (JC114) and FS Sonne (SO238). As part of the sea-going experiment, three types of controlled-source seismic acquisition were undertaken contemporaneously along three ∼330-km-long flow-line profiles, running from 40 km north of the CRR to 40 km south of ODP borehole 504B (Fig. 1b). In this study, we focus on the central of these profiles, SAP_B, which passes through the previously imaged axial magma lens at the CRR (Zhang et al.2017), and borehole 504B (Fig. 1b), where lithological samples and downhole geophysical logs (Alt et al.1993; Carlson 2010) provide ground-truth for the seismic images (Gregory et al.2018). Swath bathymetry, gravity and magnetic data were acquired throughout seismic acquisition.

To address the goals of this study, optimizing resolution, imaging depth and lateral propagation for wide-angle velocity modelling and inversion, three different seismic source arrays were utilized, with an interleaved firing pattern over a 60 s repeating window. Two arrays were towed astern of RRS James Cook: a pair of GI airguns (generator and injector volume each 105 in³, 1.72 l), and an array of six Bolt airguns (1320 in³, 21.63 l); while the FS Sonne towed an array of eight Sercel G-guns (4280 in³, 70.14 l). The three seismic data sets acquired were as follows:

• 2-D multichannel seismic reflection data, using the GI airgun array and a 4.5-km-active-length streamer with 360 channels spaced at 12.5 m intervals, towed at a depth of 10 m. Data were recorded at a sampling rate of 500 Hz, over trace record lengths of 41 s. This data set was designed to image the depth to the top igneous basement and the structure and layering of the overlying sediment cover, and will be referred to as MCS data;

• 2-D synthetic aperture multichannel seismic reflection data, using the same streamer but with airgun shots fired by both the Cook (Bolt airgun array) and the Sonne (G airgun array) following ∼8.5 km astern. This data set was designed to enhance lateral and vertical resolution imaging of layer 2 by recording both wide-angle (WA) reflections and emergent layer 2 refracted arrivals, which will be referred to as SAP MCS and SAP WA data, respectively;

• wide-angle refraction data, using 35 ocean-bottom seismographs (OBSs) deployed at 5–10 km intervals. These four-channel (three-component geophone and hydrophone) OBSs recorded all three source arrays continuously at sampling rates of either 250 Hz or 500 Hz, depending on instrument specification. This data set was designed to provide the background upper crustal velocity–depth structure for appraisal of the initial model used for detailed layer 2 modelling, and will be referred to as OBS WA data.

4 DATA PROCESSING AND MODELLING

4.1 MCS data

The aims of MCS imaging were to determine the depth to the top of the igneous basement, provide an estimate of sediment thickness, and determine basement roughness along profile. Streamer feather for the central flow-line profile was <3° throughout acquisition, equivalent to a maximum lateral offset, from in-line, of 250 m at the farthest-offset channel. Given the size of the Fresnel zone at the seabed, and that the dominant out-of-plane structures are orientated perpendicular to the profile (ridge-parallel), this in-line offset was deemed within lateral positioning error and the MCS reflection data were assigned a 2-D geometry assuming the source and receiver towed directly astern. Following this, a simple processing sequence was applied, comprising constant velocity (1500 m s⁻¹) normal move-out correction, stacking and Stolt migration, again at 1500 m s⁻¹. Fig. 2 shows the resulting stacked GI-source MCS data, in which the primary reflections associated with the seabed, top basement surface and intrasedimentary reflections, together with evidence for pervasive basement and sediment column faulting, are observed.

Figure 2.

Stacked MCS reflection data along Profile SAP_B. The mean thickness of sedimentary cover over the high amplitude basement reflection increases with increasing crustal age, although this cover is punctuated by the footwalls of normal faults in areas of rough bathymetry between 90 and 190 km. Vertical exaggeration is ∼12 at the seabed. Inverted black triangles mark OBS locations, with red triangles for those OBSs shown in Fig. 3 (cf. Fig. 1b). Ages assume a constant half-spreading rate of 33.8 mm yr⁻¹ from the CRR to ODP borehole 504B (6.9 Ma). Annotations highlight key features that are discussed in the text.

4.2 WA data

The aim of WA data modelling was to produce a P-wave velocity–depth model that could be used to describe changes in layer 2 structure along Profile SAP_B, with the results of modelling geologically ground-truthed against the direct lithological sampling and borehole geophysical logs at ODP borehole 504B to inform interpretation. Modelling was conducted using a two-staged approach. Initially, a first-pass forward model was constructed using rayinvr (Zelt & Smith 1992), to determine the longer-wavelength P-wave velocity structure used to inform construction and testing of a range of initial models prior to the selection of the starting point for inversion. The chosen inversion starting model used a generic one dimensional velocity–depth structure that contained the minimum pre-conceived characteristics of oceanic crustal structure, based on the MCS subseabed image (Fig. 2). Subsequent inverse traveltime modelling, conducted using FAST (Zelt & Barton 1998), was then used to produce the final velocity–depth model by iterative updates until the best-fitting solution within pick errors was reached.

4.2.1 Traveltime picking and initial model preparation

Wide-angle arrivals were picked preferentially on the OBS hydrophone record section, as this had the best signal-to-noise ratio (SNR). Each OBS was relocated to its actual position on the seabed by matching the observed traveltimes of the direct water wave arrival (P_w), with traveltimes calculated by forward ray-tracing using rayinvr (Zelt & Smith 1992). The along-profile seabed depth was constructed using ship-acquired swath bathymetry data, and a water column model was constructed using the average sound velocity profile calculated from temperature and salinity measurements made during eight conductivity–temperature–depth (CTD) casts (Banyte et al.2018). This water column model was used in, and remained fixed for, all subsequent subseabed modelling.

Example OBS record sections from instruments located near the ridge axis (OBS 04) and borehole 504B (OBS 32) are shown in Fig. 3. Refracted arrivals turning in the upper crust, P_g1, and mid-to-lower crust, P_g2, were picked and assigned an uncertainty of ±30 ms to represent the combined location and pick errors. Around 5300 crustal (P_g1 and P_g2) traveltime picks were made, with shot-receiver offsets ranging from 1.7 to 31.6 km.

Figure 3.

Example OBS hydrophone record sections from OBSs 04, located at the CRR axis, and 32, located at borehole 504B. (a) and (c) Record sections displayed after minimum phase bandpass filtering (2–4–24–36 Hz). (b) and (d) Indicative traveltime picks for refracted phases used in modelling are annotated, with symbol width corresponding to the assigned pick uncertainty and colour to phase identification. Record sections are plotted with a reduction velocity of 6 km s⁻¹. See Fig. 1 for location of instrument relative to profile bathymetry. P_w—water waves; P_g1–upper crustal refracted arrivals; and P_g2–middle-to-lower crustal refracted arrivals.

For the SAP WA data, as the water depth is greater than 2400 m, most of the refracted energy for the shallowest subsurface arrives at shot-receiver offsets <4.2 km, behind the primary seabed reflection. The refracted arrival turning in the upper-to-mid crust (P_g) becomes an emergent phase, arriving before the seabed reflection, between 4.2 and 5.5 km shot-receiver offset (Fig. 4). More than 500 000 traveltime picks were made of these emergent P_g phases, with offsets ranging from 4.2 to 8.8 km. Each pick was assigned an uncertainty of ±20 ms. The SAP WA data set significantly increases the sampling density of the upper 1–2 km of the crust; for comparison only ∼2000 picks from the OBS data set fall within a comparable shot-receiver offset range of <8.5 km.

Figure 4.

Example SAP WA data recorded by the MCS streamer, displayed with a minimum phase bandpass filter (2–4–24–36 Hz). Phase identification are annotated: P_sbP, seabed reflection; P_g, crustal arrivals. Record sections are plotted with a reduction velocity of 6 km s⁻¹. Indicative traveltime picks are annotated and have a ± 20 ms assigned uncertainty.

4.2.2 Forward modelling

First-pass forward modelling of the OBS WA traveltimes using rayinvr (Zelt & Smith 1992) aimed to determine the longer-wavelength P-wave velocity structure to act as a starting point for and inform subsequent inversion. A simple starting model was constructed with the water column as described above and the sediment thickness calculated from the MCS data using a mean interval velocity of 1630 m s⁻¹ (determined from stacking velocity analysis of the MCS data, and wireline logs run in borehole 504B). Beneath the top basement surface, the igneous crust was constructed as a series of layers, interpolating between two 1-D velocity–depth profiles (Fig. 5a), based on the Main Endeavour segment of the Juan de Fuca Ridge (van Ark et al.2007) for the ridge axis and a 3-D tomographic inversion around borehole 504B (Detrick et al.1998).

Figure 5.

Model initialization. (a) Representative 1-D velocity–depth profiles used to construct initial velocity models. Black lines show the layered model, created for forward modelling by interpolating between the velocity–depth structures of the Juan de Fuca Ridge (van Ark et al.2007) for the ridge axis, and borehole 504B (Detrick et al.1998). Red lines show the smooth initial inversion model. (b) Ray-trace forward model for Profile SAP_B. Model has a rms misfit of ∼56 ms, corresponding to a χ² fit of ∼3.4, fitting ∼5300 P_g1 and P_g2 traveltime picks. Inverted black triangles indicate OBS locations. Solid black line marks the model layer 2A/2B boundary. Dotted line shows the 4.3 km s⁻¹ contour. Note how this contour shallows with distance from the ridge axis until it reaches the depth of the top basement surface.

Iterative forward modelling was conducted using a top-down approach. Lateral velocity variations were introduced, over a scale of 50–100 km, to improve the fit between observed and calculated traveltimes, until a χ² fit of ∼3.4 was achieved for the ∼5300 traveltime picks. For the 30 ms traveltime pick error, this corresponds to a traveltime root-mean-square (rms) residual of ∼56 ms. The resulting model (Fig. 5b), showing the longer-wavelength, background P-wave velocity–depth structure will, henceforth, be referred to as the forward model.

The 4.3 km s⁻¹ contour, representing the base of layer 2A at the ridge axis, can be used as a proxy for the progressive change in velocity with crustal ageing due to infilling of porosity (e.g. Carlson 1998). In the forward model, this contour shallows from the base of layer 2A at the ridge axis to the top of the basement at ∼170 km distance along profile (Fig. 5b). The forward model shows only limited variability within the upper part of layer 2B, although this might reflect the smaller number of velocity nodes included at greater depth, and that forward modelling only continued until the rms misfit decreased to 56 ms. However, there is a resolvable increase in velocity at the top of this layer, from 5.6 km s⁻¹ at the ridge axis to 6.0 km s⁻¹ at 504B.

Resolution of the shallowest crust using OBS-acquired WA data sets alone is limited. However, the forward model does constrain the general along flow-line velocity–depth characteristics, and represents a first pass solution to a limited subset of the total traveltime pick set, allowing us to obtain an indicative view of crustal structure variation along profile. Consequently, we use the forward model to provide constraint on the starting point for inversion of the MCS WA traveltime picks. Since the forward model has been developed from a no preconceived view stand point, we regard such an inversion starting point as inheriting this origin.

4.2.3 Initial inversion modelling

A 2-D tomographic inversion was conducted using the widely adopted FAST, following the approach of Zelt & Barton (1998). The inversion algorithm iteratively updates a smooth velocity field, aiming to minimize traveltime residuals along ray paths, to reach a χ² fit of 1, which represents a fit to the error bounds. The inverse model was parametrized on a uniform square grid, with the seabed and top basement surface depths taken from the forward model and fixed as invariant interfaces throughout inversion. To represent the igneous crust, a 1-D velocity gradient was added below the basement surface along profile, with velocity increasing from 4.0 km s⁻¹ at the basement to 8.2 km s⁻¹ at the base of the model at 12 km depth. Fig. 5a) shows a comparison of 1-D velocity–depth profiles from the initial models used as part of the forward and inversion modelling, at various locations along the profile.

As its name implies, FAST (First Arrival Seismic Tomography) is an approach that inverts first arrival traveltime picks. However, for the SAP WA data, the first arrival by definition is the water wave travelling directly between shot and receiver through the water column. Due to geometry complications caused by the non-constant distance between the two ships, and aliasing caused by the shot spacing, it was not possible to apply downward continuation to effectively relocate the streamer to the seabed. Therefore, to invert the SAP WA traveltime picks, it was instead necessary to adapt the FAST inversion procedure. For each ray path calculated, the water layer was temporarily modified to introduce an artificially slow region (0.001 km s⁻¹) at the midpoint between shot and receiver to, in effect, delay the direct water wave to arrive later than any subseabed phase. This slow region was applied using a five-model node-wide column, or band, that extended through the entire water layer of the model.

The OBS WA and SAP WA traveltime data sets were first independently inverted over a series of five non-linear iterations, with ray paths computed and tested over five potential models to obtain the smoothest velocity model that reduced the misfit. Using the forward model as a reference for comparison, the observed similarity in structure with the inversion results for both data subsets (OBS WA and SAP WA), accounting for inherent differences between the two modelling approaches, allows us to be confident that the inversion models are a true function of the data and are not biased by, for example, the parametrization of the inversion or the starting model.

The resulting best-fitting inversion models for the OBS WA and SAP WA data subsets are shown in Figs 6a) & c), and, henceforth, are referred to as the OBS inversion and SAP inversion models, respectively. Comparing these two models, a clear difference is observed in the shallowest 600 m subseabed, where the SAP inversion model shows that the P-wave velocity is lower than that recovered in the OBS inversion model, by up to 1 km s⁻¹ for young off-axis crust (<70 km south of the ridge axis–Fig. 6e). This shallow difference tends to reduce with increasing distance from the ridge axis. Any remaining difference between the models is concentrated at deeper (>600 m) depths and higher (≥5 km s⁻¹) velocities. However, this feature is also an artefact of the shallower imaging depth of the SAP WA data compared to that of the OBS WA data.

Figure 6.

P-wave velocity models resulting from independent inversions of OBS WA and SAP WA traveltime picks and corresponding smallest resolvable checkerboard patterns. (a) Velocity model resulting from inversion of the ∼5300 P_g traveltimes picked from the OBS WA data, masked to show extent of ray coverage. Black dashed line shows 600 m depth region below seabed (see text), smoothed over a 10 km window. (b) Recovery of a ±5% velocity anomaly checkerboard (10 km x 1.5 km) using the same shot-receiver geometry and picks as the OBS WA inversion. (c) Velocity model resulting from inversion of the ∼500 000 P_g traveltimes picked from the SAP WA data, plotted as for (a). (d) Recovery of a ±5% velocity anomaly checkerboard (5 km x 1 km) using the same shot-receiver geometry and picks as the SAP WA inversion. (e) Sampled 1-D velocity–depth profiles from the OBS (black) and SAP (red) models at locations between the ridge axis and borehole 504B. Note that the SAP WA data generally results in lower velocities in the upper 600 m of the model.

4.2.4 Resolution

The minimum resolvable velocity anomalies in each of the OBS and SAP inversion models were determined using checkerboard testing (Zelt & Barton 1998). A ±5% sinusoidal velocity perturbation pattern was added to each inversion model and then synthetic traveltimes were calculated using the shot-receiver offsets of the observed traveltime picks. Random Gaussian noise was then added to the synthetic traveltimes, at a level corresponding to the pick uncertainties. These synthetic picks were then inverted with the same inversion parameters used to generate the inverse models and the degree of recovery determined.

For the OBS WA traveltimes in isolation, the smallest recoverable checkerboard has cell dimensions of 10 × 1.5 km (Fig. 6b). Poorer recovery is observed between 30–40 km and 200–210 km distance along profile, which coincides with an instrument that failed to record. Additionally, in several places sub-seabed, the amplitude of the applied velocity anomaly is over-recovered in the shallowest 300 m, reflecting the limitations of OBS data for shallowest imaging due to the generally larger intra-instrument spatial separation, and the sparser short-offset traveltime picks resulting from the larger shot spacing (∼150 m).

The SAP WA traveltime data, on the other hand, provide an even and dense coverage of the upper 2 km of the crust, as a result of the smaller shot-receiver spacing (∼50 m). The smallest recoverable checkerboard in this case is 5 × 1 km (Fig. 6d). Most notably, recovery is enhanced beneath the ridge axis, and only worsens due to gaps in data coverage (e.g. 215–220 km distance along profile, Fig. 6d). Over-recovered checkerboard amplitudes can still be observed within the shallowest 300 m subseabed, but this is much reduced compared to the solely OBS data-derived resolution.

4.2.5 Final inversion

Having determined the characteristics that each traveltime pick subset imposes on a resulting inversion model, a further, final, inversion was conducted using both the OBS WA and SAP WA data jointly, with the goal of achieving the best constraint of the shallow subsurface (<1 km) velocity structure, whilst also ensuring the deeper part of layer 2 and below are also well resolved. Due to the large mismatch in the number of picks from each data set (∼5300 versus ∼500 000), the SAP WA data were reduced in number by selecting an evenly spaced 25% subset of shots and receivers. At ∼30 000 picks, the subset of SAP WA picks still outweighed the OBS WA data set, so the traveltimes were weighted towards the top of the model (i.e. more rays turning in the upper 1–2 km of the model than deeper down). To ensure that the resolution of the model was not artificially improved beyond that achievable from OBS WA data alone (i.e. to avoid imposing bias), a horizontal 10 km-wide Gaussian filter was applied to the resulting joint inversion model, to smooth the velocity field to the feature-scale resolution of the OBS inversion model. The final joint inversion model is shown in Fig. 7b) and is, henceforth, referred to as the final inversion model.

Figure 7.

Final P-wave velocity–depth models showing the evolution of layer 2 from the ridge axis to borehole 504B. (a) Forward model (cf. Fig. 5b). (b) Preferred final inversion model from combined inversion of OBS WA and SAP WA traveltimes, with a rms misfit of 30 ms. In both parts inverted black triangles mark OBS locations, the 4.3 km s⁻¹ velocity contour is plotted as dashed black line and the forward model layer 2A/2B boundary is shown as a solid grey line for reference. Ages assume a linear half-spreading rate that matches the 6.9 Ma crustal age at borehole 504B. (c) Relative velocity model showing the difference in final inversion model P-wave velocity compared to a reference 1-D velocity–depth profile from the ridge axis. A lateral variation in layer 2 velocity is observed along profile. The forward model layer boundaries are shown as solid grey lines for reference. Dashed black line shows the depth limit of confidence in the model at ∼1.8 km subseabed. (d) Sampled 1-D velocity–depth profiles from the forward (black) and final inverse (red) models at locations between the ridge axis and borehole 504B. For profiles other than at 0 Ma, the blue line is the 1-D velocity–depth profile at ridge axis shifted to match top of layer 2 of the others.

5 RESULTS

5.1 Surface expression of crustal formation

Morphological characterization of the bathymetry of the south flank of the CRR (Fig. 1b) reveals that the majority of topographic features are elongate, axis-parallel ridges, with steeper northward-dipping sides that have been inferred to be normal faults (Haughton et al.2017). The stacked MCS section (Fig. 2) displays clear, discrete reflection events representing the seabed and top basement surface along the majority of the profile. Exceptions to this are close to the ridge axis (±30–40 km), where the top of the oceanic crust is the seabed by definition, and on surfaces dipping at greater than ∼10° (e.g. 94 km, 115 km and 122 km distance along profile). Over time, the irregularly undulating basement surface has become progressively infilled by sediment, giving the appearance of more subdued seabed topography at distances of >80 km from the ridge axis. Sediment cover becomes continuous south of ∼190 km from the ridge axis, however intrasediment reflectors show deformation patterns which suggest on-going active fault movement and/or fluid venting along the whole profile.

To investigate the extent of tectonism, and negate the effect of progressively thickening sediment cover, we quantify crustal roughness using the top basement surface as imaged by the MCS data, calculating the rms residual of its topography relative to a predicted depth dependent on the wavelength of roughness being investigated (Figs 8b and c). To determine the longer-wavelength roughness, we apply this calculation using a 100-km-wide sliding window and obtain a range of rms roughness of 110–180 m (Fig. 8c) relative to the predicted subsidence associated with plate cooling with distance away from the ridge (Hasterok 2013), described in the following section. This range of roughness agrees well with that for crust formed at similar intermediate spreading rates (e.g. Goff 1991; Ehlers & Jokat 2009).

Figure 8.

Crustal roughness and subsidence analysis. (a) Bathymetry along Profile SAP_B, extended southward towards the Carnegie Ridge. Shipboard bathymetry shown as a black line and GEBCO (2008) bathymetry as a grey line. (b) Top basement surface depth (black line) calculated using sediment two-way traveltime ‘thickness’ and an interval velocity of 1630 m s⁻¹, filtered using a 25-km-wide Gaussian window (green line). Predicted depth (red line) from the Hasterok (2013) subsidence curve, using d0 = 2.7 km and a half-spreading rate of 33.8 mm yr⁻¹ shown for comparison. (c) Long-wavelength rms roughness of the top basement (red line), calculated using a 100-km-wide sliding window, shown as a residual relative to the predicted depth in (b). Short-wavelength rms roughness (green line) calculated using a 25-km-wide sliding window as a residual relative to the filtered top-basement depth. (d) Expected crustal subsidence (Hasterok 2013) for a half-spreading rate of 33.8 mm yr⁻¹ and initial ridge axis depths (d0) of 2.7 km (red line), 3.0 km (purple line) and 3.3 km (blue line). Regions of rough bathymetry lie shallower than the expected subsidence. (e) Half-spreading rate for the south flank of the Costa Rica Rift, calculated from the full rates and degree of asymmetry of Wilson & Hey (1995).

We also calculate the basement roughness attributable to shorter wavelength features (<25 km; henceforth the short-wavelength basement roughness), following the method of Hayes & Kane (1991). The basement depth, measured by converting seabed-to-basement MCS traveltimes using a sediment average velocity of 1630 m s⁻¹, was filtered with a 25 km-wide Gaussian filter, and the raw (unfiltered) depths were subtracted to give the short-wavelength rms roughness at 5 km intervals along profile. To avoid influence in the filter window of measurements from the opposite ridge flank, and those associated with the youngest crust that may still be tectonically evolving, the results of this analysis are only considered valid at distances >12.5 km (∼0.4 Myr) from the ridge axis. The short-wavelength basement roughness ranges from 35 m to 136 m, and alternates between rough (which we define here as >86 m) and smooth (<86 m) on length scales of 50–100 km (Fig. 8c). This measure is considered a good representation of the qualitative observations based on the MCS stacked section (Fig. 2).

Rough basement is observed along two sections of the profile: (i) near the ridge axis (<40 km); and (ii) between 90 and 190 km from the ridge axis. These sections are dominated by a series of irregularly spaced normal faults, dipping towards the ridge axis at ∼30°, that have throws of >200 m (Fig. 2). Sediment accumulation is controlled by the basement topography, predominantly occurring on the hanging walls between faults and increasing in thickness towards the south (footwall). Where significant cumulative fault throw is greater than the thickness of accumulated sediment, the basement remains exposed either on the fault scarp itself (e.g. 114 and 170 km from the ridge axis), or directly at the seabed where the slope prevents sediment accumulating (e.g. 148 and 189 km). Smooth basement is observed between 40 and 90 km, and at a distance of >190 km from the ridge axis. The subdued top basement topography has an average fluctuation in depth of ∼50 m. Faulting is still observed, although only in a few instances are throws larger than 100 m seen (Fig. 2, 203 km, 208 km, 222 km), and nowhere south of 190 km is the basement the sea floor. Instead, the sediment cover is laterally continuous, and has a thickness of ∼275 m around borehole 504B.

5.2 Bathymetric subsidence

Assuming constant thermal conditions at the ridge axis and its initial depth, the plate cooling model of Hasterok (2013) can be applied to predict the amount of subsidence expected following crustal formation (Fig. 8d), where subsidence rate is proportional to t^1/2, and t is the time since accretion. The ridge depth is correlated to its thermal structure, where a hotter ridge axis corresponds to a shallower initial depth (e.g. Mendel et al.2003; Dannowski et al.2010).

The predicted subsidence was calculated using a ridge depth (d0) of 2.7 km, which represents a midpoint between the present-day depth of ∼2.9 km observed at the CRR (which may not represent the spreading-averaged depth for the past ∼7 Myr) and the average depth of 2.5 km from the global depth-age study of Crosby et al. (2006). A constant half-spreading rate of 34 mm yr⁻¹ was used to match the measured age of 6.9 Ma at borehole 504B, located 233 km off-axis.

Overall, the predicted subsidence provides a good first-order approximation of the top basement depth to the south of the CRR (Fig. 8d). In areas where crustal fabric is rough, the topography created by normal faulting leads to a high variance about the mean predicted depth. However, in two regions the top basement is significantly shallower than predicted over along-profile lengths of >10 km. The first region, located between 32 and 40 km south of the ridge axis, has a steep (∼30°) axis-facing fault scarp ∼400 m in height, a gently sloping (< 5°) top surface and an overall square, block-like appearance (Fig. 2), although this may have only a limited along-axis extent (Fig. 1b). A similar feature is located between 20 and 30 km north of the ridge axis, and together these mark the edges of the zone of rough bathymetry proximal to the ridge axis. The second location of anomalously shallow bathymetry, between 163 and 190 km south of the ridge axis, also displays multiple fault scarps which are superimposed on an overall dome shape that reaches ∼200–300 m shallower than predicted for a constant d0 and half-spreading rate (Fig. 8d), making it the largest sustained, anomalously shallow feature along profile. This coincidence of rough and anomalously shallow basement appears, therefore, to contradict the expectation that, for magma-dominated, tectonically enhanced spreading which should produce rougher seafloor, the ridge axis would be colder and, hence, deeper. If so, this would make the uplift of this region even greater as indicated by the possible curves based on a larger d0, corresponding to a cooler axis (Fig. 8d). Alternatively, it may also be possible that the observed ‘uplift’ is simply related solely to greater fault throw as noted, for example in Fig. 2, by comparing fault throw between 40–90 km and 155–190 km from the ridge axis. However, this interpretation may require spreading rates lower than those observed at the CRR (Wilson & Hey 1995), and, consequently, values of M < 0.5. Consequently, we favour the former interpretation because it provides a better fit to the determined spreading rate characteristics.

5.3 Layer 2 velocity structure

In describing our models, we follow the traditional nomenclature for layers 2A and 2B based on the initial layering structure formed at the ridge axis. Rather than using absolute velocity, we define the transition between these two layers as a region of high velocity gradient (>2 s⁻¹; Gregory et al.2018), which may represent changes in lithology, alteration, porosity and fracturing (e.g. Vera et al.1990; Christeson et al.1992, 2007; Harding et al.1993; Carlson 2010). FAST inversion results in a smooth velocity model rather than one containing discrete layer boundaries. We, therefore, overlay the layer 2A/2B boundary from the forward model onto the final inversion model (Fig. 7b), to allow descriptive comparisons of layers 2A and 2B to be made.

In the final inversion model, layer 2A retains the lower P-wave velocity resolved by the SAP inversion (Fig. 6c), whereas layer 2B generally retains the higher P-wave velocity resolved by the OBS inversion (Fig. 6a), reflecting the turning depth ranges of the two individual data subsets. Overall, we find that the P-wave velocity for zero-age crust at the seabed at the CRR is higher (∼3.2–3.4 km s⁻¹) than at the Juan de Fuca Ridge (∼2.5–3.2 km s⁻¹, Nedimović et al.2008), although this may, in part, be modelling approach dependent. Within layer 2A the P-wave velocity at the basement ranges between 3.2 and 4.5 km s⁻¹, where the lowest velocities are observed in the near-axis region and progressively increase towards the oldest crust to the south. This is a characteristic of all models (Figs 5b, 6a & c), and is highlighted by the 4.3 km s⁻¹ contour. Layer 2B in the final inversion model also shows variability in velocity structure along-profile, especially between –30 to 25 km and 100 to 195 km distance from the ridge axis, which are areas associated with rough oceanic basement. Additional short-wavelength lateral velocity variations present throughout the profile likely arise from the deeper-turning arrivals within the SAP WA data.

It is conceivable that the apparent high variability in P-wave velocity is a manifestation of the juxtaposition of adjacent crustal blocks across a series of normal faults. In order to rule out this possibility, and study the variability in more detail, the final inversion model was converted into a relative velocity model (Fig. 7c) by normalizing to a 1-D reference velocity–depth profile from the CRR; thus, the relative velocity model shows lateral change in P-wave velocity with respect to zero-age crust. If normal faulting were the cause of the observed variability, then across each fault the P-wave velocity throughout layer 2 should be relatively lower on one side and higher on the other. However, this pattern is not observed. Furthermore, the locations of relative variations in layers 2A and 2B are not correlated, indicating different mechanisms, or effects, of the crustal ageing process.

The primary features revealed by the relative velocity model are:

1. relatively lower layer 2A and 2B velocity, <0.5 km s⁻¹ faster than at the ridge axis, is observed from 30 km to the north to 25 km to the south of the ridge axis, encompassing crust up to ∼1 Ma (Fig. 8d, 0 km);

2. relatively higher velocity, ∼0.5 km s⁻¹ faster in layer 2A and ∼0.8 km s⁻¹ faster in layer 2B (Fig. 8d, 70 km), is observed from 25 to 110 km south of the ridge axis;

3. compared to that expected for ageing of crust, anomalously higher velocity in layer 2A and anomalously lower velocity in layer 2B is observed from 110  to 135 km and from 150  to 190 km south of the ridge axis (Fig. 8d, 165 km); and

4. significantly higher velocity throughout the whole of layer 2, 1.0–1.4 km s⁻¹ faster than at the CRR (Fig. 8d, 233 km), is observed >190 km south of the ridge axis.

5.4 Summary comments

The results described above reveal the following trends:

• the bathymetry can be partitioned into zones matching the expected subsidence curve (10–35 km, 45–160 km and 190–250 km distance along profile) and those where it does not (0–10 km, 35–45 km and 160–190 km);

• basement roughness can be classified into zones where the perturbations around the local mean depth are high (0–45 km, 95–190 km and >250 km), and the crust is defined as rough, and low (45–95 km and 190–250 km), where it is defined as smooth; and

• there is an overall linear increase in P-wave velocity with increasing distance from the ridge axis, suggesting a gradual evolution as the crust ages. However, there are also regions with relatively slower (0–25 km, 100–190 km), faster (25–100 km) and much faster (>190 km) seismic velocities when compared to the current ridge axis.

6 SYNTHESIS CRUSTAL MODEL FOR THE SOUTHERN FLANK OF THE CRR

From the observations of the oceanic basement roughness, subsidence, and crustal velocity structure relative to the ridge axis described in earlier sections, we have developed a synthesis model (Fig. 9) of the crust comprising the southern flank of the CRR, in which the crust can be partitioned into zones, from the ridge axis southwards, comprising:

Figure 9.

Synthesis model of the crustal formation history of the southern flank of the CRR. (a) Relative velocity model (Fig. 7c). (b) Summary of crustal subsidence, roughness and layer 2A and 2B P-wave velocity anomaly, showing relationships between characteristics. (c) Half-spreading rate for the south flank of the Costa Rica Rift, from Wilson & Hey (1995; black line), and this study (red and blue lines, cf. Fig. 10). Grey-shaded region shows the proposed threshold half-spreading rate for the transition between magma-dominated, 0.5 < M < 1 (<∼33–35 mm yr⁻¹ HSR) and magmatic, M ≈ 1 (>35 mm yr⁻¹) spreading. (d) Proposed transitions between 0.5 < M < 1 and M ≈ 1 spreading, and where there appears to be a marked increase in faulting (F), over past 7 Myr.

Distance from ridge axis (km)	Descriptors	Interpretation
0–10	deep bathymetry, low P-wave velocity, rough basement	axial valley characteristics of slow-to-intermediate spreading indicative of magma-dominated (0.5 < M < 1) crustal formation
10–35	typical bathymetry, low P-wave velocity, rough basement	young (<1 Ma) ocean crust may not be ‘mature’ with respect to low-temperature, shallow hydrothermal alteration, although the rough bathymetry is indicative of magma-dominated spreading with some tectonic extension
35–45	shallow bathymetry, low 2A and high 2B P-wave velocity, rough basement	transition to more mature crust, with roughness indicative of magma-dominated, increasingly faulting-enhanced spreading—generally shallow bathymetry suggests a more buoyant crust which could be either thicker or of lower density
45–100	typical bathymetry, low 2A and high 2B P-wave velocity, smooth basement	magmatic accretion (M ≈ 1)
100–160	typical or deeper bathymetry, variable P-wave velocity, rough basement	increasing roughness with distance south of ridge axis implies some tectonic extension, but as the crust lies at its predicted depth it suggests it has a bulk density expected for magma-dominated crust—has spatially limited, localized velocity and, therefore, density anomalies
160–190	shallow bathymetry, fast 2A and slow 2B P-wave velocity, rough basement	period of magma-dominated accretion, with significantly increased faulting (defined here by the term F—cf. Fig. 9), and low density at depth accounting for the observed uplift
190–250	typical bathymetry, typical/fast(?) velocity in both 2A and 2B, smooth basement	a phase of magmatic accretion (M ≈ 1), corresponding to the lithological ground-truth at borehole 504B at 233 km along profile

It is evident that there is significant diversity in the crust formed at the CRR over the past 7 Myr. Furthermore, there appear to be distinct correlations between variability in the crustal velocity and basement roughness. Thus, the variable velocity pattern we observe may reflect the balance in the mode of formation, between magmatic and magma-dominated accretion vs an increasing prevalence of tectonic extension, the effect of which influences the permeability of the crust and the pattern of fluid flow it facilitates. We therefore suggest that the classification of the crust here as simply ‘intermediate-spreading’ is too simplistic, and that there is significant variability in its mode of formation and subsequent evolution. Hence, we consider that the term ‘intermediate’ refers instead to a spectrum of spreading styles between faster, magmatic (M ≈ 1) and slower, magma-dominated, faulting-enhanced (0.5 < M < 1) end-members, within which the actual mode of crustal growth at any time may fluctuate along this spectrum between them as summarized in Fig. 9. In order to test this hypothesis, we compare the aforementioned crustal ‘zones’ to existing (Wilson & Hey 1995) and new determinations of spreading rate and asymmetry.

7 SPREADING RATE, SYMMETRY AND CONTINUITY AT THE CRR

When considering the influence of spreading rate on resulting crustal structure, it is also necessary to include across-axis symmetry and along-axis continuity of its formation processes. During magmatic accretion the focus of magma production, at least at depth, may not be directly below the ridge axis, instead being offset to either flank (e.g. Scheirer et al.1998; Vargas et al.2018) as a result of asymmetry in plate motion vectors (influenced by far-field tectonic stresses), or due to ridge geometry re-organizations. Alternatively, during increasingly magma-dominated, tectonic-enhanced spreading, the orientation of normal faulting, and hence which flank represents the hanging wall and which the footwall, may change over time (e.g. Sauter et al.2013). Finally, where one flank is accommodating spreading via a large-scale detachment fault, magmatic accretion may occur in the hanging wall on the opposite flank (Allerton et al.2000), rather than it simply playing a passive role.

7.1 Magnetic anomaly modelling

To quantify both the spreading rate and degree of asymmetry at the CRR, modelling of ship-acquired and EMAG2v2 global compilations of magnetic anomaly data (Maus et al.2009) was undertaken. During JC114, two north–south oriented (ridge orthogonal) profiles (SAP_A and SAP_C–Figs 1, 10c and d), which form a corridor around SAP_B, were acquired which have been modelled using MODMAG (Mendel et al.2005). Fig. 9c) shows that our modelled half-spreading rates for the south flank of the CRR demonstrate good general agreement with the estimates of Wilson & Hey (1995), for distances >40 km from the ridge axis, corresponding to ages >1 Ma (Fig. 9c).

Figure 10.

Magnetic anomaly modelling showing CRR spreading asymmetry. (a) GEBCO (2008) bathymetry for the crust 250 km to north and south of the CRR. (b) CRR magnetic anomaly data. Foreground comprises available marine track-line data (see text for source) combined with that acquired during JC114, overlain on the EMAG2v2 (Maus et al.2009—faded background) downward continued to sea level for comparison, showing a clear correlation between ship-acquired anomalies and the pattern within EMAG2v2 (c) SAP_A ship (black), EMAG2v2 (red) (both solid lines) and magnetic anomalies modelled using MODMAG (Mendel et al. 2005) (dashed). Vertical dashed lines labelled at 1 Ma intervals. Shipboard data are modelled with a sea level datum, while EMAG2v2 is modelled for its datum of 4 km above sea level. (d) As (c), for Profile SAP_C with EMAG2v2 anomaly (blue). (e) Calculated half-spreading rates from modelling of Profiles SAP_A (red) and SAP_C (blue) in (c) and (d), respectively. (f) Calculated full-spreading rates for Profiles SAP_A (red) and SAP_C (blue), the grey-shaded region shows the proposed threshold FSR between 0.5 < M < 1 and M ≈ 1 spreading modes.

Following this, the degree of spreading asymmetry was determined by modelling the EMAG2v2 magnetic anomaly grid for both the north and south flanks of the CRR, having first downward continued that anomaly to sea level, and ground-truthed it against JC114 and all other available ship-acquired track-line data available from the National Center for Environmental Information (https://www.ngdc.noaa.gov; Fig. 10b). EMAG2v2 is produced by combining a range of types of track-line magnetic measurements, with gaps between tracks interpolated based on modelled oceanic crustal ages. Using the entire magnetic anomaly data set up to and including JC114 we show that, at the CRR, EMAG2v2 mirrors the anomaly pattern observed in ship-acquired magnetic measurements, sufficient to allow spreading rate modelling along flow-lines for a time frame of at least 6–8 Myr, for both the north and south flanks.

We find that spreading has been consistently asymmetric for the past 4 Myr, occurring at half-spreading rates of >33 mm yr⁻¹ to the south and ≤30 mm yr⁻¹ to the north (Fig. 10e). Throughout this period the FSR (Fig. 10f), which ultimately controls the degree of upwelling, has been consistently ≥60 mm yr⁻¹. Between 4 and 5 Ma there is a significant reduction in the FSR, to 50 mm yr⁻¹, that is partitioned relatively equally between the two flanks. The fastest FSR at the CRR of >70 mm yr⁻¹, occurs at ages >6 Ma and corresponds to crust sampled at borehole 504B that is interpreted to have resulted from magmatic accretion.

7.2 Crustal depth anomalies

Hasterok's (2013) plate subsidence cooling model shows that the oceanic basement depth at any point is related to the initial axial depth and the age of the crust and, thus, the spreading rate. Fig. 8d) shows that there is overall good agreement between the depth of the southern flank of the CRR and the subsidence curve for an HSR of 33.8 mm yr⁻¹ and an initial axial depth of 2.7 km. However, where the basement is most heavily faulted, between 155 and 190 km distance along profile, we observe uplift of ∼200–300 m with respect to a simple half-plate cooling model. This may imply that the crust formed during a phase of magmatic accretion, contrary to expectation based on the seismic velocity and bathymetry anomalies. Similar instances of anomalously shallow bathymetry are observed at ∼110 and 150 km to the north of the CRR axis (Fig. 11a), and indicate that the subsidence model does not adequately explain observed basement depths along profile. Consequently, other mechanisms need to be considered that could result in the observed uplift.

Figure 11.

RMBA gravity reduction. (a) GEBCO (2008) bathymetry the crust 250 km to north and south of the CRR. (b) Satellite-derived free-air gravity anomaly (FAA, Sandwell et al. 2014). (c) Calculated RMBA using a crustal thickness and density contrast of the crust relative to the mantle of 6 km and 700 kg m⁻³ respectively, and an intermediate half-spreading rate for the north and south flanks of 30 mm yr⁻¹. A datum shift has been applied to distribute anomalies around zero for relative polarity analysis. Dashed grey lines indicate the window in which profiles are sampled at 1 km-along-ridge intervals and averaged in (d). Sampled and averaged RMBA at 1 mm yr⁻¹ increments of HSR, assuming symmetrical spreading (grey lines), with slowest (red), intermediate (green) and fastest (blue) CRR half-spreading rates highlighted. Black dotted line shows the 4th order polynomial used for regional detrending. (e) Detrended RMBA (rRMBA) for slowest (red), intermediate (green) and fastest (blue) CRR half-spreading rates.

In the following analysis, we only consider crustal depth anomalies at distances >30–40 km from the ridge axis. The CRR axis currently lies at a depth of ∼2.9 km, ∼200 m deeper than would be expected if present accretion was magmatic, and indicating that the ridge axis may either be in a cooler thermal or reduced magma flux state, or that magmatic accretion is occurring predominately on the north flank. Alternatively, the anomalous axial bathymetry may represent a dynamic magma-induced deformation associated with dyke intrusion at the ridge axis (Carbotte et al.2006), or that in near-ridge settings the relative thermal immaturity of the newly formed (<1 Myr) crust makes comparison with the expected ridge axis depth determined from older crust invalid.

7.2.1 Crustal density and buoyancy models

In order to appraise the potential source of the observed bathymetric uplift, it is necessary to consider the density (or buoyancy) structure of the lithosphere when modelled against the observed gravity anomaly. To reduce the number variables in the density model, we calculate the residual mantle Bouguer anomaly (RMBA, Fig. 11c) from the free-air anomaly (FAA, Fig. 11b, Sandwell et al.2014) using the method of Parker (1972) to correct for the thickness of the water (density of 1035 kg m⁻³), and the gravitational effect of an average crust of density 2600 kg m⁻³ and thickness 6 km (White et al.1992), overlying a mantle of density 3300 kg m⁻³, followed by removal of the thermal effect of ridge cooling at a specified half-rate for the north and south flanks individually, consistent with the results of magnetic modelling. As the thermal effect has a sensitivity of ∼4 mGal for a range of HSRs between 25 and 36 mm yr⁻¹ (Fig. 11d), hereafter we show the calculated RMBA as an envelope, with the limits being those associated with the fastest (36 mm yr⁻¹) and slowest (25 mm yr⁻¹) CRR half-spreadng rates (Fig. 11d,e).

A long (>100 km) wavelength trend (Fig. 11d), which may represent gradual changes in crustal thickness and density with age and the progressive accumulation of sediment cover, is removed by fitting a 4th order polynomial, which preserves localized anomalies. Finally, a datum shift is applied to distribute the remaining anomalies about zero. The RMBA residuals (rRMBA; Fig. 11e) reflect either localized lateral variations in crustal thickness and/or crustal and upper mantle densities, relative to the values specified in the reduction, with positive values representing either thinner crust or higher density crust or mantle, and negative values the opposite.

Along the CRR ridge segment from west to east, there is a variation from lower to higher RMBA values (Fig. 11c), suggesting the crust may be thinner and/or the mantle colder toward the east, and that this area may be undergoing magma-dominated, increasingly faulting-enhanced spreading. Asymmetry in the RMBA is also distributed across-axis, with the largest positive anomaly present on the north flank, adjacent to the PFZ, and extending ∼100 km from the ridge axis, in contrast to the southern flank where there is less overall variation both along and across axis. The observed locations of anomalously shallow bathymetry, at 110 and 150 km on the northern flank and between 150 and 190 km on the southern flank, correlate to RMBA/rRMBA lows indicating a mass deficit at depth.

7.2.2 2-D gravity and buoyancy modelling

Indicative modelling of the rRMBA was conducted using grav2d, which applies the method of Talwani et al. (1959). The model is constructed with planar horizontal interfaces representing the top of the oceanic crust and the Moho. Here, we consider gravity anomalies associated with three possible mechanisms that would generate buoyant uplift, which we test for plausibility and consistency with our other observations using end-member cases for each. To prevent local features with limited along-segment extent biasing the resultant models, the observed rRMBA was averaged over a swath between Profiles SAP_A and SAP_C (see Fig 11a).

Three end-member models were tested: a Pratt-type isostatic model with lateral density perturbations across a constant thickness crust; an Airy-type isostatic model with a variable thickness crust with constant density contrast at the Moho; and finally, two models based on the interpretation that either crustal or mantle densities have been locally altered as a result of metamorphism associated with fluid ingress down fault pathways. For the Pratt-type model, relative density anomalies of ±60–70 kg m⁻³ in the upper 2 km (layer 2), above a constant density layer 3, were sufficient to reproduce the observed rRMBA. Using the velocity–density relationship of Carlson & Herrick (1990), this corresponds to velocity perturbations of ∼0.5–1.0 km s⁻¹, which are consistent with the observed lateral velocity contrasts in the upper crust (Fig. 7c). However, the modelled density contrasts are less than half that required to produce the buoyancy necessary to support the observed uplifted zones. For the Airy-type model, a two-layer crust with constant densities was used with a relative density contrast of 300 kg m⁻³ at the base of the crust, which then requires about ±1 km depth perturbation to fit the rRMBA. Although this model is consistent with the degree of buoyancy required to support the uplifted zones, it is inconsistent with the FSR model of Bown & White (1994), which suggests that crust formed at slower FSRs is thicker than that formed at faster rates, as the Airy-type model predicts thicker crust is formed during periods of slowest spreading. However, these periods correspond to those that we interpret as being formed during periods of magma-dominated accretion, implying a thinner rather than thicker crust.

Finally, we test models where the density differences are the result of alteration of the lower crust and/or mantle by fluid ingress. The magnitude of the associated gravity anomaly and buoyant uplift are both strongly dependent on the density difference assumed between the alteration products and the surrounding host rock. We test two end-members, in which all of the metamorphism is accommodated either within the lower crust, from the bottom up (Fig. 12b), or the upper mantle, from the Moho down (Fig. 12c). The density of the altered material in both cases is set to 2900 kg m⁻³. If metamorphism occurs solely in the lower crust, the rRMBA can be matched using columns of material ∼3.0–3.5 km in height, which have a calculated positive buoyancy effect (i.e. uplift) of ∼0.17–0.18 km, which is within the observed error bounds. If, instead, metamorphism occurs in the upper mantle, for example with peridotite being converted to serpentinite, an ∼1 km thickness is required to reproduce the rRMBA, generating a corresponding uplift of ∼0.2 km, consistent with observations. Alternatively, it is possible that lower density material in the mantle may reflect frozen basaltic melts (e.g. Lizarralde et al.2004). However, their location coincides with zones interpreted as having formed during periods of magma-dominated accretion, from which we conclude that this explanation is unlikely.

Figure 12.

End-member models of density anomalies that fit the residual mantle Bouguer anomaly (rRMBA) after removal of the long wavelength (>100 km) trend. (a) Observed and modelled rRMBA. Grey line shows the observed rRMBA, overlain with the 2-D-modelled curves for the two models below, coloured black where the density anomalies are located entirely within layer 3 of the crust, and red for within the mantle. (b) and (c) End-member models where alteration occurs solely in the crust and mantle respectively, which produce the rRMBA shown in (a). Shaded zones centred on the present-day spreading axis and at the edges of the model are not considered reliable due to either distal or near-ridge tectonic complexities. Note that a 200 m thicker crust is modelled between 175 and 260 km to better fit the ‘background’ RMBA in this region. See text for discussion of model parametrization and modelling approach.

Gravity modelling is inherently non-unique such that a large number of potential models can fit an observed anomaly. It was not, therefore, possible to produce a single model of the lower crustal and upper mantle density structure that satisfies all of the rRMBA, variation in basement depth, and crustal velocity variability in layer 2 observations, which we have shown to be related in their manifestation and origin, and coupled both to each other and variations in full-spreading rate and asymmetry. However, each of the models outlined above variously support different characteristics of a model of crustal formation in which increased normal faulting at the seabed is a result of variation in magma supply and phases of increased tectonic extension, and that decreases in density result from increased porosity and/or permeability that allows fluid-driven alteration to occur within the crust and mantle. Thinner crust would be observed during phases of magma-dominated accretion where the crust undergoes a higher degree of tectonic extension as part of spreading.

8 DISCUSSION

8.1 Crustal velocity evolution

For magmatically accreted crust (M ≈ 1), layer 2A comprises extrusive pillow lavas and high porosity at zero-age. Large-scale seismic studies (e.g. Houtz & Ewing 1976; Grevemeyer & Weigel 1996, Carlson 1998) have demonstrated that during the crustal ageing process P-wave velocity increases with age, doubling from its initial value within ∼10 Myr. This progression is also observed in the results of our study, by the shallowing of the 4.3 km s⁻¹ contour in the final inverse model (Fig. 7b), from approximately the depth of the layer 2A/2B transition at zero-age, to the top basement surface depth by 170–220 km from the ridge axis. The velocity increase is principally attributed to the infilling of porosity with leached mineral deposits (e.g. Christeson et al.2007), from the base of layer 2A upwards, as a result of hydrothermal fluid flow.

In contrast, velocities in layer 2B tend to increase at a much higher rate, reaching ‘maturity’ within ∼0.5 Myr of crustal formation (Newman et al.2011), which may indicate that layer 2B velocity evolution is rapid and dissociated from that of layer 2A, consistent with models of hydrothermal fluid circulation where high-temperature systems penetrate into layer 2B in the proximity of an axial magma lens at the ridge axis, while low-temperature hydrothermal systems are restricted to layer 2A and continue only as long as there is permeability. However, for crust that is formed during magma-dominated periods with a higher degree of tectonic extension (0.5 < M < 1), the normal faults that facilitate that extension overprint and predominate the primary porosity and permeability, and act as open conduits for the ingress of fluids and subsequent chemical stripping of minerals, possibly allowing fluids to penetrate deeper into the lithosphere and remain active over significantly longer periods of time.

Within the oceanic crust, hydrothermal systems are strongly influenced by the degree of sediment cover (e.g. Nedimović et al.2008), with sealing of the basement to fluid flow occurring once complete cover by low permeability pelagic ooze is established. South of the CRR, complete cover occurs at distances of >190 km for crust older than ∼5.9 Myr (Fig. 2). This system sealing results in changes in the degree and pattern of hydrothermal circulation, and coincides with the highest observed layer 2A P-wave velocity anomalies in the final inversion model (Fig. 7c). Where the sediment cover is discontinuous (Nedimović et al.2008), the rate of layer 2A velocity increase since formation is observed to be slower. We infer that where the system remains open, heat is efficiently advected from the crust and results in a lower temperature in layer 2A that, in turn, preserves the permeability. Once the system is sealed though, the basement temperature increases by an order of magnitude (Kolandaivelu et al.in review) and results in the mobilization and precipitation of minerals that reduce the porosity and increase velocity. The ability of the sediment cover to overwhelm the underlying crustal topography is primarily a function of sedimentation rate, but also the topography inherited at the time of crust formation at the ridge. Where the crust is most heavily faulted, between 100 and 190 km distance from the ridge axis, the basement of fault footwalls is still exposed to seawater, as is young crust between ±30 km of the ridge axis where sediment cover is yet to accumulate (Fig. 2). These are locations that coincide with observations of reduced layer 2B P-wave velocity (Fig. 7c), indicating that these exposed-to-seawater fault pathways remain open to fluid circulation, facilitating extensive focussed flow and, thus, the alteration of this and the deeper crust.

8.2 Spreading rate and crustal formation style

Spreading rate can be used as a proxy for the thermal structure of a ridge axis. At fast spreading rates, magma supply is effectively constant and accretion is effectively entirely magmatic (M ≈ 1). At slower rates, magma supply may be intermittent or sparse, with spreading accommodated increasingly by faulting (0.5 < M < 1). At the slowest spreading rates, the ridge system may even approach an amagmatic state. The effect of spreading rate on crustal formation depends on both the full-spreading rate, which effectively controls the amount of magmatic upwelling, and the degree of spreading asymmetry which partitions the full-spreading rate into half-spreading rates for each ridge flank. Over the past 7 Myr, the CRR has spread at full-spreading rates of 50–72 mm yr⁻¹, with the south flank HSR ranging from 25 to 36 mm yr⁻¹ (Fig. 9c; Wilson & Hey 1995), covering the full range expected for intermediate rate ridges (Dick et al.2003).

Lithological analysis of samples from borehole 504B shows that ∼6.9 Ma crust, now located 233 km south of the ridge axis, was accreted magmatically (Alt et al.1996). Between 190–250 km south of the ridge axis, the bathymetry matches the d0 = 2.7 km subsidence curve (Fig. 8d), the crustal topography is smooth (Fig. 8c) and the layer 2B relative P-wave velocity is high (Figs 7c and 9b). This region of crust corresponds with the fastest full-spreading rates at the CRR of over 70 mm yr⁻¹ (Fig. 10f). The ground-truth provided by borehole 504B allows the observations of smooth basement topography and high relative layer 2B P-wave velocity to be used as an indicator of magmatic accretion of M ≈ 1. Between 50 and 100 km south of the CRR, although the full-spreading rate is slower by 5–10 mm yr⁻¹ and there is a significant asymmetry in the HSR, the crustal morphological and velocity characteristics are still consistent with a predominantly magma-dominated origin, yielding crust similar to that at borehole 504B.

In contrast, the crust located between 140 and 190 km south of the ridge axis (4.0–5.8 Myr) was formed at full-spreading rates of 50–55 mm yr⁻¹, corresponding to the slower end of the intermediate range (Fig. 10f). In this region there is evidence for increasingly faulting-enhanced spreading in the form of normal faulting (Fig. 2) and high basement roughness (Fig. 8c). The coincident lower velocity of layer 2B (Fig. 9a) indicates that the fault pathways here remain open and facilitate on-going hydrothermal circulation. The conjugate to this region on the north flank would also be expected to show similar features if crustal creation was symmetric. Although we have not modelled the velocity–depth structure of the crust to the north, there is evidence from the rRMBA (Fig. 11e) of two lows at ∼–110 and –140 km which occur within an equivalent time frame (∼4–6 Myr). The coincident observations of reduced crustal P-wave velocity, rRMBA lows, and bathymetric uplift, together suggest that increased faulting during periods of magma-dominated, faulting-enhanced spreading increases crustal porosity, potentially allowing fluid ingress to depth within the crust, or even into the upper mantle. Such fluid flow would produce metamorphism and result in density anomalies.

8.3 Transitions between modes of spreading

The rate of transition between magmatic and magma-dominated phases of crustal formation, and the time frame of complete change from one to the other, provides a means to estimate the temporal variability in longer-term thermal state of the ridge axis and deeper magmatic source. In the absence of any other controls, the full-spreading rate is expected to be directly correlated to the degree of mantle upwelling and, hence, melt supply, and should be partitioned equally between the two ridge flanks. Under these circumstances, the increase in tectonism may not occur instantaneously (in geological timescales), instead taking place as a transition as magma flux wanes and the ridge axis cools. However, the response to increasing magma flux might be expected to be more rapid, effectively reducing faulting as spreading becomes supported by crustal growth.

Our observations show that for the CRR, this ‘passive’ progression (or cycle) does not actually appear to be the case. Instead, the magmatic accretion documented and ground-truthed around borehole 504B is followed not by gradual slowing of the spreading rate, but by a rapid ∼20 mm yr⁻¹ decrease in FSR to ∼50–55 mm yr⁻¹, partitioned roughly equally between the two ridge flanks, and occurring at ∼6 Ma. The following period during which faulting is enhanced is relatively short, with the magma supply beginning to increase, together with a matching increase in the FSR, over the subsequent 2 Myr, with magmatic accretion preferentially occurring to the south of the ridge axis due to asymmetry in the spreading (Fig. 10e).

Regardless of the specific driving/controlling force, the structure of the oceanic crust formed at the CRR appears to be highly sensitive to relatively small changes in spreading rate. Patterns in the observed top basement roughness and layer 2 velocity structure appear to follow the variation in FSR, with 62–65 mm yr⁻¹ here acting as the threshold for the transition between spreading states.

If the magma supply sufficiently decreases such that faulting-enhanced spreading progresses to a tectonic-dominated mode, eventually a large-scale detachment will form that could result in the exhumation of the lower crust and the uppermost mantle at the seabed (e.g. Cannat 1993; Ranero & Reston 1999; MacLeod et al.2002; Escartin et al.2008; Dannowski et al.2010; Peirce et al.2019). Modelling suggests that such large-scale crustal detachments form if M becomes less than ∼0.3–0.5 (Buck et al.2005; Tucholke et al.2008; Olive et al.2010). Instead, at the CRR the magmatic system operates, at a minimum, in an episodic state of ‘just in time’ magma delivery, such that M does not fall to values indicative of tectonic-dominant spreading.

8.4 Possible controls on spreading rate and asymmetry variations

Over the past 6 Myr, multiple plate motion changes have occurred that may have influenced the stress regime at the CRR and, hence, contribute to the changes in the character of spreading described in this study. The rapid reduction in spreading rate at the CRR at ∼6 Ma (Fig. 9c) corresponds to a clockwise shift in the absolute motion of the Pacific Plate at 5.8–5.9 Ma (Krijgsman et al.1999). Plate reconstructions (e.g. Morell 2015) also show that the separation of the Ecuador Ridge from the Galapagos Ridge occurred between 4 and 6 Ma, which resulted in redistribution of transform motion along the ridge system. Around 4 Ma, spreading at the Malpelo Ridge is believed to have stagnated resulting in the reorientation of the Panama Fracture Zone (formerly the Coiba Fracture Zone; Lonsdale & Klitgord 1978) and, potentially, an increase in extension elsewhere within the basin. A possible manifestation of this is observed at the CRR, which shows an overall gradual increase in spreading rate, focussed to the south, and transition in magma supply state over several million years. The Cocos and Carnegie aseismic ridges are also proposed to have encountered the subduction zones of Central and South America, between 1 and 2 Myr (Lonsdale & Klitgord 1978; Gutscher et al.1999; Meschede & Barckhausen 2001). It is possible that this may have led to a suppression, impediment or slowing of northward and eastward plate motion. This timing coincides with the transition in crustal structure observed at ∼35–45 km to the south of the CRR, where the roughness and uplift increase. There is also a concurrent decrease in the FSR and increase in the asymmetry at this time (Fig 10e).

Apparent spreading asymmetry may, alternatively, result from migration of ridge-axis discontinuities (e.g. Cormier & Macdonald 1994; Grevemeyer et al.2002), which may transfer lithospheric material from one flank of the ridge to the other. Bathymetric data along the CRR axis identify a left-stepping non-transform ridge discontinuity at 3°20′N, 83°44′W, which displays an overlap and offset between the two segments of ∼2.5 and ∼1.5 km, respectively (Haughton et al.2017). However, this feature is much smaller than the well-developed overlapping spreading centres on which the plate transfer hypothesis is based, and such overlapping spreading centres are noted to manifest low upper crustal velocity anomalies beneath them (e.g. Lonsdale 1983; Christeson et al.1997; Canales et al.2005), to an extent that is not observed in the final inversion model. Consequently, this mechanism is considered unlikely as an explanation of the observations at the CRR.

Ridge jumps represent another mechanism which may manifest as apparent spreading asymmetry. This process has been observed in a number of locations (e.g. Brozena & White 1990; Small 1995) and is often associated with proximity to a hotspot. However, there is no evidence of repeated paired chrons within the coverage and resolution of the total ship-acquired magnetic anomaly data (Fig. 10b), which indicates that this mechanism is unlikely to have occurred at the CRR.

Morphological analysis of the bathymetry in the vicinity of the CRR shows evidence for extensive volcanism in the recent geological past, in the form of seamounts and lava domes and flows (Haughton et al.2017). These features indicate that the CRR experiences intermittent and spatially discrete periods of focussed enhanced magmatism. The asymmetric distribution of these features, being located predominantly on the north flank of the CRR, is consistent with previous observations which suggest that asymmetric distribution of seamounts across spreading ridge flanks may be an indication that ridge migration has occurred (e.g. Davis & Karsten 1986) or that the deep heat source at the CRR is displaced to the north, as suggested by seismic attenuation studies (e.g. Vargas et al.2018).

Finally, an alternative to the spreading rate control on magma supply and, hence, resultant crustal structure would be variability in the mantle composition providing the source melts. Correlations between melt composition, ridge-axis topography and crustal thickness have been shown to occur across a range of localities (Klein & Langmuir 1987). However, the results of geochemical analyses of recovered samples from Deep Sea Drilling Project sites south of the CRR, indicate a near steady-state in magma composition between 3.9–5.9 Ma (Natland et al.1983), which may not, therefore, support this mechanism.

Consequently, none of these possible alternatives to the spreading rate control on magma supply is considered likely to explain all of the observed features of the CRR, at least for the last ∼6–8 Myr.

9 GLOBAL CONTEXT

Our results from the CRR reveal variability in both the crustal morphology and layer 2 P-wave velocity, which manifest over a range of spatial and temporal scales. Previous studies of layer 2A at intermediate-spreading ridges have, to date, predominantly focused on mapping the thickness of this layer and the depth to the axial magma lens, where present, in order to determine ridge thermal structure (e.g. Blacic et al.2004; Baran et al.2005; Canales et al.2006; van Ark et al.2007), while resolution of the observed layer 2B variable velocity structure has been relatively limited, due to experiments consisting of ridge‐parallel profiles (e.g. Grevemeyer et al.1998) or profiles not extending sufficiently far-enough off-axis to sample the necessary time frames over which variation occurs (Baran et al.2010). However, layer 2B velocity anomalies of up to ∼1 km s⁻¹ have been observed to be associated with areas of several hundred metres of basement topography across the Northern Symmetric and Cleft segments of the JdFR (Newman et al.2011), consistent with the observations of layer 2B characteristics at the CRR described here. Furthermore, based on observations of variable axial morphology along the Cleft and Vance segments of the JdFR (Canales et al.2005) and along the GSC (Blacic et al.2004), it has been proposed that along-axis spatial variability in morphology of the crust formed at the ridges can be explained by relatively subtle changes in extent of magma supply, as demonstrated by the models of Kappel & Ryan (1986) and Phipps Morgan & Chen (1993). Our study, therefore, suggests that these same models can be applied generically to temporal variations in crustal structure variability observed across axis, and driven by the same processes.

Thus, the CRR, and by extension all intermediate spreading ridges with spreading rates between the fast end of Atlantic ‘slow spreading’, and the slow end of EPR ‘fast spreading’, is likely to express variations in the mode of crustal formation. The threshold for the transition between magmatic and magma-dominated, faulting-enhanced spreading modes appears to be controlled by the full-spreading rate and degree of spreading asymmetry between the two flanks. For the southern flank of the CRR, this threshold full-spreading rate appears to be ∼62–65 mm yr⁻¹. There appears to be several different scales of magma supply occurring at the CRR and, by extension, intermediate-rate spreading centres in general. The first is a longer-period (several Myr) variability in the overall degree to which spreading is supported by magmatism, and which may reflect either waxing and waning of a deeper mantle source, or the effects of changing regional tectonic stresses. Secondly, during periods of magmatic spreading, potentially including at the present-day, there is evidence for a small axial magma lens (or lenses) at the ridge axis (Buck et al.1997; Zhang et al.2017), the relatively small volume of which indicates that it, or they, must be subject to regular short-period episodic replenishment, otherwise ‘freezing’ quickly occurs (Liu & Lowell 2009). Finally, there is evidence for intermittent and spatially discrete periods of focussed enhanced magmatism (Haughton et al.2017).

The evolution of upper crustal properties with time has been shown to be strongly dependent on the mode of crustal formation, and the extent and continuity of hydrothermal alteration. How long hydrothermal circulation can continue is governed not only by temperature (to act as a driving force), but also by the availability of potential fluid pathways. Thus, in addition to the significant role played by the crustal fabric imposed at accretion, the rate of crustal ageing is also controlled by the continuity of sediment cover (Nedimović et al.2008). Furthermore, there may be an important relationship between the locations of significant basement faulting at the time of formation, sediment insulation and the evolution of heat flow from the crust to the ocean (Kolandaivelu et al.in review).

10 CONCLUSIONS

Joint analyses of MCS, SAP WA and OBS WA data, with additional constraint provided by potential field data, has revealed the detailed structure of the upper crust formed at the Costa Rica Rift over the last ∼6–8 Myr. Overall, the results of our study indicate the following:

• The upper crustal formation process at the intermediate spreading Costa Rica Rift is sensitive to relatively small changes in full-spreading rate (50–70 mm yr⁻¹) with corresponding changes in half-spreading rate within the range of 25–36 mm yr⁻¹ which are asymmetric about the ridge axis. These fluctuations correspond to different characteristics of the spreading processes occurring at the ridge axis that imprint strongly on the resulting primary fabric of layer 2 and influence its subsequent evolution through secondary alteration. Even at intermediate spreading ridges there appear to be two end-member modes of crustal formation: one magmatic (M ≈ 1) which occurs at the faster end of the spectrum, and the other magma-dominated and accompanied by enhanced tectonic extension (0.5 < M < 1), that occurs at the slower end. Generally, there is a fluctuating transition back and forth between the two, finely balanced at a threshold spreading rate that may be ridge location specific, and dictated by the broader regional plate tectonic setting.

• During times of magmatic spreading, corresponding to full-spreading-rates >∼62–65 mm yr⁻¹ at the CRR, the resulting top basement surface is smooth and the crust subsides over time at the rate predicted by a simple half-plate cooling model. The across-axis P-wave velocity structure of both layer 2A and layer 2B displays an overall trend of increasing velocity with age, with the former, shallower layer experiencing greater change than the latter.

• During periods of slower spreading, with a reduced magma supply, a rough top basement is observed together with pervasive faulting, and the depth to the top basement surface is shallower than expected based on cooling models. This observation is similar to that recorded at slow- and ultraslow-spreading ridges, where oceanic core complexes are ultimately exhumed along low-angle detachment faults. Deeper fluid ingress along faults and through layer 2B, characterized by a lower than expected P-wave velocity, leads to alteration of olivine-rich mineral assemblages within the lower crust or upper mantle, lowering their density and resulting in bathymetric uplift.

• The potential for, and resulting effect of, hydrothermal alteration is directly related to the permeability structure of the crust and is, therefore, strongly influenced by the initial structure inherited at the time of accretion. For smooth crust, once a continuous sediment cover is established, hydrothermal circulation occurs in a closed system and the P-wave velocity and density of the crust increase due to secondary mineralization. For rough crust with discontinuous sediment cover, significant hydrothermal circulation exchange with the overlaying ocean may continue off axis, despite the full sealing and healing of primary porosity, by exploiting fault pathways created by increased tectonic extension at the time of formation.

• It remains unclear what drives the observed changes in spreading at the CRR. However, the timing and sequence of transitions between different modes of accretion suggest a several-Myr cyclicity in broad-scale magma supply may not be the controlling factor and, instead, these changes may be influenced by tectonic events associated with the plate boundary, or changes in the spreading rate or patterns of spreading of adjacent plates.

• Magmatic and magma-dominated, faulting-enhanced modes of spreading are not mutually exclusive and, instead, the oceanic crust at any location can be considered as forming via a mix-and-match of these two mechanisms, with varying predominance and/or durations of each and where the formation style may be different for crust on opposite ridge flanks. At ‘intermediate’ spreading rates there is, thus, a fine balance between these two modes. Therefore, the term ‘crust formed at an intermediate-spreading ridge’ is not an accurate description of the underlying crustal formation process(es). Instead, it represents a temporally and spatially integrated view of the crustal structure over a much longer time span, where variations in formation mode have occurred on timescales of the order of a few Myr.

ACKNOWLEDGEMENTS

This research was funded by the Natural Environmental Research Council (NERC) grant NE/I027010/1. We would like to thank all those involved in the planning and acquisition of data during research cruise JC114 and SO238 including the officers, engineers and crew of the RRS James Cook and the FS Sonne, the scientific party, and all seagoing NERC facility technicians and engineers. The NERC Ocean-Bottom Instrumentation Facility (Minshull et al.2005) provided the OBSs and their technical support at sea. The MCS data were processed using Claritas™ and Seismic Unix (Cohen & Stockwell 2010). All figures were prepared using the Generic Mapping Tools (GMT–Wessel & Smith 1998). All data from cruise JC114 are archived at the NERC's British Oceanographic Data Centre and are available on request from the OSCAR PIs, and the final accepted version of this manuscript is available through Durham Research Online (dro.dur.ac.uk). We thank the reviewers for their helpful comments that improved the clarity of this paper.