19 August 2012
Oceanic lithosphere is rapidly recycled into the mantle through subduction, an important part of the dynamic evolution of the Earth. Cratonic continental lithosphere, however, can exist for billions of years, moving coherently with the tectonic plates. At the Caribbean–South American Plate margin, a complex subduction system and continental transform fault is adjacent to the South American cratonic keel. Parallel to the transform fault plate boundary, an anomalous region of seismic anisotropy1—created when minerals become aligned during mantle flow—is observed2, 3, 4, 5. This region of anisotropy has been attributed to stirring of the mantle by subducting slabs2, 3. Here we use seismological measurements and global geodynamic models adapted to this unique region to investigate how mantle flow, induced by subduction beneath the Antilles volcanic arc, is influenced by the stiff, deep continental craton. We find that three components—a stiff cratonic keel, a weak asthenospheric layer beneath the oceans and an accurate representation of the subducted slabs globally—are required in the models to match the unusual observed seismic anisotropy in the southeast Caribbean region. We conclude that mantle flow near the plate boundary is deflected and enhanced by the keel of the South American craton, rather than by slab stirring.
The tectonic setting of the southeastern Caribbean involves subduction of oceanic lithosphere in opposite polarity beneath the Antilles arc and the Americas, which are linked by a transform plate boundary between continental South America and the Caribbean Plate (Fig. 1). The Cenozoic tectonic history has been primarily influenced by the Caribbean Plate migrating east relative to the Americas. About 55 million years (Myr) ago, the Caribbean Plate collided with both the Bahamas bank and northern South America6 as it moved eastwards. The Caribbean Plate rotated clockwise as a result of to this collision, which initiated orogenesis in northern South America and development of the San Sebastian–El Pilar right lateral strike-slip system along the margin of northern South America to accommodate the sustained, east–west South American–Caribbean Plate motion7. Just a few hundred kilometres south of this complex plate boundary system lies the Guyana Shield, a ~ 1.7-billion-year-old craton, inferred to have a relatively rigid keel extending into the upper mantle (Fig. 1).
Data from broadband seismic instrumentation from temporary deployments and permanent national networks in Venezuela and in the surrounding region allow for the investigation of lithospheric and upper-mantle structure8, 9, 10 beneath this complex plate boundary (Fig. 1). Rayleigh wave tomography9 reveals a linear shear-wave velocity change that parallels the dextral strike-slip fault system along the northern coast of Venezuela, imaging the differences between the South American continental lithosphere, the Venezuelan archipelago and the Caribbean oceanic lithosphere. The location of the Paria cluster earthquakes (
Waveform-splitting analysis of core shear wave phases (SKS/SKKS) used to study seismic anisotropy1 in the southeastern Caribbean2, 3, 4, 5 has found anomalously large delay times with fast polarization orientations parallel to the transform plate boundary. These 1.5–2.5 s shear-wave splits near the El Pilar–San Sebastian Fault, which marks the boundary between the Caribbean and South American plates (Fig. 1), are much larger in magnitude than the global average of ~ 1.0 s (ref. 13). However, the delay times decrease considerably (
The complex mantle structure and tectonic history of the region have been previously studied, but the interactions of the subduction zones, transform plate boundary and stiff continental keel are still poorly understood. We use high-resolution, global mantle flow models adapted for the seismically and geologically derived mantle and lithospheric structure of the Caribbean and South America, quantitatively predict shear-wave splitting from flow and compare the results with dense seismic observations to evaluate the influence of subducted slabs and cratonic keels on mantle convection and regional plate motions to test some of the dynamical questions of how the unusually strong seismic anisotropy has developed.
We ran 176 global geodynamic models to explore a range of viscosity (for example, weak asthenosphere, stiff keels) and density structures (for example, various slab morphologies, lithospheric structures) to evaluate which geodynamic scenarios are compatible with the observed anisotropy (Supplementary Figs S1 and S2). This approach provides for new quantitative testing of the tectonic relationships between the El Pilar–San Sebastian strike-slip system bounding the South American and Caribbean plates, the influence of the subducting slab beneath the Antilles arc and its possible detachment of the oceanic portion of the South American lithosphere from continental lithosphere. We model mantle circulation by solving the infinite Prandtl number, Stokes equation for incompressible fluid flow in a global, spherical shell14, 15 (see Methods and Supplementary Information). Density anomalies are inferred by scaling velocity anomalies from seismic tomography to temperature and all parameters, other than viscosity, are the same as those in ref. 16. Our reference model uses density anomalies derived from the SMEAN composite tomography model17, but we also replace upper-mantle (cold/fast velocity) structure with subducting slabs (globally) inferred from Wadati–Benioff zone seismicity18 and other seismicity-defined slab geometries (Supplementary Fig. S2). By computing global geodynamic models of mantle flow and inferred lattice preferred orientation (LPO) of olivine we are able to test and provide an explanation for the observed seismic anisotropy in the southeastern Caribbean. Although this modelling step from flow to LPO involves several uncertainties19, 20, 21, 22, we expect that the effect of density and viscosity variations will be dominant for the models considered here. Model performance is evaluated by comparing the real and synthetic splitting measures, for both station-averaged and back-azimuth specific splits, by computing the mean angular deviation of apparent fast polarization orientations (fast azimuths) Δα (0≤Δα≤90°) and mean delay-time misfits Δ(δt) = δtmodel−δtdata.
Although we computed 176 global models, which are described in the Supplementary Information, we found three primary factors contribute to the best fit to the seismic anisotropy data. The reference model, based on SMEAN (ref. 17), predicts the splitting orientations well at some stations near the plate boundary and in northwestern Venezuela, Mérida Andes and the Maracaibo block, but underpredicts delay-time magnitudes slightly (Fig. 2a). However, overall this model, without any adjustment for regional lateral viscosity variations, does not predict the decrease in magnitude or the orientation of the splitting measurements away from the plate boundary and into the interior of the South American continent or onto the Caribbean Plate. The second model (Fig. 2b) replaces the upper-mantle structure from tomography with slabs inferred from Benioff zone seismicity18 because such models have been shown to lead to improvements in plate motion fits23, 24. The third model (Fig. 2c) adds a weak, suboceanic asthenospheric layer to the global model (100–300 km depth, Supplementary Fig. S1). This model overpredicts the delay times for the shear-wave splits and degrades the fit for the stations in northeastern Venezuela (Fig. 2c), however it is an improvement from the SMEAN and slab model. Our preferred model is shown in Fig. 2d, which includes 300-km-deep keels based on cratonic geometry from Nataf and Ricard25 to the global model with an asthenospheric layer at 100–300 km depths underneath all the oceans. Supplementary Fig. S2 illustrates the performance of all 176 models and the range of viscosity and density structures tested in terms of global velocity correlation, angular fast-orientation misfit and delay-time misfit. A range of robust misfit systematics exists for the various models. For example, the combination of cratonic keels and a weak asthenosphere beneath the oceanic plates always improves the match between the predicted and observed shear-wave splitting, as does implementation of lateral viscosity variations, compared with models without such effects. Specifically the addition of the cratonic keels globally, including the Guyana Shield to 300 km (Fig. 1), produces a stirring effect in the global model that overall markedly improves both the fast azimuth and delay-time misfit. The magnitude and orientation of the splitting predictions match the east–west-oriented measurements along both the plate boundary, the northeast–southwest-trending splits in the Maracaibo block and, more importantly, the smaller delay times that are measured in the continental interior and the craton, which none of the other ~ 160 models were able to do (see Supplementary Figs S2–S4).
Investigation of the calculated flow patterns at asthenospheric depths (165 km) illustrates the changes due to the lateral viscosity variations imposed in each of the four models described in Fig. 2. Adding the inferred slab structure produces only a very minor increase in velocity along the northern portion of the Antilles subduction zone and a slight decrease in velocities beneath the Caribbean Plate (Fig. 3a,b). However, the weak asthenosphere markedly increases the downward velocity of the subducting slab and flow around the edges of the plate. Yet, neither the orientation of flow along the South America–Caribbean Plate boundary nor beneath continental South America changes. With the addition of the stiff, deep continental keels into the global model the flow is deflected around and below the craton instead of into the interior of the continent (Fig. 3c,d). This deflection of flow around the craton improves the predicted splits beneath the cratonic continental lithosphere and improves the fit of the splits in northeastern Venezuela, while preserving the match to the observations in northwestern Venezuela. At greater mantle depths the flow returns to a generally north–south orientation. Although the Caribbean Plate motion is somewhat overpredicted, the surface-velocity orientations broadly match estimates based on geodetic data, providing insight into internal plate deformation and how plate velocities are affected by the stiff continental keel (Supplementary Figs S4 and S5).
When adapted to the expected lateral viscosity variations of the Caribbean and South America slab/craton system, our global models newly indicate a possible link to channelling of flow by the deep South American craton, where mantle is escaping around the southern end of the Antilles slab facilitated by a weak mantle wedge. Although our models are based on instantaneous-flow calculations, the anisotropic fabric is thought to have evolved into its preferred orientation over ~ 10 Myr (see Supplementary Information). Therefore, we can assume that the present east–west velocity of the Caribbean Plate relative to stationary South America7 has remained mostly constant during this timeframe and been a major factor in the development of the strong plate boundary parallel anisotropy. This supports the notion that the oceanic portion of the lithosphere is tearing away from continental South America. Intriguingly, it is not slab stirring3 but the deep cratonic keel associated with the Guyana Shield that deflects and enhances mantle flow into a wide channel near the transform plate boundary from around the southern edge of the subducting slab. The keel acts as an anchor that decreases the overall surface velocity of South America, leads to increased internal deformation in the plate and contributes to the strong seismic anisotropy observed in the southeastern Caribbean.
We model mantle circulation by solving the infinite Prandtl number, Stokes equation for incompressible fluid flow (Boussinesq approximation) in a global, spherical shell26, 27. We use a variant of the finite element code CitcomS14, 15 to arrive at steady-state solutions for density-driven flow using an approach similar to ref. 15 by prescribing plate boundaries as weak zones and solving for the resulting surface plate motions self-consistently. Our numerical resolution of ~ 25 km laterally and vertically in the upper mantle allows us to resolve regional mantle flow at the required level of detail in a global model16. Density anomalies in the mantle are inferred by scaling the isotropic anomalies of a shear-wave velocity tomography model by a constant conversion factor dlnρ/dlnvS = 0.2, or by assigning constant density anomalies to upper-mantle slabs as inferred from the Wadati–Benioff zone geometry. All parameters, other than viscosity, are the same and discussed in ref. 16.
Seismic anisotropy, such as that measured most directly by shear-wave (SKS/SKKS) splitting, can be used to infer flow in the upper mantle because mantle rocks record deformation by means of the formation of LPO of olivine13, 28. By computing geodynamic models of mantle flow and inferred LPO we can test possible explanations of the observed seismic anisotropy patterns. Both laboratory and field evidence show that the crystallographic axes of olivine align with shear and LPO development can be estimated quantitatively20, 28. Synthetic LPO calculated for mantle flow models using kinematic LPO texturing theory19 was benchmarked in ref. 21 to show that the synthetic LPO heterogeneity matches that observed in mantle xenoliths, implying that mineral physics methods and lab experiments may be valid when used for global applications. This is corroborated by a number of studies that have previously explored the link between flow and SKS splitting in simpler tectonic settings22, 29.
In our models, LPO is computed along streamlines as inferred from the instantaneous mantle flow solutions using the DREX method with parameters as discussed in ref. 19. We follow tracers in the mantle on forward advection until a logarithmic, finite strain of 0.75 is reached (typically within ~ 10 Myr; ref. 21), assuming no feedback between seismic and mechanical anisotropy. Then from the LPO, seismic anisotropy is estimated by Voigt averaging single-crystal elastic tensors taking into account the depth-dependence of moduli21. Synthetic splitting is computed using a cross-correlation method from full waveforms that incorporate full anisotropy along the path and finite frequency effects22, assuming an average 5 ° incidence angle. Model performance is then evaluated by comparison of the real and synthetic splitting measures by computing the absolute, angular deviation of fast apparent polarization direction (fast azimuths), (0≤Δα≤90°), where Δα = 45°equal to the random mean and the delay-time misfit (Δ(δt) = δtmodel−δtdata), where Δδt>0 indicates overprediction of anisotropy. Model performance is evaluated by comparing station-averaged splits with averages computed from full back-azimuth scans of synthetic waveforms, as well as split by split for the actual event back-azimuths if this information is available.