Probing layered arc crust in the Lesser Antilles using receiver functions

Oceanic arcs can provide insight into the processes of crustal growth and crustal structure. In this work, changes in crustal thickness and composition along the Lesser Antilles Arc (LAA) are analysed at 10 islands using receiver function (RF) inversions that combine seismological data with v_P/v_S ratios estimated based on crustal lithology. We collected seismic data from various regional networks to ensure station coverage for every major island in the LAA from Saba in the north to Grenada in the south. RFs show the subsurface response of an incoming signal assuming horizontal layering, where phase conversions highlight discontinuities beneath a station. In most regions of the Earth, the Mohorovičić discontinuity (Moho) is seismically stronger than other crustal discontinuities. However, in the LAA we observe an unusually strong along-arc variation in depth of the strongest discontinuity, which is difficult to explain by variations in crustal thickness. Instead, these results suggest that in layered crust, especially where other discontinuities have a stronger seismic contrast than the Moho, H–kstacking results can be easily misinterpreted. To circumvent this problem, an inversion modelling approach is introduced to investigate the crustal structure in more detail by building a one-dimensional velocity–depth profile for each island. Using this method, it is possible to identify any mid-crustal discontinuity in addition to the Moho. Our results show a mid-crustal discontinuity at about 10–25 km depth along the arc, with slightly deeper values in the north (Montserrat to Saba). In general, the depth of the Moho shows the same pattern with values of around 25 km (Grenada) to 35 km in the north. The results suggest differences in magmatic H₂O content and differentiation history of each island.

1. Introduction

1.1. Overview

Subduction zones are regions on Earth where new continental crust is thought to have formed; however, even though the origin of continental crust has been studied for a long time, major details, such as the discrepancy in composition between average continental crust and that beneath many island arcs, remain unclear [1,2]. A better understanding of crustal structure provides insight into the link between subduction processes and the formation of continental crust through arc volcanism.

The Mohorovičić discontinuity (Moho), the boundary between the crust and the mantle, marks a sharp change in seismic velocities that is thought to be due to chemical composition and/or rheology changes. In addition to the Moho, a mid-crustal discontinuity (MCD), in some areas referred to as the Conrad discontinuity [3], can be observed in many subduction environments (e.g. [4,5] and references therein). Even though it is normally found to be the dominant crustal discontinuity, the Moho is sometimes weak and difficult to resolve, especially beneath volcanic arcs [2,4].

Ideally an active source along-arc seismic experiment would be carried out to provide a comprehensive investigation of crustal structure, but long offsets and large sources are needed. Furthermore, this option is expensive and long offsets, which are required to guarantee that the Moho will be visible, might not be easy to acquire in many regions such as curved island arcs (e.g. [4,6]). Passive–seismic observations in a well-monitored arc setting provide an alternative approach. Receiver functions (RFs) and related H–k stacking are now a common method for studying crustal structures (e.g. [2,7–15]).

The H–k stacking method makes use of the difference between P- and S-wave velocities to estimate the depth of discontinuities at which strong changes in seismic velocities occur (termed H) and the average v_P/v_S ratio between the receiver and the discontinuity (termed k; [8]). The ratio of the P- and S-wave velocities can be used to better constrain the average material properties that are present in the crust between the surface and the discontinuity [16]. To investigate crustal thickness, H–k stacking is normally used with the assumption that the largest P-to-S conversions occur at the Moho.

In this work, crustal structure variation along the Lesser Antilles Arc (LAA; figure 1) is studied to investigate potential influences of subduction on the overlying crust. Our approach integrates seismology and petrological observations with a specific emphasis on the LAA. We use extensive seismic data from 26 stations on 10 islands and use RFs to explore crustal discontinuities along the arc. This is complemented by published work on reconstructed crustal structure and compositions of fossil and currently active arcs (e.g. [2,21–23]). The results are compared with structural features and Moho depth estimates from previous works, to propose hypotheses about the link between subduction-related processes and the crustal structure beneath the LAA.

