Research ArticleGEOPHYSICS

Imaging paleoslabs in the D″ layer beneath Central America and the Caribbean using seismic waveform inversion

See allHide authors and affiliations

Science Advances  29 Nov 2017:
Vol. 3, no. 11, e1602700
DOI: 10.1126/sciadv.1602700

Abstract

D″ (Dee double prime), the lowermost layer of the Earth’s mantle, is the thermal boundary layer (TBL) of mantle convection immediately above the Earth’s liquid outer core. As the origin of upwelling of hot material and the destination of paleoslabs (downwelling cold slab remnants), D″ plays a major role in the Earth’s evolution. D″ beneath Central America and the Caribbean is of particular geodynamical interest, because the paleo- and present Pacific plates have been subducting beneath the western margin of Pangaea since ~250 million years ago, which implies that paleoslabs could have reached the lowermost mantle. We conduct waveform inversion using a data set of ~7700 transverse component records to infer the detailed three-dimensional S-velocity structure in the lowermost 400 km of the mantle in the study region so that we can investigate how cold paleoslabs interact with the hot TBL above the core-mantle boundary (CMB). We can obtain high-resolution images because the lowermost mantle here is densely sampled by seismic waves due to the full deployment of the USArray broadband seismic stations during 2004–2015. We find two distinct strong high-velocity anomalies, which we interpret as paleoslabs, just above the CMB beneath Central America and Venezuela, respectively, surrounded by low-velocity regions. Strong low-velocity anomalies concentrated in the lowermost 100 km of the mantle suggest the existence of chemically distinct denser material connected to low-velocity anomalies in the lower mantle inferred by previous studies, suggesting that plate tectonics on the Earth’s surface might control the modality of convection in the lower mantle.

INTRODUCTION

The purpose of this study is to obtain high-resolution three-dimensional (3D) images of the S-velocity structure in the D″ (Dee double prime) layer beneath Central America and the Caribbean to search for evidence of paleoslabs above the core-mantle boundary (CMB) and for evidence of small-scale low-velocity anomalies that might suggest chemical heterogeneity (1).

Recently, dense seismic arrays such as the USArray, which includes many portable stations that have steadily been moved eastward to cover the entire conterminous area of the United States, are providing excellent data for high-resolution imaging of localized regions of D″ using waveform inversion. Our group has recently conducted two small-scale feasibility tests of this method to invert for the 3D S-velocity structure in the D″ layer beneath Central America (2) and the western Pacific (3), and applied the same method to a much larger data set to invert for the 3D S-velocity structure in D″ beneath the northern Pacific (4). Here, we use the full USArray to obtain dense coverage of the D″ layer beneath Central America and the Caribbean. The use of short-period (up to 8 s) waveforms makes it possible to image small-scale structure with finer resolution than travel-time tomography or global waveform inversion studies.

The D″ layer at the base of the mantle is, after the Earth’s crust and uppermost mantle, the second most seismically laterally heterogeneous region of the Earth’s mantle (5, 6). Strong large low-velocity provinces (LLVPs) beneath Africa and the South Pacific, and high-velocity regions beneath the circum-Pacific, are the large-scale features found ubiquitously by travel-time tomography or global waveform inversion studies (5, 7).

Low–seismic velocity regions in the lowermost mantle can be explained by high temperatures, chemically distinct material, or a combination of the two. Pyrolite is widely thought to be the average composition of the lower mantle (8, 9), but the details of the bulk composition of the lower mantle remain controversial (10, 11). Chemical compositions with increased amounts of impurities, such as Fe and Al, have lower shear velocities than pyrolite (12, 13). This chemical heterogeneity, resulting from exchanges between the core and the mantle, from partial melting in the thermal boundary layer (TBL), from basalt entrained to the base of the mantle by past subduction, or as long-lived remnants of chemical differentiation in the early Earth, is expected at the CMB (1). However, the extent to which chemical anomalies contribute to the lowermost mantle seismic structure is still unclear; the LLVPs, which could possibly be large-scale chemically distinct regions, have been variously suggested to be due to temperature anomalies only (3, 1416) or to be chemically distinct from the rest of the lower mantle (17, 18). High-velocity anomalies in the lowermost mantle can be explained by the combined effect of low temperatures and the bridgmanite (abbreviated as MgPv below) to Mg-postperovskite (abbreviated as MgPPv below) phase transition (19, 20).

It is generally thought that high-velocity anomalies inferred from seismic tomography in the upper mantle and in the upper part of the lower mantle correspond to remnants of past subduction (21, 22). Seismic tomography shows high-velocity anomalies continuous from the upper mantle down to the 660-km discontinuity, where the negative Clapeyron slope of spinel decomposition and a jump in viscosity make some slabs stagnate, whereas others can be seen to descend to the mid-mantle (~1800 km depth). In the lowermost ~500 km of the mantle, strong and broad high-velocity anomalies are seen beneath slabs that extend to the mid-mantle (5, 7, 23).

Evidence supporting the existence of paleoslabs at the CMB includes reports of continuous high-velocity anomalies from the transition zone down to the CMB beneath North America (Farallon slab) and beneath the southern Indian Ocean (24) and the fact that high-velocity anomalies in the lowermost mantle are generally consistent with the location of ancient slabs (25). Estimates of the global average sinking rate of slabs in the lower mantle based on seismic tomography vary from 1.1 to 1.9 cm/year (26, 27). Subducted material older than 160 to 260 million years ago (Ma) might thus have reached the CMB. Slabs at the CMB are strongly heated and probably disintegrate in ~100 million years (My) (28); thus, material older than ~260 to 360 Ma is unlikely to cause high-velocity anomalies at the CMB.

