Mobilization and thinning of cratonic lithosphere by a lower mantle slab

Mobilization and thinning of cratonic lithosphere by a lower mantle slab | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Mobilization and thinning of cratonic lithosphere by a lower mantle slab Junlin Hua, Steve Grand, Thorsten Becker, Helen Janiszewski, Chujie Liu, and 2 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3254038/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 28 Mar, 2025 Read the published version in Nature Geoscience → Version 1 posted You are reading this latest preprint version Abstract Continental cratons are characterized by thick lithospheric roots that remain intact for billions of years. However, some cratonic roots appear to have been thinned or completely removed, with the reasons for such thinning being debated. In this study, we obtain a high-resolution full-waveform seismic tomographic model for North America which newly illuminates ongoing craton-thinning. Extensive drip-like transport of lithosphere is imaged from the base of the craton beneath the central United States to the mantle transition zone. Geodynamical modeling suggests that such dripping may be mobilized by the sinking of the deep Farallon slab, whose associated mantle flow can drag material at the base of the craton from afar to the dripping location. There, lithospheric material can descend within the ambient downward mantle flow, even though the slab is presently in the lower mantle. Dripping lithosphere could be further facilitated by prior lithospheric weakening such as due to volatiles released from the slab. Our findings show how cratonic lithosphere can be altered by external forces, and that subduction can play a key role in craton mobilization and thinning even when slabs are at great depths in the mantle. Earth and environmental sciences/Solid Earth sciences/Seismology Earth and environmental sciences/Solid Earth sciences/Geodynamics Earth and environmental sciences/Solid Earth sciences/Tectonics Figures Figure 1 Figure 2 Figure 3 Figure 4 Main Cratons are old (no younger than Proterozoic) regions which make up ~60% of continental lithosphere 1 and are often underlain by thick lithosphere which can extend to over 200 km depth 1,2 . Cratons are believed to have been overall stable over billions of years due to their near neutral effective buoyancy and high viscosity 3,4 . However, the North China Craton lost its root during the Mesozoic 5 , for example, indicating that the stability of cratons is not universal. Different mechanisms have been proposed for the loss of cratonic lithosphere, including convective removal 6 , flat-slab rollback 7 , and erosion with melting 8,9 . Other studies suggest the cratonic lithosphere is not stable through history but subject to reworking 10 . Given that most of these examples are from processes in the geological past, it is not straightforward to distinguish between different mechanisms, and new observations of an ongoing craton-thinning event would offer a unique perspective to the puzzle. North America is an ideal place to study such present-day behavior of cratons. Core continental North America consists of multiple Archean cratons and volcanic arcs accreted during the Proterozoic. Although most cratonic regions in the U.S. were stabilized during the Proterozoic from 1.35-1.8 Ga 11 , they still share key features with their Archean counterparts, including very thick lithosphere 2 and surrounding kimberlite volcanism 12 . Meanwhile, thanks to dense seismometer deployments, the Central and Eastern United States cratons can be studied in detail seismically. In this study, we obtain a full-waveform tomography model for North America (Fig. 1), named Seismic Adjoint Tomography of North America (SATONA). Compared with existing North America tomographic models, SATONA focuses more on fitting the detailed shape of waveforms at relatively high frequencies to better reveal fine-scale structures. SATONA reveals clear lithosphere-asthenosphere boundaries (LAB) for the craton, and high seismic velocity materials beneath the craton from the LAB to the mantle transition zone (MTZ). These high-velocity, drip-shape features indicate thermo-chemical differences relative to the surrounding mantle and likely originate from the cratonic lithosphere. To find the cause for such dripping, we test the influence of the lower mantle Farallon slab on upper mantle convective flow through geodynamic modeling, and find that the dripping may be a consequence of the slab sinking. These findings lead us to propose a new mechanism for craton thinning —large-scale downwelling induced by a deep-seated slab that mobilizes the craton, possibly mediated by fluid release, leading to partial removal of the cratonic lithosphere. Dripping cratonic lithosphere Model SATONA is obtained through full-waveform adjoint tomography. During inversion, we maximize the correlation coefficient between synthetic and observed waveforms 13 . Hence, besides arrival times of seismic phases, the detailed shape of phases including secondary arrivals are required to match the data. In this study, both synthetic seismograms and adjoint kernels were calculated with SpecFem3D Globe 14,15 , and all available stations in the study region were utilized (Extended Data Fig. 1). We used multiple body and surface wave phases to constrain S- and P-wave velocities 16 . To optimize computational efficiency and avoid local minima, the inversion consists of two stages 17 : a coarse-grid stage using 206 earthquakes, and a fine-grid stage using 110 earthquakes with higher quality (Extended Data Fig. 1). We incorporate the mini-batch algorithm 18,19 by using a subset of earthquakes in each iteration, which, compared with common full-waveform tomography practices using only a few tens of iterations on supercomputers 16 , enables a total of 216 iterations with 4 GPUs. With these, the improvement on the correlation coefficient and the arrival time difference between synthetic and observed waveforms is significant (Extended Data Fig. 2). Body waves show higher accuracy in timing than surface waves, while surface waves have the correlation better recovered. Nonetheless, the vast majority of both wave types show correlation coefficients over 0.9, suggesting an excellent fitting of the waveform shape. Detailed mantle structures are newly imaged in SATONA. S-wave velocity (V S , Extended Data Fig. 3) at 100 km depth correlates well with tectonic regions 20 , and at greater depths, the model captures features such as the Cascadia slab in fine detail. P-wave velocity (V P , Extended Data Fig. 4), which is independently inverted for, shows similar patterns. A resolution test with only 16 iterations (versus 206 for the full model) shows that V S has good resolution down to 1000 km, especially within the upper mantle (Extended Data Fig. 5). V P also has good resolution down to the MTZ, but deeper resolution is reduced (Extended Data Fig. 5). We only discuss regions with reliable resolution. The velocity perturbation model also shows some horizontally aligned anomalies near reference model discontinuities (e.g., near 660 km depth in Fig. 2a, b). However, as discontinuity depths are taken from previous models upon initialization of the inversion, these anomalies might be introduced to reconcile inaccuracies about discontinuity depth in those starting models, and are therefore not interpreted much in this study. Cratonic lithosphere is well-imaged in the model. At 200 km depth, a broad region beneath the Central and Eastern U.S. shows high V S (Fig. 1), marking the geographic extent of the deep cratonic root. Cratonic LABs are clearly characterized as a velocity decrease with depth around 200 km in both V S and V P (Figs. 2a,b and Extended Data Fig. 6). Beneath the cratonic LAB, extensive high V S anomalies are observed down to the MTZ (Fig. 2, Extended Data Fig. 3). A specific model resolution test (Extended Data Fig. 7) suggests these anomalies can be well imaged with our dataset and inversion strategy. Similar high velocity anomalies are also found in previous models 21-25 (Extended Data Fig. 8). However, without the help of multiple seismic phases, their cratonic lithosphere is not well constrained (Extended Data Fig. 8), and the deep structures sometimes extend into the lithosphere depth range, making them harder to interpret. To distinguish whether these high V S anomalies originate from the cratonic lithosphere, we convert V S and V P perturbations of the dripping bodies to temperature and compositional perturbations. Compared with V S , V P shows much weaker anomalies at 300-500 km depths (Extended Data Fig. 6), distinct from the strong V P anomaly in the Cascadia slab, suggesting a potentially different composition. To constrain the compositional perturbation, two endmembers are considered. One composition is assumed to be pyrolite 26 and the other to be cratonic lithosphere. We estimate a representative major element composition for the craton from kimberlite xenoliths (Extended Data Fig. 9). The V S and V P for these two endmember compositions at different temperatures and pressures are then calculated using Perple_X 27,28 , accounting for anelastic effects 29 . Derivatives of V S and V P with respect to temperature and composition are obtained, and velocity anomalies are then converted to temperature and compositional perturbations. We performed the conversion for all geographic locations based on a smoothed version of SATONA (Extended Data Fig. 9) to eliminate bias due to the misalignment of V S and V P anomalies. Taking the compositional perturbation between the endmembers as unity, at all depths, the high V S anomalies show 20-50% more cratonic composition than the ambient mantle (Extended Data Fig. 9), which indicates the entrainment of lithospheric material. Especially below 300 km, due to the phase transition of orthopyroxene at ~320 km 28 , the dependence of V S -V P discrepancy on composition is significant, making the compositional perturbation better constrained. At 350 km depth, the median of the dripping high V S bodies has 45% more cratonic composition than normal (Fig. 2d), which is close to the 75% percentile of compositional perturbations for regions excluding the dripping area (Extended Data Fig. 9m), indicating the high V S anomalies are substantially more cratonic than normal. Despite having higher V S , when converted to temperature perturbation, we found these anomalies are only ~60ºC colder than normal (Fig. 2e). Even for less constrained depths, the thermal anomaly is less than 200ºC (Extended Data Fig. 9j). This test demonstrates these compositionally distinct anomalies could be lithospheric drips, and they are hereafter referred to as dripping anomalies. Deep-seated slab induced dripping Similarly imaged dripping anomalies in the region have previously been interpreted as remnant slabs 22 or dripping due to Rayleigh-Taylor instability 30 , but both mechanisms have some limitations. For slab relics, relative values of V S and V P anomalies do not match those of the Cascadia slab whose V P anomalies are much stronger (Fig. 2c and Extended Data Fig. 6). There is also no evidence of arc volcanism in the study region since at least the Mesozoic 31 , which is a long time for stalling thermal anomalies. In particular, large scale downward flow likely occurred in this region for 70 Ma and thus it is not likely for slab relics to persist. If the anomalies are purely from a Rayleigh-Taylor instability due to a dense craton root 30 , it is hard to explain the relatively weak V P anomaly for the drips (Extended Data Fig. 6), since with the inclusion of a large percentage of dense minerals like garnet, V P would also be increased 32 . Meanwhile, if a uniform dense root is common among cratons, similar dripping would likely be observed beneath every craton, and likely have occurred for over a billion years and hence might have consumed all dense material. We propose a new mechanism that the deep Farallon slab is the main driver of dripping, though it is in the lower mantle and not connected to the surface at present. To evaluate the dynamic influence of the Farallon slab, we calculated mantle flow 33 for two different density structures (Fig. 3): mantle with or without the Farallon slab. These structures are modified from a global mantle tomography model 34 . Results show that large-scale mantle flow is strongly controlled by the Farallon slab, which pulls shallow mantle from the East, West, and North to the dripping area where these mantle materials flow downward (Fig. 3a), in line with previous work 35 . Such flow is only predicted when the Farallon slab is present (Fig. 3a, b), and is absent without it (Fig. 3c, d). The Farallon slab may induce horizontal mantle flow over a domain wider than ~3000 km from the sinker (Fig. 3a versus 3c), and associated shear flow may serve to thin and entrain the base of the craton (Fig. 3b). That is, even though the imaged dripping area has relatively limited extent as in Fig. 2c, the material involved in the sinker may come from a much wider region. Such basal erosion by slab-induced flow leading to dripping removal of continental lithosphere may be assisted by prior weakening of the lithosphere, e.g. due to volatile influx from the slab. Associated slab-induced dripping could have started at ~70 Ma, as the Farallon slab has stayed at a similar relative location with respect to North America since then 36 (Fig. 1). Based on the tomographic model, by assuming the V S anomaly is proportional to the amount of cratonic materials, with the consideration of horizontal flow concentrating material from the whole imaged cratonic region (Fig. 1), the total thinning associated with the currently imaged anomaly may be of order ~20 km. This number might be an overestimate if the slab-induced horizontal flow could influence an even larger area (Fig. 3), for example, than the cratonic region defined in Fig. 1 that is covered by this model, and the edge of the craton (Fig. 2) at shallower depth could extend wider (Extended Data Fig. 3b) than the cratonic region at 200 km here. Additionally, we chose to use the stronger V S anomaly instead of V P anomaly for this estimation. Because the dripping material appears more cratonic than mantle in composition, its effective negative buoyancy may be reduced 3 . To evaluate whether the flow caused by the Farallon slab is strong enough to pull down these materials, we merged our V S anomalies from SATONA with global tomography 34 to estimate the effect of different density anomalies (Extended Data Fig. 10). To account for the range of potential densities of the dripping body, we tested different scaling factors between V S and density (Extended Data Fig. 10), and in one case we make density negatively correlated with velocity for the dripping material, so that they are positively buoyant. While the sinking speeds are expectedly reduced for this case, the mantle still flows downward at the dripping locations (Extended Data Fig. 10d). This is partly because strong horizontal inflow driven by the Farallon slab always converges at the dripping location. This test shows that the downward dragging force by the slab could overcome any neutral and even small positive buoyancy of lithospheric material once mobilized. Compound process of craton thinning To enable dripping of cratonic materials, the base of the craton cannot be very rigid, so that downward tractions applied by the deep slab and shear at the LAB could indeed mobilize the lithosphere and entrain material. Such a relatively weak lower part of the craton has previously been inferred. For instance, cratonic kimberlite xenoliths from larger depths often show more evidence for deformation than shallower ones 37,38 . Seismically, cratons often consist of two layers, separated by the mid-lithospheric discontinuity 39 , as also evident in some locations in SATONA (Fig. 2a,b), and the lower half can have strong seismic anisotropy (beneath ~110 km in Extended Data Fig. 6) consistent with a higher degree of deformation. If such deformation is recent, this may suggest a lower strength for the deep craton, and/or discontinuities may be mechanically weak. It has also been argued that the lower half of the craton is not always the original lithosphere, and could have experienced reworking 40 . Hence, it is likely that the base of the lithosphere could be weakened 41 (Fig. 4a), dynamically only marginally stable, and subject to deformation when external forces exist. Moreover, with the inclusion of a small amount of high-density minerals like garnet in the deep part of the craton 10 , the dripping could be easier to develop. Besides driving the flow to mobilize the lithosphere mechanically, past subduction could also facilitate sub-lithospheric convection and lithospheric weakening 42 (Fig. 4a). When subducting slabs enter the MTZ, volatiles can be released due to dehydration and decarbonization 43-46 . Some of these volatiles would then ascend due to mechanisms such as subduction-induced poloidal flow 47 and the breakdown and re-oxidation of hydrous phases 9 , reduce the asthenospheric bulk viscosity, and possibly form hydrous-carbonated melts ponded beneath the LAB 45,46 , which could lower the strength of the bottom of the lithosphere. There is evidence that such subduction-induced weakening has occurred. Relative to North America, the Farallon slab was located further west before ~100 Ma 36 . Around that time, kimberlite, lamproite, and carbonatite magmatism occurred at the western edge of the craton (Fig. 1), and such magmatism types favor a carbon-enriched mantle source 45,48 . After 70 Ma, the deep slab moved eastward 36 , and kimberlite magmatism appeared near the eastern edge of the craton in Kentucky (Fig. 1). Such findings are consistent with ponded hydrous-carbonated melts (Fig. 4), causing low V S anomalies observed at the eastern and western edges of the craton and along the east coast (Figs. 1 and 2a, b). Past plume activity could also contribute to the weakening. Previous studies suggest hotspot tracks across the continental US 49,50 , which could weaken the lithosphere along their trajectories. When a plume rises up from the lower mantle, the 660-km phase change with a negative Clapeyron slope could temporally stall plume ascent, causing hot materials to pond beneath the MTZ and lead to a transient boundary layer 51 (Fig. 4a), and plume instabilities might then ascend when encountering the Farallon slab, leading to further weakening of the craton. We propose that a deep, compositionally distinct craton is being pulled down by slab-induced flow. The small ~60°C thermal difference of the drips (Fig. 2e) suggests they are likely from the bottom of the craton, which could have been weakened by the processes discussed above. These bottom materials are continuously entrained and mobilized by the slab-induced horizontal flow which concentrates them in the dripping area (Fig. 3a), where they then follow the downward flow to larger depths. Therefore, compared to delamination, where a large volume of lithosphere is lost locally, the thinning in this study is achieved by constantly mobilizing and removing the bottom of the lithosphere from areas affected by the horizontal flow, which could be larger than the surface projection of the dripping area (Fig. 4b). Combining all evidence, the lower part of a craton may be weaker and prone to deformation (Fig. 4a). Normally, these weakened portions remain part of the lithosphere due to lower density, or convective compression effects due to steep LAB gradients 52 . However, if strong tractions from below are present, such as when the slab transitions into the lower mantle, the bottom of the weakened craton may be mobilized and drip into the deeper mantle (Fig. 4b). Though a lower mantle slab anomaly beneath a thick craton is presently only found in North America 34,53 , it could have contributed to previous craton thinning events. This suggests that slabs may not only erode cratons from the sides 4 but craton mobilization and thinning may arise from deep mantle effects of subduction. Declarations Acknowledgements We thank R. W. Clayton, X. Pérez Campos, and C. Cardenas Monroy for acquiring data from the Mexican National Network. We thank M. Wiederspahn for managing the computation resources and data. We thank Z. Zhao and X. Li for the discussions regarding the inversion method. We thank E. Sandvol, C. Sun, D. B. Rowley, and E. K. Heilman for constructive discussions on the interpretation of structures. Mexican National Network data was obtained by the Servicio Sismológico Nacional (México), station maintenance, data acquisition and distribution is thanks to its personnel. J.H. and S.P.G. were partially supported by NSF EAR-1902400 and the Jackson School of Geosciences. T.W.B. was partially supported by NSF EAR-2045292. Author Contributions J.H. conducted the tomography, data analysis, composition conversion, and modeling in the paper. S.P.G. advised on the seismological aspects of the paper. T.W.B. advised on the geodynamical modeling aspects. J.H., S.P.G., and T.W.B. put together the main conclusions of the paper. H.A.J. processed seismic data from the ocean-bottom seismometers. C.L. helped with the inversion approach. D.T.T. helped with the computational environment. H.Z. advised on the inversion approach. Detailed interpretation of the results reflects discussions among the authors. The manuscript was written by J.H. with contributions from S.P.G., T.W.B and other authors. Competing Interests The authors declare no competing interests. Data availability Seismograms from networks other than the Mexican National Network were downloaded from the IRIS Data Management Center (http://ds.iris.edu/ds/nodes/dmc/). Seismograms from the Mexican National Network (https://doi.org/10.21766/SSNMX/SN/MX) were downloaded via the SSNstp client, which is free (http://www2.ssn.unam.mx:8080/getData/SSNdata_UseAndPolicy.pdf) upon request without requirements for authorization. The SATONA model is deposited to the IRIS archive (https://ds.iris.edu/ds/products/emc-earthmodels/), which also hosts US-SL-2014 22 , CAP22 93 , BBNAP19 21 , and NA13 94 that were compared. For other compared models, US22 57 is available on H. Zhu’s website (https://labs.utdallas.edu/seismic-imaging-lab/download/), and models in ref. 23 and ref. 25 are available as a supplement in their publications. All figures and maps are generated by the Generic Mapping Tools 95 . Code availability Computer codes used for data processing, inversion, composition analysis, and plotting are available upon request. SpecFem3D Globe and the geodynamic tool to calculate the flow field are openly available on GitHub (https://github.com/SPECFEM/specfem3d_globe, https://github.com/geodynamics/hc). References Pearson, D. G. et al. Deep continental roots and cratons. Nature 596 , 199-210 (2021). Yuan, H. & Romanowicz, B. Lithospheric layering in the North American craton. Nature 466 , 1063-1068 (2010). Jordan, T. H. Structure and formation of the continental tectosphere. Journal of Petrology , 11-37 (1988). Lenardic, A. & Moresi, L. N. Some thoughts on the stability of cratonic lithosphere: Effects of buoyancy and viscosity. Journal of Geophysical Research: Solid Earth 104 , 12747-12758 (1999). Gao, S., Rudnick, R. L., Carlson, R. W., McDonough, W. F. & Liu, Y.-S. Re–Os evidence for replacement of ancient mantle lithosphere beneath the North China craton. Earth and Planetary Science Letters 198 , 307-322 (2002). Lee, C.-T. A., Luffi, P. & Chin, E. J. Building and destroying continental mantle. Annual Review of Earth and Planetary Sciences 39 , 59-90 (2011). Wu, F.-Y., Yang, J.-H., Xu, Y.-G., Wilde, S. A. & Walker, R. J. Destruction of the North China craton in the Mesozoic. Annual Review of Earth and Planetary Sciences 47 , 173-195 (2019). Foley, S. F. Rejuvenation and erosion of the cratonic lithosphere. Nature geoscience 1 , 503-510 (2008). Liu, L. et al. Mantle dynamics of the North China Craton: new insights from mantle transition zone imaging constrained by P-to-S receiver functions. Geophysical Journal International 231 , 629-637 (2022). Wang, Y. et al. Secular craton evolution due to cyclic deformation of underlying dense mantle lithosphere. Nature Geoscience , 1-9 (2023). Whitmeyer, S. J. & Karlstrom, K. E. Tectonic model for the Proterozoic growth of North America. Geosphere 3 , 220-259 (2007). Tappe, S., Smart, K., Torsvik, T., Massuyeau, M. & de Wit, M. Geodynamics of kimberlites on a cooling Earth: Clues to plate tectonic evolution and deep volatile cycles. Earth and Planetary Science Letters 484 , 1-14 (2018). Tao, K., Grand, S. P. & Niu, F. Full-waveform inversion of triplicated data using a normalized-correlation-coefficient-based misfit function. Geophysical Journal International 210 , 1517-1524 (2017). Komatitsch, D., Erlebacher, G., Göddeke, D. & Michéa, D. High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster. Journal of Computational Physics 229 , 7692-7714 (2010). Tromp, J., Komatitsch, D. & Liu, Q. Spectral-element and adjoint methods in seismology. Communications in Computational Physics 3 , 1-32 (2008). Tao, K., Grand, S. P. & Niu, F. Seismic structure of the upper mantle beneath eastern Asia from full waveform seismic tomography. Geochemistry, Geophysics, Geosystems 19 , 2732-2763 (2018). Bunks, C., Saleck, F. M., Zaleski, S. & Chavent, G. Multiscale seismic waveform inversion. Geophysics 60 , 1457-1473 (1995). Bottou, L. in Proceedings of COMPSTAT'2010: 19th International Conference on Computational StatisticsParis France, August 22-27, 2010 Keynote, Invited and Contributed Papers. 177-186 (Springer). van Herwaarden, D. P. et al. Accelerated full-waveform inversion using dynamic mini-batches. Geophysical Journal International 221 , 1427-1438 (2020). Hasterok, D. et al. New maps of global geological provinces and tectonic plates. Earth-Science Reviews 231 , 104069 (2022). Boyce, A. et al. Variable modification of continental lithosphere during the Proterozoic Grenville orogeny: Evidence from teleseismic P-wave tomography. Earth and Planetary Science Letters 525 , 115763 (2019). Schmandt, B. & Lin, F. C. P and S wave tomography of the mantle beneath the United States. Geophysical Research Letters 41 , 6342-6349 (2014). Sigloch, K. Mantle provinces under North America from multifrequency P wave tomography. Geochemistry, Geophysics, Geosystems 12 (2011). Liang, X., Zhao, D., Hua, Y. & Xu, Y. G. Seismic anisotropy tomography and mantle dynamics of central‐eastern USA. Journal of Geophysical Research: Solid Earth 127 , e2022JB025484 (2022). Wang, H., Zhao, D., Huang, Z. & Wang, L. Tomography, seismotectonics, and mantle dynamics of central and eastern United States. Journal of Geophysical Research: Solid Earth 124 , 8890-8907 (2019). McDonough, W. F. & Sun, S.-S. The composition of the Earth. Chemical geology 120 , 223-253 (1995). Connolly, J. The geodynamic equation of state: what and how. Geochemistry, geophysics, geosystems 10 (2009). Stixrude, L. & Lithgow-Bertelloni, C. Thermal expansivity, heat capacity and bulk modulus of the mantle. Geophysical Journal International 228 , 1119-1149 (2022). Yamauchi, H. & Takei, Y. Polycrystal anelasticity at near‐solidus temperatures. Journal of Geophysical Research: Solid Earth 121 , 7790-7820 (2016). Biryol, C. B., Wagner, L. S., Fischer, K. M. & Hawman, R. B. Relationship between observed upper mantle structures and recent tectonic activity across the Southeastern United States. Journal of Geophysical Research: Solid Earth 121 , 3393-3414, doi:https://doi.org/10.1002/2015JB012698 (2016). Garrity, C. P. & Soller, D. R. Database of the geologic map of North America— Adapted from the map by J.C. Reed, Jr. and others (2005). Report No. 424, (Reston, VA, 2009). Schutt, D. & Lesher, C. Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite. Journal of Geophysical Research: Solid Earth 111 (2006). Hager, B. H. & O'Connell, R. J. A simple global model of plate dynamics and mantle convection. Journal of Geophysical Research: Solid Earth 86 , 4843-4867, doi:https://doi.org/10.1029/JB086iB06p04843 (1981). Lu, C., Grand, S. P., Lai, H. & Garnero, E. J. TX2019slab: A New P and S Tomography Model Incorporating Subducting Slabs. Journal of Geophysical Research: Solid Earth 124 , 11549-11567, doi:https://doi.org/10.1029/2019JB017448 (2019). Forte, A. M., Mitrovica, J. X., Moucha, R., Simmons, N. A. & Grand, S. P. Descent of the ancient Farallon slab drives localized mantle flow below the New Madrid seismic zone. Geophysical Research Letters 34 , doi:https://doi.org/10.1029/2006GL027895 (2007). Straume, E. O., Steinberger, B., Becker, T. W. & Faccenna, C. The impact of mantle convection and dynamic topography on the Cenozoic paleogeography of Central Eurasia and the West Siberian Seaway (accepted). Earth and Planetary Science Letters (2023). Bascou, J. et al. Seismic velocities, anisotropy and deformation in Siberian cratonic mantle: EBSD data on xenoliths from the Udachnaya kimberlite. Earth and Planetary Science Letters 304 , 71-84, doi:https://doi.org/10.1016/j.epsl.2011.01.016 (2011). Kopylova, M. G., Russell, J. K. & Cookenboo, H. Petrology of Peridotite and Pyroxenite Xenoliths from the Jericho Kimberlite: Implications for the Thermal State of the Mantle beneath the Slave Craton, Northern Canada. Journal of Petrology 40 , 79-104, doi:10.1093/petroj/40.1.79 (1999). Abt, D. L. et al. North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions. Journal of Geophysical Research: Solid Earth 115 (2010). Hu, J. et al. Modification of the Western Gondwana craton by plume–lithosphereinteraction. Nature Geoscience 11 , 203-210, doi:10.1038/s41561-018-0064-1 (2018). Liu, L. & Gao, S. S. Lithospheric layering beneath the contiguous United States constrained by S-to-P receiver functions. Earth and Planetary Science Letters 495 , 79-86 (2018). Wang, Z., Kusky, T. M. & Capitanio, F. A. Lithosphere thinning induced by slab penetration into a hydrous mantle transition zone. Geophysical Research Letters 43 , 11,567-511,577 (2016). Komabayashi, T., Omori, S. & Maruyama, S. Petrogenetic grid in the system MgO‐SiO2‐H2O up to 30 GPa, 1600 C: Applications to hydrous peridotite subducting into the Earth’s deep interior. Journal of Geophysical Research: Solid Earth 109 (2004). Richard, G., Bercovici, D. & Karato, S.-I. Slab dehydration in the Earth's mantle transition zone. Earth and Planetary Science Letters 251 , 156-167 (2006). Pokhilenko, N., Agashev, A., Litasov, K. & Pokhilenko, L. Carbonatite metasomatism of peridotite lithospheric mantle: implications for diamond formation and carbonatite-kimberlite magmatism. Russian Geology and Geophysics 56 , 280-295 (2015). Zhang, Y., Wang, C. & Jin, Z. Decarbonation of stagnant slab in the mantle transition zone. Journal of Geophysical Research: Solid Earth 125 , e2020JB019533 (2020). Yang, J. & Faccenda, M. Intraplate volcanism originating from upwelling hydrous mantle transition zone. Nature 579 , 88-91 (2020). Sun, C. & Dasgupta, R. Slab–mantle interaction, carbon transport, and kimberlite generation in the deep upper mantle. Earth and Planetary Science Letters 506 , 38-52 (2019). Chu, R., Leng, W., Helmberger, D. V. & Gurnis, M. Hidden hotspot track beneath the eastern United States. Nature Geoscience 6 , 963-966 (2013). Cox, R. T. & Van Arsdale, R. B. The Mississippi Embayment, North America: a first order continental structure generated by the Cretaceous superplume mantle event. Journal of Geodynamics 34 , 163-176 (2002). Dongmo Wamba, M., Montagner, J.-P. & Romanowicz, B. Imaging deep-mantle plumbing beneath La Réunion and Comores hot spots: Vertical plume conduits and horizontal ponding zones. Science Advances 9 , eade3723, doi:doi:10.1126/sciadv.ade3723 (2023). Paul, J., Conrad, C. P., Becker, T. W. & Ghosh, A. Convective Self‐Compression of Cratons and the Stabilization of Old Lithosphere. Geophysical Research Letters 50 , e2022GL101842 (2023). French, S. W. & Romanowicz, B. Broad plumes rooted at the base of the Earth's mantle beneath major hotspots. Nature 525 , 95-99 (2015). Methods Full-waveform adjoint seismic tomography. Seismograms were obtained for 206 earthquakes occurred between 2004 and 2021, at 6099 individual stations (Extended Data Fig. 1). Magnitudes are mainly between 5.5 and 6.8, to both avoid low-quality data and satisfy the point source assumption which is often not valid for larger earthquakes 16 . All 206 earthquakes are used during the coarse-grid stage which took 126 iterations; among those, 110 earthquakes are selected for the following 90-iteration fine-grid stage (Extended Data Fig. 1). Due to strong noise on horizontal components, only vertical components were used for ocean-bottom seismometers, and both tilt and compliance noise were removed 54 . We used multiple body wave phases including P , S , PP , SS , PPP , SSS and their depth phases, as well as surface waves. For each body wave phase, we cut a waveform segment that starts 30 s before and ends 50 s after the arrival time for waveform fitting. During the coarse-grid stage, body waves were filtered between 17-100 s, and for the fine-grid stage, the filter is 14-100 s. Rayleigh and Love waves were filtered in multiple period bands, including 20-40 s, 40-60 s, 60-80 s, and 80-150 s. Waveform segments were also cut for surface waves to contain their full envelopes. In total, we used 1.5 × 10 6 waveform segments. The inversion scheme is similar to the one in ref. 16 . During the inversion, we tried to maximize the correlation coefficient between synthetic and observed waveforms, which guides the velocity structure to fit phase shapes, while not being strongly affected by less constrained factors like attenuation, and compared to methods based purely on travel times, this method can exploit more information in the detailed shape of waveforms. A preconditioned conjugate gradient method is used to update the model 16 , and a diagonal Hessian is approximated by calculating the sensitivity kernel difference after perturbing the model 55 . Weights for different segments are based on four metrics: correlation coefficients, maximum cross-correlation values, time differences between synthetic and observed arrivals, and the signal-to-noise ratio of observations 16 . In this study, a segment is weighted as unity when these four criteria are met: correlation coefficient > 0.8, maximum cross-correlation value > 0.95, time shift 12; in contrast, it is weighted as zero, if one of these criteria is met: correlation coefficient < 0.55, maximum cross-correlation value 8 s, or signal-to-noise ratio < 8. For metrics in between these thresholds, the weight is determined by the cosine-shape function in ref. 16 . To overcome the cost of wavefield simulations that often limits adjoint tomography to a few tens of iterations 16,55 , and to reduce the possibility of being trapped in local minima, we included the mini-batch algorithm by using only a portion of earthquakes for each iteration. This algorithm is commonly used for stochastic gradient descent in machine learning 18 . Here, the algorithm is similar to the one in ref. 19 , where mini-batches greatly enhance the convergence rate by reducing the redundancy among sensitivity kernels from different earthquakes. However, unlike ref. 19 , the Hessian is not estimated via the L-BFGS algorithm 56 , so we do not need to worry about the positive definiteness of the Hessian during mini-batch selection. Hence, there only remain two key procedures for mini-batch selection: 1) choosing earthquakes to be excluded for the next iteration; 2) adding earthquakes not used in the current iteration. We used a similar algorithm to ref. 19 to select earthquakes to be excluded. The angular distance between the summed kernel using all earthquakes in the current iteration and the summed kernel for all preserved earthquakes with one earthquake removed is calculated, and the earthquake with the smallest angular distance is removed. Such earthquake exclusion is performed iteratively 19 until the angular distance reaches 30º or the total number of excluded stations exceeds 55% of the earthquakes used for the current iteration. The way to add new earthquakes for the next iteration is more stochastic to avoid some earthquakes being constantly chosen. Specifically, we give each candidate earthquake a probability based on: 1) the average distance from a candidate earthquake to its three nearest earthquakes that have secured a spot for the next iteration, 2) the number of iterations that this earthquake was used in. Hence, earthquakes either far from other earthquakes or not frequently used have a better chance of being selected. We used model US22 57 as our starting structure. To account for attenuation, we assume the 3D attenuation model QRFSI12 58 , and while SpecFem3D Globe handles attenuation with a constant Q across frequencies 59 , velocities shown in this study are for 1 s. We assume CRUST1.0 60 for Moho topography, and S362ANI 61 for the topography of the 410- and 660-discontinuities. Accurate earthquake sources are crucial for tomography. In this study, we invert all source parameters 16 within the gCMT solution 62 before inverting for structure, or once the structure inversion has produced substantial updates from the model used for the last source inversion. Sources were updated eight times, and each time, the number of iterations for different earthquakes depends on their convergence rates. Meanwhile, to overcome the uneven distribution of stations, in addition to the original weight based on the four metrics, the weight for a certain station is further normalized by the sum of weights for all stations around it with azimuth difference < 7.5º and epicentral distance difference < 6º. A water level of 2‰ the summed weight of all stations is applied during normalization to not overweight isolated stations. The forward and adjoint simulations for both structural and source inversions were performed using SpecFem3D Globe v8.0.0 14 on a server with 4 Nvidia A100 GPUs. For each earthquake, the duration of simulated wavefields is determined by the end time of the 20-40 s Rayleigh wave segment at the farthest station. The simulation domain has a dimension of 57º × 79º (Extended Data Fig. 1), and during the more influential fine-grid stage, horizontal element size is ~44 km on the free surface, making the spacing between Gauss-Lobatto-Legendre interpolation points ~11 km, and the minimum resolved period 11 s. In this study, we mainly focus on isotropic velocity anomalies. During simulations, we make the structure isotropic below the 410-discontinuity, and transversely isotropic above it. Hence, both isotropic velocity and radial anisotropy were solved for during inversion, and here we used the Voigt average 63 to approximate the isotropic velocity. To compare velocity anomalies at different depths, seismic velocities are shown as velocity perturbations with respect to a reference model. A 1D mantle reference model for V S and V P is obtained by averaging the final velocity model for all regions with good resolution (Extended Data Fig. 6), while in the crust, 3D velocities from CRUST1.0 60 are assumed as the reference. Resolution tests of the model. Two resolution tests are performed separately for V S and V P (Extended Data Fig. 5). In each test, the same point-shape velocity perturbations were added to the final ΔV S /V S and ΔV P /V P models. These perturbations are horizontally 4º and vertically 200 km apart from each other and have reversed polarization between neighbors (Extended Data Fig. 5). The maximum amplitude of the perturbations is 0.01, and their shapes are characterized by a Gaussian with a standard deviation of 40 km. Such spatial extent is comparable with the wavelength of body waves used in this study, making the test challenging, and as set up, preferable to checkerboard tests 64 . During the test, we generated waveforms for the perturbed model as “observations”, and then started from the unperturbed model to invert for the perturbations. Meanwhile, to resemble the resolution of real data, we use the same set of weighting for waveform segments as in the last tomography iteration. Due to computation resources, we did not replicate the 206 iterations for tomography, but 16 iterations were still performed for each test with the mini-batch strategy (Extended Data Fig. 5). Compared to commonly used one-iteration tests for full-waveform tomography 16,55 , this multi-iteration test refines structures at places with lower data coverage to better indicate the actual resolution. A specific test is also designed for the resolution on the high V S dripping bodies (Extended Data Fig. 7). Here, due to computation limits, we use the point-spread function method 16,65 with one iteration. We first made a model with all high V S anomalies around the dripping region removed, and then calculate the V S sensitivity kernel for this modified model. After that, we obtain the difference in preconditioned V S sensitivity kernel between the original model and the modified one, which represents the V S change that the data would favor to compensate for the removed dripping body. Based on results from these tests, the dripping bodies are not likely to be artifacts based on their geometries. Though the general shape of these bodies is sub-vertical, they also contain secondary features that are horizontally aligned (Fig. 2a, b), which cannot be produced by smearing along ray paths, and the resolution test shows no distortion of structures around the region (Extended Data Fig. 5). Also, given the earthquake-station distribution (Extended Data Fig. 1), there are no ray paths following the eastward dripping direction of high V S bodies on the east side (Fig. 2b). Composition and temperature estimation. Because there are two observables (V S and V P ), and temperature is one variable to be constrained, only one compositional variable can be solved for, so we model compositions between two endmembers. Endmember compositions are expressed as oxide weight percentages for the SiO 2 -Al 2 O 3 -MgO-FeO-CaO-Na 2 O system. We take one endmember to be pyrolite following ref. 26 with SiO 2 = 45%, Al 2 O 3 = 4.45%, MgO = 37.8%, FeO = 8.05%, CaO = 3.55%, and Na 2 O = 0.36%. For the endmember composition of a depleted cratonic mantle lithosphere, we compiled whole rock major element compositions for peridotite xenolith samples from four cratons (Extended Data Fig. 9): Greenland 66-68 , Slave 69 , Kaapvaal 70-72 and Siberia 73-75 . Based on these samples, clear negative correlations are seen between the Mg number (Mg# = 100×Mg/(Mg+Fe)) and the FeO content; MgO and SiO 2 content; MgO and Al 2 O 3 content; MgO and CaO content, and all these relationships could be generalized through linear regression (Extended Data Fig. 9). Therefore, by setting the Mg# to 93.7 for the craton endmember, we can get a corresponding FeO from the first relationship, and a MgO content is obtained based on that Mg# and FeO. With the MgO, based on other relationships, other major element contents are determined. The amount of Na 2 O is very small and can be ignored. Together, the craton endmember has SiO 2 = 46.55%, Al 2 O 3 = 1.1%, MgO = 45.99%, FeO = 5.7%, CaO = 0.62%. With these endmember compositions, we estimated the dependence of ΔV S /V S and ΔV P /V P on composition and temperature at different depths. In this study, we performed the conversion between 200 and 380 km depths, because dripping bodies are clearly observed there (Extended Data Fig. 3); the depth is not too deep to be affected by potential inaccuracy from the assumed 410-discontinuity topography; and ΔV P /V P is best resolved at ~300 km depth (Extended Data Fig. 5). To estimate the temperature dependence, we calculated seismic velocities at different temperature and pressure conditions. We assumed an ambient mantle potential temperature of 1350ºC 76 , and an adiabatic temperature gradient of 0.4ºC/km 77 . At each depth, we estimate the elastic V S and V P 50ºC above and below the adiabat using Perple_X 27 with thermodynamic data and solutions in ref. 28 , and pressures from PREM 78 . Then, we estimated the anelasticity effect for these two temperatures at the period of 20 s based on the relationship in ref. 29 through the very-broadband rheology calculator 79 with the solidus from ref. 80 assumed for calculation in ref. 29 . We assume the attenuation effect on bulk modulus to be negligible, and how complex anelastic compliance in ref. 29 would affect V P is derived in the same way as in ref. 81 for V S . The period of 20 s is chosen as it is close to the frequency of body waves used in this study. After that, gradients of V S and V P with respect to temperature are obtained by calculating the differences in anelastic V S and V P between these two temperatures and dividing them by the temperature difference of 100ºC. These gradients are further divided by the anelastic V S and V P at adiabatic temperatures for further conversion. The compositional dependence was estimated similarly. By assuming mechanical mixing 82 , at different depths along the assumed 1350ºC adiabat, we calculated the difference in anelastic V S and V P between the two endmember compositions, and these differences are divided by unity to represent the compositional gradients. These gradients of V S and V P are then divided by the 1D reference V S and V P of SATONA (Extended Data Fig. 6) at these depths for conversion. V S and V P perturbations are then converted to temperature and compositional perturbations. Here, we convert velocity perturbations rather than absolute velocities, because absolute velocities could be strongly influenced by the assumed reference anelasticity model and the exact composition of mantle rocks, while perturbations are less dependent on these reference conditions. Before conversion, to avoid misalignment of independently constrained V S and V P features, ΔV S /V S and ΔV P /V P are smoothed by convolving with a 3D Gaussian, whose horizontal and vertical standard deviations are 120 km and 20 km. Velocities are also converted to 20 s based on the attenuation model 58 for consistency. Velocity perturbations are expressed by temperature perturbation ( Δ T ) and compositional perturbation obtained by solving these two equations (Fig. 2d, e, and Extended Data Fig. 9). When converting for the dripping anomalies (e.g., Fig. 2d, e), only geographic locations within the dripping area with positive ΔV S /V S and ΔV P /V P are considered. The range of converted compositional perturbation appears to be large for the whole model (Extended Data Fig. 9i,m). However, considering these endmember compositions might not be suitable for other areas; the anelasticity model could be more complicated for hot environments 29 such as in the western US and Mexico (Extended Data Fig. 3). Since our goal is mainly to show whether overall V S and V P patterns for dripping bodies could correspond to an increased amount of cratonic material, we focus more on the relative compositional difference between the dripping area and regions excluding the dripping area rather than the exact values. Estimation of thinning. An approximate estimate of the amount of lithospheric thinning follows from assuming the amount of cratonic material is proportional to ΔV S /V S , and all previously dripped lithosphere is still present in the tomographic model as fast anomalies. We first set ΔV S /V S at 180 km as a reference for craton material, since it is within the craton while not too shallow to be strongly influenced by the cold temperature. Then we calculated the sum of positive ΔV S /V S in the dripping area (Fig. 2c) between 250 and 550 km depths and obtained the ratio between them and the reference summed ΔV S /V S at 180 km depth, which represents the volume proportion of cratonic materials at those depths and is on average ~20%. Hence, the thinning would be (550 km – 250 km) × 20% = 60 km if dripping locally. Then, since horizontal flow induced by the slab (Fig. 3a versus c) could focus material from the whole cratonic region covered by this model (Fig. 1) with an area of ~4×10 6 km 2 , and the dripping area only has an area of ~1.5×10 6 km 2 (~38% of the cratonic area), the average amount of thinning is calculated to be 60 km × 38% ~ 23 km. However, some regions may experience more thinning than others, for instance, flanks at the western edge of the craton in Fig. 2a, b could partly be caused by the thinning. Meanwhile, although we have estimated compositional perturbation, and it acts as a powerful way to show the significant distinction in composition, we did not use it for the thinning estimation as what is obtained there are perturbations with large uncertainties and are also subject to many assumptions. Geodynamic modeling of mantle convection. Mantle flow fields driven by different density structures were modeled. The calculation is based on ref. 33,83 up to spherical harmonic degree 63. The 1D viscosity distribution from ref. 84 and the surface MORVEL 85 plate motion are assumed, and we only consider vertical variations in viscosity, for simplicity, which leads us to trust patterns more than amplitudes of flow, e.g. ref. 86 . We first estimated flow fields from two density models based on the global tomographic model TX2019slab 34 . Starting from its V S model, we removed all velocity perturbations around the dripping area at 280-600 km depths (Fig. 3b,d) to eliminate local buoyancy flows. Then for one model, at each depth slice below 600 km, we outlined the geometry of the high-velocity Farallon slab and removed all perturbations within it (Fig. 3d). After that, two V S models with or without the Farallon slab are converted to density based on a V S -density relationship 84 , and their corresponding flow fields were calculated (Fig. 3). To evaluate whether the downward flow induced by the Farallon slab is strong enough to drag down the dripping materials, we calculated flow fields after merging SATONA with TX2019slab. For SATONA, we only used reliable regions based on resolution tests (Extended Data Fig. 5). For locations that are >200 km within the resolution boundary (Extended Data Fig. 5), ΔV S /V S in model SATONA is used, and for locations >200 km outside the boundary, TX2019slab 34 is used. For locations within 200 km of the boundary, we weighted averaged ΔV S /V S between these two models to make the merged model continuous at ± 200 km from the boundary. A similar merging strategy is designed for the bottom of SATONA, with depths 900 km, and in between, weighted averaged values are used. After smoothing, ΔV S /V S of the dripping area is about 1% (Extended Data Fig. 9e), in agreement with some longer wavelength previous models (Extended Data Fig. 7). Then, to understand how dripping body densities could affect flows, we tested four scenarios by making the V S -density scaling factor at 260-600 km depths for positive V S anomalies within the dripping area to be 1, 0.5, 0, -0.5 times the factor for the rest of the region 84 (Extended Data Fig. 10), which in the last scenario, makes dripping bodies positively buoyant (Extended Data Fig. 10d). We also used past slab locations to help our interpretation (Fig. 1). Mantle structures in history were predicted in ref. 36 by advecting velocity structures in TX2019slab 34 back in time. Because the present-day Farallon slab does not always appear above the 660-discontinuity, we chose to use the contour of 0.12% density anomaly at 800 km depth to outline the slab for those paleo-density models. Then, to show the relative location of the Farallon slab to the North American plate, these outlines were translated to their present location based on the assumed plate motion 87 in ref. 36 . These past slab locations were analyzed with dated magmatism (Fig. 1) in Kansas 88 , Arkansas 89 , Monroe uplift 90 , Kentucky 91 , and Virginia 92 . Methods References 54 Janiszewski, H. A., Gaherty, J. B., Abers, G. A., Gao, H. & Eilon, Z. C. Amphibious surface-wave phase-velocity measurements of the Cascadia subduction zone. Geophysical Journal International 217 , 1929-1948 (2019). 55 Zhu, H., Bozdağ, E. & Tromp, J. Seismic structure of the European upper mantle based on adjoint tomography. Geophysical Journal International 201 , 18-52 (2015). 56 Liu, D. C. & Nocedal, J. On the limited memory BFGS method for large scale optimization. Mathematical programming 45 , 503-528 (1989). 57 Zhu, H., Komatitsch, D. & Tromp, J. Radial anisotropy of the North American upper mantle based on adjoint tomography with USArray. Geophysical Journal International 211 , 349-377 (2017). 58 Dalton, C. A., Ekström, G. & Dziewoński, A. M. The global attenuation structure of the upper mantle. Journal of Geophysical Research: Solid Earth 113 (2008). 59 Komatitsch, D. & Tromp, J. Spectral-element simulations of global seismic wave propagation—I. Validation. Geophysical Journal International 149 , 390-412 (2002). 60 Laske, G., Masters, G., Ma, Z. & Pasyanos, M. in Geophysical research abstracts. 2658 (EGU General Assembly 2013, Vienna, Austria). 61 Kustowski, B., Ekström, G. & Dziewoński, A. Anisotropic shear‐wave velocity structure of the Earth's mantle: A global model. Journal of Geophysical Research: Solid Earth 113 (2008). 62 Ekström, G., Nettles, M. & Dziewoński, A. The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. Physics of the Earth and Planetary Interiors 200 , 1-9 (2012). 63 Panning, M. & Romanowicz, B. A three-dimensional radially anisotropic model of shear velocity in the whole mantle. Geophysical Journal International 167 , 361-379 (2006). 64 Rawlinson, N. & Spakman, W. On the use of sensitivity tests in seismic tomography. Geophysical Journal International 205 , 1221-1243 (2016). 65 Fichtner, A. & Leeuwen, T. v. Resolution analysis by random probing. Journal of Geophysical Research: Solid Earth 120 , 5549-5573 (2015). 66 Bernstein, S., Kelemen, P. B. & Brooks, C. K. Depleted spinel harzburgite xenoliths in Tertiary dykes from East Greenland: restites from high degree melting. Earth and Planetary Science Letters 154 , 221-235 (1998). 67 Hanghøj, K., Kelemen, P., Bernstein, S., Blusztajn, J. & Frei, R. Osmium isotopes in the Wiedemann Fjord mantle xenoliths: a unique record of cratonic mantle formation by melt depletion in the Archaean. Geochemistry, Geophysics, Geosystems 2 (2001). 68 Wittig, N. et al. Origin of cratonic lithospheric mantle roots: A geochemical study of peridotites from the North Atlantic Craton, West Greenland. Earth and Planetary Science Letters 274 , 24-33 (2008). 69 Liu, J. et al. Plume-driven recratonization of deep continental lithospheric mantle. Nature 592 , 732-736 (2021). 70 Griffin, W., Graham, S., O'Reilly, S. Y. & Pearson, N. Lithosphere evolution beneath the Kaapvaal Craton: Re–Os systematics of sulfides in mantle-derived peridotites. Chemical Geology 208 , 89-118 (2004). 71 Simon, N. S., Carlson, R. W., Pearson, D. G. & Davies, G. R. The origin and evolution of the Kaapvaal cratonic lithospheric mantle. Journal of Petrology 48 , 589-625 (2007). 72 Simon, N. S., Irvine, G. J., Davies, G. R., Pearson, D. G. & Carlson, R. W. The origin of garnet and clinopyroxene in “depleted” Kaapvaal peridotites. Lithos 71 , 289-322 (2003). 73 Boyd, F. et al. Composition of the Siberian cratonic mantle: evidence from Udachnaya peridotite xenoliths. Contributions to Mineralogy and Petrology 128 , 228-246 (1997). 74 Ionov, D. A., Doucet, L. S. & Ashchepkov, I. V. Composition of the lithospheric mantle in the Siberian craton: new constraints from fresh peridotites in the Udachnaya-East kimberlite. Journal of petrology 51 , 2177-2210 (2010). 75 Ionov, D. A., Doucet, L. S., Carlson, R. W., Golovin, A. V. & Korsakov, A. V. Post-Archean formation of the lithospheric mantle in the central Siberian craton: Re–Os and PGE study of peridotite xenoliths from the Udachnaya kimberlite. Geochimica et Cosmochimica Acta 165 , 466-483 (2015). 76 Herzberg, C. & Gazel, E. Petrological evidence for secular cooling in mantle plumes. Nature 458 , 619-622 (2009). 77 Katsura, T. A revised adiabatic temperature profile for the mantle. Journal of Geophysical Research: Solid Earth 127 , e2021JB023562 (2022). 78 Dziewonski, A. M. & Anderson, D. L. Preliminary reference Earth model. Physics of the earth and planetary interiors 25 , 297-356 (1981). 79 Havlin, C., Holtzman, B. K. & Hopper, E. Inference of thermodynamic state in the asthenosphere from anelastic properties, with applications to North American upper mantle. Physics of the Earth and Planetary Interiors 314 , 106639 (2021). 80 Hirschmann, M. M. Mantle solidus: Experimental constraints and the effects of peridotite composition. Geochemistry, Geophysics, Geosystems 1 (2000). 81 McCarthy, C., Takei, Y. & Hiraga, T. Experimental study of attenuation and dispersion over a broad frequency range: 2. The universal scaling of polycrystalline materials. Journal of Geophysical Research: Solid Earth 116 (2011). 82 Xu, W., Lithgow-Bertelloni, C., Stixrude, L. & Ritsema, J. The effect of bulk composition and temperature on mantle seismic structure. Earth and Planetary Science Letters 275 , 70-79 (2008). 83 HC, a global mantle circulation solver (2009). 84 Steinberger, B. & Calderwood, A. R. Models of large-scale viscous flow in the Earth's mantle with constraints from mineral physics and surface observations. Geophysical Journal International 167 , 1461-1481 (2006). 85 DeMets, C., Gordon, R. G. & Argus, D. F. Geologically current plate motions. Geophysical journal international 181 , 1-80 (2010). 86 Ghosh, A., Becker, T. & Humphreys, E. Dynamics of the North American continent. Geophysical Journal International 194, 651-669 (2013). 87 Torsvik, T. H. et al. Pacific‐Panthalassic reconstructions: Overview, errata and the way forward. Geochemistry, Geophysics, Geosystems 20 , 3659-3689 (2019). 88 Blackburn, T. J., Stockli, D. F., Carlson, R. W. & Berendsen, P. (U–Th)/He dating of kimberlites—A case study from north-eastern Kansas. Earth and Planetary Science Letters 275 , 111-120 (2008). 89 Eby, G. N. & Vasconcelos, P. Geochronology of the Arkansas alkaline province, southeastern United States. The Journal of Geology 117 , 615-626 (2009). 90 Baksi, A. K. The timing of Late Cretaceous alkalic igneous activity in the northern Gulf of Mexico basin, southeastern USA. The Journal of Geology 105 , 629-644 (1997). 91 Heaman, L. M., Kjarsgaard, B. A. & Creaser, R. A. The temporal evolution of North American kimberlites. Lithos 76 , 377-397 (2004). 92 Zartman, R. E. Geochronology of some alkalic rock provinces in eastern and central United States. Annual Review of Earth and Planetary Sciences 5 , 257-286 (1977). 93 Boyce, A. et al. A New P‐Wave Tomographic Model (CAP22) for North America: Implications for the Subduction and Cratonic Metasomatic Modification History of Western Canada and Alaska. Journal of Geophysical Research: Solid Earth 128 , e2022JB025745 (2023). 94 Bedle, H., Lou, X. & van der Lee, S. High-resolution imaging of continental tectonics in the mantle beneath the United States, through the combination of USArray data sets. Authorea Preprints (2022). 95 Wessel, P. et al. The generic mapping tools version 6. Geochemistry, Geophysics, Geosystems 20 , 5556-5564 (2019). 96 Kennett, B. L., Engdahl, E. & Buland, R. Constraints on seismic velocities in the Earth from traveltimes. Geophysical Journal International 122 , 108-124 (1995). Additional Declarations There is NO Competing Interest. Supplementary Files ExtendedData.docx Cite Share Download PDF Status: Published Journal Publication published 28 Mar, 2025 Read the published version in Nature Geoscience → Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-3254038","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":273858570,"identity":"b6f52861-5129-4168-815f-4cf6540813b7","order_by":0,"name":"Junlin Hua","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA1klEQVRIiWNgGAWjYBACxhlAIqFCgocfxEsoIFbLhzMWcpINIC0GxFgjAdQ2s6XC2OAAiEeMFubZPYafeRskEjefX5344YEBgzy/2AECDptzxliad4dE4rYbbzdLAB1mOHN2AgEtM3IMpHnPgLSc3QDSkmBwm7AW49+8bUCHzTi7+QexWswkZ7ZJGBvw924j1pa0MosPZyTkJG7wbrNIMJAg7BfDGcmbbyRU1PHw95/dfPNHhY08vzQhLQ0wlgRYpQR+5SAgD2fxHyCsehSMglEwCkYmAACjEUZ751tiBAAAAABJRU5ErkJggg==","orcid":"https://orcid.org/0000-0002-6137-4319","institution":"The University of Texas at Austin","correspondingAuthor":true,"prefix":"","firstName":"Junlin","middleName":"","lastName":"Hua","suffix":""},{"id":273858571,"identity":"3f76cba7-cc0f-4c66-a9fb-1f75d55340b6","order_by":1,"name":"Steve Grand","email":"","orcid":"","institution":"University of Texas","correspondingAuthor":false,"prefix":"","firstName":"Steve","middleName":"","lastName":"Grand","suffix":""},{"id":273858572,"identity":"0faaa7a4-4413-4e62-a565-1738c4e2ba05","order_by":2,"name":"Thorsten Becker","email":"","orcid":"https://orcid.org/0000-0002-5656-4564","institution":"Uni Texas Austin","correspondingAuthor":false,"prefix":"","firstName":"Thorsten","middleName":"","lastName":"Becker","suffix":""},{"id":273858573,"identity":"7506caef-05db-401f-b4c5-fda88ac09082","order_by":3,"name":"Helen Janiszewski","email":"","orcid":"","institution":"University of Hawaii","correspondingAuthor":false,"prefix":"","firstName":"Helen","middleName":"","lastName":"Janiszewski","suffix":""},{"id":273858574,"identity":"cc9caada-13d0-4bb8-a637-688f218fee9c","order_by":4,"name":"Chujie Liu","email":"","orcid":"","institution":"The University of Texas at Austin","correspondingAuthor":false,"prefix":"","firstName":"Chujie","middleName":"","lastName":"Liu","suffix":""},{"id":273858575,"identity":"6340da85-fe43-4144-9d99-e0c3236a15dc","order_by":5,"name":"Daniel Trugman","email":"","orcid":"","institution":"University of Nevada, Reno","correspondingAuthor":false,"prefix":"","firstName":"Daniel","middleName":"","lastName":"Trugman","suffix":""},{"id":273858576,"identity":"f14b42ad-baf7-4b03-86d3-ea3484247bea","order_by":6,"name":"Hejun Zhu","email":"","orcid":"https://orcid.org/0000-0002-7452-075X","institution":"The University of Texas at Dallas","correspondingAuthor":false,"prefix":"","firstName":"Hejun","middleName":"","lastName":"Zhu","suffix":""}],"badges":[],"createdAt":"2023-08-11 03:40:48","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3254038/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3254038/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1038/s41561-025-01671-x","type":"published","date":"2025-03-28T04:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":52483511,"identity":"ffea4745-71e1-4a5c-8481-47f2511f915d","added_by":"auto","created_at":"2024-03-12 06:58:47","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1766279,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIsotropic V\u003c/strong\u003e\u003csub\u003e\u003cstrong\u003eS\u003c/strong\u003e\u003c/sub\u003e\u003cstrong\u003e distribution at 200 km depth, revealing the craton root.\u003c/strong\u003e Background color shows the V\u003csub\u003eS\u003c/sub\u003e perturbation. The craton root is characterized by its high V\u003csub\u003eS\u003c/sub\u003e and outlined by the black dashed line. Low V\u003csub\u003eS\u003c/sub\u003e anomalies along the east coast are outlined by the light brown dashed line. Red triangles show past magmatism around the craton with their ages labeled. Solid lines with different colors show surface projections of current and past locations of the Farallon slab at 800 km depth with respect to the North American plate\u003csup\u003e36\u003c/sup\u003e.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-3254038/v1/ef3bb55f5654ca6e57c539be.png"},{"id":52483512,"identity":"4d74c0ab-b554-4b05-9309-d66e2e65c183","added_by":"auto","created_at":"2024-03-12 06:58:47","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1224297,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe dripping cratonic lithosphere. a \u003c/strong\u003eand\u003cstrong\u003e b.\u003c/strong\u003e V\u003csub\u003eS\u003c/sub\u003e perturbation profiles across High velocity dripping anomalies at different latitudes. Locations of these profiles are shown in \u003cstrong\u003ec\u003c/strong\u003e. Green dots at the top are geographic markers. Important features in the cross-sections are labeled. \u003cstrong\u003ec. \u003c/strong\u003eMap for the average V\u003csub\u003eS\u003c/sub\u003e perturbation between 300 and 500 km depths. The gray dashed line outlines the craton root as in Fig. 1, and the black dashed line outlines where the dripping high V\u003csub\u003eS\u003c/sub\u003e body is evident at this depth range. White lines show locations for profiles in \u003cstrong\u003ea \u003c/strong\u003eand\u003cstrong\u003e b\u003c/strong\u003e, and green geographic markers correspond to those at the top of the cross-sections. \u003cstrong\u003ed.\u003c/strong\u003e Histogram for the converted compositional perturbation distribution of dripping bodies (within black dashed line in \u003cstrong\u003ec\u003c/strong\u003e) at 350 km depth. Each sample represents one geographic location. Zero in this plot represents the average composition for all regions with good resolution, and positive or negative values mean more cratonic or more pyrolytic, respectively. The red solid line shows the median value, while the dashed lines show 25% and 75% percentiles. \u003cstrong\u003ee. \u003c/strong\u003eLike \u003cstrong\u003ed\u003c/strong\u003e, but for the temperature perturbation.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-3254038/v1/a89a3e9640e568adbc8a9971.png"},{"id":52483513,"identity":"b4cfb358-aa15-4da1-8049-e577721d8e09","added_by":"auto","created_at":"2024-03-12 06:58:47","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1170402,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMantle flow induced by the subducting Farallon slab.\u003c/strong\u003e \u003cstrong\u003ea.\u003c/strong\u003e The radial flow speed between 300 and 500 km depths when the Farallon slab is present. Colors show radial flow, so a negative value represents flow to deeper depths. Arrows show the average horizontal flow fields at the same depth range. The white line shows the location of the cross-section in \u003cstrong\u003ec\u003c/strong\u003e and \u003cstrong\u003ed\u003c/strong\u003e, and green geographic markers correspond to those at the top of the cross-section. The black dashed line is the same as in Fig. 2c. \u003cstrong\u003eb.\u003c/strong\u003e Cross-section for the input density anomaly that drives the mantle flow. Black arrows show the projection of the modeled 3D flow field on the cross-section. \u003cstrong\u003ec\u003c/strong\u003e and \u003cstrong\u003ed. \u003c/strong\u003eLike \u003cstrong\u003ea\u003c/strong\u003e and \u003cstrong\u003eb\u003c/strong\u003e, but based on an input model without the Farallon slab.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-3254038/v1/3085564cda1d10732c00b464.png"},{"id":52483514,"identity":"92a3041e-ebf2-4c67-a20f-1be6bba2924d","added_by":"auto","created_at":"2024-03-12 06:58:47","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":942611,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSchematics for craton mobilization and thinning.\u003c/strong\u003e \u003cstrong\u003ea.\u003c/strong\u003eInferred process of successive, slab induced craton modification. The intrinsic layering of cratonic lithosphere is depicted. The sequence may involve the deep slab motion towards the craton; volatile enrichment in the MTZ; hydration of the asthenosphere; partial melt ponded beneath the craton; magmatism at craton edges; and potential plume-originated materials spreading beneath the MTZ that could later ascend. \u003cstrong\u003eb.\u003c/strong\u003e Lithospheric dripping mobilized by the deep slab. Gray arrows mark large-scale flow lines induced by the slab. The lithosphere drips and mobilized lithospheric bottoms are also depicted.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-3254038/v1/8205e92f681d4b1e99478ba1.png"},{"id":79482155,"identity":"a226ea04-ed48-41b6-8529-2e2346e127f3","added_by":"auto","created_at":"2025-03-29 07:11:23","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":7539663,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3254038/v1/9ba3e020-78d3-450a-9455-e504942b7ed7.pdf"},{"id":52483516,"identity":"f6f753ef-549b-48be-abe9-62bf9f3873a2","added_by":"auto","created_at":"2024-03-12 06:58:49","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":30664866,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"ExtendedData.docx","url":"https://assets-eu.researchsquare.com/files/rs-3254038/v1/98b81378c927dd588d98d0b3.docx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Mobilization and thinning of cratonic lithosphere by a lower mantle slab","fulltext":[{"header":"Main","content":"\u003cp\u003eCratons are old (no younger than Proterozoic) regions which make up ~60% of continental lithosphere\u003csup\u003e1\u003c/sup\u003e and are often underlain by thick lithosphere which can extend to over 200 km depth\u003csup\u003e1,2\u003c/sup\u003e. Cratons are believed to have been overall stable over billions of years due to their near neutral effective buoyancy and high viscosity\u003csup\u003e3,4\u003c/sup\u003e. However, the North China Craton lost its root during the Mesozoic\u003csup\u003e5\u003c/sup\u003e, for example, indicating that the stability of cratons is not universal. Different mechanisms have been proposed for the loss of cratonic lithosphere, including convective removal\u003csup\u003e6\u003c/sup\u003e, flat-slab rollback\u003csup\u003e7\u003c/sup\u003e, and erosion with melting\u003csup\u003e8,9\u003c/sup\u003e. Other studies suggest the cratonic lithosphere is not stable through history but subject to reworking\u003csup\u003e10\u003c/sup\u003e. Given that most of these examples are from processes in the geological past, it is not straightforward to distinguish between different mechanisms, and new observations of an ongoing craton-thinning event would offer a unique perspective to the puzzle.\u003c/p\u003e\n\u003cp\u003eNorth America is an ideal place to study such present-day behavior of cratons. Core continental North America consists of multiple Archean cratons and volcanic arcs accreted during the Proterozoic. Although most cratonic regions in the U.S. were stabilized during the Proterozoic from 1.35-1.8 Ga\u003csup\u003e11\u003c/sup\u003e, they still share key features with their Archean counterparts, including very thick lithosphere\u003csup\u003e2\u003c/sup\u003e and surrounding kimberlite volcanism\u003csup\u003e12\u003c/sup\u003e. Meanwhile, thanks to dense seismometer deployments, the Central and Eastern United States cratons can be studied in detail seismically.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eIn this study, we obtain a full-waveform tomography model for North America (Fig. 1), named Seismic Adjoint Tomography of North America (SATONA). Compared with existing North America tomographic models, SATONA focuses more on fitting the detailed shape of waveforms at relatively high frequencies to better reveal fine-scale structures. SATONA reveals clear lithosphere-asthenosphere boundaries (LAB) for the craton, and high seismic velocity materials beneath the craton from the LAB to the mantle transition zone (MTZ). These high-velocity, drip-shape features indicate thermo-chemical differences relative to the surrounding mantle and likely originate from the cratonic lithosphere. To find the cause for such dripping, we test the influence of the lower mantle Farallon slab on upper mantle convective flow through geodynamic modeling, and find that the dripping may be a consequence of the slab sinking. These findings lead us to propose a new mechanism for craton thinning \u0026mdash;large-scale downwelling induced by a deep-seated slab that mobilizes the craton, possibly mediated by fluid release, leading to partial removal of the cratonic lithosphere.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDripping cratonic lithosphere\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eModel SATONA is obtained through full-waveform adjoint tomography. During inversion, we maximize the correlation coefficient between synthetic and observed waveforms\u003csup\u003e13\u003c/sup\u003e. Hence, besides arrival times of seismic phases, the detailed shape of phases including secondary arrivals are required to match the data. In this study, both synthetic seismograms and adjoint kernels were calculated with SpecFem3D Globe\u003csup\u003e14,15\u003c/sup\u003e, and all available stations in the study region were utilized (Extended Data Fig. 1). We used multiple body and surface wave phases to constrain S- and P-wave velocities\u003csup\u003e16\u003c/sup\u003e. To optimize computational efficiency and avoid local minima, the inversion consists of two stages\u003csup\u003e17\u003c/sup\u003e: a coarse-grid stage using 206 earthquakes, and a fine-grid stage using 110 earthquakes with higher quality (Extended Data Fig. 1). We incorporate the mini-batch algorithm\u003csup\u003e18,19\u003c/sup\u003e by using a subset of earthquakes in each iteration, which, compared with common full-waveform tomography practices using only a few tens of iterations on supercomputers\u003csup\u003e16\u003c/sup\u003e, enables a total of 216 iterations with 4 GPUs.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eWith these, the improvement on\u0026nbsp;the correlation coefficient and the arrival time difference between synthetic and observed waveforms\u0026nbsp;is significant (Extended Data Fig. 2).\u0026nbsp;Body waves show higher accuracy in timing than surface waves, while surface waves have the correlation better recovered. Nonetheless, the vast majority of both wave types show correlation coefficients over 0.9, suggesting an excellent fitting of the waveform shape.\u003c/p\u003e\n\u003cp\u003eDetailed mantle structures are newly imaged in SATONA. S-wave velocity (V\u003csub\u003eS\u003c/sub\u003e, Extended Data Fig. 3) at 100 km depth correlates well with tectonic regions\u003csup\u003e20\u003c/sup\u003e, and at greater depths, the model captures features such as the Cascadia slab in fine detail. P-wave velocity (V\u003csub\u003eP\u003c/sub\u003e, Extended Data Fig. 4), which is independently inverted for, shows similar patterns. A resolution test with only 16 iterations (versus 206 for the full model) shows that V\u003csub\u003eS\u003c/sub\u003e has good resolution down to 1000 km, especially within the upper mantle (Extended Data Fig. 5). V\u003csub\u003eP\u003c/sub\u003e also has good resolution down to the MTZ, but deeper resolution is reduced (Extended Data Fig. 5). We only discuss regions with reliable resolution.\u0026nbsp;The velocity perturbation model also shows some horizontally aligned anomalies near reference model discontinuities (e.g., near 660 km depth in Fig. 2a, b). However, as discontinuity depths are taken from previous models upon initialization of the inversion, these anomalies might be introduced to reconcile inaccuracies about discontinuity depth in those starting models, and are therefore not interpreted much in this study.\u003c/p\u003e\n\u003cp\u003eCratonic lithosphere is well-imaged in the model. At 200 km depth, a broad region beneath the Central and Eastern U.S. shows high V\u003csub\u003eS\u003c/sub\u003e (Fig. 1), marking the geographic extent of the deep cratonic root. Cratonic LABs are clearly characterized as a velocity decrease with depth around 200 km in both V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e (Figs. 2a,b and Extended Data Fig. 6).\u003c/p\u003e\n\u003cp\u003eBeneath the cratonic LAB, extensive high V\u003csub\u003eS\u003c/sub\u003e anomalies are observed down to the MTZ (Fig. 2, Extended Data Fig. 3). A specific model resolution test (Extended Data Fig. 7) suggests these anomalies can be well imaged with our dataset and inversion strategy. Similar high velocity anomalies are also found in previous models\u003csup\u003e21-25\u003c/sup\u003e (Extended Data Fig. 8). However, without the help of multiple seismic phases, their cratonic lithosphere is not well constrained (Extended Data Fig. 8), and the deep structures sometimes extend into the lithosphere depth range, making them harder to interpret.\u003c/p\u003e\n\u003cp\u003eTo distinguish whether these high V\u003csub\u003eS\u003c/sub\u003e anomalies originate from the cratonic lithosphere, we convert V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e perturbations of the dripping bodies to temperature and compositional perturbations. Compared with V\u003csub\u003eS\u003c/sub\u003e, V\u003csub\u003eP\u003c/sub\u003e shows much weaker anomalies at 300-500 km depths (Extended Data Fig. 6), distinct from the strong V\u003csub\u003eP\u003c/sub\u003e anomaly in the Cascadia slab, suggesting a potentially different composition. To constrain the compositional perturbation, two endmembers are considered. One composition is assumed to be pyrolite\u003csup\u003e26\u003c/sup\u003e and the other to be cratonic lithosphere. We estimate a representative major element composition for the craton from kimberlite xenoliths (Extended Data Fig. 9). The V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e for these two endmember compositions at different temperatures and pressures are then calculated using Perple_X\u003csup\u003e27,28\u003c/sup\u003e, accounting for anelastic effects\u003csup\u003e29\u003c/sup\u003e. Derivatives of V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e with respect to temperature and composition are obtained, and velocity anomalies are then converted to temperature and compositional perturbations. We performed the conversion for all geographic locations based on a smoothed version of SATONA (Extended Data Fig. 9) to eliminate bias due to the misalignment of V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e anomalies.\u003c/p\u003e\n\u003cp\u003eTaking the compositional perturbation between the endmembers as unity, at all depths, the high V\u003csub\u003eS\u003c/sub\u003e anomalies show 20-50% more cratonic composition than the ambient mantle (Extended Data Fig. 9), which indicates the entrainment of lithospheric material. Especially below 300 km, due to the phase transition of orthopyroxene at ~320 km\u003csup\u003e28\u003c/sup\u003e, the dependence of V\u003csub\u003eS\u003c/sub\u003e-V\u003csub\u003eP\u003c/sub\u003e discrepancy on composition is significant, making the compositional perturbation better constrained. At 350 km depth, the median of the dripping high V\u003csub\u003eS\u003c/sub\u003e bodies has 45% more cratonic composition than normal (Fig. 2d), which is close to the 75% percentile of compositional perturbations for regions excluding the dripping area (Extended Data Fig. 9m), indicating the high V\u003csub\u003eS\u003c/sub\u003e anomalies are substantially more cratonic than normal. Despite having higher V\u003csub\u003eS\u003c/sub\u003e, when converted to temperature perturbation, we found these anomalies are only ~60\u0026ordm;C colder than normal (Fig. 2e). Even for less constrained depths, the thermal anomaly is less than 200\u0026ordm;C (Extended Data Fig. 9j). This test demonstrates these compositionally distinct anomalies could be lithospheric drips, and they are hereafter referred to as dripping anomalies.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDeep-seated slab induced dripping\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eSimilarly imaged dripping anomalies in the region have previously been interpreted as remnant slabs\u003csup\u003e22\u003c/sup\u003e or dripping due to Rayleigh-Taylor instability\u003csup\u003e30\u003c/sup\u003e, but both mechanisms have some limitations. For slab relics, relative values of\u0026nbsp;V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e anomalies do not match those of the Cascadia slab whose V\u003csub\u003eP\u003c/sub\u003e anomalies are much stronger (Fig. 2c and Extended Data Fig. 6). There is also no evidence of arc volcanism in the study region since at least the Mesozoic\u003csup\u003e31\u003c/sup\u003e, which is a long time for stalling thermal anomalies. In particular, large scale downward flow likely occurred in this region for 70 Ma and thus it is not likely for slab relics to persist. If the anomalies are purely from a\u0026nbsp;Rayleigh-Taylor instability\u0026nbsp;due to a dense craton root\u003csup\u003e30\u003c/sup\u003e, it is hard to explain the relatively weak V\u003csub\u003eP\u003c/sub\u003e anomaly for the drips (Extended Data Fig. 6), since with the inclusion of a large percentage of dense minerals like garnet, V\u003csub\u003eP\u003c/sub\u003e would also be increased\u003csup\u003e32\u003c/sup\u003e. Meanwhile, if a uniform dense root is common among cratons, similar dripping would likely be observed beneath every craton, and likely have occurred for over a billion years and hence might have consumed all dense material.\u003c/p\u003e\n\u003cp\u003eWe propose a new mechanism that the deep Farallon slab is the main driver of dripping, though it is in the lower mantle and not connected to the surface at present. To evaluate the dynamic influence of the Farallon slab, we calculated mantle flow\u003csup\u003e33\u003c/sup\u003e for two different density structures (Fig. 3): mantle with or without the Farallon slab. These structures are modified from a global mantle tomography model\u003csup\u003e34\u003c/sup\u003e. Results show that large-scale mantle flow is strongly controlled by the Farallon slab, which pulls shallow mantle from the East, West, and North to the dripping area where these mantle materials flow downward (Fig. 3a), in line with previous work\u003csup\u003e35\u003c/sup\u003e. Such flow is only predicted when the Farallon slab is present (Fig. 3a, b), and is absent without it (Fig. 3c, d). The Farallon slab may induce horizontal mantle flow over a domain wider than ~3000 km from the sinker (Fig. 3a versus 3c), and associated shear flow may serve to thin and entrain the base of the craton (Fig. 3b). That is, even though the imaged dripping area has relatively limited extent as in Fig. 2c, the material involved in the sinker may come from a much wider region. Such basal erosion by slab-induced flow leading to dripping removal of continental lithosphere may be assisted by prior weakening of the lithosphere, e.g. due to volatile influx from the slab. Associated slab-induced dripping could have started at ~70 Ma, as the Farallon slab has stayed at a similar relative location with respect to North America since then\u003csup\u003e36\u003c/sup\u003e (Fig. 1).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eBased on the tomographic model, by assuming the\u0026nbsp;V\u003csub\u003eS\u003c/sub\u003e anomaly is proportional to the amount of cratonic materials, with the consideration of horizontal flow concentrating material\u0026nbsp;from the whole imaged cratonic region (Fig. 1), the total thinning associated with the currently imaged anomaly may be of order ~20 km.\u0026nbsp;This number might be an overestimate if the slab-induced horizontal flow could influence an even larger area (Fig.\u0026nbsp;3), for example, than the cratonic region defined in Fig. 1 that is covered by this model, and the edge of the craton (Fig. 2) at shallower depth could extend wider (Extended Data Fig. 3b) than the cratonic region at 200 km here. Additionally, we chose to use the stronger V\u003csub\u003eS\u003c/sub\u003e anomaly\u0026nbsp;instead of\u0026nbsp;V\u003csub\u003eP\u003c/sub\u003e anomaly\u0026nbsp;for this estimation.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eBecause the dripping material appears more cratonic than mantle in composition, its effective negative buoyancy may be reduced\u003csup\u003e3\u003c/sup\u003e. To evaluate whether the flow caused by the Farallon slab is strong enough to pull down these materials, we merged our\u0026nbsp;V\u003csub\u003eS\u003c/sub\u003e anomalies from SATONA with global tomography\u003csup\u003e34\u003c/sup\u003e to estimate the effect of different density anomalies (Extended Data Fig. 10). To account for the range of potential densities of the dripping body, we tested different scaling factors between V\u003csub\u003eS\u003c/sub\u003e and density (Extended Data Fig. 10), and in one case we make density negatively correlated with velocity for the dripping material, so that they are positively buoyant. While the sinking speeds are expectedly reduced for this case, the mantle still flows downward at the dripping locations (Extended Data Fig. 10d). This is partly because strong horizontal inflow driven by the Farallon slab always converges at the dripping location. This test shows that the downward dragging force by the slab could overcome any neutral and even small positive buoyancy of lithospheric material once mobilized.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompound process of craton thinning\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;To enable dripping of cratonic materials, the base of the craton cannot be very rigid, so that downward tractions applied by the deep slab and shear at the LAB could indeed mobilize the lithosphere and entrain material. Such a relatively weak lower part of the craton has previously been inferred. For instance, cratonic kimberlite xenoliths from larger depths often show more evidence for deformation than shallower ones\u003csup\u003e37,38\u003c/sup\u003e. Seismically, cratons often consist of two layers, separated by the mid-lithospheric discontinuity\u003csup\u003e39\u003c/sup\u003e, as also evident in some locations in SATONA (Fig. 2a,b), and the lower half can have strong seismic anisotropy (beneath ~110 km in Extended Data Fig. 6) consistent with a higher degree of deformation. If such deformation is recent, this may suggest a lower strength for the deep craton, and/or discontinuities may be mechanically weak. It has also been argued that the lower half of the craton is not always the original lithosphere, and could have experienced reworking\u003csup\u003e40\u003c/sup\u003e. Hence, it is likely that the base of the lithosphere could be weakened\u003csup\u003e41\u003c/sup\u003e (Fig. 4a), dynamically only marginally stable, and subject to deformation when external forces exist. Moreover, with the inclusion of a small amount of high-density minerals like garnet in the deep part of the craton\u003csup\u003e10\u003c/sup\u003e, the dripping could be easier to develop.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Besides driving the flow to mobilize the lithosphere mechanically, past subduction could also facilitate sub-lithospheric convection and lithospheric weakening\u003csup\u003e42\u003c/sup\u003e (Fig. 4a). When subducting slabs enter the MTZ, volatiles can be released due to dehydration and decarbonization\u003csup\u003e43-46\u003c/sup\u003e. Some of these volatiles would then ascend due to mechanisms such as subduction-induced poloidal flow\u003csup\u003e47\u003c/sup\u003e and the breakdown and re-oxidation of hydrous phases\u003csup\u003e9\u003c/sup\u003e, reduce the asthenospheric bulk viscosity, and possibly form hydrous-carbonated melts ponded beneath the LAB\u003csup\u003e45,46\u003c/sup\u003e, which could lower the strength of the bottom of the lithosphere.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eThere is evidence that such subduction-induced weakening has occurred. Relative to North America, the Farallon slab was located further west before ~100 Ma\u003csup\u003e36\u003c/sup\u003e. Around that time, kimberlite, lamproite, and carbonatite magmatism occurred at the western edge of the craton (Fig. 1), and such magmatism types favor a carbon-enriched mantle source\u003csup\u003e45,48\u003c/sup\u003e. After 70 Ma, the deep slab moved eastward\u003csup\u003e36\u003c/sup\u003e, and kimberlite magmatism appeared near the eastern edge of the craton in Kentucky (Fig. 1). Such findings are consistent with ponded hydrous-carbonated melts (Fig. 4), causing low\u0026nbsp;V\u003csub\u003eS\u003c/sub\u003e anomalies observed at the eastern and western edges of the craton and along the east coast (Figs. 1 and 2a, b).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Past plume activity could also contribute to the weakening. Previous studies suggest hotspot tracks across the continental US\u003csup\u003e49,50\u003c/sup\u003e, which could weaken the lithosphere along their trajectories. When a plume rises up from the lower mantle, the 660-km phase change with a negative Clapeyron slope could temporally stall plume ascent, causing hot materials to pond beneath the MTZ and lead to a transient boundary layer\u003csup\u003e51\u003c/sup\u003e (Fig. 4a), and plume instabilities might then ascend when encountering the Farallon slab, leading to further weakening of the craton.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;We propose that a deep, compositionally distinct craton is being pulled down by slab-induced flow. The small ~60\u0026deg;C thermal difference of the drips (Fig. 2e) suggests\u0026nbsp;they are likely from the bottom of the craton, which could have been weakened by the processes discussed above. These bottom materials are continuously entrained and mobilized by the slab-induced horizontal flow which concentrates them in the dripping area (Fig. 3a), where they then follow the downward flow to larger depths. Therefore, compared to delamination, where a large volume of lithosphere is lost locally, the thinning in this study is achieved by constantly mobilizing and removing the bottom of the lithosphere from areas affected by the horizontal flow, which could be larger than the surface projection of the dripping area (Fig. 4b).\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Combining all evidence, the lower part of a craton may be weaker and prone to deformation (Fig. 4a). Normally, these weakened portions remain part of the lithosphere due to lower density, or convective compression effects due to steep LAB gradients\u003csup\u003e52\u003c/sup\u003e. However, if strong tractions from below are present, such as when the slab transitions into the lower mantle, the bottom of the weakened craton may be mobilized and drip into the deeper mantle (Fig. 4b). Though a lower mantle slab anomaly beneath a thick craton is presently only found in North America\u003csup\u003e34,53\u003c/sup\u003e, it could have contributed to previous craton thinning events. This suggests that slabs may not only erode cratons from the sides\u003csup\u003e4\u003c/sup\u003e but craton mobilization and thinning may arise from deep mantle effects of subduction.\u0026nbsp;\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe thank R. W. Clayton, X. Pérez Campos, and C. Cardenas Monroy for acquiring data from the Mexican National Network. We thank M. Wiederspahn for managing the computation resources and data. We thank Z. Zhao and X. Li for the discussions regarding the inversion method. We thank E. Sandvol, C. Sun, D. B. Rowley, and E. K. Heilman for constructive discussions on the interpretation of structures. Mexican National Network data was obtained by the Servicio Sismológico Nacional (México), station maintenance, data acquisition and distribution is thanks to its personnel. J.H. and S.P.G. were partially supported by NSF EAR-1902400 and the Jackson School of Geosciences.\u0026nbsp;T.W.B. was partially supported by NSF EAR-2045292.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor Contributions\u003c/strong\u003e\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eJ.H. conducted the tomography, data analysis, composition conversion, and modeling in the paper. S.P.G. advised on the seismological aspects of the paper. T.W.B. advised on the geodynamical modeling aspects. J.H., S.P.G., and T.W.B. put together the main conclusions of the paper. H.A.J. processed seismic data from the ocean-bottom seismometers. C.L. helped with the inversion approach. D.T.T. helped with the computational environment. H.Z. advised on the inversion approach. Detailed interpretation of the results reflects discussions among the authors. The manuscript was written by J.H. with contributions from S.P.G., T.W.B and other authors.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting Interests\u003c/strong\u003e\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eThe authors declare no competing interests.\u003cstrong\u003e\u003cbr\u003e\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eSeismograms from networks other than the Mexican National Network were downloaded from the IRIS Data Management Center (http://ds.iris.edu/ds/nodes/dmc/). Seismograms from the Mexican National Network (https://doi.org/10.21766/SSNMX/SN/MX) were downloaded via the SSNstp client, which is free (http://www2.ssn.unam.mx:8080/getData/SSNdata_UseAndPolicy.pdf) upon request without requirements for authorization. The SATONA model is deposited to the IRIS archive (https://ds.iris.edu/ds/products/emc-earthmodels/), which also hosts US-SL-2014\u003csup\u003e22\u003c/sup\u003e, CAP22\u003csup\u003e93\u003c/sup\u003e, BBNAP19\u003csup\u003e21\u003c/sup\u003e, and NA13\u003csup\u003e94\u003c/sup\u003e that were compared. For other compared models, US22\u003csup\u003e57\u003c/sup\u003e is available on H. Zhu’s website (https://labs.utdallas.edu/seismic-imaging-lab/download/), and models in ref.\u003csup\u003e23\u003c/sup\u003e and ref.\u003csup\u003e25\u003c/sup\u003e are available as a supplement in their publications. All figures and maps are generated by the Generic Mapping Tools\u003csup\u003e95\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCode availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eComputer codes used for data processing, inversion, composition analysis, and plotting are available upon request. SpecFem3D Globe and the geodynamic tool to calculate the flow field are openly available on GitHub (https://github.com/SPECFEM/specfem3d_globe, https://github.com/geodynamics/hc).\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003ePearson, D. G.\u003cem\u003e et al.\u003c/em\u003e Deep continental roots and cratons. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e596\u003c/strong\u003e, 199-210 (2021).\u003c/li\u003e\n\u003cli\u003eYuan, H. \u0026amp; Romanowicz, B. Lithospheric layering in the North American craton. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e466\u003c/strong\u003e, 1063-1068 (2010).\u003c/li\u003e\n\u003cli\u003eJordan, T. H. Structure and formation of the continental tectosphere. \u003cem\u003eJournal of Petrology\u003c/em\u003e, 11-37 (1988).\u003c/li\u003e\n\u003cli\u003eLenardic, A. \u0026amp; Moresi, L. N. Some thoughts on the stability of cratonic lithosphere: Effects of buoyancy and viscosity. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e104\u003c/strong\u003e, 12747-12758 (1999).\u003c/li\u003e\n\u003cli\u003eGao, S., Rudnick, R. L., Carlson, R. W., McDonough, W. F. \u0026amp; Liu, Y.-S. Re\u0026ndash;Os evidence for replacement of ancient mantle lithosphere beneath the North China craton. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e198\u003c/strong\u003e, 307-322 (2002).\u003c/li\u003e\n\u003cli\u003eLee, C.-T. A., Luffi, P. \u0026amp; Chin, E. J. Building and destroying continental mantle. \u003cem\u003eAnnual Review of Earth and Planetary Sciences\u003c/em\u003e \u003cstrong\u003e39\u003c/strong\u003e, 59-90 (2011).\u003c/li\u003e\n\u003cli\u003eWu, F.-Y., Yang, J.-H., Xu, Y.-G., Wilde, S. A. \u0026amp; Walker, R. J. Destruction of the North China craton in the Mesozoic. \u003cem\u003eAnnual Review of Earth and Planetary Sciences\u003c/em\u003e \u003cstrong\u003e47\u003c/strong\u003e, 173-195 (2019).\u003c/li\u003e\n\u003cli\u003eFoley, S. F. Rejuvenation and erosion of the cratonic lithosphere. \u003cem\u003eNature geoscience\u003c/em\u003e \u003cstrong\u003e1\u003c/strong\u003e, 503-510 (2008).\u003c/li\u003e\n\u003cli\u003eLiu, L.\u003cem\u003e et al.\u003c/em\u003e Mantle dynamics of the North China Craton: new insights from mantle transition zone imaging constrained by P-to-S receiver functions. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e231\u003c/strong\u003e, 629-637 (2022).\u003c/li\u003e\n\u003cli\u003eWang, Y.\u003cem\u003e et al.\u003c/em\u003e Secular craton evolution due to cyclic deformation of underlying dense mantle lithosphere. \u003cem\u003eNature Geoscience\u003c/em\u003e, 1-9 (2023).\u003c/li\u003e\n\u003cli\u003eWhitmeyer, S. J. \u0026amp; Karlstrom, K. E. Tectonic model for the Proterozoic growth of North America. \u003cem\u003eGeosphere\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 220-259 (2007).\u003c/li\u003e\n\u003cli\u003eTappe, S., Smart, K., Torsvik, T., Massuyeau, M. \u0026amp; de Wit, M. Geodynamics of kimberlites on a cooling Earth: Clues to plate tectonic evolution and deep volatile cycles. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e484\u003c/strong\u003e, 1-14 (2018).\u003c/li\u003e\n\u003cli\u003eTao, K., Grand, S. P. \u0026amp; Niu, F. Full-waveform inversion of triplicated data using a normalized-correlation-coefficient-based misfit function. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e210\u003c/strong\u003e, 1517-1524 (2017).\u003c/li\u003e\n\u003cli\u003eKomatitsch, D., Erlebacher, G., G\u0026ouml;ddeke, D. \u0026amp; Mich\u0026eacute;a, D. High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster. \u003cem\u003eJournal of Computational Physics\u003c/em\u003e \u003cstrong\u003e229\u003c/strong\u003e, 7692-7714 (2010).\u003c/li\u003e\n\u003cli\u003eTromp, J., Komatitsch, D. \u0026amp; Liu, Q. Spectral-element and adjoint methods in seismology. \u003cem\u003eCommunications in Computational Physics\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 1-32 (2008).\u003c/li\u003e\n\u003cli\u003eTao, K., Grand, S. P. \u0026amp; Niu, F. Seismic structure of the upper mantle beneath eastern Asia from full waveform seismic tomography. \u003cem\u003eGeochemistry, Geophysics, Geosystems\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 2732-2763 (2018).\u003c/li\u003e\n\u003cli\u003eBunks, C., Saleck, F. M., Zaleski, S. \u0026amp; Chavent, G. Multiscale seismic waveform inversion. \u003cem\u003eGeophysics\u003c/em\u003e \u003cstrong\u003e60\u003c/strong\u003e, 1457-1473 (1995).\u003c/li\u003e\n\u003cli\u003eBottou, L. in \u003cem\u003eProceedings of COMPSTAT\u0026apos;2010: 19th International Conference on Computational StatisticsParis France, August 22-27, 2010 Keynote, Invited and Contributed Papers.\u003c/em\u003e 177-186 (Springer).\u003c/li\u003e\n\u003cli\u003evan Herwaarden, D. P.\u003cem\u003e et al.\u003c/em\u003e Accelerated full-waveform inversion using dynamic mini-batches. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e221\u003c/strong\u003e, 1427-1438 (2020).\u003c/li\u003e\n\u003cli\u003eHasterok, D.\u003cem\u003e et al.\u003c/em\u003e New maps of global geological provinces and tectonic plates. \u003cem\u003eEarth-Science Reviews\u003c/em\u003e \u003cstrong\u003e231\u003c/strong\u003e, 104069 (2022).\u003c/li\u003e\n\u003cli\u003eBoyce, A.\u003cem\u003e et al.\u003c/em\u003e Variable modification of continental lithosphere during the Proterozoic Grenville orogeny: Evidence from teleseismic P-wave tomography. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e525\u003c/strong\u003e, 115763 (2019).\u003c/li\u003e\n\u003cli\u003eSchmandt, B. \u0026amp; Lin, F. C. P and S wave tomography of the mantle beneath the United States. \u003cem\u003eGeophysical Research Letters\u003c/em\u003e \u003cstrong\u003e41\u003c/strong\u003e, 6342-6349 (2014).\u003c/li\u003e\n\u003cli\u003eSigloch, K. Mantle provinces under North America from multifrequency P wave tomography. \u003cem\u003eGeochemistry, Geophysics, Geosystems\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e (2011).\u003c/li\u003e\n\u003cli\u003eLiang, X., Zhao, D., Hua, Y. \u0026amp; Xu, Y. G. Seismic anisotropy tomography and mantle dynamics of central‐eastern USA. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e127\u003c/strong\u003e, e2022JB025484 (2022).\u003c/li\u003e\n\u003cli\u003eWang, H., Zhao, D., Huang, Z. \u0026amp; Wang, L. Tomography, seismotectonics, and mantle dynamics of central and eastern United States. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e124\u003c/strong\u003e, 8890-8907 (2019).\u003c/li\u003e\n\u003cli\u003eMcDonough, W. F. \u0026amp; Sun, S.-S. The composition of the Earth. \u003cem\u003eChemical geology\u003c/em\u003e \u003cstrong\u003e120\u003c/strong\u003e, 223-253 (1995).\u003c/li\u003e\n\u003cli\u003eConnolly, J. The geodynamic equation of state: what and how. \u003cem\u003eGeochemistry, geophysics, geosystems\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e (2009).\u003c/li\u003e\n\u003cli\u003eStixrude, L. \u0026amp; Lithgow-Bertelloni, C. Thermal expansivity, heat capacity and bulk modulus of the mantle. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e228\u003c/strong\u003e, 1119-1149 (2022).\u003c/li\u003e\n\u003cli\u003eYamauchi, H. \u0026amp; Takei, Y. Polycrystal anelasticity at near‐solidus temperatures. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e121\u003c/strong\u003e, 7790-7820 (2016).\u003c/li\u003e\n\u003cli\u003eBiryol, C. B., Wagner, L. S., Fischer, K. M. \u0026amp; Hawman, R. B. Relationship between observed upper mantle structures and recent tectonic activity across the Southeastern United States. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e121\u003c/strong\u003e, 3393-3414, doi:https://doi.org/10.1002/2015JB012698 (2016).\u003c/li\u003e\n\u003cli\u003eGarrity, C. P. \u0026amp; Soller, D. R. Database of the geologic map of North America\u0026mdash; Adapted from the map by J.C. Reed, Jr. and others (2005). Report No. 424, (Reston, VA, 2009).\u003c/li\u003e\n\u003cli\u003eSchutt, D. \u0026amp; Lesher, C. Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e111\u003c/strong\u003e (2006).\u003c/li\u003e\n\u003cli\u003eHager, B. H. \u0026amp; O\u0026apos;Connell, R. J. A simple global model of plate dynamics and mantle convection. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e86\u003c/strong\u003e, 4843-4867, doi:https://doi.org/10.1029/JB086iB06p04843 (1981).\u003c/li\u003e\n\u003cli\u003eLu, C., Grand, S. P., Lai, H. \u0026amp; Garnero, E. J. TX2019slab: A New P and S Tomography Model Incorporating Subducting Slabs. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e124\u003c/strong\u003e, 11549-11567, doi:https://doi.org/10.1029/2019JB017448 (2019).\u003c/li\u003e\n\u003cli\u003eForte, A. M., Mitrovica, J. X., Moucha, R., Simmons, N. A. \u0026amp; Grand, S. P. Descent of the ancient Farallon slab drives localized mantle flow below the New Madrid seismic zone. \u003cem\u003eGeophysical Research Letters\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, doi:https://doi.org/10.1029/2006GL027895 (2007).\u003c/li\u003e\n\u003cli\u003eStraume, E. O., Steinberger, B., Becker, T. W. \u0026amp; Faccenna, C. The impact of mantle convection and dynamic topography on the Cenozoic paleogeography of Central Eurasia and the West Siberian Seaway (accepted). \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e (2023).\u003c/li\u003e\n\u003cli\u003eBascou, J.\u003cem\u003e et al.\u003c/em\u003e Seismic velocities, anisotropy and deformation in Siberian cratonic mantle: EBSD data on xenoliths from the Udachnaya kimberlite. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e304\u003c/strong\u003e, 71-84, doi:https://doi.org/10.1016/j.epsl.2011.01.016 (2011).\u003c/li\u003e\n\u003cli\u003eKopylova, M. G., Russell, J. K. \u0026amp; Cookenboo, H. Petrology of Peridotite and Pyroxenite Xenoliths from the Jericho Kimberlite: Implications for the Thermal State of the Mantle beneath the Slave Craton, Northern Canada. \u003cem\u003eJournal of Petrology\u003c/em\u003e \u003cstrong\u003e40\u003c/strong\u003e, 79-104, doi:10.1093/petroj/40.1.79 (1999).\u003c/li\u003e\n\u003cli\u003eAbt, D. L.\u003cem\u003e et al.\u003c/em\u003e North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e115\u003c/strong\u003e (2010).\u003c/li\u003e\n\u003cli\u003eHu, J.\u003cem\u003e et al.\u003c/em\u003e Modification of the Western Gondwana craton by plume\u0026ndash;lithosphereinteraction. \u003cem\u003eNature Geoscience\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 203-210, doi:10.1038/s41561-018-0064-1 (2018).\u003c/li\u003e\n\u003cli\u003eLiu, L. \u0026amp; Gao, S. S. Lithospheric layering beneath the contiguous United States constrained by S-to-P receiver functions. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e495\u003c/strong\u003e, 79-86 (2018).\u003c/li\u003e\n\u003cli\u003eWang, Z., Kusky, T. M. \u0026amp; Capitanio, F. A. Lithosphere thinning induced by slab penetration into a hydrous mantle transition zone. \u003cem\u003eGeophysical Research Letters\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, 11,567-511,577 (2016).\u003c/li\u003e\n\u003cli\u003eKomabayashi, T., Omori, S. \u0026amp; Maruyama, S. Petrogenetic grid in the system MgO‐SiO2‐H2O up to 30 GPa, 1600 C: Applications to hydrous peridotite subducting into the Earth\u0026rsquo;s deep interior. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e109\u003c/strong\u003e (2004).\u003c/li\u003e\n\u003cli\u003eRichard, G., Bercovici, D. \u0026amp; Karato, S.-I. Slab dehydration in the Earth\u0026apos;s mantle transition zone. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e251\u003c/strong\u003e, 156-167 (2006).\u003c/li\u003e\n\u003cli\u003ePokhilenko, N., Agashev, A., Litasov, K. \u0026amp; Pokhilenko, L. Carbonatite metasomatism of peridotite lithospheric mantle: implications for diamond formation and carbonatite-kimberlite magmatism. \u003cem\u003eRussian Geology and Geophysics\u003c/em\u003e \u003cstrong\u003e56\u003c/strong\u003e, 280-295 (2015).\u003c/li\u003e\n\u003cli\u003eZhang, Y., Wang, C. \u0026amp; Jin, Z. Decarbonation of stagnant slab in the mantle transition zone. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e125\u003c/strong\u003e, e2020JB019533 (2020).\u003c/li\u003e\n\u003cli\u003eYang, J. \u0026amp; Faccenda, M. Intraplate volcanism originating from upwelling hydrous mantle transition zone. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e579\u003c/strong\u003e, 88-91 (2020).\u003c/li\u003e\n\u003cli\u003eSun, C. \u0026amp; Dasgupta, R. Slab\u0026ndash;mantle interaction, carbon transport, and kimberlite generation in the deep upper mantle. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e506\u003c/strong\u003e, 38-52 (2019).\u003c/li\u003e\n\u003cli\u003eChu, R., Leng, W., Helmberger, D. V. \u0026amp; Gurnis, M. Hidden hotspot track beneath the eastern United States. \u003cem\u003eNature Geoscience\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 963-966 (2013).\u003c/li\u003e\n\u003cli\u003eCox, R. T. \u0026amp; Van Arsdale, R. B. The Mississippi Embayment, North America: a first order continental structure generated by the Cretaceous superplume mantle event. \u003cem\u003eJournal of Geodynamics\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, 163-176 (2002).\u003c/li\u003e\n\u003cli\u003eDongmo Wamba, M., Montagner, J.-P. \u0026amp; Romanowicz, B. Imaging deep-mantle plumbing beneath La R\u0026eacute;union and Comores hot spots: Vertical plume conduits and horizontal ponding zones. \u003cem\u003eScience Advances\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, eade3723, doi:doi:10.1126/sciadv.ade3723 (2023).\u003c/li\u003e\n\u003cli\u003ePaul, J., Conrad, C. P., Becker, T. W. \u0026amp; Ghosh, A. Convective Self‐Compression of Cratons and the Stabilization of Old Lithosphere. \u003cem\u003eGeophysical Research Letters\u003c/em\u003e \u003cstrong\u003e50\u003c/strong\u003e, e2022GL101842 (2023).\u003c/li\u003e\n\u003cli\u003eFrench, S. W. \u0026amp; Romanowicz, B. Broad plumes rooted at the base of the Earth\u0026apos;s mantle beneath major hotspots. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e525\u003c/strong\u003e, 95-99 (2015).\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Methods","content":"\u003cp\u003e\u003cstrong\u003eFull-waveform adjoint seismic tomography.\u003c/strong\u003e Seismograms were obtained for 206 earthquakes occurred between 2004 and 2021, at 6099 individual stations (Extended Data Fig. 1). Magnitudes are mainly between 5.5 and 6.8, to both avoid low-quality data and satisfy the point source assumption which is often not valid for larger earthquakes\u003csup\u003e16\u003c/sup\u003e. All 206 earthquakes are used during the coarse-grid stage which took 126 iterations; among those, 110 earthquakes are selected for the following 90-iteration fine-grid stage (Extended Data Fig. 1). Due to strong noise on horizontal components, only vertical components were used for ocean-bottom seismometers, and both tilt and compliance noise were removed\u003csup\u003e54\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eWe used multiple body wave phases including \u003cem\u003eP\u003c/em\u003e, \u003cem\u003eS\u003c/em\u003e, \u003cem\u003ePP\u003c/em\u003e, \u003cem\u003eSS\u003c/em\u003e, \u003cem\u003ePPP\u003c/em\u003e, \u003cem\u003eSSS\u003c/em\u003e and their depth phases, as well as surface waves. For each body wave phase, we cut a waveform segment that starts 30 s before and ends 50 s after the arrival time for waveform fitting. During the coarse-grid stage, body waves were filtered between 17-100 s, and for the fine-grid stage, the filter is 14-100 s. Rayleigh and Love waves were filtered in multiple period bands, including 20-40 s, 40-60 s, 60-80 s, and 80-150 s. Waveform segments were also cut for surface waves to contain their full envelopes. In total, we used 1.5 \u0026times; 10\u003csup\u003e6\u003c/sup\u003e waveform segments.\u003c/p\u003e\n\u003cp\u003eThe inversion scheme is similar to the one in ref.\u003csup\u003e16\u003c/sup\u003e. During the inversion, we tried to maximize the correlation coefficient between synthetic and observed waveforms, which guides the velocity structure to fit phase shapes, while not being strongly affected by less constrained factors like attenuation, and compared to methods based purely on travel times, this method can exploit more information in the detailed shape of waveforms. A preconditioned conjugate gradient method is used to update the model\u003csup\u003e16\u003c/sup\u003e, and a diagonal Hessian is approximated by calculating the sensitivity kernel difference after perturbing the model\u003csup\u003e55\u003c/sup\u003e. Weights for different segments are based on four metrics: correlation coefficients, maximum cross-correlation values, time differences between synthetic and observed arrivals, and the signal-to-noise ratio of observations\u003csup\u003e16\u003c/sup\u003e. In this study, a segment is weighted as unity when these four criteria are met: correlation coefficient \u0026gt; 0.8, maximum cross-correlation value \u0026gt; 0.95, time shift \u0026lt; 4 s, and signal-to-noise ratio \u0026gt; 12; in contrast, it is weighted as zero, if one of these criteria is met: correlation coefficient \u0026lt; 0.55, maximum cross-correlation value \u0026lt; 0.7, time shift \u0026gt; 8 s, or signal-to-noise ratio \u0026lt; 8. For metrics in between these thresholds, the weight is determined by the cosine-shape function in ref.\u003csup\u003e16\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eTo overcome the cost of wavefield simulations that often limits adjoint tomography to a few tens of iterations\u003csup\u003e16,55\u003c/sup\u003e, and to reduce the possibility of being trapped in local minima, we included the mini-batch algorithm by using only a portion of earthquakes for each iteration. This algorithm is commonly used for stochastic gradient descent in machine learning\u003csup\u003e18\u003c/sup\u003e. Here, the algorithm is similar to the one in ref.\u003csup\u003e19\u003c/sup\u003e, where mini-batches greatly enhance the convergence rate by reducing the redundancy among sensitivity kernels from different earthquakes. However, unlike ref.\u003csup\u003e19\u003c/sup\u003e, the Hessian is not estimated via the L-BFGS algorithm\u003csup\u003e56\u003c/sup\u003e, so we do not need to worry about the positive definiteness of the Hessian during mini-batch selection. Hence, there only remain two key procedures for mini-batch selection: 1) choosing earthquakes to be excluded for the next iteration; 2) adding earthquakes not used in the current iteration.\u003c/p\u003e\n\u003cp\u003eWe used a similar algorithm to ref.\u003csup\u003e19\u003c/sup\u003e to select earthquakes to be excluded. The angular distance between the summed kernel using all earthquakes in the current iteration and the summed kernel for all preserved earthquakes with one earthquake removed is calculated, and the earthquake with the smallest angular distance is removed. Such earthquake exclusion is performed iteratively\u003csup\u003e19\u003c/sup\u003e until the angular distance reaches 30\u0026ordm; or the total number of excluded stations exceeds 55% of the earthquakes used for the current iteration.\u003c/p\u003e\n\u003cp\u003eThe way to add new earthquakes for the next iteration is more stochastic to avoid some earthquakes being constantly chosen. Specifically, we give each candidate earthquake a probability based on: 1) the average distance from a candidate earthquake to its three nearest earthquakes that have secured a spot for the next iteration, 2) the number of iterations that this earthquake was used in. Hence, earthquakes either far from other earthquakes or not frequently used have a better chance of being selected.\u003c/p\u003e\n\u003cp\u003eWe used model US22\u003csup\u003e57\u003c/sup\u003e as our starting structure. To account for attenuation, we assume the 3D attenuation model QRFSI12\u003csup\u003e58\u003c/sup\u003e, and while SpecFem3D Globe handles attenuation with a constant \u003cem\u003eQ\u003c/em\u003e across frequencies\u003csup\u003e59\u003c/sup\u003e, velocities shown in this study are for 1 s. We assume CRUST1.0\u003csup\u003e60\u003c/sup\u003e for Moho topography, and S362ANI\u003csup\u003e61\u003c/sup\u003e for the topography of the 410- and 660-discontinuities.\u003c/p\u003e\n\u003cp\u003eAccurate earthquake sources are crucial for tomography. In this study, we invert all source parameters\u003csup\u003e16\u003c/sup\u003e within the gCMT solution\u003csup\u003e62\u003c/sup\u003e before inverting for structure, or once the structure inversion has produced substantial updates from the model used for the last source inversion. Sources were updated eight times, and each time, the number of iterations for different earthquakes depends on their convergence rates. Meanwhile, to overcome the uneven distribution of stations, in addition to the original weight based on the four metrics, the weight for a certain station is further normalized by the sum of weights for all stations around it with azimuth difference \u0026lt; 7.5\u0026ordm; and epicentral distance difference \u0026lt; 6\u0026ordm;. A water level of 2\u0026permil; the summed weight of all stations is applied during normalization to not overweight isolated stations.\u003c/p\u003e\n\u003cp\u003eThe forward and adjoint simulations for both structural and source inversions were performed using SpecFem3D Globe v8.0.0\u003csup\u003e14\u003c/sup\u003e on a server with 4 Nvidia A100 GPUs. For each earthquake, the duration of simulated wavefields is determined by the end time of the 20-40 s Rayleigh wave segment at the farthest station. The simulation domain has a dimension of 57\u0026ordm; \u0026times; 79\u0026ordm; (Extended Data Fig. 1), and during the more influential fine-grid stage, horizontal element size is ~44 km on the free surface, making the spacing between Gauss-Lobatto-Legendre interpolation points ~11 km, and the minimum resolved period 11 s.\u003c/p\u003e\n\u003cp\u003eIn this study, we mainly focus on isotropic velocity anomalies. During simulations, we make the structure isotropic below the 410-discontinuity, and transversely isotropic above it. Hence, both isotropic velocity and radial anisotropy were solved for during inversion, and here we used the Voigt average\u003csup\u003e63\u003c/sup\u003e to approximate the isotropic velocity. To compare velocity anomalies at different depths, seismic velocities are shown as velocity perturbations with respect to a reference model. A 1D mantle reference model for V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e is obtained by averaging the final velocity model for all regions with good resolution (Extended Data Fig. 6), while in the crust, 3D velocities from CRUST1.0\u003csup\u003e60\u003c/sup\u003e are assumed as the reference.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eResolution tests of the model.\u003c/strong\u003e Two resolution tests are performed separately for V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e (Extended Data Fig. 5). In each test, the same point-shape velocity perturbations were added to the final \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e and \u0026Delta;V\u003csub\u003eP\u003c/sub\u003e /V\u003csub\u003eP\u003c/sub\u003e models. These perturbations are horizontally 4\u0026ordm; and vertically 200 km apart from each other and have reversed polarization between neighbors (Extended Data Fig. 5). The maximum amplitude of the perturbations is 0.01, and their shapes are characterized by a Gaussian with a standard deviation of 40 km. Such spatial extent is comparable with the wavelength of body waves used in this study, making the test challenging, and as set up, preferable to checkerboard tests\u003csup\u003e64\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eDuring the test, we generated waveforms for the perturbed model as \u0026ldquo;observations\u0026rdquo;, and then started from the unperturbed model to invert for the perturbations. Meanwhile, to resemble the resolution of real data, we use the same set of weighting for waveform segments as in the last tomography iteration. Due to computation resources, we did not replicate the 206 iterations for tomography, but 16 iterations were still performed for each test with the mini-batch strategy (Extended Data Fig. 5). Compared to commonly used one-iteration tests for full-waveform tomography\u003csup\u003e16,55\u003c/sup\u003e, this multi-iteration test refines structures at places with lower data coverage to better indicate the actual resolution.\u003c/p\u003e\n\u003cp\u003eA specific test is also designed for the resolution on the high V\u003csub\u003eS\u003c/sub\u003e dripping bodies (Extended Data Fig. 7). Here, due to computation limits, we use the point-spread function method\u003csup\u003e16,65\u003c/sup\u003e with one iteration. We first made a model with all high V\u003csub\u003eS\u003c/sub\u003e anomalies around the dripping region removed, and then calculate the V\u003csub\u003eS\u003c/sub\u003e sensitivity kernel for this modified model. After that, we obtain the difference in preconditioned V\u003csub\u003eS\u003c/sub\u003e sensitivity kernel between the original model and the modified one, which represents the V\u003csub\u003eS\u003c/sub\u003e change that the data would favor to compensate for the removed dripping body.\u003c/p\u003e\n\u003cp\u003eBased on results from these tests, the dripping bodies are not likely to be artifacts based on their geometries. Though the general shape of these bodies is sub-vertical, they also contain secondary features that are horizontally aligned (Fig. 2a, b), which cannot be produced by smearing along ray paths, and the resolution test shows no distortion of structures around the region (Extended Data Fig. 5). Also, given the earthquake-station distribution (Extended Data Fig. 1), there are no ray paths following the eastward dripping direction of high V\u003csub\u003eS\u003c/sub\u003e bodies on the east side (Fig. 2b).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eComposition and temperature estimation.\u003c/strong\u003e Because there are two observables (V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e), and temperature is one variable to be constrained, only one compositional variable can be solved for, so we model compositions between two endmembers. Endmember compositions are expressed as oxide weight percentages for the SiO\u003csub\u003e2\u003c/sub\u003e-Al\u003csub\u003e2\u003c/sub\u003eO\u003csub\u003e3\u003c/sub\u003e-MgO-FeO-CaO-Na\u003csub\u003e2\u003c/sub\u003eO system. We take one endmember to be pyrolite following ref.\u003csup\u003e26\u003c/sup\u003e with SiO\u003csub\u003e2\u0026nbsp;\u003c/sub\u003e= 45%, Al\u003csub\u003e2\u003c/sub\u003eO\u003csub\u003e3\u0026nbsp;\u003c/sub\u003e= 4.45%, MgO = 37.8%, FeO = 8.05%, CaO = 3.55%, and Na\u003csub\u003e2\u003c/sub\u003eO = 0.36%. For the endmember composition of a depleted cratonic mantle lithosphere, we compiled whole rock major element compositions for peridotite xenolith samples from four cratons (Extended Data Fig. 9): Greenland\u003csup\u003e66-68\u003c/sup\u003e, Slave\u003csup\u003e69\u003c/sup\u003e, Kaapvaal\u003csup\u003e70-72\u003c/sup\u003e and Siberia\u003csup\u003e73-75\u003c/sup\u003e. Based on these samples, clear negative correlations are seen between the Mg number (Mg# = 100\u0026times;Mg/(Mg+Fe)) and the FeO content; MgO and SiO\u003csub\u003e2\u003c/sub\u003e content; MgO and Al\u003csub\u003e2\u003c/sub\u003eO\u003csub\u003e3\u003c/sub\u003e content; MgO and CaO content, and all these relationships could be generalized through linear regression (Extended Data Fig. 9). Therefore, by setting the Mg# to 93.7 for the craton endmember, we can get a corresponding FeO from the first relationship, and a MgO content is obtained based on that Mg# and FeO. With the MgO, based on other relationships, other major element contents are determined. The amount of Na\u003csub\u003e2\u003c/sub\u003eO is very small and can be ignored. Together, the craton endmember has SiO\u003csub\u003e2\u0026nbsp;\u003c/sub\u003e= 46.55%, Al\u003csub\u003e2\u003c/sub\u003eO\u003csub\u003e3\u0026nbsp;\u003c/sub\u003e= 1.1%, MgO = 45.99%, FeO = 5.7%, CaO = 0.62%.\u003c/p\u003e\n\u003cp\u003eWith these endmember compositions, we estimated the dependence of \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e and \u0026Delta;V\u003csub\u003eP\u003c/sub\u003e /V\u003csub\u003eP\u003c/sub\u003e on composition and temperature at different depths. In this study, we performed the conversion between 200 and 380 km depths, because dripping bodies are clearly observed there (Extended Data Fig. 3); the depth is not too deep to be affected by potential inaccuracy from the assumed 410-discontinuity topography; and \u0026Delta;V\u003csub\u003eP\u003c/sub\u003e /V\u003csub\u003eP\u003c/sub\u003e is best resolved at ~300 km depth (Extended Data Fig. 5).\u003c/p\u003e\n\u003cp\u003eTo estimate the temperature dependence, we calculated seismic velocities at different temperature and pressure conditions. We assumed an ambient mantle potential temperature of 1350\u0026ordm;C\u003csup\u003e76\u003c/sup\u003e, and an adiabatic temperature gradient of 0.4\u0026ordm;C/km\u003csup\u003e77\u003c/sup\u003e. At each depth, we estimate the elastic V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e 50\u0026ordm;C above and below the adiabat using Perple_X\u003csup\u003e27\u003c/sup\u003e with thermodynamic data and solutions in ref.\u003csup\u003e28\u003c/sup\u003e, and pressures from PREM\u003csup\u003e78\u003c/sup\u003e. Then, we estimated the anelasticity effect for these two temperatures at the period of 20 s based on the relationship in ref.\u003csup\u003e29\u003c/sup\u003e through the very-broadband rheology calculator\u003csup\u003e79\u003c/sup\u003e with the solidus from ref.\u003csup\u003e80\u003c/sup\u003e assumed for calculation in ref.\u003csup\u003e29\u003c/sup\u003e. We assume the attenuation effect on bulk modulus to be negligible, and how complex anelastic compliance in ref.\u003csup\u003e29\u003c/sup\u003e would affect V\u003csub\u003eP\u003c/sub\u003e is derived in the same way as in ref.\u003csup\u003e81\u003c/sup\u003e for V\u003csub\u003eS\u003c/sub\u003e. The period of 20 s is chosen as it is close to the frequency of body waves used in this study. After that, gradients of V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e with respect to temperature are obtained by calculating the differences in anelastic V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e between these two temperatures and dividing them by the temperature difference of 100\u0026ordm;C. These gradients are further divided by the anelastic V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e at adiabatic temperatures \u003cimg src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAKAAAAAcCAYAAADmx7QjAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAABfaVRYdFNuaXBNZXRhZGF0YQAAAAAAeyJjbGlwUG9pbnRzIjpbeyJ4IjowLCJ5IjowfSx7IngiOjE2MCwieSI6MH0seyJ4IjoxNjAsInkiOjI4fSx7IngiOjAsInkiOjI4fV19dBTZ1gAABjtJREFUeF7tm09LVU8Yx0dfgFDuRFpkexdKgZagYApuhQpatIhEX4BRtglS0BdgvgN1IbhRtEULE8FWtRA3KSLiSukl2PnM73yvz53Ovfd4NK/9mg8c5u+ZeeaZZ2aeOVcbzhJcJFInGtMwEqkL0QAjdSUa4DVzenrqGhoafnsGBwfTGv8W0QCvmd3dXTc6Oupwvdva2ty3b9/cly9fXEdHR1rj72RiYiJzYW1ubqY1KsAlBPb29riM+FDx+fn5tLQ2iTLPEsWmqeJcRTvJhHr5CW8aJycn/kHHt27dKuWRrgdXpasnT56cJYvJ2wxxGBgY8GMD9RNS2gEfP37skpfd3bt305zIn6C5udk/X79+dU+fPi3l/e16n5qacu3t7W5jY8P19vb6vOnpaT82ePjwoU+HroY3wIWFBZ+QQiJ/nuXlZdfT05Om/n60gD59+uT6+/t9HIO0jI+Pu/X19bJj2Rvgu3fv3KtXr3xGJTDS8HzP4t69e6XysbGxNPecmZmZsjZk/NWw9Wk/hAHZOltbW2lJdaysPMhmQX5bnjWeInARYaJaW1vTnHxYWXis7jQ/oS5CH6yorvLw/ft3H1bbzfF/P3z4kKYSdDZbHyT0AQnDOpzv+GuCOHWSbTbN8R+4y/w5ymw7YT8Q+oBhG8Rtv6EPozZtXhZhP2pHskhWkSVrUdAdbRHmhfpWt5oTjVFpqxvGR54oqqu8IF+tMUkG2QDncpmQECo7nCwIjYk6YedSSqmzJG6VCKSt0mxfWbKBlY0+K8lWSamhXII8tUWYpczwnSxwxnHEaY/xIM9FjC0kSw+aI8pAY7JjVh3lFdGVZWVl5ayzs9PXJ+Tdt2/f+jIZVq22JJPmr/Hg4MB/DqjGjx8/3OzsrI/rWHr9+rVPW/r6+tLYf9y/f9+HONw6CoaHh30ourq6XCKU29/fT3PO+fz5s9+yQxJFemcX8ClCX4o2q4Gvm4zdHxX0q+PIQpu0Tb6VrdZlgbo44bdv33bJDdBtb297XV3mMwu+E/KCjtBKc9bS0pLGfqeIrsTq6qobGhpyL1++9LK8f//ej6u7u9uXc8kgn4d4JaS/w8NDH+b6Dmj9v48fP/pOEutPS/NxdHTkQxSntngePXrk84+Pj31owfDpz9bnQZFZBnsR5P8hD2PRBAuMlIUBVuZaJDuCX3gsWN14QRNVFPWPvpLdoyRbHqT7y/D8+XOvp5GREZ9uamry4YMHD3xYlFwG+OzZs9IkhROVB+tsozi1Y59Kq4YdMKv+2tpaWuN8NeWFy4aVgx0mC1ar6mjCsy5BlsXFRffixYs0de6YhxPFotYi4CNuNTh1WASS5TJfKy6qK2DX/fnzZ9nptbOz45JjuLTAitLIllxtNWmnuXPnjg8FR3cIR6aFoxcwLnscW7S7Zu1ofJvkthjCxOlGynEc9lvrZpfldtTaUTFGFmE1XcnNsAuOYyqcKPpC/qWlJX9MY7Qy1CyoHxp+1olRiyK6AtXR8ckt/s2bN35+Lor0XDr65TxaxzFRcpmjmExWmRMth5eHukAd+47aSCbNpwEH2L4DpK1jTDtKqw1bTnu2Dckf9stTyRmW/LZc8qsvxkueJdRDiPqem5vzvwDgoHMZsfIDaTnvQLtWJyHhmIE0j/SrMdk6kke6KaIrUNtcrngH2ZFZ7QDjZVzU44JCe4Qhakv4GI1ZQwkFt4LySGHE9Z7aIFQ926YgT+VZdXg/nDBbnydEitWjQVZTqpRl31GewNhsnWrGJ5gc1WXCGI+dKGBiKNOkhePNIpSFsZFH+6AxMy8inEcooivkVP/IiuzEbV+UaxzogDTvhVCHMuG1jRAaSKQYfKKwuxpocu1EENdvwJSF79w02CCsAQNGZHc3LTRBOfrIIjT20nLPWqmR/LCq9SM8oGQMLTQwJkb1MEYdVzcVDAaXQhAnzxoY8ssgsSGNKTRCjNnuflAyQF6gYbutRvLDLoDi0SEPC7qSC6KjCl3f9IWPbMiocTHG0LBYSNRhwaEHFp3iopJ9xf8JuWb4axC+byYT4ZJJ8p+ZJicn09J/j2iA1wyfnJId49Lfz/4v5PoQHbka+E6YHGHR+AxxB4zUlbgDRupKNMBIXYkGGKkr0QAjdcS5X6PDoMZR0oVDAAAAAElFTkSuQmCC\"\u003e\u0026nbsp;for further conversion.\u003c/p\u003e\n\u003cp\u003eThe compositional dependence was estimated similarly. By assuming mechanical mixing\u003csup\u003e82\u003c/sup\u003e, at different depths along the assumed 1350\u0026ordm;C adiabat, we calculated the difference in anelastic V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e between the two endmember compositions, and these differences are divided by unity to represent the compositional gradients. These gradients of V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e are then divided by the 1D reference V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e of SATONA (Extended Data Fig. 6) at these depths \u003cimg src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAKEAAAAcCAYAAAAJBd8dAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAABfaVRYdFNuaXBNZXRhZGF0YQAAAAAAeyJjbGlwUG9pbnRzIjpbeyJ4IjowLCJ5IjowfSx7IngiOjE2MSwieSI6MH0seyJ4IjoxNjEsInkiOjI4fSx7IngiOjAsInkiOjI4fV19MjD7rwAABbxJREFUeF7tm79OFU8UxwffQGipwDeQSAJCooXyACYIFRWJPoAml5oCYw021mikNVEKGwgJVNADMYZYKa+A+5k735tzD3vvXfYHXOA3n2Syd2Znz8z5szNnljBwVhAymT5yL10zmb6RgzDTd3IQ9pG/f/+G169fh4GBgTAzM5Nabz8X1SsHYR/5+PFjOD09DaTl3759S623nwvrxcEkc/00Gg0OhLE8f/48td5+6ujVCsKtra34oP3NtSqjo6Nnr169SrX6MO7Kykqq1aPO/PvB/fv3z46OjlKtPzA+tlpfX08t/51OehEjZeO0tuOpqalQOC3VMlfNwcFBGBwcDCMjI6nlbtBNr83NzTA3NxeOj49TS5MYhO/evQvF0hkeP34cGzNXz/fv38OzZ89S7e7QTS8Cs9gtw/v371NLkxiEb9++DQsLC7GhEwQqpx2VBw8epDvt0K4+nJA8Xs6nT5/SnXJ4a2z/stPW9vZ2W5+dnZ10pzv2GYqfi054KmX61OXLly9heno61aphbUvBlkI28LbwOnl7bmxspDuXQy+95ufnw9raWvtqyB5dXNLu3MTnVORovg/7u008qdPH5nPUbZ4oOcoXyvIRK0P3rUzGtOP6ueoZ21aGlys7dNK5bK512d/fj7Jkhyr4nFt6az6q2zlXtbdvq0tVvdDF2p7lMTZapJAcwm/7EPAc7cIHJcixmlSZHOp2fNunbG4ypObGmNY5wPO2j0f3LZJrx/b6QC8DA86YnZ2N8mRwycImJO4Xcbq3o6BNustnXq5tK9NJsqvM57L04hk7j3ssi522VlH0C2/evIm/2Q5ZxllSPU+fPk2/mjx69Che9/b24jYBL168iFcxMTERCuOeS1aBRHZxcTHVmpBXFAZobbnkIH75R2Y30AWdQNsXMi3IRDb37Nx6HSTo++TJk5ic//nzJ+zu7sZ05+HDh/H+y5cv4zc0rlWhL/NlbORrKy1jeHg4/ToP9uzko15cpl7ocXh4mGoVP1bbPI6HMUjxVqW71Tg5OYlXnC1ZFE7l8Pv373i1EJwoavtTaP/582fqVQ/JYvzi7Y0yLRhTbXbOvWg0GtGxq6urYWhoKBaYnJyM17ooH2QuxQrUeomq8OvXr/SrPlelF/QMQt4AAgFHobiN4KrYtxPHIseXTidzGdwXjCEuamQOGDhTsjq9vbyx6qOA7LVrfP78ue2QxycLGB8fj9c6sAhYu2lXqkPdl/cq9BL3WJ67BZZWKL/M2y1K/PjxI/1qwjYMBJjdmi2c3njDy+QVecM5mUB/nQzL+vQ6HTOWD6ayldhCQPJC+BXTopTD2ooXeGxsrLVyCPTW6ra0tJRayyFwfLpQZq9e8OmELdni/VHGRfSqAnNv+4xDIllc2pLesoOJTf5JSGmjiMJIsa7EFHnU6SuQQZsdi7qVbZ/RPMpkCJ+Qa1yK5u/R/P08KBqLxBmdLNRtQu3R2B8+fDgr8qb4JywSeasf0I8knkSffsjldyfkI6uP7C3ZsoPX2eqk+Vl7UqfIfmVU1asqfrzoTd/oFVJdBUeoTc9hFJSTcShWWUGb7pf18W0ygAryPX5+ZU7zoIN9hr4+8HyfbgEo9LdT+hJYyPMOxnn0E/RBz27o5VNBptpANvA609bNnrKVn6Onl14Ep+ZTrJBxHlw9mqfVN2rAJKsYONOZr1+/tgUWyME4yIJztAriOMpNpapexI/0oD91rzeU6dva1xDq36JMdTA6W5TAlmy53oE4hnbA5v7+TaOKXloZBUFG8Hq0CvtVvxWEiu5MPXAEKxw2pOAUn2oAzpFTCUhtXTeVKnoxf229xJF08oHYySY56q4ZnKDtiBXB51a3EV4m9NBhi1VSv6uQ/9vumuEvTnwSKwIwFI6KH/2Xl5fT3f8nOQivGb4LFitHre9rd5VKf7bLXA589C3ypRyAjrwSZvpOXgkzfScHYabPhPAPIwGpcgk81eAAAAAASUVORK5CYII=\"\u003e\u0026nbsp;for conversion.\u003c/p\u003e\n\u003cp\u003eV\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e perturbations are then converted to temperature and compositional perturbations. Here, we convert velocity perturbations rather than absolute velocities, because absolute velocities could be strongly influenced by the assumed reference anelasticity model and the exact composition of mantle rocks, while perturbations are less dependent on these reference conditions. Before conversion, to avoid misalignment of independently constrained V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e features, \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e and \u0026Delta;V\u003csub\u003eP\u003c/sub\u003e/V\u003csub\u003eP\u003c/sub\u003e are smoothed by convolving with a 3D Gaussian, whose horizontal and vertical standard deviations are 120 km and 20 km. Velocities are also converted to 20 s based on the attenuation model\u003csup\u003e58\u003c/sup\u003e for consistency. Velocity perturbations are expressed by temperature perturbation (\u003cem\u003e\u0026Delta;\u003c/em\u003e\u003cem\u003eT\u003c/em\u003e) and compositional perturbation \u003cimg src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAksAAAAfCAYAAAD3PYnKAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAABfaVRYdFNuaXBNZXRhZGF0YQAAAAAAeyJjbGlwUG9pbnRzIjpbeyJ4IjowLCJ5IjowfSx7IngiOjU4NywieSI6MH0seyJ4Ijo1ODcsInkiOjMxfSx7IngiOjAsInkiOjMxfV19SmzjqAAAERlJREFUeF7tncmLJcUTx2t+d5cZPYmIqAcVRdFRwQ0UXFAZPIijIuJBXAcv7tvNBdeD4IqCiNgqiiLugoIbbsgIetIREfE0bn+Bv/fJqW9PdHRmVdZ7r7fX8YGcepWVlZURGRkZVZ1Vs+G/EU0QBEEQBEGQ5X/tNgiCIAiCIMgQwVKw7nnnnXeaQw45pNmwYUPz2WeftblBEMADDzzQbNq0KY2PIFivRLAULILgAcfoE05z1vjzzz+bSy+9tHnttdca/iJ98sknt0eCIODm4emnn25++umnND6CYL0SwVKwiB9++KGZm5trduzY0WzcuDE5yfvvv7854ogj2hKzwffff9+ccMIJzd9//90cffTR8VQpqOKOO+7I3kzMmv289NJLzZYtW5If2HfffdvcIFifdAZLPEk4++yz273FXHvttclJsM3B+XIkBxxwQJsbrHYuuOCC5qKLLmq+/vrr5swzz5zPI7CYJY466qjmpptuaq655pqZeqrEmO0at8FkEDxs37493VBs3bo12c5ZZ53VHHbYYW2J2QAfwPhHzll+qsT8xZ/hh8D4Yl7LPW3Xn/S70nI+pR9HviFze6nMrFEMlrhLuuWWW5r33nuvzVnIL7/80jzxxBPNwQcfnH7n4HyeSDDQDj/88GU1kJWEP2PxN/61ykEHHZS2n3zySXP++een3+Tts88+6fcs8fHHHzennnpquzc5dv3Tcccd1+YGsNbHhbj33ntToM34OO2001Iefm4Wx8cHH3zQHH/88e3e5DAHaP0TT+jWIsyN77//fpr7fv311zZ3Nz///HOa85RA86DSzTffnPJXI0Pn9scff7zNnW2KwdLdd9+dlFHioYceSndTV155ZTIOj5R84oknpu2dd96ZFLzSMFBzkbCiZJIep7OexeazX8Pzzz+fnlaA7kBINlhkIlU+k8hqhMfwNY6SO9DcnYu9+8glzquBepYq0H755ZerZCzZjQW7Offcc9P6p08//XRNTJ5PPfVU6gvvFG3f2adU9k9QQyc7Oy6Wov7lQjcTBBJnnHFG+k3wVKJkO9Y35BJ9UwP1LMWfAKmTtkveEvw5m/biL7qgz7X+ibllr732ao+sLZgbsWOeupWCCaF+0Ty4Fhg6t68bRpHhIkaOnnD4vx07drQ5C9Hxubm5lHzZkaJTntJoYKT8UaQ6/7uE6rPJwnX88VI7Pap748aNbc5u3n777XTMt4998pG5BrVPbdq5c2faRycW6bBPHyuF2teH7Y/t27e3ubvYunVrSuhA+pVe0Eet7JRbCj3RJmyyjy67sSDTyIm2e/1QVrojeRvReLHjyZcBf1ypBq7BeU8++WSbs5vNmzcnmek/wW/y1K+1+HEB06x/ucHWJ7Ud8mTXt99++3x90kGtz6Gva8sOgTbV2LPsmD4rIZnQx7igK66j5PVfO17suOMc9mv6EtSf6Fvt6UJlrN2XUN02WZCFtvr2eyaRT36ftqg9tu1WtyTZr8e2geT7QefquOxX11eintVCtqdpYM7IBMekfGs8FnWsBeV0dVqpc+w5HLcdNMQQcM5MCtTBROnJ5VN/ySByUNbrjn0cjwVZfbnVBLrwus6BXJTBEXoZ0bcmPDsZAOd4mylB2b52jAPOvWYw9tkNMHlynH6tIWe3nG/bw3FbpwIOqwtsyNoR51OmxraQhXLoARk9XMfn0wb6emggQ12+TdOsf7nJyZOjZDvIac+nnO37mroFZWvHUi3oH/vrs2cFQegDGUv9Jh1Y3z4E1W+hfVZPNePFjzvVa/O6sNfUuV0y5cZ5Dtrs6+I69lz2vTzsW7uZVD57TbWpZm635GT27WSfZFGgZG2OerqutZxkg6WuBnqBZJBWQCDPdirkDMKSu641SF3Ld16pPou9EyzJ59vMdXKTSBfU7XWBcZGEHFFNuy3SfV/y+smBPpgkKU9bkNu2sRYcJXKgT+k3h58MhkDbvC3VwDm0CRkVGGnCIh95+yblGrvhOlb/JM4rITvO2YnVodpo4frKkz14OyKvpi/RCW0ojUvJZaHeUsDYBbJ4eadZP+T6IZdqoA3YLOXZUrduBuw47BprNbYD2CB1ef3Ugs5qxryFa9Iexi8J2ZCTfNmnv/nJQRBEPSWbBtonfZG43lA4jz6wcF3yRd94URu9rjhP/dSFxonOlx106b6r3y25crJnjUvk8O2kjPImlU/yqA9LfUqe7wtROse3nTJeXsr4/lObpIOVJOs5upSBwF4gX77UaRLc5+fQQFCSsrg++zUGaKG8OpC25gYs9UoOBTRdk54HB5szSm8ETFI4mZUCXSI/OkFOOesa52hBn+oH9W1OX6p/kkmwZI8lkAUZZWvs0wbaMoQauwHKjDMJgGyaZO2H397O2Zct0Z4ae8uBDehcdELbvU1SP20SHMd2h1IaF9Oqf9rQXtolfWh/qP3W2o7qH2qbgr6u8akWAiOS/Kr2h0K/qg7aUeo/dDGtvuU66EtJ9I0X+sKWF5TJ2aeHc239yE2e+tij4/T9ELgG5ylZ/UoWYX3ApPJRxtfv2y+Zau2NOiWHbYOvF8jL6bKUv9wM+s4SC/hGykqLtS0jJTQfffRRu9ekV85hv/32S9taqF+LG1mNP2rfokXmLDgbdX46rrJ9i39ZmE3dWozJa/B8W8cvSBw5i/ZX0zzyyCNpgVvXwk0PC1g5x3Pssce2v3a9EfTPP/80V111VZuz/IwCh7SombcYWISshcgnnXRS2tby8MMPN5dcckn6zWv32EFuUepXX32VtrWfHvALX3kxgGTz7KJgDwtOeWPpzTffnP8cAItJ6d8hi65r7QZ4M2rIW0Ms+JUsLI7H1keOqj3aTW7R5VCeffbZeVtFJyy2f/DBB9O+sN/VYlEnxx977LE2p57SuJhW/dOGj5TidzRG99xzz7Qd8umMIbbz+eefD7JNOw5IvJl1yimnLMjr8omMUfw4bzpr8TbX1mdCatGbn6rj8ssvTy9M5BY9sxj+mGOOafeGY18I0NtmzAM1TGO8SJ833nhj2oLkZuzn+OOPP9K2ZiE0diH5SnNfCfpyUrj+tOb2cXybbObiiy+eP1cJfvvtt7RdSQYFSygSof33aFCINUgEQ8l9b1F4UBQGgnJJJZjkVYYBw0Sac0KCiYHJQM6IduGc3njjjbQvdJw3GBjc9vVOOQbeYmPr32Cjs3EUOEWP3vrAgV5//fULXrUkD+Oizk2bNnW+IUa7vCHlUt+bMbQTxyYILsBPBuiUdlGnfzOJa+D8rS2M7hyz/fDWW28Nmgxw4upfEjZh7YJU+qQF6LVe27bvvvtu8GRQazeAvZx++untXjfoDoc4ujvrlaWPcRwlNsf1ra2ed955qS7ZAihIAGz00UcfXdCH9o1ObDfX913jYhr1W+yE2pW6kF3b9v7444+DA+0htoN+htimHQckfLJsSanr1fTXX399QduASXDozRL9Zf2IAsMPP/wwbQU2gG2N+/YU5+Pj59rvPU0a/GhiHgLXZ67xcxp+plTfF198kbZ+vsxRO/fVMI5805rbJ/Vt6mOfVsOnFrLBEkrz34/AEeUiT2CSsE6bSBQle37//fe0zRmPOth/vDL3HQuLAo+uyJPXVf2TnCuuuCI5KSYOy7///pscwDPPPNPm7EL/JcY333yTggLr6AHnSH7JiOC6665LH0G0ZV599dXmr7/+SgbJK7U8dSqB3nKG5FPX4FQgtf/++6ct4Aj8ZEB/MIEh886dO5Ou7ET64osvJlksF154YZpofCBJIDE0UJkEb3/ITPv9ZEDfM2kwgVLeyge1dsN5QyYDjQN/dzZ0EtD11KeCYLELJjP6w9rhOeeck77W/sorr7Q5u+FJxN57753KWG644YbkK7C5F154ITl8T824mKR+Cw7VjoNS6kITnNpLP992222D7bfWdhhn2M7QQGUSsI8DDzyw3dvl2xm3Q56c0W4CLHtzh/+gr73vVPBUEzTk0BMa67NgaFCgJ796OiLwT13gB8E+VRLYRWnc4ocILPoYd+7zjCvfNOZ2Ma5vY7yhK/+UDt3gn/tulJaDbLBEw72CFVnnDF6dLKeNYlCyB0Xg/HJIWc8991ybs+vJBlGqoH6vOP0uTVQcp3P9n9NKd0H8+abvz2/33HPPIj3Qzssuu6zdWwwOqu/PbzibSZ4y1CAj5m4Zp80TIyYq/wFFvrVBf6MHBVF77LFH2mLA6NU/LaAsfciTJEGfMdiW85sq2BI2iHxc/9133035hx56aNqKu+66Kz2xYAIl8LvvvvvaI8Ps5ssvv0zb2slATo1AWeCQrVOqgesxnuzdvRx7FzjFbdu2tXu7YeJjMvcQLOT+PMZ41jgvTbR94wImqX/ayJcRAGPn/DkeGznyyCNTPmBX6BlfxLjBxuz4GWI7sgF/87WUMEZ50goEqkzK5NmbpT7wD7feemu7txs+Yvvtt98uCGT48GvJ79egccUNmmCC77sp8OAX8Gk26FagUAI5sGGe+nC+h6CT8628Ah9UE2TXzH01jCMfTGNuF5P4Nr5fhczILjgX3digfMUYTRSLYPEWh0YCpv2Roaf9vsQiLM7Rb89I6M7FbjpXifLK03lqm03klRjd6SwqbxPHBdfYXFjkyCJMFmjmFkEjK20tQfs4V/q0sKiTNnD+kMXkk6DFzvSr3tjx/YUeOEb7RgNpwcJGFr56PdpkF7La/K5+6oJ+6bIbj2SiHZyX6x/ZlRbVUsbKOMRu+I0uh8D1bH1cW3nSE222bQL2vSx2fHKc/VJ70I3KlpK1Q/ZLC5tpB22Vjfhr5vTumaT+pYBrSZ9cU/qyY5fjHAPGEvuyI6i1HWSz+ePC9YeMLfqFsYF+0T2ySB5BPsdpF7+5hu0nzrdt98m+LKBxWEJ2n/OPwutKMvOb84H2ejnY9zZInq0nV0ZIB33J61/+RW3rQ+WVuK7ypDvaSrJwjDKWIfJxXGW7EnLUyqT+VOL6ypOe+J2ziVw/rxaKI5RGdhn4UKSEtQwOEUdnHQFMQ1foxwYZSwHOzgd7MmLr7PmttnAsFyAuJ+i2Vr/0hXe6OAovA7qwkxZlbJAwBHTlbWLW8UEXurQ2BJOMi5r6pw1t9RMBjt7eQNEuO/Fw3AYRKwF61iTUBeOCshZNgFYG+SJklS/omyBLqB/ZdoEe/bgNgtVEMXqZdnAzjYBiJcBZWAfjgyU5m3EGOhO4nBzOZKmDJeRQgAByil2BBPIzWdQ445VGNivHTNtxwsjo+weZ6UdscpJAiTpy9c866K30FBYmGRfQV/9SQHvt2OY3eT6QULsIIDQ2VjpgqkG2rqCTvqH9Xs+MfQVHKjMu1FVzvg1Ag2A10hkNMdH4O5Fx0CBdizCJyqEggw8s2LcByBA4lzrRMdtxJ+xaJAsTAIlr0jce8uh7wFlSbtw7y+UGnRK8SEb6JqdX9KB8ZBvHztVva2GinCZMttJvKYieZFzU1L8UYAf0p66Njfi+pW2UwcawH9nbUo/daUDb6RPJR7sZ5wqeBHaN3Ix95GMf+YYGvlxD53ZBH+f8UBCsJjbwz8iog2AevmHEQr6Rc0xvSI0calrUPkuwOBf5WBTJIkjeXkLmIFjvsOj76quvTuOeN8BGQWNaYLte/nf5IMgRwVKwCAKJ0d3moLdj1hK82bFly5b0m1emmQx43bnrDcggCIJg/TLoo5TB7EMgQfAwq4ES8C0d7pr5vhX3Cnw7KwKlIAiCoEQES8EC+NYGwcMsw4fVlvMjgEEQBMHaJv4MFwRBEARB0EE8WQqCIAiCIOgggqUgCIIgCIIOIlgKgiAIgiDoIIKlIAiCIAiCDiJYCoIgCIIg6CCCpSAIgiAIgiJN83+gl5mjlEZa6gAAAABJRU5ErkJggg==\"\u003e obtained by solving these two equations (Fig. 2d, e, and Extended Data Fig. 9). When converting for the dripping anomalies (e.g., Fig. 2d, e), only geographic locations within the dripping area with positive \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e and \u0026Delta;V\u003csub\u003eP\u003c/sub\u003e/V\u003csub\u003eP\u003c/sub\u003e are considered.\u003c/p\u003e\n\u003cp\u003eThe range of converted compositional perturbation appears to be large for the whole model (Extended Data Fig. 9i,m). However, considering these endmember compositions might not be suitable for other areas; the anelasticity model could be more complicated for hot environments\u003csup\u003e29\u003c/sup\u003e such as in the western US and Mexico (Extended Data Fig. 3). Since our goal is mainly to show whether overall V\u003csub\u003eS\u003c/sub\u003e and V\u003csub\u003eP\u003c/sub\u003e patterns for dripping bodies could correspond to an increased amount of cratonic material, we focus more on the relative compositional difference between the dripping area and regions excluding the dripping area rather than the exact values.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEstimation of thinning.\u0026nbsp;\u003c/strong\u003eAn approximate estimate of the amount of lithospheric thinning follows from assuming the amount of cratonic material is proportional to \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e, and all previously dripped lithosphere is still present in the tomographic model as fast anomalies. We first set \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e at 180 km as a reference for craton material, since it is within the craton while not too shallow to be strongly influenced by the cold temperature. Then we calculated the sum of positive \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e in the dripping area (Fig. 2c) between 250 and 550 km depths and obtained the ratio between them and the reference summed \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e at 180 km depth, which represents the volume proportion of cratonic materials at those depths and is on average ~20%. Hence, the thinning would be (550 km \u0026ndash; 250 km) \u0026times; 20% = 60 km if dripping locally. Then, since horizontal flow induced by the slab (Fig. 3a versus c) could focus material from the whole cratonic region covered by this model (Fig. 1) with an area of ~4\u0026times;10\u003csup\u003e6\u003c/sup\u003e km\u003csup\u003e2\u003c/sup\u003e, and the dripping area only has an area of ~1.5\u0026times;10\u003csup\u003e6\u003c/sup\u003e km\u003csup\u003e2\u003c/sup\u003e (~38% of the cratonic area), the average amount of thinning is calculated to be 60 km \u0026times; 38% ~ 23 km. However, some regions may experience more thinning than others, for instance, flanks at the western edge of the craton in Fig. 2a, b could partly be caused by the thinning. Meanwhile, although we have estimated compositional perturbation, and it acts as a powerful way to show the significant distinction in composition, we did not use it for the thinning estimation as what is obtained there are perturbations with large uncertainties and are also subject to many assumptions.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGeodynamic modeling of mantle convection.\u003c/strong\u003e Mantle flow fields driven by different density structures were modeled. The calculation is based on ref.\u003csup\u003e33,83\u003c/sup\u003e up to spherical harmonic degree 63. The 1D viscosity distribution from ref.\u003csup\u003e84\u003c/sup\u003e and the surface MORVEL\u003csup\u003e85\u003c/sup\u003e plate motion are assumed, and we only consider vertical variations in viscosity, for simplicity, which leads us to trust patterns more than amplitudes of flow, e.g. ref.\u003csup\u003e86\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eWe first estimated flow fields from two density models based on the global tomographic model TX2019slab\u003csup\u003e34\u003c/sup\u003e. Starting from its V\u003csub\u003eS\u003c/sub\u003e model, we removed all velocity perturbations around the dripping area at 280-600 km depths (Fig. 3b,d) to eliminate local buoyancy flows. Then for one model, at each depth slice below 600 km, we outlined the geometry of the high-velocity Farallon slab and removed all perturbations within it (Fig. 3d). After that, two V\u003csub\u003eS\u003c/sub\u003e models with or without the Farallon slab are converted to density based on a V\u003csub\u003eS\u003c/sub\u003e-density relationship\u003csup\u003e84\u003c/sup\u003e, and their corresponding flow fields were calculated (Fig. 3).\u003c/p\u003e\n\u003cp\u003eTo evaluate whether the downward flow induced by the Farallon slab is strong enough to drag down the dripping materials, we calculated flow fields after merging SATONA with TX2019slab. For SATONA, we only used reliable regions based on resolution tests (Extended Data Fig. 5). For locations that are \u0026gt;200 km within the resolution boundary (Extended Data Fig. 5), \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e in model SATONA is used, and for locations \u0026gt;200 km outside the boundary, TX2019slab\u003csup\u003e34\u003c/sup\u003e is used. For locations within 200 km of the boundary, we weighted averaged \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e between these two models to make the merged model continuous at \u0026plusmn; 200 km from the boundary. A similar merging strategy is designed for the bottom of SATONA, with depths \u0026lt; 700 km using SATONA, while using TX2019slab\u003csup\u003e34\u003c/sup\u003e velocity for depths \u0026gt; 900 km, and in between, weighted averaged values are used. After smoothing, \u0026Delta;V\u003csub\u003eS\u003c/sub\u003e /V\u003csub\u003eS\u003c/sub\u003e of the dripping area is about 1% (Extended Data Fig. 9e), in agreement with some longer wavelength previous models (Extended Data Fig. 7). Then, to understand how dripping body densities could affect flows, we tested four scenarios by making the V\u003csub\u003eS\u003c/sub\u003e-density scaling factor at 260-600 km depths for positive V\u003csub\u003eS\u003c/sub\u003e anomalies within the dripping area to be 1, 0.5, 0, -0.5 times the factor for the rest of the region\u003csup\u003e84\u003c/sup\u003e (Extended Data Fig. 10), which in the last scenario, makes dripping bodies positively buoyant (Extended Data Fig. 10d).\u003c/p\u003e\n\u003cp\u003eWe also used past slab locations to help our interpretation (Fig. 1). Mantle structures in history were predicted in ref.\u003csup\u003e36\u003c/sup\u003e by advecting velocity structures in TX2019slab\u003csup\u003e34\u003c/sup\u003e back in time. Because the present-day Farallon slab does not always appear above the 660-discontinuity, we chose to use the contour of 0.12% density anomaly at 800 km depth to outline the slab for those paleo-density models. Then, to show the relative location of the Farallon slab to the North American plate, these outlines were translated to their present location based on the assumed plate motion\u003csup\u003e87\u003c/sup\u003e in ref.\u003csup\u003e36\u003c/sup\u003e. These past slab locations were analyzed with dated magmatism (Fig. 1) in Kansas\u003csup\u003e88\u003c/sup\u003e, Arkansas\u003csup\u003e89\u003c/sup\u003e, Monroe uplift\u003csup\u003e90\u003c/sup\u003e, Kentucky\u003csup\u003e91\u003c/sup\u003e, and Virginia\u003csup\u003e92\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eMethods References\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e54 Janiszewski, H. A., Gaherty, J. B., Abers, G. A., Gao, H. \u0026amp; Eilon, Z. C. Amphibious surface-wave phase-velocity measurements of the Cascadia subduction zone. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e217\u003c/strong\u003e, 1929-1948 (2019).\u003c/p\u003e\n\u003cp\u003e55 Zhu, H., Bozdağ, E. \u0026amp; Tromp, J. Seismic structure of the European upper mantle based on adjoint tomography. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e201\u003c/strong\u003e, 18-52 (2015).\u003c/p\u003e\n\u003cp\u003e56 Liu, D. C. \u0026amp; Nocedal, J. On the limited memory BFGS method for large scale optimization. \u003cem\u003eMathematical programming\u003c/em\u003e \u003cstrong\u003e45\u003c/strong\u003e, 503-528 (1989).\u003c/p\u003e\n\u003cp\u003e57 Zhu, H., Komatitsch, D. \u0026amp; Tromp, J. Radial anisotropy of the North American upper mantle based on adjoint tomography with USArray. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e211\u003c/strong\u003e, 349-377 (2017).\u003c/p\u003e\n\u003cp\u003e58 Dalton, C. A., Ekstr\u0026ouml;m, G. \u0026amp; Dziewoński, A. M. The global attenuation structure of the upper mantle. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e113\u003c/strong\u003e (2008).\u003c/p\u003e\n\u003cp\u003e59 Komatitsch, D. \u0026amp; Tromp, J. Spectral-element simulations of global seismic wave propagation\u0026mdash;I. Validation. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e149\u003c/strong\u003e, 390-412 (2002).\u003c/p\u003e\n\u003cp\u003e60 Laske, G., Masters, G., Ma, Z. \u0026amp; Pasyanos, M. in \u003cem\u003eGeophysical research abstracts.\u003c/em\u003e 2658 (EGU General Assembly 2013, Vienna, Austria).\u003c/p\u003e\n\u003cp\u003e61 Kustowski, B., Ekstr\u0026ouml;m, G. \u0026amp; Dziewoński, A. Anisotropic shear‐wave velocity structure of the Earth\u0026apos;s mantle: A global model. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e113\u003c/strong\u003e (2008).\u003c/p\u003e\n\u003cp\u003e62 Ekstr\u0026ouml;m, G., Nettles, M. \u0026amp; Dziewoński, A. The global CMT project 2004\u0026ndash;2010: Centroid-moment tensors for 13,017 earthquakes. \u003cem\u003ePhysics of the Earth and Planetary Interiors\u003c/em\u003e \u003cstrong\u003e200\u003c/strong\u003e, 1-9 (2012).\u003c/p\u003e\n\u003cp\u003e63 Panning, M. \u0026amp; Romanowicz, B. A three-dimensional radially anisotropic model of shear velocity in the whole mantle. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e167\u003c/strong\u003e, 361-379 (2006).\u003c/p\u003e\n\u003cp\u003e64 Rawlinson, N. \u0026amp; Spakman, W. On the use of sensitivity tests in seismic tomography. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e205\u003c/strong\u003e, 1221-1243 (2016).\u003c/p\u003e\n\u003cp\u003e65 Fichtner, A. \u0026amp; Leeuwen, T. v. Resolution analysis by random probing. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e120\u003c/strong\u003e, 5549-5573 (2015).\u003c/p\u003e\n\u003cp\u003e66 Bernstein, S., Kelemen, P. B. \u0026amp; Brooks, C. K. Depleted spinel harzburgite xenoliths in Tertiary dykes from East Greenland: restites from high degree melting. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e154\u003c/strong\u003e, 221-235 (1998).\u003c/p\u003e\n\u003cp\u003e67 Hangh\u0026oslash;j, K., Kelemen, P., Bernstein, S., Blusztajn, J. \u0026amp; Frei, R. Osmium isotopes in the Wiedemann Fjord mantle xenoliths: a unique record of cratonic mantle formation by melt depletion in the Archaean. \u003cem\u003eGeochemistry, Geophysics, Geosystems\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e (2001).\u003c/p\u003e\n\u003cp\u003e68 Wittig, N.\u003cem\u003e\u0026nbsp;et al.\u003c/em\u003e Origin of cratonic lithospheric mantle roots: A geochemical study of peridotites from the North Atlantic Craton, West Greenland. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e274\u003c/strong\u003e, 24-33 (2008).\u003c/p\u003e\n\u003cp\u003e69 Liu, J.\u003cem\u003e\u0026nbsp;et al.\u003c/em\u003e Plume-driven recratonization of deep continental lithospheric mantle. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e592\u003c/strong\u003e, 732-736 (2021).\u003c/p\u003e\n\u003cp\u003e70 Griffin, W., Graham, S., O\u0026apos;Reilly, S. Y. \u0026amp; Pearson, N. Lithosphere evolution beneath the Kaapvaal Craton: Re\u0026ndash;Os systematics of sulfides in mantle-derived peridotites. \u003cem\u003eChemical Geology\u003c/em\u003e \u003cstrong\u003e208\u003c/strong\u003e, 89-118 (2004).\u003c/p\u003e\n\u003cp\u003e71 Simon, N. S., Carlson, R. W., Pearson, D. G. \u0026amp; Davies, G. R. The origin and evolution of the Kaapvaal cratonic lithospheric mantle. \u003cem\u003eJournal of Petrology\u003c/em\u003e \u003cstrong\u003e48\u003c/strong\u003e, 589-625 (2007).\u003c/p\u003e\n\u003cp\u003e72 Simon, N. S., Irvine, G. J., Davies, G. R., Pearson, D. G. \u0026amp; Carlson, R. W. The origin of garnet and clinopyroxene in \u0026ldquo;depleted\u0026rdquo; Kaapvaal peridotites. \u003cem\u003eLithos\u003c/em\u003e \u003cstrong\u003e71\u003c/strong\u003e, 289-322 (2003).\u003c/p\u003e\n\u003cp\u003e73 Boyd, F.\u003cem\u003e\u0026nbsp;et al.\u003c/em\u003e Composition of the Siberian cratonic mantle: evidence from Udachnaya peridotite xenoliths. \u003cem\u003eContributions to Mineralogy and Petrology\u003c/em\u003e \u003cstrong\u003e128\u003c/strong\u003e, 228-246 (1997).\u003c/p\u003e\n\u003cp\u003e74 Ionov, D. A., Doucet, L. S. \u0026amp; Ashchepkov, I. V. Composition of the lithospheric mantle in the Siberian craton: new constraints from fresh peridotites in the Udachnaya-East kimberlite. \u003cem\u003eJournal of petrology\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 2177-2210 (2010).\u003c/p\u003e\n\u003cp\u003e75 Ionov, D. A., Doucet, L. S., Carlson, R. W., Golovin, A. V. \u0026amp; Korsakov, A. V. Post-Archean formation of the lithospheric mantle in the central Siberian craton: Re\u0026ndash;Os and PGE study of peridotite xenoliths from the Udachnaya kimberlite. \u003cem\u003eGeochimica et Cosmochimica Acta\u003c/em\u003e \u003cstrong\u003e165\u003c/strong\u003e, 466-483 (2015).\u003c/p\u003e\n\u003cp\u003e76 Herzberg, C. \u0026amp; Gazel, E. Petrological evidence for secular cooling in mantle plumes. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e458\u003c/strong\u003e, 619-622 (2009).\u003c/p\u003e\n\u003cp\u003e77 Katsura, T. A revised adiabatic temperature profile for the mantle. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e127\u003c/strong\u003e, e2021JB023562 (2022).\u003c/p\u003e\n\u003cp\u003e78 Dziewonski, A. M. \u0026amp; Anderson, D. L. Preliminary reference Earth model. \u003cem\u003ePhysics of the earth and planetary interiors\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 297-356 (1981).\u003c/p\u003e\n\u003cp\u003e79 Havlin, C., Holtzman, B. K. \u0026amp; Hopper, E. Inference of thermodynamic state in the asthenosphere from anelastic properties, with applications to North American upper mantle. \u003cem\u003ePhysics of the Earth and Planetary Interiors\u003c/em\u003e \u003cstrong\u003e314\u003c/strong\u003e, 106639 (2021).\u003c/p\u003e\n\u003cp\u003e80 Hirschmann, M. M. Mantle solidus: Experimental constraints and the effects of peridotite composition. \u003cem\u003eGeochemistry, Geophysics, Geosystems\u003c/em\u003e \u003cstrong\u003e1\u003c/strong\u003e (2000).\u003c/p\u003e\n\u003cp\u003e81 McCarthy, C., Takei, Y. \u0026amp; Hiraga, T. Experimental study of attenuation and dispersion over a broad frequency range: 2. The universal scaling of polycrystalline materials. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e116\u003c/strong\u003e (2011).\u003c/p\u003e\n\u003cp\u003e82 Xu, W., Lithgow-Bertelloni, C., Stixrude, L. \u0026amp; Ritsema, J. The effect of bulk composition and temperature on mantle seismic structure. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e275\u003c/strong\u003e, 70-79 (2008).\u003c/p\u003e\n\u003cp\u003e83 HC, a global mantle circulation solver (2009).\u003c/p\u003e\n\u003cp\u003e84 Steinberger, B. \u0026amp; Calderwood, A. R. Models of large-scale viscous flow in the Earth\u0026apos;s mantle with constraints from mineral physics and surface observations. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e167\u003c/strong\u003e, 1461-1481 (2006).\u003c/p\u003e\n\u003cp\u003e85 DeMets, C., Gordon, R. G. \u0026amp; Argus, D. F. Geologically current plate motions. \u003cem\u003eGeophysical journal international\u003c/em\u003e \u003cstrong\u003e181\u003c/strong\u003e, 1-80 (2010).\u003c/p\u003e\n\u003cp\u003e86 Ghosh, A., Becker, T. \u0026amp; Humphreys, E. Dynamics of the North American continent. Geophysical Journal International 194, 651-669 (2013).\u003c/p\u003e\n\u003cp\u003e87 Torsvik, T. H.\u003cem\u003e\u0026nbsp;et al.\u003c/em\u003e Pacific‐Panthalassic reconstructions: Overview, errata and the way forward. \u003cem\u003eGeochemistry, Geophysics, Geosystems\u003c/em\u003e \u003cstrong\u003e20\u003c/strong\u003e, 3659-3689 (2019).\u003c/p\u003e\n\u003cp\u003e88 Blackburn, T. J., Stockli, D. F., Carlson, R. W. \u0026amp; Berendsen, P. (U\u0026ndash;Th)/He dating of kimberlites\u0026mdash;A case study from north-eastern Kansas. \u003cem\u003eEarth and Planetary Science Letters\u003c/em\u003e \u003cstrong\u003e275\u003c/strong\u003e, 111-120 (2008).\u003c/p\u003e\n\u003cp\u003e89 Eby, G. N. \u0026amp; Vasconcelos, P. Geochronology of the Arkansas alkaline province, southeastern United States. \u003cem\u003eThe Journal of Geology\u003c/em\u003e \u003cstrong\u003e117\u003c/strong\u003e, 615-626 (2009).\u003c/p\u003e\n\u003cp\u003e90 Baksi, A. K. The timing of Late Cretaceous alkalic igneous activity in the northern Gulf of Mexico basin, southeastern USA. \u003cem\u003eThe Journal of Geology\u003c/em\u003e \u003cstrong\u003e105\u003c/strong\u003e, 629-644 (1997).\u003c/p\u003e\n\u003cp\u003e91 Heaman, L. M., Kjarsgaard, B. A. \u0026amp; Creaser, R. A. The temporal evolution of North American kimberlites. \u003cem\u003eLithos\u003c/em\u003e \u003cstrong\u003e76\u003c/strong\u003e, 377-397 (2004).\u003c/p\u003e\n\u003cp\u003e92 Zartman, R. E. Geochronology of some alkalic rock provinces in eastern and central United States. \u003cem\u003eAnnual Review of Earth and Planetary Sciences\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, 257-286 (1977).\u003c/p\u003e\n\u003cp\u003e93 Boyce, A.\u003cem\u003e\u0026nbsp;et al.\u003c/em\u003e A New P‐Wave Tomographic Model (CAP22) for North America: Implications for the Subduction and Cratonic Metasomatic Modification History of Western Canada and Alaska. \u003cem\u003eJournal of Geophysical Research: Solid Earth\u003c/em\u003e \u003cstrong\u003e128\u003c/strong\u003e, e2022JB025745 (2023).\u003c/p\u003e\n\u003cp\u003e94 Bedle, H., Lou, X. \u0026amp; van der Lee, S. High-resolution imaging of continental tectonics in the mantle beneath the United States, through the combination of USArray data sets. \u003cem\u003eAuthorea Preprints\u003c/em\u003e (2022).\u003c/p\u003e\n\u003cp\u003e95 Wessel, P.\u003cem\u003e\u0026nbsp;et al.\u003c/em\u003e The generic mapping tools version 6. \u003cem\u003eGeochemistry, Geophysics, Geosystems\u003c/em\u003e \u003cstrong\u003e20\u003c/strong\u003e, 5556-5564 (2019).\u003c/p\u003e\n\u003cp\u003e96 Kennett, B. L., Engdahl, E. \u0026amp; Buland, R. Constraints on seismic velocities in the Earth from traveltimes. \u003cem\u003eGeophysical Journal International\u003c/em\u003e \u003cstrong\u003e122\u003c/strong\u003e, 108-124 (1995).\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-3254038/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3254038/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eContinental cratons are characterized by thick lithospheric roots that remain intact for billions of years. However, some cratonic roots appear to have been thinned or completely removed, with the reasons for such thinning being debated. In this study, we obtain a high-resolution full-waveform seismic tomographic model for North America which newly illuminates ongoing craton-thinning. Extensive drip-like transport of lithosphere is imaged from the base of the craton beneath the central United States to the mantle transition zone. Geodynamical modeling suggests that such dripping may be mobilized by the sinking of the deep Farallon slab, whose associated mantle flow can drag material at the base of the craton from afar to the dripping location. There, lithospheric material can descend within the ambient downward mantle flow, even though the slab is presently in the lower mantle. Dripping lithosphere could be further facilitated by prior lithospheric weakening such as due to volatiles released from the slab. Our findings show how cratonic lithosphere can be altered by external forces, and that subduction can play a key role in craton mobilization and thinning even when slabs are at great depths in the mantle.\u003c/p\u003e","manuscriptTitle":"Mobilization and thinning of cratonic lithosphere by a lower mantle slab","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-03-12 06:58:42","doi":"10.21203/rs.3.rs-3254038/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"nature-geoscience","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"ngeo","sideBox":"Learn more about [Nature Geoscience](http://www.nature.com/ngeo/)","snPcode":"","submissionUrl":"","title":"Nature Geoscience","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Nature Research","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"ab86b613-207f-4ae1-9bb3-2413723a14da","owner":[],"postedDate":"March 12th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":28862415,"name":"Earth and environmental sciences/Solid Earth sciences/Seismology"},{"id":28862416,"name":"Earth and environmental sciences/Solid Earth sciences/Geodynamics"},{"id":28862417,"name":"Earth and environmental sciences/Solid Earth sciences/Tectonics"}],"tags":[],"updatedAt":"2025-03-29T07:11:11+00:00","versionOfRecord":{"articleIdentity":"rs-3254038","link":"https://doi.org/10.1038/s41561-025-01671-x","journal":{"identity":"nature-geoscience","isVorOnly":false,"title":"Nature Geoscience"},"publishedOn":"2025-03-28 04:00:00","publishedOnDateReadable":"March 28th, 2025"},"versionCreatedAt":"2024-03-12 06:58:42","video":"","vorDoi":"10.1038/s41561-025-01671-x","vorDoiUrl":"https://doi.org/10.1038/s41561-025-01671-x","workflowStages":[]},"version":"v1","identity":"rs-3254038","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3254038","identity":"rs-3254038","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}