1.2. The Lesser Antilles Arc

The LAA (figure 1), extending some 800 km northwards from the South American continent to the Greater Antilles, is an expression of slow (18–20 mm yr⁻¹), westward subduction of Atlantic oceanic crust (North and South American Plates) beneath the Caribbean Plate (e.g. [24]). Subduction is sub-orthogonal in the vicinity of Martinique, with sinistral obliquity to the north and dextral to the south [17]. A comprehensive review of the geological and tectonic setting of the LAA is provided by Smith et al. [25] and summarized briefly here.

Magmatism along the LAA dates from the Eocene [26]. The present-day arc consists of 11 major volcanic islands; a string of 19 small islands (the Grenadines) lies between St. Vincent and Grenada. Several volcanic centres in the LAA are currently or recently active, including those on Montserrat, Martinique, Dominica and Guadeloupe and the submarine Kick-'em Jenny [27]. The variation in size and spacing of the volcanic islands reflects spatial and temporal variation in magmatic output. Magma production rates in the LAA are at the low end of intra-oceanic arcs worldwide (162 km³ km⁻¹ Myr⁻¹ [28]), possibly a consequence of slow convergence. The LAA lies on the eastern margin of thickened oceanic crust of the Caribbean Plate [29], although the extent to which vestigial Caribbean Plate material, attenuated or otherwise, is present beneath the LAA is not known.

An unusual feature of the LAA is a marked bifurcation north of Martinique into an inactive eastern limb and active western limb; there is evidence for an abandoned, Mesozoic volcanic arc to the west of the Lesser Antilles (the Aves Ridge), separated from the active arc by the Grenada Basin [30]. The westward jump in the northern LAA accounts for a hiatus in volcanism there between mid- and late Miocene times. The cause of the jump is not known. However, there is a marked change in the dip and orientation of the Wadati–Benioff zone along the arc, interpreted by Wadge & Shepherd [24] to indicate that either a single American plate was torn and deformed during subduction or that the North and South American plates were subducted with different velocities. In either case, bifurcation of the LAA just north of Martinique would correspond to the triple junction where Caribbean, North American and South American plates meet.

The division between the northern and southern parts of the LAA is also reflected in the presence and character of the sedimentary cover. The incoming plate in the south is rich in clastic detritus from the South American continent, partially scraped off to form the Barbados accretionary prism (e.g. [31]). To the north, sediment supply is limited by the presence of submarine highs, such as the Tiburon Ridge, and the incoming plate is blanketed by pelagic marine sediments. The spatial variation in sediments plays a role in changing magma chemistry along the LAA [32].

The LAA transects five major fracture zones on the down-going plate (figure 1). The down-going plate displays the topographic expression of strong tectonic extension including normal faults with large amounts of rotation and dome-shaped faulted detachment surfaces, or core complexes, at the edge of the inner valley floor. The presence of serpentine in the down-going plate, associated with fracture zones and/or core complexes, could introduce significant H₂O to the mantle wedge, perhaps accounting for along-strike variation in magma productivity [33] and subduction zone seismicity [20].

1.3. Crustal structure of the Lesser Antilles Arc

The LAA has been the subject of three major geophysical experiments designed to elucidate crustal structure [4,34,35], as well as an attempt to map an along-arc transect of crustal thickness using RFs and H–k stacking [15]. In [12], they have used the same method on the island of Montserrat. Estimates of crustal thickness (depth to Moho) beneath the arc from these studies range from 22 to 37 km. In [4], from their along-arc survey of the southern and central part of the arc (Grenada to Guadeloupe), they reveal the presence of two refractors that split the crust into discrete layers. Their upper layer, with an average velocity of 6.2 km s⁻¹, has significant along-strike variation in depth and velocity. The average upper layer thickness is 10 km, but varies from 2 to 20 km [4]. In [36], they interpreted this layer as being built of dense, solidified volcanic rock and plutons of intermediate composition. The uppermost portion of the upper layer has significantly lower seismic velocities and densities, and is likely to be composed of volcaniclastic and sedimentary rocks with abundant fractures [35]. Gravity data from Guadeloupe [37] show that this uppermost layer (v_P < 6 km s⁻¹) is approximately 4 km thick. The lower crustal layer of [4] that immediately overlies the mantle has average v_P = 6.9 km s⁻¹ and is thought to represent dense, more mafic igneous rocks, including cumulates.