An increase in the reported amplitudes of the high-velocity anomalies in cold regions in the lowermost mantle, and thus in the visibility of hypothetical paleoslabs, seems most likely to be due to the MgPv to MgPPv phase transition (19, 20). The detailed structure of broad high-velocity anomalies in the lowermost ~500 km of the mantle cannot be resolved in current tomographic models because of their coarse parametrization in the lowermost mantle. Thus, it was heretofore not possible to say whether high-velocity anomalies near the CMB result from spreading of accumulated slabs or from separate paleoslabs that followed different paths to the lower mantle, possibly originating at different subduction zones.

Study area

The western margin of North and South America (previously the western margin of the supercontinent Pangaea) is a region of long-lived subduction, where the oceanic Farallon plate is believed to have initiated eastward subduction around 180 to 207 Ma (26); tectonic studies suggest subduction of the Pacific Ocean beneath the west coast of Pangaea since ~250 Ma (29). The Farallon plate is seen in tomographic models as an eastward-dipping high-velocity feature reaching depths of ~2000 to 2500 km (21, 26, 30), and its location at those depths is in general agreement with the reconstructed plate boundary ~180 Ma (31).

Continuity of the Farallon plate (as an eastward-dipping high-velocity anomaly) from depths of ~2000 to 2500 km to the CMB beneath Central America is not a strong feature of previous models. However, previous models show a broad high-velocity anomaly beneath and westward of the location of the Farallon slab at a depth of ~2000 km. Tearing and breaking of the slab at the 660-km discontinuity or around a depth of ~1000 km could explain the apparent discontinuity of the Farallon slab (22). This may be supported by geological evidence for voluminous igneous activity ~200 Ma (29). Thus, high-velocity anomalies deeper than ~2500 km might correspond to subduction older than 200 Ma at the western margin of Pangaea. In addition to subduction at the western margin of Pangaea, tectonic studies suggest that there were subduction zones within the Pacific Ocean, where the ocean floor subducted beneath active volcanic arcs (32). This intra-oceanic subduction (33) has been suggested as a possible explanation for the high-velocity seismic structure in the deep mantle beneath the western margin of North America (34) and in the lowermost mantle beneath Central America (26).

RESULTS

We conduct waveform inversion for the detailed 3D S-velocity structure of the lowermost 400 km of the mantle beneath Central America and the Caribbean using ~7700 transverse component records cut 20 s before the arrival of the direct S wave and 60 s after the arrival of the ScS phase (see Materials and Methods). The events are deep- and intermediate-focus events recorded at epicentral distances 70° < Δ < 100° at broadband seismic stations of the USArray, Canadian Northwest Experiment (CANOE), Global Seismographic Network (GSN-IRIS/USGS), Southern California Seismic Network (SCSN), Pacific Northwest Seismic Network (PNSN), Berkeley Digital Seismic Network (BDSN), and Canadian National Seismograph Network (CNSN) (Fig. 1). The azimuthal- and epicentral-distance distribution of the stations is shown in fig. S1. The data are filtered in the period range of 8 to 200 s using a Butterworth bandpass filter. The 3D model is obtained by linearized inversion with respect to a spherically symmetric initial model.

Fig. 1 Target region.

Waveforms from deep- and intermediate-focus earthquakes beneath South America (red stars) recorded at stations of the USArray, CNSN, CANOE, and other seismic networks (see text) (blue inverted triangles) provide dense raypath coverage of the target region 0 to 400 km above the CMB (yellow squares). Red curves show ScS raypaths that sample the target region, and black crosses show ScS bounce points at the CMB. The pink solid circle at 30°N and 110°W shows the location for the shallow structure trade-off test (fig. S3). The inset shows the location of the cross sections presented in Fig. 4 and the location of the Farallon plate boundary at 180 Ma (31).

We first conduct 1D waveform inversion (35, 36), with respect to the global reference model PREM (37), to infer a regional 1D model (hereafter called PREM′) of the lowermost mantle beneath Central America (fig. S2). This model was derived by optimizing the fit of the synthetics to our full data set over a wide range of azimuths and thus differs somewhat from previously published 1D models of the lowermost mantle beneath Central America that used a significantly smaller number of waveforms from a narrower range of azimuths (35, 38).

The 3D model is parametrized with 744 3D cells (voxels) of dimensions 5° × 5° × 50 km (equivalent to approximately 300 × 300 × 50 km at the CMB) in eight horizontal layers of 50 km thickness from the CMB to 400 km above the CMB (Fig. 1). To study the dependence of the 3D models on the initial 1D model, we conducted two inversions, with respect to PREM and PREM′, respectively. We call the resulting 3D models CACAR (an abbreviation of Central America–Caribbean) and CACAR′, respectively. Map views of these 3D models (with the average perturbation in each layer set to zero) are shown in Figs. 2 and 3, respectively. The two 3D models are in good general agreement.

Fig. 2 Model CACAR.

The eight panels show the results of the inversion for the eight depth layers from 400 km above the CMB (upper left) to the CMB (lower right), with the lateral average of the 3D perturbation set to zero in each layer. The perturbation is relative to the initial 1D model PREM (37). Two distinct high-velocity regions at the CMB (lower right) suggest two distinct cold paleoslabs. A 3% velocity decrease beneath Mexico concentrated within 100 km of the CMB suggests the possible existence of chemically distinct material with enriched iron content (for example, basaltic composition). The location of high- and low-velocity anomalies is generally consistent with recently inferred topography of the D″ discontinuity (42). The Farallon plate boundary at 180 Ma (31) is shown in red in the lower right panel. Its location is consistent with the high-velocity anomaly beneath Venezuela but is ~1000 km away from the high-velocity anomaly beneath Central America. This might indicate past intra-oceanic subduction, or breaking or tearing of an ancient paleoslab in the upper mantle (22).

