Subduction-driven recycling of continental margin lithosphere

LETTER doi:10.1038/nature13878 Subduction-driven recycling of continental margin lithosphere A. Levander1, M. J. Bezada2,3, F. Niu1,4, E. D. Humphreys2, I. Palomeras1, S. M. Thurner1, J. Masy1, M. Schmitz5, J. Gallart6, R. Carbonell6 & M. S. Miller7 Body wave tomography shows that the fast ATL anomaly is larger and lies twice as far south of the plate boundary under northern South America6 than predicted by North American-South American plate motion models for the Cenozoic21,22. The anomaly, well resolved in an east–west direction, extends into the mantle transition zone as a steeply dipping, fast slab (Fig. 2). Lithosphere thickness determined from surface wave tomography and Ps and Sp receiver functions in northern South America8 (Figs 2, 3) indicate that the continental lithosphere is ,50–100% thicker adjacent to and east of the subduction zone than to the west: east of the descending ATL plate at about 64u W, the South America lithosphere is almost uniformly ,110–120 km thick from the Guayana shield to the plate boundary (Fig. 3). To the west, along the coast that the subduction zone has traversed, South American lithosphere is 55–90 km thick between the plate boundary and the Guayana shield (64–68u W, Fig. 3). Restoring the fast ATL anomaly to the surface, the Caribbean-South America CAR OF EPF SSF SdI CdlC Maturin en 10º N Antilles trench Antilles arc 15º N ATL basin pi no gr ab BF SA Es Whereas subduction recycling of oceanic lithosphere is one of the central themes of plate tectonics, the recycling of continental lithosphere appears to be far more complicated and less well understood1. Delamination and convective downwelling are two widely recognized processes invoked to explain the removal of lithospheric mantle under or adjacent to orogenic belts2–5. Here we relate oceanic plate subduction to removal of adjacent continental lithosphere in certain plate tectonic settings. We have developed teleseismic body wave images from dense broadband seismic experiments that show higher than expected volumes of anomalously fast mantle associated with the subducted Atlantic slab under northeastern South America and the Alboran slab beneath the Gibraltar arc region6,7; the anomalies are under, and are aligned with, the continental margins at depths greater than 200 kilometres. Rayleigh wave analysis8,9 finds that the lithospheric mantle under the continental margins is significantly thinner than expected, and that thin lithosphere extends from the orogens adjacent to the subduction zones inland to the edges of nearby cratonic cores. Taking these data together, here we describe a process that can lead to the loss of continental lithosphere adjacent to a subduction zone. Subducting oceanic plates can viscously entrain and remove the bottom of the continental thermal boundary layer lithosphere from adjacent continental margins. This drives surface tectonics and pre-conditions the margins for further deformation by creating topography along the lithosphere–asthenosphere boundary. This can lead to development of secondary downwellings under the continental interior, probably under both South America and the Gibraltar arc8,10, and to delamination of the entire lithospheric mantle, as around the Gibraltar arc11. This process reconciles numerous, sometimes mutually exclusive, geodynamic models proposed to explain the complex oceanic-continental tectonics of these subduction zones12–17. In both the southern Caribbean (CAR) and the Gibraltar arc the subduction zones propagate towards the Atlantic along Mesozoic passive margins, with the oceanic plates tearing away from the adjacent continents as they descend into the mantle. We start with the simpler subduction zone in the southeastern CAR (Fig. 1), where the southern edge of the Atlantic slab, ATL, is interacting with the northern edge of the South American continental margin. Since the Late Palaeocene, the Antilles subduction zone has traversed northern South America as the Americas have moved westward15,18 (Fig. 1), leaving a steeply dipping ATL slab near the southeastern corner of the Antilles subduction zone6. This subduction zone terminates at the El Pilar-San Sebastian strike-slip fault system in northern Venezuela, characterized as a STEP fault19. A zone of intermediate depth seismicity (less than ,150 km), the Paria cluster, is hypothesized to be the mechanical expression of the descending ATL slab tearing from the South America continental margin20. Here we are concerned with what is occurring on the underside of the lithosphere where it is nearly aseismic and viscous, rather than in the brittle regime near the plate surface. Guayana shield 5º N 70º W 65º W 4 3 2 1 60º W 0 1 2 3 4 km Figure 1 | Colour topography and bathymetry of the southeastern Caribbean showing plate boundaries and significant tectonic features. The Antilles subduction zone has traversed northern South America from west to east, creating the coastal mountain belts. The Antilles subduction zone and the limits of orogenic belts are shown in heavy and light white lines with teeth indicating the upper plate. Plate and block-bounding right-lateral strike-slip faults are shown as heavy black and white dashed lines. Solid black line is an inactive strike-slip fault. Colour key shows topography and bathymetry relative to sea level. Red circle indicates the Paria seismic cluster. Black dots are locations of broadband seismograph stations. CAR, Caribbean plate; SA, South American plate; ATL, Atlantic plate attached to South America; SSF, San Sebastian fault; EPF, El Pilar fault; BF, Bocono fault; OF, Oca fault; SdI, Serrania del Interior; CdlC, Cordillera de la Costa. 1 Earth Science Department, Rice University, Houston, Texas 77005-1892, USA. 2Department of Geological Sciences, University of Oregon, Eugene, Oregon 97043, USA. 3Department of Earth Sciences, University of Minnesota, Minneapolis, Minnesota 55455-0231, USA. 4State Key Laboratory of Petroleum Resources and Prospecting, and Unconventional Natural Gas Institute, China, University of Petroleum, Beijing 102249, China. 5Fundación Venezolana de Investigaciones Sismológicas, Caracas 1073, Venezuela. 6Institut de Cièncias de la Terra Jaume Almera, CSIC, Barcelona 08028, Spain. 7 Department of Earth Sciences, University of Southern California, Los Angeles, California 90089-0740, USA. 1 3 N O V E M B E R 2 0 1 4 | VO L 5 1 5 | N AT U R E | 2 5 3 ©2014 Macmillan Publishers Limited. All rights reserved RESEARCH LETTER a a CAR-SA viewed from east Guayana shield Latitude (ºN) 7 Maturin basin 8 Guayana shield 9 SdI EPF 10 SA lithospheric thickness 11º N SSF 11 VS ≥ 4.5 km s–1 Cariaco EPF basin Serrania del Interior Cordillera de la Costa SA lithosphere 100 Base of pre-subduction lithosphere Secondary downwelling Espino graben 9º N 200 Zone of viscous removal 300 dlnVP ≥ 1.5% 400 500 Depth (km) ATL Nose of subducted CM lithosphere Guayana shield 7º N 68º W 600 66º W 64º W 60 80 100 Depth (km) 40 b 120 ATL slab restored to surface CAR SA viewed from north along SA-CAR plate boundary –66 SSF –67 Base of post-subduction lithosphere ATL-SA EPF SdI 10º N 60 80 100 120 140 160 180 Figure 2 | Composite three-dimensional seismic images of the Atlantic plate and South American lithosphere. a, The composite seismic image viewed from the east showing the CAR-SA (Caribbean-South American plate boundary) and SA lithosphere. An iso-surface from a teleseismic P-wave tomography model encloses dlnVP anomalies faster than 1.5% (purple). Another iso-velocity surface from a Rayleigh wave tomography model encloses VS $ 4.50 km s21 (blue), which we take as a proxy for the SA lithospheric mantle. The P-wave tomography images the subducted Atlantic (ATL) slab and subducted continental mantle (CM) lithosphere. The anomaly’s southern end is south of the plate-bounding El Pilar fault (EPF) and extends to the west and south as a nose of high velocity. The lithosphere thickness is ,110–120 km in the region between the Guayana Shield and the ATL slab (dashed red line). A secondary downwelling is imaged by both seismic methods under the Maturin Basin. Seismicity is shown as small orange dots, the red circle indicates the Paria seismic cluster. b, Same surface wave volume as a, viewing the SA-CAR boundary from the north along the Paria seismic cluster. Lithosphere thickness decreases from ,100 km east of the subducting slab to ,75 km west of this slab, as indicated by the red dotted lines. Subduction is ongoing against the SA margin east of 263u and has progressed beyond the SA margin to the west. slab extends several hundred kilometres underneath northern South America, suggesting that the southern edge of the slab anomaly includes South American mantle lithosphere removed from at least as far west as ,69u W (Fig. 3). We hypothesize that the excess mass south of the plate boundary is the base of the continental margin lithosphere that was viscously removed with the ATL plate as it subducted. The Espino graben (Fig. 1) is part of the rifted continental margin that formed along northern and eastern South America at the opening of the central Atlantic ,160 Myr ago. Sufficient time has elapsed since the Atlantic opening for thermal boundary layer (TBL) mantle to form to a depth of at least 100 km under the entire passive margin23,24, comparable to the lithosphere thickness we observe in northeastern South America. The Espino graben lithosphere is now only 55–70 km thick (Figs 2 and 3), 30–45 km less than expected from TBL growth, and thinner than the lithosphere elsewhere to the west. Because rifting left the graben lithosphere en –65 BF ab –64 Depth (km) Base of pre-subduction lithosphere –63 gr –62 Maturin basin ATL no Longitude (ºW) Cordillera de la Costa pi Serrania del Interior SA Es b Antilles trenc h 700 62º W Guayana shield 200 km 5º N Slab restored to the surface 70º W 65º W 60º W Figure 3 | Map of lithospheric thickness in northeastern South America and restoration of the Atlantic slab P-wave anomaly to the surface. a, The lithosphere thickness is less than 100 km everywhere west of about 264u, and north of the Guayana shield, and thicker than 100 km to the east. Red circle indicates the Paria seismic cluster, white circle the location of the secondary downwelling. Abbreviations as Fig. 1, black dashed lines are active plate and block bounding strike-slip faults. Triangles are broadband seismograph stations. b, The ATL slab anomaly restored to the surface. The restored anomaly underlies twice as much of the South American continent as North America–South America plate convergence predicts. The thick dashed line is the South America-Caribbean plate boundary, the thin dashed line is location of the modern passive margin off eastern South America and the nominal location of the palaeo-passive margin off northern South America. Arrow indicates the direction of subduction zone motion relative to South America. Arrowhead terminates at the Paria cluster, where the ATL is currently tearing from South America. thinner than along other parts of the margin, it grew a greater thickness of TBL lithosphere than elsewhere. This lithosphere has now been removed partly or entirely. To us this suggests that the continental lithospheric mantle removed by the descending ATL plate is TBL lithosphere, rather than chemically depleted mantle. An additional deep high-velocity anomaly ,250 km south of the subducting ATL (near 8.75u N) has been imaged by body and surface wave tomography (Fig. 2). It extends to ,200 km depth beneath the southern Maturin basin, the foreland basin of the coastal fold belt. We suggest that this anomaly is a convective instability that developed in response to mantle flow closer to the subduction zone. Instabilities in the continental interior can result from hydration and weakening of the continental lithosphere from exposure to mantle wedge asthenosphere flowing from the subduction zone, combined with continental margin loading by the foreland basin, orogenic belt and descending ATL. Flexure due to loading can drive lower crustal mafic lithologies into eclogite facies, 2 5 4 | N AT U R E | VO L 5 1 5 | 1 3 NO V E M B E R 2 0 1 4 ©2014 Macmillan Publishers Limited. All rights reserved LETTER RESEARCH creating a high-density anomaly in the lower crust. We suggest that viscous removal of TBL mantle during subduction has helped to trigger convective instability under the continental interior. The Gibraltar arc is part of a diffuse plate boundary that extends north across the Atlas Mountains, Alboran Sea and Betic Mountains (Fig. 4). Slab rollback, slab breakoff, lithospheric downwelling and lithospheric delamination have been invoked to explain different aspects of the Gibraltar arc structure12,14,17,25–29. Prior to the initiation of western Mediterranean subduction tectonics in the Oligocene, both southern Iberia and northwestern Africa were passive margins. About 18 Myr of subduction rollback of east-dipping Alpine Tethys seafloor between the two continents (Fig. 4) has left an oceanic slab suspended in the upper mantle beneath the Gibraltar arc7,29 (Extended Data Fig. 1), the remnant of Neogene slab rollback25,26. Although the details are still debated, mantle tomography supports a history in which rollback swept across the Alboran Sea region with rapid east-to-west propagation along the African margin while slowing along the Iberian margin, so that the subduction zone rotated clockwise as it propagated westward towards and then under the Betics, the Rif Mountains, and Gibraltar (Fig. 3). The Alboran crust was drawn westward above the retreating subduction zone7. The simultaneous occurrence of rollback and Eurasian–African convergence led to highly oblique thrusting of Alboran crust over the African and Iberian margins. Below the crustal overthrusting, slab rollback requires that the subducting ocean lithosphere detach from the passive margins along the edges of the Iberian and African plates. The almost vertically hanging Alboran slab is an arcuate feature, concave southeast, underlying both the Iberian and African continental margins at 70–100 km depth, and extending through the 410-km discontinuity7 (Extended Data Fig. 1). Interpretation of GPS data28 and a variety of Rif seismic images9,11 (Extended Data Fig. 2) shows that the Alboran slab is attached to and delaminating from the Rif lithosphere under northern Morocco, the African element of the Gibraltar arc. The seismic measurements detect thicker Rif crust than is predicted from local elevation. a b Gibraltar arc and Atlas 40º N Receiver function images show the top of the Alboran slab merging with the Rif Moho at ,50–55 km depth, evidence that the descending Alboran plate is attached to the Rif lithosphere at the Moho. Under southern Iberia the Alboran slab lies close to the base of the Betic crust on an east–west axis. Here geodetic data, receiver functions, surface wave tomography, and seismicity all suggest that the Alboran slab is still detaching or has recently detached from the base of the Betics9,11,30 (Extended Data Fig. 3). Thus the continental margins are in places still attached to the subducting slab on either side of the Alboran Sea. Rayleigh wave tomography shows thin continental margin lithosphere (50–70 km) around the margins of the Alboran Sea surrounding the descending slab9 (Fig. 4, Extended Data Fig. 1), with thin Iberian lithosphere under the Betics parallel to the direction of slab retreat. Similar to the descending ATL slab, the Alboran slab under southern Iberia is imaged as a high-velocity curtain extending through the transition zone. Restored to the surface, this curtain occupies the Alboran Sea and the adjacent continental margins7 (Fig. 4). We suggest that in addition to oceanic lithosphere, the high-velocity anomaly includes continental margin lithosphere drawn from under Spain and Africa. The structure and tectonics of North Africa are complex, as different tectonic domains interact. The thin lithosphere along the Rif margin extends into the Middle and High Atlas Mountains, where the lithosphere is only slightly thicker than the crust (Fig. 4). The High Atlas, an inverted Triassic-Jurassic graben adjacent to the Sahara craton, has anomalously thin crust for its elevation31. Plio-Pleistocene uplift of the High Atlas is attributed to piecemeal convective downwellings enabled by mantle flow from the Canary plume and/or around the Alboran slab30 that undermined the mantle lithosphere of the Middle and High Atlas and produced intermittent basaltic volcanism. Sites of thin lithosphere in Spain also have associated basaltic volcanism, with a range of chemistries indicative of a hydrated mantle wedge, or a dry asthenospheric source. We suggest that the Alboran subduction has removed continental margin lithosphere from the base of Iberia and North Africa, creating Figure 4 | Gibraltar arc region tectonic features, lithosphere thickness, positive P-wave tomography anomalies, and P-wave anomaly restored to the surface. a, The Gibraltar arc region is defined by the arcuate Rif-Betic mountain system. The diffuse plate boundary separating Eurasia (Iberia) from Africa (AF) extends from the Betics across Gibraltar (G), the Alboran Sea and Rif Mountains to the southern edge of the Atlas Mountains. Slab rollback proceeded from approximately northeast to west (see d). Limits of Betic and Rif thrust system are shown in white. ATL, Atlantic plate. TASZ, left-lateral TransAlboran shear zone. Black dots are broadband seismograph stations. Solid black lines are faults and geologic province boundaries b, Map of lithospheric thickness in the Gibraltar arc region. The lithosphere is thin north, east and south of the Alboran slab, and thinnest to the northeast and east of the slab beneath the Betics and the Alboran slab along the direction of trench retreat. c, The slab anomaly is shown as an iso-velocity surface enclosing dlnVP $ 1.5%. The anomaly extends downward beneath the Rif (A) and Betics (B) and dips east from Gibraltar under the Alboran Sea. A curtain of high velocity underlies southern Spain (C). d, The Alboran slab anomaly restored to the surface. Dashed lines show Spain 30 Myr ago with nominal Iberia and African passive margins. Arrow indicates the direction of subduction zone motion relative to Spain’s current position. LAB depth 40º N Iberia Iberia tics SZ 35º N ATL At la s Rif ATL M id dl e id dl e At la s 35º N TA Alboran Sea G tics Be G Alboran Sea Rif AF M as Atl 30º N 10º W as Atl Sahara craton Sahara craton 5º W 0º 30º N 10º W AF 5º W 0º 30 40 50 60 70 80 90 100 110 120 130 Depth (km) 4 3 2 1 0 1 2 3 4 km c TA tics Be SZ Be d Alboran slab Present Slab restored to the surface Depth (km) 0 A B 200 C 400 600 45 800 10 5 Longit u 0 de (º) 35 –5 30 Latit 40 (º) ude 200 km ~30 Myr ago 1 3 N O V E M B E R 2 0 1 4 | VO L 5 1 5 | N AT U R E | 2 5 5 ©2014 Macmillan Publishers Limited. All rights reserved RESEARCH LETTER lithosphere–asthenosphere boundary (LAB) topography that helped generate secondary downwellings interior to both continental margins. Both the southeastern Caribbean and the Gibraltar arc subduction zones are adjacent to former Mesozoic passive margins that were largely unaffected by tectonic activity between the time of rifting and the arrival of the modern subduction zones. Following rifting, like the developing ocean basin lithosphere, the passive margins form a TBL lithosphere of thickness comparable to oceanic lithosphere, that is, ,100 km, within ,70 Myr (refs 23, 24). The seismic images from northeastern South America and the Gibraltar arc demonstrate that the adjacent continental margin lithosphere is thinner than expected along the margins traversed by the subduction zones. We attribute this to convective removal of the continental margin mantle lithosphere by the adjacent subducting oceanic plate. In the Gibraltar arc this process extends to levels as shallow as the Moho, resulting in lithospheric delamination beneath both adjacent continental margins. We speculate that removal of lithospheric mantle to within a plate thickness of the continental margin creates gradients in the LAB that result in additional downwellings under the continental interior. We have imaged one such downwelling beneath the South America continent, and infer that previous downwellings in the South America interior removed lithosphere north of, but close to the Guayana shield. Slab rollback and lithosphere removal from the Gibraltar arc margins occurred shortly before large-scale downwellings beneath the Middle and High Atlas, suggesting a relation between the two. Potentially hotter asthenosphere and asthenospheric flow around the narrow Alboran slab resulted in delamination and widespread volcanism. Loss of mantle lithosphere to downwellings and delamination beneath growing and collapsing mountain belts is well recognized locally in the Andes and in the western United States as an indirect consequence of subduction and its aftermath2–5. What we describe is more similar to ‘ablative subduction’ predicted by two-dimensional geodynamic modelling32, in which a subducting oceanic plate viscously removes the TBL mantle on an overriding continental plate. Here, oceanic plate subduction, the primary convection system, viscously removes the TBL from beneath an adjacent continental margin, creating LAB topography under the continent that can give rise to secondary, edge convection33. The mantle lithosphere that remains beneath the margins after removal of the basal TBL is likely to be the buoyant, basalt-depleted peridotite that forms the chemical boundary layer of the continental mantle lithosphere. Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper. Received 6 September 2013; accepted 19 September 2014. 1. Gao, S. et al. Recycling deep cratonic lithosphere and generation of intraplate magmatism in the North China Craton. Earth Planet. Sci. Lett. 270, 41–53 (2008). 2. Houseman, G. & Molnar, P. in Continental Reactivation and Reworking (eds Miller, J. A., Holdsworth, R. E., Buick, I. S. & Hand, M.) 13–38 (Spec. Publ. 184, Geological Society of London, 2001). 3. Zandt, G. et al. Active foundering of a continental arc root beneath the southern Sierra Nevada in California. Nature 431, 41–46 (2004). 4. Kay, R. W. & Kay, S. M. Delamination and delamination magmatism. Tectonophysics 219, 177–189 (1993). 5. Levander, A. et al. Continuing Colorado plateau uplift by delamination-style convective lithospheric downwelling. Nature 472, 461–465 (2011). 6. Bezada, M. J., Levander, A. & Schmandt, B. Subduction in the southern Caribbean: Images from finite-frequency P wave tomography. J. Geophys. Res. 115, B12333 http://dx.doi.org/10.1029/2010JB007682 (2010). 7. Bezada, M. J. et al. Evidence for slab rollback in westernmost Mediterranean from improved upper mantle imaging. Earth Planet. Sci. Lett. 368, 51–60 (2013). 8. Masy, J., Levander, A., Niu, F. & Schmitz, M. Lithospheric expressions of Cenozoic subduction, Mesozoic rifting and the Precambrian Shield in Venezuela. Earth Planet. Sci. Lett. (in the press). 9. Palomeras, I., Levander, A., Liu, K., Thurner, V. A. & Gallart, J. Finite-frequency Rayleigh wave tomography of the western Mediterranean. Geochem. Geophys. Geosyst. 15, 140–160 (2014). 10. Bezada, M. J., Humphreys, E. D., Palomeras, I., Levander, A. & Carbonell, R. Piecewise delamination of Moroccan lithosphere from beneath the Atlas Mountains. Geochem. Geophys. Geosyst. 15, 975–985 (2014). 11. Thurner, S. M., Palomeras, I., Levander, A., Carbonell, R. & Lee, C. T. A. Ongoing lithospheric removal in the Western Mediterranean: Ps receiver function results from the PICASSO project. Geochem. Geophys. Geosyst. 15, 1113–1127 (2014). 12. Lonergan, L. & White, N. Origin of the Betic-Rif mountain belt. Tectonics 16, 504–522 (1997). 13. Platt, J. P. & Vissers, R. L. M. Extensional collapse of thickened continental lithosphere: a working hypothesis for the Alboran Sea and Gibraltar arc. Geology 17, 540–543 (1989). 14. Duggen, S., Hoernle, K., Bogaard, P. d. & Harris, C. Magmatic evolution of the Alboran region: the role of subduction in forming the western Mediterranean and causing the Messinian Salinity Crisis. Earth Planet. Sci. Lett. 218, 91–108 (2004). 15. Pindell, J. & Kennan, L. in The Geology and Evolution of the Region Between North and South America (eds James, K., Lorente, M. A. & Pindell, J.) 1–55 (Spec. Publ. 328, Geological Society of London, 2009). 16. VanDecar, J. C., Russo, R. M., James, D. E., Ambeh, W. B. & Franke, M. Aseismic continuation of the Lesser Antilles slab beneath continental South America. J. Geophys. Res. 108, 2043 http://dx.doi.org/10.1029/2001JB000884 (2003). 17. Calvert, A., Sandvol, E., Seber, D. & Barazangi, M. Geodynamic evolution of the lithosphere and upper mantle beneath the Alboran region of the western Mediterranean: constraints from travel time tomography. J. Geophys. Res. 105, 10871–10898 (2000). 18. Burke, K. Tectonic evolution of the Caribbean. Annu. Rev. Earth Planet. Sci. 16, 201–230 (1988). 19. Govers, R. & Wortel, M. J. R. Lithosphere tearing at STEP faults: response to edges of subduction zones. Earth Planet. Sci. Lett. 236, 505–523 (2005). 20. Clark, S. et al. Eastern Venezuelan tectonics driven by lithospheric tear geodynamics, not oblique collision. Geochem. Geophys. Geosyst. 9, Q11004 http://dx.doi.org/10.1029/2008GC002084 (2008). 21. Pindell, J., Kennan, L., Stanek, K.-P., Maresch, W. & Draper, G. Foundations of Gulf of Mexico and Caribbean evolution: eight controversies resolved. Geol. Acta 4, 303–341 (2006). 22. Müller, R. D., Royer, J.-Y., Cande, S. C., Roest, W. R. & Maschenkov, S. in Sedimentary Basins of the World Vol. 4 (ed. Mann, P.) 33–59 (Elsevier, 1999). 23. Sleep, N. H. Thermal effects of the formation of Atlantic continental margins by continental break up. Geophys. J. R. Astron. Soc. 24, 325–350 (1971). 24. Sclater, J. G., Parsons, B. & Jaupart, C. Oceans and continents: similarities and differences in the mechanisms of heat loss. J. Geophys. Res. 86, 11535–11552 (1981). 25. Rosenbaum, G. Tectonic Reconstruction of the Alpine Orogen in the Western Mediterranean Region PhD thesis, Monash Univ. (2003). 26. Royden, L. H. Evolution of retreating subduction boundaries formed during continental collision. Tectonics 12, 629–638 (1993). 27. Seber, D., Barazangi, M., Ibembrahim, A. & Demnati, A. Geophysical evidence for lithospheric delamination beneath the Alboran Sea and Rif-Betic mountains. Nature 379, 785–790 (1996). 28. Pérouse, E., Vernant, P., Chéry, J., Reilinger, R. & McClusky, S. Active surface deformation and sub-lithospheric processes in the western Mediterranean constrained by numerical models. Geology 38, 823–826 (2010). 29. Gutscher, M.-A. et al. Evidence for active subduction beneath Gibraltar. Geology 30, 1071–1074 (2002). 30. Mancilla, F. L. et al. Delamination in the Betic Range: deep structure, seismicity, and GPS motion. Geology 41, 307–310 (2013). 31. Ayarza, P. et al. Crustal thickness and velocity structure across the Moroccan Atlas from long offset wide-angle reflection seismic data: the SIMA experiment. Geochem. Geophys. Geosyst. 15, 1698–1717 (2014). 32. Tao, W. C. & O’Connell, R. J. Ablative subduction: a two-sided alternative to convectional subduction model. J. Geophys. Res. 97, 8877–8904 (1992). 33. Kaislaniemi, L. & van Hunen, J. Dynamics of lithospheric thinning and mantle melting by edge-driven convection: application to Moroccan Atlas mountains. Geochem. Geophys. Geosyst. 15, 3175–3189 (2014). Acknowledgements We thank R. Govers for suggestions that improved the clarity and quality of the manuscript, and E. Engquist for aid in using the Rice DAVinCI Visualization Laboratory. We especially thank M. Harnafi and the Scientific Institute of Rabat for their contributions to the project. This research was supported by US National Science Foundation grants EAR 0003572, 0607801 and 0808939 (A.L.), EAR 0808931 (E.D.H.), EAR 0809023 and 1054638 (M.S.M.), the Venezuelan National Fund for Science, Technology and Innovation grant G-2002000478 and PDVSA-INTEVEP-FUNVISIS cooperative agreement 2004-141 (M.S.), the Spanish Ministry of Science and Innovation grants CSD2006-00041, CGL2009-09727 and CGL2010-15146 (J.G. and R.C.), and by an A. v. Humboldt Foundation Research Prize (A.L.). Author Contributions M.S., J.G., R.C., E.D.H., F.N., M.S.M. and A.L. oversaw different aspects of the field data acquisition. A.L., F.N., S.M.T., J.M., M.S. and R.C. contributed to the receiver function data analysis. I.P., F.N., J.M., M.S.M., J.G. and A.L. contributed to the Rayleigh tomography analysis. M.J.B., E.D.H. and A.L. contributed to the body wave tomography analysis. M.S., R.C. and J.G. provided geologic, tectonic and geophysical background, which allowed A.L. and E.D.H. to pose the lithosphere removal hypothesis for testing. M.J.B., I.P., S.M.T., J.M. and A.L. constructed and interpreted the 3D images. A.L primarily wrote the manuscript, with substantive input from E.D.H., M.J.B. and F.N. and with additional input from all of the co-authors. Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to A.L. ([email protected]). 2 5 6 | N AT U R E | VO L 5 1 5 | 1 3 NO V E M B E R 2 0 1 4 ©2014 Macmillan Publishers Limited. All rights reserved LETTER RESEARCH METHODS Tomography. Regional finite-frequency seismic body wave tomography is a means of including Fresnel zone effects in constructing three-dimensional Earth images from teleseismic travel time residuals. Nonlinear tomography incorporates threedimensional ray tracing through the imaged structure7. The finite-frequency Rayleigh wave tomography34,35 images are from regional studies of the southeastern Caribbean plate boundary and the Gibraltar arc8,9. This form of Rayleigh wave tomography gives absolute radially averaged VS measurements, and, because the waves travel horizontally, provides somewhat better vertical resolution of upper mantle VS structure than teleseismic body wave tomography. Rayleigh waves in the band we examined (0.0067–0.05 Hz) are sensitive to VS structure from the mid-crust, ,20 km depth, into the upper mantle, ,250 km depth. Receiver functions. Receiver function imaging makes a scattered wave image of subsurface seismic impedance boundaries using P to S or S to P converted waves from teleseismic earthquakes36. Ps receiver functions are made by deconvolving direct P wave on the vertical component of motion from the radial motion to remove the earthquake source function, and replace it with a known shaping pulse. Using the incidence angle of the P wave and an estimated velocity model, the receiver functions are individually back-projected to the conversion points. The partial images made from many earthquakes recorded by many stations are summed, providing a three-dimensional image and improving signal to noise ratio37. Fresnel zone effects can be included in the stacking procedure38. The resulting three-dimensional image volume is a common conversion point stacked image. The process for making Sp receiver functions is similar but includes rotation of the seismograms into the direction of the S wave at the surface. Direct teleseismic S waves have lower frequency (0.2–0.05 Hz) content than the teleseismic P waves (1.0–2.0 Hz), and therefore produce a lower resolution image of the subsurface, but have other advantages, notably, Sp receiver functions are not contaminated by families of multiple S and P wave reverberations in the crust. Composite seismic images. Composite seismic images are made by superimposing elements of images made from different seismic measurements. In three dimensions, volumes for superposition are constructed by choosing threshold values for isosurfaces. Since different seismic probes identify different earth structures, the isosurface thresholds are chosen to represent some aspect of seismic structure that is either unique or common to the different measurements. Regional body wave tomography is good at resolving lateral velocity variations such as those due to a subducting slab. For slabs identified in body wave tomography we chose an isosurface enclosing dlnVP $ 11.5%. Surface wave tomography, which measures absolute velocity, is better for identifying the high and low velocities in the upper mantle lithosphere and asthenosphere, respectively, and when used in conjunction with receiver functions provides a robust means of determining lithosphere thickness. We chose an isosurface enclosing VS $ 4.5 km s21 as indicative of the lithospheric mantle. Combined the two volumes show both lithosphere and slab structure, which neither seismic probe determines completely by itself. 34. 35. 36. 37. 38. Forsyth, D. W. & Li, A. in Seismic Earth: Array Analysis of Broadband Seismograms (eds Levander, A. & Nolet, G.) 81–97 (Geophys. Monogr. Ser. Vol. 157, American Geophysical Union, 2005). Yang, Y. & Forsyth, D. W. Rayleigh wave phase velocities, small-scale convection and azimuthal anisotropy beneath southern California. J. Geophys. Res. 111, B07306 http://dx.doi.org/doi:10.1029/2005JB004180 (2006). Rondenay, S. Upper mantle imaging with array recordings of converted and scattered teleseismic waves. Surv. Geophys. 30, 377–405 (2009). Dueker, K. G. & Sheehan, A. F. Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track. J. Geophys. Res. 102, 8313–8327 (1997). Levander, A. & Miller, M. S. Evolutionary aspects of lithospheric discontinuity structure in the western U.S. Geochem. Geophys. Geosyst. 13, Q0AK07 http://dx.doi.org/10.1029/2012GC004056 (2012). ©2014 Macmillan Publishers Limited. All rights reserved RESEARCH LETTER Extended Data Figure 1 | Composite seismic image showing the top of the Alboran slab and the lithosphere beneath the Gibraltar arc. a, Top panel is viewed from above from the east-northeast. Topography is shown at the top of the panel. The bottom of the panel is a composite of a P-body wave tomography image showing the slab (magenta, with the isosurface enclosing dlnVP $ 1.5%), and a Rayleigh wave tomography image showing the top of the slab and the lithosphere (blue, with the isosurface enclosing VS . 4.5 km s21). The dashed black line outlines the bottom of the lithosphere. Note that these lines do not represent depth in the perspective view. b, Same azimuthal view as a but viewed from below. The black and white dashed lines outline the bottom of the lithosphere. ©2014 Macmillan Publishers Limited. All rights reserved LETTER RESEARCH Extended Data Figure 2 | Surface wave tomography model and receiver function images from northern Morocco. Top panel, Rayleigh wave tomography model along 35u N. Middle and bottom panels, 2 Hz Ps receiver function CCP stacks along 35u N (middle) and 34.75u N (bottom) showing the top of the lower crust (dashed black lines), the Moho (solid black line) and the top of the Alboran slab (dashed white line) beneath the Moroccan Rif. In the two receiver function images the Moho and the top of the Alboran slab merge at ,50 km depth at 24.5u and diverge to the east. Moho depth from unpublished refraction profiles is shown by heavy grey line. Seismicity, shown as white diamonds, is concentrated at the Trans-Alboran shear zone (TASZ). The seismic images are shown with no vertical exaggeration. ©2014 Macmillan Publishers Limited. All rights reserved RESEARCH LETTER Extended Data Figure 3 | Surface wave tomography model and receiver function images from southern Spain. Top panel, Rayleigh wave tomography model along 37u N. Middle and bottom panels, 2 Hz Ps receiver function CCP stacks along 37u N (middle) and 36.75u N (bottom) showing the top of the lower crust (dashed black lines), the Moho (black solid lines) and the top of the Alboran slab (heavy dashed white line) beneath the Betics. In the two receiver function images the Moho and the top of the Alboran slab merge at ,50–55 km depth at 24u and diverge in either direction. Seismicity, shown as white diamonds, occurs in the upper crust and in the zone of detachment near the base of the Iberian crust and the top of the Alboran slab. ©2014 Macmillan Publishers Limited. All rights reserved