Two cross-arc seismic surveys, between Dominica and Guadeloupe [35] and south of Grenada [34], provide a more detailed picture of crustal structure. The layering persists, although the overall vertical velocity gradient is smoother than that suggested by Boynton et al. [4]. The crust between Dominica and Guadeloupe is 26 km thick and south of Grenada it is 24 km thick. The Moho is not well resolved in either location; mantle v_P varies from 7.7 km s⁻¹ in the south to 8 km s⁻¹ in the centre. Neither survey shows any significant deepening of the Moho beneath the active arc. West of Grenada, the Moho shallows beneath the Grenada Basin (to less than 20 km in places), thickening to 27 km beneath the Aves Ridge. Seismic surveys have been unable to identify unequivocally any vestiges of Caribbean Plate in the sub-arc crust.

Here, we refine the along-arc crustal image of [4] using seismic data collected from 29 remote seismic stations along the active LAA, in combination with insights from experimental and igneous petrology, to develop a method for constraining crustal structure using RF. Our approach, which refines a recent investigation of Moho depths along the LAA using conventional RF analysis [15], has widespread applicability to volcanic arcs and layered crust more generally.

2. Receiver functions and H–k stacking

We use seismic broadband stations from different networks, located on most of the major islands of the LAA (figure 1 and table 1); for Montserrat, Guadeloupe and Martinique, more than one station is available. For this study, we limit our catalogue to events greater than magnitude 5.5. Teleseismic events are required (30° to 90° distance) to ensure subvertical incidence angles. Events are filtered using a second-order Butterworth bandpass filter from 0.4 Hz to 3 Hz (after [13]). Only events with a clear P-phase are then selected for this study.

The method uses the coda of an arriving signal, which contains mode-converted energy due to the structure beneath the receiver [38,39]. A large velocity contrast at a seismic discontinuity leads to a strong P-to-S-converted phase [8]. The signal at the receiver is a convolution of the initial signal with the subsurface structure. Therefore, assuming horizontal layering, a deconvolution can be carried out to remove the source effects and produce a sequence of pulses representing this structure by isolating the P-to-S conversions [8,13,39,40]. This resulting sequence is called a ‘receiver function’ (RF).

In this study, the extended-time multitaper frequency-domain cross-correlation receiver function (ETMTRF [40]) is used to create the RFs. ETMTRF, based on the work of [41], includes later arriving multiple converted phases and has the advantage of being less sensitive to noise. We use a high-frequency cut-off at 1.5 Hz, and for the purpose of this study three overlapping tapers are sufficient. The radial RFs are then stacked by jackknifing [42] from −10 s to + 30 s relative to the P-peak. This produces a standard variation that can be used as a pointwise uncertainty for the RF. Here, the 2σ level is used. The H–k stacking method follows the work of [8–10], which involves applying a bootstrapping algorithm [43] to determine the uncertainties of the model parameters.

Based on theoretical arrival times of converted phases, the method derives values for the depth of the discontinuity (H) and the average P-wave to S-wave (v_P/v_S = k) ratio of the crust between that point and the surface. The amplitudes at the theoretical arrival times are summed as follows:

2.1

where N is the number of RFs used in the stack, r_n(t) is the RF amplitude at time t, which is the predicted arrival time for the individual phases (the indices are 1 for Ps, 2 for PpPs and 3 for PsPs/PpSs), and w₁, w₂ and w₃ are the weighting factors with .