Fig. 3 Model CACAR′.

Same as Fig. 2 but using PREM′ (fig. S2) as the initial 1D model. The velocity perturbations are shown with respect to PREM′, with the lateral average of the 3D perturbation set to zero in each layer. CACAR (Fig. 2) and CACAR′ are in good general agreement, which suggests that the inversion is robust. The strong low-velocity anomaly beneath Mexico observed in CACAR is still present, but the strongest low-velocity anomaly is beneath Ecuador. This suggests that strong low-velocity anomalies are present at several sites around the high-velocity anomalies.

The key features of our models are the following: (i) two distinct high-velocity anomalies just above the CMB, one beneath Central America and another beneath Venezuela; (ii) strong low-velocity anomalies concentrated in the lowermost 100 km of the mantle at the edge of the high-velocity anomalies; and (iii) vertically continuous low-velocity structures above the strong low-velocity anomalies at the CMB. The past Farallon plate boundary at 180 Ma (31), indicated with a red line in the lowermost right panel (0 to 50 km above the CMB) in Figs. 2 and 3, is consistent with the location of the high-velocity anomaly beneath Venezuela, whereas it is ~1000 km from the high-velocity anomaly beneath Central America. In the depth range of 0 to 50 km above the CMB, the high-velocity anomaly beneath Venezuela in CACAR does not extend to the north along the past Farallon plate boundary (Fig. 2). The same high-velocity anomaly beneath Venezuela in CACAR′ extends up to a few degrees north of the Caribbean islands, but a strong high-velocity anomaly is not present to the north of the Caribbean islands (Fig. 3). However, as the altitude above the CMB increases, the high-velocity anomaly beneath Venezuela gradually extends along the past Farallon plate boundary to the north of the Caribbean islands and is strongest 250 to 300 km above the CMB (Figs. 2 and 3).

Cross sections and comparison to previous studies

Several works in the previous decade considered this study area using finite-frequency travel-time tomography with a model parametrization with a scale similar to ours (39) or migration of phases refracted by the D″ discontinuity (40, 41); a more recent study conducted forward modeling for the topography of the D″ discontinuity beneath Central America and the Caribbean (42). Our results can also be compared to the structure in this region obtained by global waveform inversion (7). From the above studies, we select the finite-frequency travel-time tomography study (39) and the global waveform inversion study (7) for a detailed comparison to our models, as shown in Fig. 4.

Fig. 4 Cross sections.

Left column: Cross sections through two previous models obtained by global waveform inversion (labeled FR2014) (7) and by regional finite-frequency travel-time tomography (labeled H+2005) (39). Middle column: Cross sections for our models CACAR and CACAR′. Right column: Whole-mantle cross section through the global model (7). The horizontal axis of each profile shows degrees along the corresponding great-circle cross section; vertical axis shows elevation above the CMB (in kilometers) for profiles in left and middle columns and depth (in kilometers) for the whole-mantle profiles in the right column. Locations of the cross sections are shown in the inset in Fig. 1. The leftmost part of the cross section for model H+2005 in the left column of (B) has been grayed out due to the lack of resolution. (Note that CACAR and CACAR′ also have little resolution in this region.) (A) The strong velocity contrasts within 100 km above the CMB in CACAR and CACAR′ (blue dashed line in the upper middle panel) suggest the presence of paleoslabs at the CMB and dense chemical heterogeneities. The paleoslab beneath Central America (labeled “CA”) is perched above a strong low-velocity anomaly (blue dashed line in the upper middle panel), which might suggest dense iron-enriched material at the CMB (for example, basaltic composition). Low-velocity vertically continuous structures (brown arrows with unfilled points in cross-section A-A′ in the middle column) at the edges of “CA” suggest upwelling from the possibly iron-enriched material at the CMB. (B) A low-velocity vertically continuous structure (red filled arrow in cross-section B-B′ in middle column) between two distinct paleoslabs beneath Central America (“CA”) and Venezuela (“VZ”), and connecting to a low-velocity region in the mid-mantle in a previous global waveform inversion model (7) (right column) suggests upsplashing of hot TBL material caused by two paleoslabs “CA” and “VZ” sinking to the CMB. Past location of the Farallon plate boundary at ~180 Ma (31) (green vertical dashed line labeled FA) is consistent with location of “VZ”, but is laterally ~1000 km away from “CA”.

Cross sections of the various models along two profiles (see inset of Fig. 1 for locations) are shown in Fig. 4. Our models CACAR and CACAR′ are shown in the middle column of Fig. 4. The left column of Fig. 4 shows the models along these cross sections obtained by the recent global waveform inversion (7), labeled FR2014, and the models from the finite-frequency travel-time tomography study (39), labeled H+2005. Whole-mantle cross sections through the global waveform inversion model (7) are shown in the right column of Fig. 4.