The weighting factors are chosen so that phases that are more apparent will be enhanced [9,13,44]. During the bootstrapping process, we compute 300 iterations, while changing the value of v_P using random, normally distributed values with 95% of the values in a range of Δv_P = ±0.3 around a mean of 6.5 km s⁻¹. At the same time, the weighting factors are chosen so that 95% fall in Δw = ±0.05 around the means of 0.6 for w₁ and 0.3 for w₂ (with ). The best estimation occurs at the location in H–k space where the three phases are stacked coherently [8] among all possible cases in terms of varying v_P and the weighting factors.

As the assumed value of v_P can also be a source of error [8,14], the mean of v_Pis changed in subsequent calculations from 6.5 to 6.3 and 6.7 km s⁻¹. These values are maximum deviations of the crustal mean v_P as determined by Boynton et al. [4]. Excluding the noisy part of data from Dominica, our test shows that the bootstrap error is either similar (Martinique, Montserrat, St. Vincent) to the uncertainty due to the change in v_P or smaller (in all other cases). In the case of the v_P/v_S ratio, the bootstrap error is always larger. The larger uncertainty between the bootstrap error estimate and the error arising from a change in the mean v_P is used in this study.

In the case of Martinique, we can see a strong, well-resolved discontinuity at 28.3 ± 1.1 km (figure 2a). On St. Lucia, the discontinuity is placed at 46.5 km, much deeper than expected in this area and the solution is very poorly constrained. Furthermore, on some islands H values are arguably too shallow (e.g. St. Vincent with 19.9 ± 0.7 km, figure 2b) to be the Moho. We note that the data are in general noisy, as they are from island stations (figure 2c,d). To help mitigate this issue, we use the RFs with the highest signal-to-noise (S/N) ratio and employ stacking to further improve the S/N ratio.

An explanation for the unexpectedly shallow or deep results is the existence of a weak Moho beneath some of the islands and a stronger MCD, in which case H–kstacking cannot resolve the Moho, returning instead the depth (H) of the MCD as the preferred result. Additionally, near-surface complexity affects H–k results, often leading to estimates of discontinuities that do not reflect any real structure. We, therefore, conclude that H–k stacking on its own may not be appropriate for mapping layered crustal structure in arc settings without additional constraints. The method is limited by the fact that it will only search for one pair of values. In the case of multiple discontinuities, however, it is likely that peaks in the RF caused by different discontinuities will overlie each other and distort the RF to a point where H–k stacking may find values that do not represent any real discontinuity depth and layer v_P/v_S ratio.

3. Inversion for a layered crust

To overcome limitations of the H–k stacking method when applied to layered crust, we adopted a grid-search inversion of a three-layer crust overlying the mantle within the following petrological framework: (i) upper crustal layer composed of loosely consolidated and fractured volcanoclastic sediments and lavas; (ii) middle crustal layer composed of plutonic rocks (solidified magma); (iii) lower crustal layer composed of mafic and ultramafic crustal cumulates; and (iv) mantle layer (figure 3). In a subsequent development of the model, we also consider the presence of vestigial crust of the over-riding, proto-Caribbean plate (layer 2a).