Figure 4A shows a northwest (A) to southeast (A′) great-circle cross section through the high-velocity anomaly beneath Central America. As shown in the inset to Fig. 1, this cross section is roughly parallel to the reconstructed Farallon plate boundary at 180 Ma (31). Our models CACAR and CACAR′ are shown along this cross section in the middle column of Fig. 4A. The two models differ slightly, but we find the following consistent features in both models. (i) Low-velocity anomalies are seen at both edges of the high-velocity anomaly (beneath Central America, labeled “CA” in the Fig. 4), marked by the brown arrows with unfilled points beneath the second panel in the middle column. (ii) The low-velocity anomalies are particularly strong in the lowermost 100 km of the mantle, where ~3% low-velocity anomalies are present. (iii) The northern part of the high-velocity anomaly label ed “CA” in the figure is perched above a 100-km-thick low-velocity anomaly just above the CMB, resulting in a strong velocity contrast (5 to 6% peak to peak over a small vertical range at the dashed blue line in the upper two panels in the middle column of Fig. 4A).

Figure 4B shows a west (B) to east (B′) cross section at 5°N, roughly perpendicular to the past Farallon plate boundary at 160 to 180 Ma (31). For this particular cross section, CACAR′ shows anomalies with smaller amplitudes than CACAR, but the pattern of high- and low-velocity anomalies for the two models is consistent. We see two high-velocity structures that correspond to the high-velocity anomalies beneath Central America (labeled “CA”) and Venezuela (labeled “VZ”), respectively. “CA” and “VZ” are separated by a vertically continuous low-velocity region (indicated by the red arrow).

We now compare our models to the recent global waveform inversion (7), labeled FR2014, and the finite-frequency travel-time tomography study (39), labeled H+2005. The high-velocity anomaly “CA” is also present in the two previous models (left column of Fig. 4A). The strong low-velocity regions on both edges of “CA” (brown arrows with unfilled points beneath the second panel in the middle column of Fig. 4A) are present in the previous models (left column of Fig. 4A) but are smaller in size and weaker by ~1.5%. As a result, the fact that near “CA” the high-velocity anomaly is perched ~100 km above a strong low-velocity zone at the CMB (dashed blue line in the upper two panels in the middle column of Fig. 4A) is not evident from the previous models. We also note in passing that the high-velocity anomaly near “CA” in Fig. 4A is generally consistent with previous regional forward modeling studies (4042). The low-velocity region to the right (south) of “CA” in the middle column of Fig. 4A (marked by the brown arrow with unfilled points) connects to a larger low-velocity region that was found to extend to the mid-mantle by a recent global waveform inversion study (right column of Fig. 4A) (7). The high-velocity anomaly “VZ” in the middle column of Fig. 4B is covered by low-velocity material at shallower depths in the global model (right column of Fig. 4B). As noted above, the location of the high-velocity anomaly beneath Venezuela agrees with the past Farallon plate boundary at 180 Ma (indicated by the vertical green dashed line labeled FA), whereas the high-velocity zone anomaly beneath Central America is ~1000 km distant from it.

The high-velocity anomalies “CA” and “VZ” are also visible in the previous regional and global models (left column of Fig. 4B). However, the amplitude of “VZ” in the previous models is somewhat weaker than that of “CA”; the improvement in the resolution of the high-velocity anomaly “CA” in our model is most likely due to a greater number of records than the previous regional study (39), made possible by the new data from the dense, transportable USArray. In addition, waveform inversion allows us to also use the phase and amplitude information in the data, which most likely explains why our model shows more consistent vertical features. Part of the discrepancy between our model and the previous regional study might also be due to the latter’s use of core phases (SKS), which can be affected by anisotropy and uncertainty in the outer-core velocity structure and thus might be affected by artifacts.

The location of high- and low-velocity anomalies in our model generally agrees with variations in the elevation of the D″ discontinuity reported on the basis of forward modeling of the recent USArray data (42). In particular, the forward modeling study reports (i) an elevated D″ discontinuity beneath Venezuela, where we observe the high-velocity anomaly “VZ”, (ii) a deeper D″ discontinuity in a corridor along Ecuador and Columbia that extends in the Caribbean Sea, where we see the low-velocity anomaly corridor separating “VZ” from “CA”, and (iii) the highest elevation of the D″ discontinuity above the CMB, where we observe the high-velocity anomaly beneath “CA”. We note that a previous study also observed this correlation (43), but our model has finer resolution.

DISCUSSION

We used waveform inversion to image the complex, small-scale D″ structure with two distinct strong high-velocity anomalies at the CMB shown in Figs. 2 to 4. Our results suggest the presence of two paleoslabs just above the CMB and dense chemical anomalies (that is, iron-enriched material) concentrated in the lowermost 100 km of the mantle. Figure 5 shows a schematic illustration of this interpretation.

Fig. 5 Possible geodynamical interpretations.

Two high-velocity anomalies (“CA” and “VZ”, see Fig. 4) just above the CMB, beneath Central America and Venezuela, respectively, suggest the existence of two distinct paleoslabs that took different subduction paths to the lowermost mantle. Strong low-velocity anomalies in the lowermost 100 km of the mantle suggest dense iron-rich material (for example, basaltic composition). This dense low-velocity material at the CMB can also explain why paleoslab “CA” is perched above the strong low-velocity anomaly beneath Mexico. The low-velocity material is upwelling from the iron-rich anomalies and between the two slabs. The Farallon plate boundary at 180 Ma (red curve) is consistent with the location of the high-velocity anomaly beneath Venezuela, suggesting that the paleoslab “VZ” subducted at the western margin of Pangaea.