Crustal cumulates are rocks formed by accumulation of near-liquidus phases from magmas undergoing chemical differentiation. Cumulates consist of assemblages that represent instantaneous solid compositions from one or more magma batches. Conversely, plutonic rocks have mineralogy and textures consistent with protracted, in situ solidification of magmatic mushes (melt + mineral phases) without attendant differentiation. Our three-layer crustal model is based on studies of currently active and exposed fossil island arcs (e.g. [22,23,34,35,47,48]). The exposed arcs of Talkeetna and Kohistan show lower crust that consists predominantly of mafic/ultramafic cumulates such as pyroxenite, hornblendites and gabbros, whereas the middle crust is composed of evolved, often felsic, plutonic rocks. The lithologies we see in fossil arcs are comparable to those we find in currently active intra-oceanic arcs. Geochemistry and thermobarometry of lavas and their igneous xenoliths along the LAA support that conclusion (e.g. [26,46,49–55]). The seismic studies of active Aleutian and Izu-Bonin island arcs also reveal layered crustal structures which have been similarly interpreted as mafic lower crust and more felsic middle-upper crust based on v_P and v_S properties of different rock types. We build our model and estimate of v_P/v_S ratios on these observations. The physical properties of the mantle (layer 4) are fixed at: v_P = 8.00 km s⁻¹, v_S = 4.53 km s⁻¹, density = 3.33 g cm⁻³. They are derived from [46] and represent putative mantle values beneath the LAA.

Multiple studies of currently active and fossil arcs demonstrated that, although v_P and v_S vary considerably for the lower-crustal layer, the v_P/v_S ratio is surprisingly constant and on average lies between 1.75 and 1.80 (e.g. [2,21,23,56,57]). A v_P/v_S ratio of 1.79 for the lower-crustal layer (3) is used in our model. The middle crust (2) in our model consists of plutonic rocks of basaltic to andesite compositions. This assumption is based on average lava composition of LAA and upper crustal xenoliths (e.g. [26,50–53]). According to [58], plutonic rocks of this composition will have a v_P/v_S ratio of 1.82 to 1.87 at mid-crustal pressures. The v_P/v_S ratios of the upper layer (1) are controlled primarily by the fracture density and degree of compaction, rather than lithology. In a refinement of our model, we also consider the presence of vestigial Caribbean Plate (layer 2a). The seismic properties of layer (2a) are estimated from rock compositions of the Caribbean oceanic plateau [59]. The v_P/v_S ratio found is 1.86.

The advantage of a fixed v_P/v_S ratio is a significant reduction in parameter space (see the electronic supplementary material for a comparison between modelling with and without petrological constraints). More importantly, v_P/v_S ratios obtained for crustal and mantle lithology allow us to reconcile the seismological and petrological interpretation of crustal structure. To test the reliability of the chosen v_P/v_S ratios, we tested the upper and lower limits of this value for layers (2), (2a) and (3) and found that a modification in the v_P/v_S ratio of ±0.05 does not change the overall discontinuity depth results significantly.

We assume that all melt has either been extracted or is isolated at a very low melt fraction along grain boundaries. Xenoliths from the LAA contain variable, but small, amounts of melt distributed along grain boundaries [49,50,53,55]. We discuss the seismological implication of the presence of melt below.

We set the thickness of layer (1) to be 5 km, consistent with the geophysical data of [34,37,51], and derive its physical properties from the RF data alone. Using the method of [39], we first invert the seismic RF (using −5 s to +15 s after the initial P-peak) for the uppermost 5 km with an initial model that consists of five 1 km thick subsidiary layers (1a–e) and two main crustal layers. We use a smoothness value of 0.1 based on visual observation to create the smoothest models that still fit the observations well. The horizontal slowness was chosen to be 0.06 s km⁻¹ (but different values have been tested for stability of the result) and the singular-value decomposition truncation fraction was chosen to be 0.001 to handle values close to zero. The thickness of each layer stays fixed during this inversion. Thicknesses of the middle (2) and lower (3) crustal layers are varied within the range of plausible values throughout different inversion runs to ensure the stability of the solution for layers (1a–e). This first step accounts for the strong effects of the highly variable structure near the surface on the RFs and can overcome the nonlinearity and non-uniqueness of this problem (e.g. [45]). Because of the nature of layer (1) this inversion does not include any petrological constraints. In the second step, we introduce a grid search to investigate the depth to the MCD and the Moho and the velocity contrast at these discontinuities, thus defining the thickness of the middle and lower crustal layers. Having already fixed the highly variable upper layer (1) in a previous step, it is possible to reduce the grid search to a reasonable number of models and computation time. In this step, we introduce v_P/v_S ratios based on the petrological considerations above. We keep the v_P/v_S ratios for individual layers fixed but allow v_P and v_S to vary. This proves to be a useful constraint, further restricting the explored parameter space. This additional step is needed to vary the thicknesses of the layers, whereas the first step only works with fixed layer thicknesses.