As mentioned in Introduction, geological evidence for subduction ~250 Ma beneath the western margin of Pangaea, together with an average subduction rate of ~1.5 cm/year in the lower mantle, suggests that remnants of past subduction should be found at the CMB beneath Central America and the Caribbean. The following features of our models suggest the presence of paleoslabs: (i) ~3% high-velocity anomalies just above the CMB, (ii) vertically continuous low-velocity structures at the edge of the inferred paleoslabs, (iii) the good agreement between the location of the past Farallon plate boundary and the high-velocity anomaly beneath Venezuela, and (iv) the correlation between the topography of the D″ discontinuity reported by previous forward modeling studies and the distribution of high- and low-velocity anomalies in our models.

To explain the 3% high-velocity anomalies in our inversion results, we note that a ~1.5% high-velocity anomaly can be explained by the MgPv to MgPPv phase transition, and the remaining 1.5% by a 390 K decrease in temperature using the temperature derivative for MgPPv under lowermost mantle conditions (44). Because the temperature of 3800 K should be homogeneous at the CMB (16), a 390 K decrease in temperature ~25 to 50 km, or less, above the CMB strongly suggests the presence of cold material. On the other hand, as noted in Introduction, it is difficult to explain a 1.5% velocity increase as the result of differences in chemical composition alone (particularly because increases in both Fe and Al content would decrease the shear velocity alone).

The presence of a low-velocity vertically continuous structure (red arrow in Fig. 4B) that separates the two high-velocity anomalies “CA” and “VZ” and connects to a low-velocity region that was found to extend to the mid-mantle by a recent global waveform inversion study (7) is consistent with the subduction of two distinct paleoslabs. This continuous low-velocity structure from the CMB to the mid-mantle suggests that the hot TBL material has been upsplashed by paleoslabs reaching the CMB, as shown in geodynamical simulations (28), and upwelled to the mid-mantle.

Finally, the positive correlation on small scales between the distribution of high- and low-velocity anomalies in our models and the lateral variations in the elevation of the D″ discontinuity reported by previous regional reflectivity and forward modeling studies (see Results) strongly suggest (assuming the MgPv to MgPPv phase transition) that the high-velocity anomalies “CA” and “VZ” are regions with lower-than-average temperatures and that the low-velocity corridor that separates “CA” and “VZ” is a region with higher-than-average temperatures. The resulting strong temperature gradients might be difficult to sustain without the presence of colder, and thus stiffer, slab remnants.

Our results thus suggest the presence of a paleoslab at the CMB beneath Venezuela, which is in good agreement with the past location of the Farallon plate and might thus be a remnant of subduction at the western margin of Pangaea. Beneath Venezuela, this paleoslab is covered by a low-velocity region ~200 to 400 km above the CMB, but it gradually extends to the north along the Farallon plate boundary as the altitude above the CMB increases (Figs. 2 and 3).

The high-velocity anomaly that we observe beneath Central America has often been imaged by seismic studies and has been interpreted as folding or spreading of the Farallon plate at the CMB (40, 41, 45). However, as noted above, the location of this paleoslab at the CMB is ~1000 km to the west of the past Farallon plate boundary, which tends to argue in favor of two separate subduction paths in the lower mantle rather than the spreading of the Farallon paleoslab at the CMB.

We present below three possible interpretations based on our models CACAR and CACAR′. We note that the velocity model for the structure above our target region should be further investigated to verify these interpretations.

As discussed in Introduction, the initiation of the subduction of the Farallon plate is estimated at 180 to 207 Ma (26) and was probably accompanied by strong igneous activity (29). The two distinct paleoslabs we observe at the CMB might possibly suggest that the strong igneous activity was related to tearing or breaking of a plate subducting beneath western Pangaea, leaving two paleoslabs that sunk to the CMB, followed by the subduction of the Farallon plate itself.

Alternatively, the intra-oceanic subduction of Pacific oceanic floor beneath a volcanic arc within the Pacific located around the current location of Central America, possibly the Stikinia-Quesnellia arc (32), could also explain the presence of two distinct slabs at the CMB and the ~1000-km discrepancy between the location of the Farallon plate boundary and the location of the paleoslab we observe beneath Central America.

We note that if the convection in the lower mantle is isolated from the convection in the upper mantle (46), the high-velocity anomalies “CA” and “VZ” (Fig. 4) that we observe at the CMB might be colder-than-average material related to downflow in the isolated lower mantle.

Chemically distinct layer

The strong velocity contrast (5 to 6% peak-to-peak over less than 100 km vertically and 300 km laterally) that we observe in the lowermost 100 km of the mantle seems too strong and sharp to be explained by temperature variations only; a 5 to 6% velocity anomaly would require a 1300 to 1560 K temperature variation (44). Chemical heterogeneities with enriched iron content (that is, basaltic composition) have a lower seismic velocity than pyrolite and can explain strong negative anomalies (12). They are also denser than pyrolite and thus could remain close to the CMB. Our model shows that strong low-velocity anomalies are concentrated in the lowermost 100 km of the mantle. Chemical anomalies at the CMB provide the most reasonable explanation for the strong velocity gradient we observe within 100 km of the CMB.

Furthermore, geodynamical simulations of a slab sinking to the CMB show that in the case where a thin basaltic layer is present just above the CMB, the center of the slab reaches the CMB, but its edges are perched above a dense basaltic layer (28). This is consistent with our results (Figs. 4A and 5) that suggest that the northwestern edge of the paleoslab “CA” is perched above a strong low-velocity anomaly.

The dynamical stability of a paleoslab perched above a low-velocity anomaly will depend on the density and viscosity contrast between the possible chemically distinct material at the CMB and the rest of the mantle. Some geodynamical studies suggest that this material beneath subduction zones will always be entrained by mantle convection and form large thermochemical piles imaged as large low–shear velocity provinces by seismic tomography (47). Other studies suggest that a denser (4.3% denser than harzburgite) primordial layer at the CMB cannot be easily entrained by mantle convection (48). The fact that in our model CA is partially perched above a possibly chemically distinct low-velocity anomaly seems to suggest that the density contrast between the ambient mantle material and the chemical anomaly is relatively strong and that that a chemical anomaly may thus not be easily entrained by mantle convection.

CONCLUSION

Our model shows that low-velocity anomalies adjacent to and below high-velocity anomalies interpreted as cold subducted paleoslabs seem to be connected to low-velocity anomalies inferred by previous tomographic studies (7). The low-velocity anomalies beneath the Caribbean in the lower mantle can be interpreted as originating at the low-velocity anomalies just above the CMB found in our model. These low-velocity anomalies have possibly been growing below the cold paleoslabs due to increased heat transfer from the core, as suggested by some geodynamical studies (28). They can be interpreted as chemical rather than thermal anomalies because of dense primordial material or iron-rich materials chemically concentrated from pyrolytic or basaltic materials due to the heat from the core or self-generated heat in basalt. Once these concentrations of dense materials are created immediately above the CMB, they are pinned there until stirred or entrained by mantle convection when they become the origin of upwelling flow in the lower mantle. Considering these pinning effects, significant upwelling flows in the lower mantle found beneath hotspots by previous studies (49) might be related to past subduction history.

The significant upwellings beneath the Pacific and Africa are located beneath past supercontinents Rodinia and Gondwana, respectively (50). Because significant downwelling flow is expected in the lowermost mantle beneath the location of accumulated subduction zones such as the present East Asian zone and past supercontinental margins, self-generated heat in subducted basaltic material is likely to be the origin of the significant upwelling (51). Our seismic velocity model is consistent with the possibility that subducted material accumulated at the TBL produces concentrations of iron-rich material due to chemical differentiation, which becomes the origin of upwelling flow from this iron-rich material pinned at the CMB. Significant downwelling from the Earth’s surface due to subduction could be responsible for transporting large amounts of material enriched in radiogenic elements, such as basalt, that produce iron-rich materials according to chemical differentiation due to self-generated heat, which are then pinned at the CMB, thus providing the origin of the hotspots. This suggests that the modality of convection in the lower mantle is controlled by plate tectonics on the Earth’s surface.

As the temperature at the CMB is isothermal because of the vigorous convection of the outer core (52), knowledge of the temperature gradient just above the CMB is essential for understanding the thermal evolution of the Earth. Because the heat flux is the product of the temperature gradient and the thermal conductivity, intermittent paleoslab subduction could make the cooling rate at the surface of the outer core both spatially and temporally heterogeneous. Because significant downwelling flow and iron-rich material pinned at the CMB are expected in the lowermost mantle beneath the location of accumulated subduction zones such as the present East Asia and past supercontinental margins, long-lasting localized significant cooling of the core might affect the geodynamo (including contributing to causing geomagnetic reversals) in the outer core.

Our results can be explained by a whole-mantle convection model where cold and dense paleoslabs sink to the base of the mantle and trigger upwelling flow of hot and less dense material. This whole-mantle convection is also supported by recent imaging of broad plumes originating at LLVPs at the CMB and continuously connecting to current hotspots at the Earth’s surface (49). We note that it is also possible that our study region is not representative of the whole Earth and that a layered-type convection with decoupling between the upper and lower mantle (or above and below a depth of ~1000-km) might be appropriate for other regions, as suggested by tomographic images of stagnating slabs (22) and geodynamical simulations (53).

MATERIALS AND METHODS

Our data set consists of ~13,000 transverse component records of ground velocity at epicentral distances 70° < Δ < 100° from 40 deep- and intermediate-focus South American earthquakes (see table S1 for details) in the period 2004–2015, recorded at broadband stations of the transportable and backbone arrays of the USArray. We augmented this data set by ~1500 records from a total of 80 South American earthquakes in the period 1993–2015 (also listed in table S1) recorded at CNSN (Canadian National Seismograph Network), CANOE, IRIS/USGS, SCSN, PNSN, and BDSN networks (Fig. 1). The records used in the inversion were selected so that (i) the amplitude ratio between the observed record and the corresponding synthetics was less than 3 and greater than 0.33 and (ii) the variance of the residual (that is, the variance of the difference between the data and synthetics) was less than 300%. After selection, 7768 and 7654 records were used in the inversions for models CACAR and CACAR′, respectively. The large number of records used in this study contributed to the stability of the inversion.

We filtered our data set between 8 and 200 s using a Butterworth bandpass filter. We cut each trace 20 s before the arrival of the direct S phase and 60 s after the arrival of the ScS phase. S and ScS arrivals were computed using the TauP Toolkit (54). ScS precursors (for example, Scd) and postcursors (for example, Sbc), which are sensitive to sharp velocity contrasts in the D″ layer (55), were included in the waveforms in our data set. We weighted each cut residual trace (observed trace minus synthetic trace) so that its maximum amplitudes were all equal to unity. We then applied a second weighting factor to the residuals to partially correct for the uneven azimuthal- and epicentral-distance distribution of the stations in our data set (fig. S1; see the Supplementary Materials for details about the computation of the weighting factors). Each record was also corrected for the effect of the 3D structure outside the target region by a time shift that aligned the onset of the direct S phase on the record with the onset of the direct S phase on the corresponding synthetic (56). We discarded the records for which the time shift was greater than 10 s. Although this did not correct for propagation effects due to strong heterogeneity in the upper mantle, using both S and ScS made our data set primarily sensitive to the lowermost mantle. To verify this, we conducted four different tests shown in figs. S3 to S5 and S12. Figure S3 shows that the partial derivative kernel for the structure in the target region is nearly completely independent from that for the structure in the shallow mantle. Figure S4 shows that (i) simultaneous inversion for the shallow upper mantle and D″ S-velocity structures left the D″ essentially unchanged compared to inversion for D″ only and (ii) inversion for the shallow upper mantle S-velocity structure only yields small amplitude S-velocity perturbations (~0.6%) and small variance reduction (~1%) compared to the inversion for the D″ model only (~5%). Figure S5 shows that a strong 10% velocity decrease in the depth range of 24.4 to 220 km can be nearly completely corrected for using the abovementioned time shifts. Finally, fig. S12 shows that three inverted models using three different data sets selected by dividing the stations into western, central, and eastern regions, where there are different upper mantle structures (57), are in good general agreement in regions where there is common raypath coverage.