A χ²-misfit is used to evaluate the match of different models with the seismological data and, thus, make them comparable. It is described by

3.1

where N is the number of data points and d(n), σ(n) and s(n) are, respectively, the data RF, its pointwise uncertainties obtained by the jackknife stacking, the model RF at point n and a weighting factor w(n), which is chosen so that it forms an envelope around the maximum at 0 s (P-arrival) and decreases exponentially to both sides (see [45] for further information). The smallest value of χ² depicts the best model; the uncertainties on this model are given by all models that reach 95% of this value, taken from a χ² distribution table.

We carry out the grid search with v_P/v_S values as high as 2.2 to explore the possible presence of melt, which is known to have a much greater effect on v_Sthan on v_P [60]. No low χ² solutions are obtained with such large v_P/v_S values, indicating that large pockets of interconnected melt are unlikely beneath the LAA, consistent with the H–k stacking results, as well as [15] and the textural evidence from xenoliths.

4. Results

The depth to the MCD and v_P for layers (2) and (3) is in excellent agreement with previous work on various segments of the arc [4,12,34,35,61]. The Moho was not observed seismically by Boynton et al. [4], but was estimated to lie at about 35 km based on their gravity data. In [15], using conventional H–k stacking, they propose that the average Moho depth beneath the LAA is 29 ± 7 km.

The best fitting models for stations FDF and SVB can be seen in figure 4. The obtained crustal structure (figure 5 and table 2) shows that the depths to the MCD and to the Moho are highly variable over surprisingly short distances of tens of kilometres. In [15], they arrived at a similar conclusion, with up to 10 km change in Moho depth across Guadeloupe alone. Furthermore, we find that the seismic velocities of layers (2) and (3) also vary laterally. Note that uncertainties in the results do not arise from these lateral variations. Although both discontinuities are present along the entire LAA, beneath some islands they compete to produce either a strong Moho and weak MCD or strong MCD and weak Moho (figure 5). For St. Eustatius and Saba, the Moho is very weak but the MCD is very strong. Beneath Grenada, Martinique, Guadeloupe and St. Kitts, the converse is true. The depth to the MCD varies between 11 and 25 km, while the depth to the Moho varies between 24 and 37 km (e.g. Grenada and St. Kitts, respectively).

Our inversion approach considers only new, magmatic arc crust. We have not considered thus far the possible presence of vestigial proto-Caribbean crust (pCc) within the arc (layer 2a). The estimated seismic properties of this layer (v_P= 7.11 km s⁻¹, v_S = 3.97 km s⁻¹, v_P/v_S = 1.79) are described above. Incorporating layer (2a) into our models does not change the depths of the discontinuities in the inversion results because the change in seismic velocity between the pCc and the adjacent crust is too small to be identified by RFs and, hence, H–kstacking and crustal inversions. This conclusion comes from models that included an 8 km layer (2a), consistent with the seismic refraction study of [34] beneath the Grenada basin where the pCc is appreciably thinned. In the unlikely scenario that layer (2a) is chosen with sufficient thickness that it takes up most of the crustal column, the MCD is suppressed, leading to a significantly different inverted crustal structure. In our models, a 20 km layer (2a) is found to cause such a change. However, the resultant misfit between the model and data RFs is considerably larger in those instances, leading to the reasonable conclusion that a pCc, if present, cannot exceed thicknesses of around 10–15 km, depending on the island, and may not be present at all.