We used the methods recently developed by our group for localized 3D waveform inversion (24). Waveform inversion uses all the information in the waveforms, that is, not only travel times of identified phases but also their shape and amplitude. Because waveform inversion compares observed waveforms to synthetic waveforms directly, phases do not have to be identified individually; that is, overlapping phases can be used. This allows the use of ScS when it partially overlaps with S at Δ > 85°. The improvements in resolving power provided by the additional use of shape and amplitude information, and overlapping phases are discussed below.

We computed partial derivatives with respect to a spherically symmetric (1D) initial model using the Born approximation, which gave the first-order perturbation to the wave equation. The 1D synthetics were computed using the direct solution method (DSM) (58). Using the Born approximation with respect to a 1D initial model significantly reduces the computation time, thus allowing the use of a large data set of relatively short period waveforms (down to 8 s).

Every inversion method used some smoothing or regularization techniques. In our case, we truncated the conjugate gradient (CG) expansion for both CACAR and CACAR′ at six CG vectors based on the criterion of minimizing Akaike’s information criterion (AIC) as defined in our recent studies (2, 3).

At present, we fixed the earthquake source parameters to the GCMT (Global Centroid Moment Tensor) catalog. Redetermination of the source parameters and a study of how this affects the 3D models obtained by the inversion is a topic for future work.

Resolution and robustness

The improvement in the fit of the synthetics for the 3D models to the data is shown in table S2. For the initial model PREM, the fit of the initial synthetics to the data (144.4%) was reduced to 79.7% after time-shifting the data to correct for 3D structures above the target region in D″. Inversion of the corrected data yielded synthetics for the final 3D model CACAR with a fit of 74.5%, improving by 5.2% the fit of the initial synthetics for PREM. Similar fit improvements were realized for CACAR′ (see table S2). We note that the relatively small variance reduction is, in part, due to the fact that it is mainly the fit of the part of the waveforms around the smaller amplitude ScS phase (rather than the larger amplitude S phase) that is improved by the inversion.

The similarity between CACAR and CACAR′ suggested that our method and data set yielded robust results. In addition, we conducted various tests (presented in the Supplementary Materials for the partial derivatives computed with respect to PREM) to evaluate and confirm the robustness and resolution of our models. We conducted two checkerboard tests showing that we could resolve an input checkerboard pattern of horizontal scale of 300 × 300 km and that the addition of artificial Gaussian noise to the input synthetics did not affect our ability to recover the input checkerboard pattern (figs. S6 and S7).

Because the Born approximation was only strictly accurate for infinitesimal perturbations, we tested our ability to recover synthetic input models with spherically symmetric perturbations (hereafter referred to as block tests) of ± 1% and 2% in 100- and 200-km-thick layers (~1.7 and 3.4 times the wavelength for 8-s waveforms); the input models were reasonably well recovered (fig. S8). The block tests also showed that (i) we could resolve a 2% low-velocity anomaly 0 to 100 km above the CMB topped by a 2% high-velocity anomaly 100 to 200 km above the CMB (fig. S8C), which suggested that the strong vertical velocity contrast beneath Mexico in CACAR and CACAR′ was not an artifact; (ii) the fact that the strongest high-velocity anomalies in our inferred model were located 0 to 100 km above the CMB was probably real and not an artifact of increase in sensitivity of the partial derivatives just above the CMB (fig. S8, A and B); (iii) CACAR might not accurately constrain the absolute amplitude of the perturbation, which might be overestimated (fig. S8, A and B). However, we note that, in fig. S8C, the recovered amplitude matched that of the input model. For this case, the radially averaged perturbation in the lowermost 400 km of the mantle of the input model, and thus the S-ScS differential travel time (for 70° < Δ < ~85°), was nearly the same as for the initial model (PREM). This is illustrated in fig. S9, in which we showed a profile of stacked waveforms for the input model of fig. S9C. Figure S9 shows that the recovered model in fig. S8C fits the difference in amplitude between the input and initial model (PREM). Because the S-ScS differential travel time was nearly zero, travel-time tomography could not resolve this input model. In addition, fig. S9 shows that the inversion is highly sensitive to data around an epicentral distance of 85°, which cannot be used by travel-time tomography because S and ScS merge at those epicentral distances.

The fact that CACAR′, which was obtained from a different initial model with ~1% higher average velocity in the lowermost 400 km of the mantle, had perturbations with amplitudes comparable to CACAR suggested that the amplitude of perturbations we inferred was relatively robust.

As an additional test of the validity of the Born approximation for inversion of S-velocity in D″, we show a “nonlinear checkerboard test” in fig. S10. The input model is shown in fig. S10A. Synthetic seismograms for this model were computed using full 3D wave calculation (SPECFEM3D GLOBE) (5961). Because of heavy computational requirements with increasing maximum frequency, we computed synthetics down to ~17 s and filtered them with a bandpass filter between 23 and 200 s (as compared to the 8- to 200-s synthetics we used in the inversion of the actual data). We also did not include anelasticity because the synthetics computed using SPECFEM3D GLOBE differed slightly from those computed using the DSM, whereas for the elastic case, which we used for this test, the synthetics for DSM and SPECFEM3D GLOBE for PREM were in nearly perfect agreement. Because we expected longer wavelengths to have less resolving power than shorter wavelengths, we increased the dimension of heterogeneities in the checkerboard pattern to 100 km in the vertical direction (we kept the same lateral dimension of 5° × 5° as in the inversion of actual data). The result (fig. S10B) shows that the inversion using the Born approximation underestimated the absolute amplitude of the perturbations by about 30% (~2% for the inversion result, compared to the 3% perturbation of the input pattern) but that the pattern of high- and low-velocity anomalies was reasonably well recovered.

To test the robustness of CACAR, we also conducted a jackknife test (fig. S11); we conducted three inversions, each of which used 50% of the data randomly picked from our total data set. The result shows that the high-velocity anomalies beneath Central America and Venezuela might be connected in the north of our target region beneath the Caribbean islands (our model was less well constrained in this region due to the smaller number of raypaths) but that a low-velocity corridor along the west coast of South America and extending to the Caribbean Sea was a robust feature of our model; thus, the two high-velocity anomalies we imaged at the CMB beneath Central America and Venezuela were most likely separated. The jackknife test also shows that the strong low-velocity anomaly beneath Mexico is a robust feature of CACAR.

To visually confirm that the final models CACAR and CACAR′ did improve the fit of the synthetics to the data, we show record sections of traces (data and synthetics for initial and final models) for two events (#35 and #49 in table S1) in fig. S13. Figure S13 shows that the synthetics for the final models CACAR and CACAR′ are closer to the data than those for PREM or PREM′. In particular, we note that the double-peaked S wave at Δ > ~95° in the observed traces, as well as an ScS precursor, which both indicate triplication of the S wave due to a strong discontinuity in the D″ layer, were better reproduced by the synthetics for the final models CACAR and CACAR′ than those for the initial models. This visual check was strictly an ancillary measure for quality control.

Although the raypaths in our data set were almost all aligned in the north-south direction, the checkerboard tests (figs. S6, S7, and S10) seemed to show nearly no smearing along the event-receiver path (north-south direction). As an additional check of the absence of smearing, we conducted two point-spread function tests (fig. S14). The results confirmed the absence of smearing in the north-south direction. This suggested that the fact that the high-velocity anomalies “CA” and “VZ” (see Fig. 4) in models CACAR and CACAR′ were elongated in the direction of the event-receiver path was probably a feature of the actual structure in D″ beneath Central America and not an artifact due to smearing.

We suggested a qualitative explanation for the absence of smearing in the north-south direction in fig. S15. We show that the large range of epicentral distances and events latitudes in our data set created a crossing raypath geometry for ScS in a vertical cross section. This crossing raypath geometry of ScS in the lowermost mantle was probably the reason why we could resolve individual voxels in the lowermost mantle using S and ScS waveforms.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/11/e1602700/DC1

Weighting factor

Resolution and robustness tests

fig. S1. Data statistics.

fig. S2. Initial models for 3D inversions.

fig. S3. Trade-off with shallow structure.

fig. S4. Simultaneous inversion for D″ and shallow upper mantle structures.

fig. S5. Possibility of contamination due to shallow structure.

fig. S6. Checkerboard test.

fig. S7. Checkerboard test with artificial noise.

fig. S8. Block tests.

fig. S9. Record section for the block test model.

fig. S10. Nonlinear checkerboard test.

fig. S11. Jackknife test.

fig. S12. Inversion of western, central, and eastern data sets.

fig. S13. Quality control stacks.

fig. S14. Point-spread function.

fig. S15. Crossing raypaths.

table S1. Earthquakes used in this study.

table S2. Variance and AIC.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank S.-H. Hung for sharing the 3D model of D″ beneath the Caribbean. We thank the USArray program, the CANOE program, and IRIS for making a large data set of high-quality data available freely and without delay. We used data from the following networks: TA (USArray Transportable Array), US (U.S. National Seismic Network), CN (CNSN), XN (CANOE), CI (SCSN), UW (PNSN), IU [Global Seismograph Network (GSN)–IRIS/USGS], and BK (BDSN). We thank the editor, G. Beroza, and the reviewers for extensive and constructive comments. We thank the Computational Infrastructure for Geodynamics (http://geodynamics.org), which is funded by the NSF under awards EAR-0949446 and EAR-1550901, for providing SPECFEM3D GLOBE 7.0.0 published under the GPL 2 license. Funding: This research was supported by grants from the Japan Society for the Promotion of Science (nos. 16K05531, 15K17744, and 15H05832). Author contributions: A.F.E.B., K. Kawai, and R.J.G. wrote the manuscript. K. Konishi, K. Kawai, and R.J.G. improved and extended the algorithms used in the inversion. A.F.E.B. and K. Konishi performed the data analysis. K. Kawai and R.J.G. designed the project. A.F.E.B. and K. Konishi helped with the design and evaluation of the project. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The seismic waveforms can all be downloaded freely from the various data centers. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article