Research ArticleGEOLOGY

Seismological evidence for the earliest global subduction network at 2 Ga ago

See allHide authors and affiliations

Science Advances  05 Aug 2020:
Vol. 6, no. 32, eabc5491
DOI: 10.1126/sciadv.abc5491


The earliest evidence for subduction, which could have been localized, does not signify when plate tectonics became a global phenomenon. To test the antiquity of global subduction, we investigated Paleoproterozoic time, for which seismic evidence is available from multiple continents. We used a new high-density seismic array in North China to image the crustal structure that exhibits a dipping Moho bearing close resemblance to that of the modern Himalaya. The relict collisional zone is Paleoproterozoic in age and implies subduction operating at least as early as ~2 billion years (Ga) ago. Seismic evidence of subduction from six continents at this age is interpreted as the oldest evidence of global plate tectonics. The sutures identified can be linked in a plate network that resulted in the assembly of Nuna, likely Earth’s first supercontinent. Global subduction by ~2 Ga ago can explain why secular planetary cooling was not appreciable until Proterozoic time.


Plate tectonics is the fundamental unifying theory in the Earth sciences; however, there is no consensus on its initiation and earliest evolution. How far back in time plate tectonics has operated as it does today and how and when it may have initiated are thus some of the largest uncertainties in Earth history (13). The debate about when plate tectonics began and how it evolved has been handicapped by terminological ambiguities and imprecision. In effect, two fundamental but separate questions of timing need to be answered: (i) When did subduction first occur? (ii) When did global subduction evolve? Many studies arguing for Archean or even Hadean plate tectonics have considered the first question, with possible starting times as old as >4 billion years (Ga) ago (4). Such earliest evidence for localized subduction, however, does not signify a fully linked plate tectonic network (2).

The Orosirian period [2050–1800 million years (Ma) ago] is named for the exceptional number of orogenic events (>20) (5) recorded globally in this particular part of the Paleoproterozoic era. Key geological signs of modern orogenic belts such as ophiolites and blueschist/cold eclogite (3) have been reported in many parts of the world in Orosirian time (69). If Orosirian underthrust structures were widely preserved and consistent with regional geological records, it would be a global indication of collision after oceanic subduction and put a robust constraint on the timing of the modern global plate tectonic network becoming fully operational.

Plate boundaries at collisional and subduction zones imaged seismically universally record an underthrust crustal structure (1013). Seismic studies available from multiple continents already suggest that relict subduction zones of Paleoproterozoic age may be widespread (14, 15). Difficulties interpreting such putatively ancient structures arise as deep seismic images represent collages of geological history and are prone to reflect the latest or strongest magmatic and tectonic events. Subduction has recycled almost all oceanic crust older than 200 Ma in age back into the mantle, and cycles of continental assembly and breakup strongly rework continental margins, which together obscure the earliest signs of such an underthrust structure. In this study, we deployed a high-density seismic array across the northern margin of the North China craton. The tectonic stability of the cratonic margin makes it ideal for testing for the presence of a relict subduction zone of Paleoproterozoic age. The purpose of the array is to provide detailed images of the crustal structure. We then take advantage of existing seismic evidence for relict subduction zones of Orosirian age on other continents, as well as recent advances in Paleoproterozoic paleogeography, to test for the existence of a globally linked network of subduction by this time.

Geological setting of the Ordos block

North China is an Archean-Paleoproterozoic craton (5), and its western part in particular retains its cratonic nature with a thick lithospheric mantle (~200 km) and crust [~42 km (16)] and a seismically slow lower crust [Vp = ~6.7 km s−1 (17)] (Fig. 1). On the northern margin of the western North China craton, a ~1.94-Ga-old ophiolitic mélange and >1.82-Ga-old high-pressure/low-temperature eclogites have been identified along an east-west trend (all orientations in the text referring to present coordinates) (7, 9). The arcuate and broadly east-west-trending Khondalite belt of 1.9- to 1.8-Ga-old metamorphic sedimentary rocks (18) is located to the south of those high-pressure rocks (Fig. 1) and represents the edge of the stable continental margin that experienced high-grade metamorphism (5). The Khondalite belt, together with the exhumed high-pressure/low-temperature rocks, was formed by the underthrusting of northern North China beneath another continent (probably Siberia) in the assembly of supercontinent Nuna at ~1.8 Ga ago (7). South of the Khondalite belt, in the interior of western North China, is the Ordos block.

Fig. 1 North China craton and high-density seismic array.

(A) Tectonic map of the North China craton and location of the high-density seismic array (green dashed line) with the Khondalite belt (yellow) and exhumed Paleoproterozoic high-pressure rocks (stars). Phanerozoic magmatism, metamorphism, and faults are from Ren et al. (19). Paleo-Asian domain is the Neoproterozoic–Early Mesozoic accretionary orogenic system that reworked the northern margin of North China. Tethyan domain is the Mesozoic to Cenozoic collisional orogenic system that affected the southwest margin of North China. Pacific domain is the Mesozoic to Cenozoic subduction system along the eastern margin of Eurasia. Globe inset depicts positions of teleseismic events in this study, with red and blue representing earthquakes from the west and east, respectively. (B) Positions of 609 short-period seismic stations (yellow) that form a north-south transect from the stable interior of the craton to its reworked margin. Piercing points for receiver functions (RFs) at a depth of 40 km produced by 16 earthquakes from west and east BAZ (red and blue, respectively).

The Ordos block contains Triassic through Cretaceous strata that were deposited on the Archean-Paleoproterozoic basement (5, 19). It exhibits little deformation and extremely low seismicity and has an absence of Phanerozoic magmatism and metamorphism (Fig. 1) (16, 19). During late Proterozoic time, North China appears to have been largely unaffected by assembly and breakup of Neoproterozoic supercontinent Rodinia, as indicated by a lack of metamorphism and the paucity of magmatism reported, with only a few dykes and sills found (20). In earlier Proterozoic time, nearly identical positions for North China, Siberia, and Laurentia in supercontinents Nuna and Rodinia imply relatively minor Mesoproterozoic plate reorganization for North China between these supercontinents (7, 21). The Ordos block thus exhibits long-term tectonic stability, and the comparatively quiescent geological history of the area makes it an ideal place to test for the presence of a relict ancient underthrust structure.


We conducted a seismic study by deploying a ~300-km-long high-density (609 stations at ~0.5 km station spacing) seismic array that ran perpendicular to regional structures (Fig. 1). The array covered the stable Ordos block and the Yinshan Mountains and operated from April to May 2019. We used the P receiver function (RF) method to image the crustal structures. The technique enhances P to S (Ps) converted waves on velocity interfaces (the Moho) beneath an array in the records of teleseismic earthquakes. In total, 16 events (>5.6 Mw; table S1) recorded by the array are windowed around the theoretical direct P arrival (20 s before and 60 s after) and band-pass–filtered (0.05 to 5.0 Hz) for the purpose of crustal RF analysis (22). Radial RFs were estimated by a time-domain iterative deconvolution of vertical and radial seismograms (23). The processing removes the earthquake’s source signature and travel path effects from the source to beneath the recording stations.

A RF cross section from a single teleseismic earthquake from the southwest backazimuth (BAZ) is shown in Fig. 2. As the dominant P-to-S conversion phase in the shallow lithosphere, the Moho can be clearly traced around 6 to 8 s along the profile: flat (around 6 s) a distance of 0 to 20 km and 160 to 300 km, north-dipping (from 6 to 8 s) between 20 and 60 km, and disrupted (around ~6 s) between 90 and 160 km (Fig. 2A). RFs from two other earthquakes from similar southwest BAZ show nearly identical Moho features across the array (see Materials and Methods; fig. S2, A and B). The RFs from the east BAZ (Fig. 2B) show a similarly disrupted Moho from 90 to 160 km and a flat Moho at 160 to 300 km distance but an unusually muffled Moho between a distance of 20 to 60 km (fig. S2, C and D). Previous studies show that for dipping interfaces, the apparent Ps conversion amplitude may be strongly altered with respect to BAZ (24, 25). We attribute the observed discrepancy in the Moho signal to a dipping Moho. To verify this, we forward-modeled synthetic RFs based on isotropic models with a dipping Moho interface and found the down-dip direction of the Moho to be 300°NW. The detailed parameters of the model are shown in Fig. 3A. The synthetics are consistent with the observed data in which the Moho is strongest from the southwest but weak from the southeast (Fig. 3B).

Fig. 2 Seismic transect across the Khondalite belt.

(A) RF section for a teleseismic event from southwest that shows the dipping Moho. (B) RF section for a teleseismic event from southeast with unclear Moho. The average RF is shown on the right, where the direct P and the averaged Moho arrival are shown at 0 and ~6 s. (C) CCP section image of the crustal structure along transect (Fig. 1). The position of the first seismic station is at 0 km. The segment between a distance of −14 and 156 km was stacked using the three teleseismic events coming from the southwest, and the segment between 156 and 280 km was stacked using all 16 teleseismic events. T-K, Triassic through Cretaceous; Q, Quaternary.

Fig. 3 Synthetic tests for the dipping Moho.

(A) Left: Seismic stations and teleseismic event distribution, with events coming from the southwest (red) and the southeast (blue). Model used in forward-modeled synthetic RFs is shown in the middle. Right: BAZ variation in the radial RFs resulting from the dipping Moho interface. Direct P arrival at 0 s and Ps conversions from the dipping interface near 5 s. The Ps conversions with azimuths between 50° and 190° (gray shading) have a smaller amplitude. (B) Synthetic RFs (left) computed from two BAZ (128° and 233°) and their corresponding RFs computed from real data (right). RFs from 128° BAZ along a north-south trend (same direction as our profile) shows weak Moho conversions, but RFs from 233° show a coherent dipping Moho.

The time delay between the converted Ps wave and the direct wave is proportional to the depth of the converting interface and the average velocity structure above it. To construct a depth-domain crustal conversion image, the RFs were migrated to depth using a common conversion point (CCP) stacking technique to produce a depth cross section of Ps converted phases along the profile (Fig. 2C) (26). We used all time-domain RFs between a distance of 156 and 280 km, but only three from the southwest between a distance of 0 and 156 km were used to enhance the dipping Moho signal from 20 to 60 km (Fig. 2C and figs. S3 and S4). We assumed a two-dimensional (2D) velocity model based on the wide-angle reflection results from the same profile (17). In the CCP image of Fig. 2C, the Moho can be continuously traced as a strong-amplitude, positive-polarity converted phase that starts flat at the south end of the array at a depth of ~35 km and dips to the north to a depth at ~50 km at a distance of 20 km to over 60 km. The off-plane westward dip is inferred from the synthetic tests. A vertical ~5-km Moho offset is seen at 60 km distance, and north of it, the Moho continues around ~40 km depth to 280 km distance, including a weak and complex Moho a distance of 90 to 160 km that may be overstepped by surface reverberations off the shallow slow-velocity sediments. One strong-amplitude phase in the shallow crust at a depth of ~20 to 30 km with a south-dipping structure between a distance of 60 and 80 km is caused by multiple reflections (fig. S6).

The observed crustal image along the array suggests a complex regional geological history. The strongest tectonic modification of the stable Ordos block is the Cenozoic graben system along its north-to-northwest edge (Fig. 1). A graben is located, for example, between 100 and 140 km along the array, where RFs clearly record its U-shaped structure at a depth of 0 to 15 km (Fig. 2B). The Moho beneath the graben zone is obscured. Seismic studies in adjacent regions to the north and south revealed a flat Moho at a depth of ~40 km (17), suggesting that the graben is not a crust-cutting feature. Therefore, the dipping and deepening Moho from a distance of 20 to 60 km along the array is not related to Cenozoic extension. North China was assembled into supercontinent Pangea by 250 Ma ago (27) and has been stable in the interior of Eurasia ever since. Only the margins (north of the Khondalite belt) of North China have been affected by Paleozoic to Cenozoic magmatism and faulting (19). The study area lacks such tectono-magmatic activity, so the abrupt Moho change is not related to those Phanerozoic events (28).

Furthermore, the North China craton rifted away from Nuna in the Mesoproterozoic along its northernmost margin ~200 km north of the Moho offset region. That means that the region where the deepening Moho was detected escaped known tectono-magmatic events after the Paleoproterozoic assembly of supercontinent Nuna. Hence, the dipping Moho should be a relict collisional structure from the Nuna assembly. Northern North China is likely a passive margin, based on the comparable geology between the Orosirian metamorphic Khondalite belt and the Cenozoic Greater-Himalaya series in northern India. A subducting oceanic slab can drag the connected passive continental margin beneath an overriding continent after closure of the ocean basin by underthrusting. A fragment of 1.94-Ga-old oceanic crust (7) to the northwest of our study region corroborates the existence of an ancient subducting oceanic slab. Exhumed 2.6 GPa and 670°C rocks at ~1.8 Ga ago in the same belt further indicate the existence of deep and cold subduction (7, 9). The relationship of observed tectonic units on the surface indicates that the study region could only represent a fold and thrust belt on a passive margin. The north-dipping Moho, together with the structure of geological units, was formed by the underthursting of northern North China by subduction in Orosirian time. To achieve continental collision and such an underthrusting structure by ~1.8 Ga ago, the subduction zone must have initiated by at least ~2 Ga ago as indicated by subduction-related magmatism in this region (29).

The Himalayan Range offers a modern analog for the Orosirian event in northern North China (Fig. 4). The Indian continent underthrusted the Eurasian continent along the Yarlung suture since 60 Ma ago to present. Because of continuous convergence of the two continents, the crust of the Indian passive margin was thickened and shortened with the development of a series of fold and thrust belts. Such a sequence, made partly of sedimentary rocks on the passive margin, underwent granulite to eclogite facies metamorphism accompanied by granitic intrusions (30, 31). Geophysical studies reveal that the Moho ~200 km south of the suture beneath the fold and thrust belt starts to dip to the north and deepens from a depth of 40 to 50 km to ~70 km beneath the suture (10). Notably, the Moho stays at a similar depth of ~70 km northward across the suture for 150 km distance beneath the Tibetan Plateau. The underthrust structure in a subduction zone exists beneath the active margin, but the structure may progressively migrate toward the deep crust in the passive continental margin in a collisional belt (10, 32). The parallels between Cenozoic India and Orosirian North China strongly suggest modern-style plate tectonic subduction, at least developing locally, by ~2 Ga ago (Fig. 4).

Fig. 4 Modern Himalayan orogen and Orosirian Khondalite orogen.

Comparison of surface tectonics (left) and deep seismic structures (right) between modern (top) and Paleoproterozoic (bottom) continent collision zones. Light blue lines on the tectonic maps represent the locations of the RF profiles. Modern Himalayan range and its RF image across the main central thrust from 50 to 105 km along the profile of Nábelek et al. (10). Tectonic map of western North China (5) and the deep RF image shown from a distance of 5 to 60 km along our array (Fig. 2C). Comparison of entire profiles is shown in fig. S7.


All continents record Orosirian orogenesis (33), and most of those that have high-resolution seismic profiles exhibit the preservation of such Orosirian underthrust structures (Fig. 5; see Materials and Methods). For example, Laurentia was one of the core constituents of supercontinent Nuna (33) and the Trans-Hudson Orogen was a major suture uniting the cratons of Laurentia (34). Various seismic studies have revealed that crustal structures are complex along the length of this orogen. In the northern part of the orogen, RF results indicate that the Moho depths are 5 to 20 km deeper than neighboring Archean cratons and exhibit a north-dipping collisional polarity. Meanwhile, in the southern part of the orogen, the Moho is more complex and recorded a double-vergent (toward the northwest and southeast) style collision. These deep seismic results are consistent with geological observations (Fig. 5) (15). Baltica is another core continent of Nuna with seismic data available. Seismic reflection results indicate that the Svecofennian Orogen has a northeast-dipping underthrust Moho, which is consistent with the occurrences of Orosirian arc magmas and ophiolites at the surface (Fig. 5) (14).

Fig. 5 Global subduction and supercontinent Nuna.

Reconstruction of supercontinent Nuna ~1.8 Ga ago (33), including the position of North China (7, 35). Convergence directions (black arrows) inferred from subduction polarity (black triangles) indicated at seismic profiles (blue). Speculated continuation of seismically confirmed sutures (red). Details of Orosirian orogens are in the Supplementary Materials. A 2.15-Ga-old blueschist immediately preceding the Orosirian period is shown (6). Green region marks configuration of core continents assembled through two to three major sutures. Putative positions in the larger supercontinent Nuna are suggested for other less studied continents (slightly transparent). Seismic evidence is consistent with previous paleomagnetic and geological configurations in the green region under the paradigm of a modern, global subduction-driven plate tectonic network.

Geological and paleomagnetic comparisons have begun to constrain the paleogeography of supercontinent Nuna, and the Laurentia, Baltica, and Siberia connection and configuration is least uncertain (5, 33). Independent correlation of high-pressure eclogite facies metamorphism and paleomagnetic data (7, 35) suggest the addition of North China juxtaposed against Siberia in this core Nuna configuration. The fact that these Nuna-forming sutures are closely associated and can be linked as longer subduction zones is reminiscent of the modern plate network. The seismic images thus support the supercontinent configuration previously constrained by paleomagnetism and geology. Future seismic work on other continents with larger paleogeographic uncertainty may help constrain their positions within Nuna. One important observation on the modern Earth is that the negative buoyancy and suction of subducting oceanic slabs are a major driving force of plate tectonics (36). Orosirian ophiolites have also been reported in those seismically studied areas, providing independent corroborating evidence that subducting oceanic slabs drove continental convergence and underthrusting (7, 37, 38). If subduction was the driving force behind the convergence of Laurentia, Baltica, and Siberia, together with North China, then at least approximately two-thirds of supercontinent Nuna was assembled by two or three major subduction zones (Fig. 5). Thus, our new result, coupled with similar images from around the world, establishes a global picture of subduction-driven continental assembly in Paleoproterozoic time. Evidence for global subduction by ~2 Ga ago is consistent with paleomagnetic evidence for relative lateral motion between cratons of Nuna at this age (39). The speculation that Nuna may represent Earth’s first supercontinent (34, 40) may further imply that our identification of global subduction during its assembly suggests the establishment of the modern global plate tectonic network at this time.

Many studies argue that plate tectonics has operated since early Archean time, with some even suggesting Hadean plate tectonics (4). However, arguing for such tectonic antiquity is difficult to reconcile with other myriad geological constraints that demonstrate that many aspects of modern-style plate tectonics (low-pressure metamorphism and ophiolites) do not appear to emerge until significantly later in Earth history (3). Large strike-slip faults documented on Superior craton are an indicator of torsionally rigid plates (41), but few Archean crustal fragments are large enough to document large faults, rendering the structural evidence for Archean tectonics uncertain and localized at best. In this paper, we demonstrate an essentially global picture of subduction-driven plate tectonics by ~2.0 Ga ago in Proterozoic time. Although this is the oldest evidence of global subduction to date, it does not rule out the possibility of future studies finding earlier evidence, if indeed it can be demonstrated to be global in scope. The earliest evidence for subduction is strictly localized (4), leaving the possibility that the development of the modern, globally interconnected plate network did not arise until billions of years later. This assertion is consistent with modeling studies that demonstrate that much of the continental crust of Archean cratons could have been generated in the absence of subduction (42, 43). The time gap between localized subduction events starting at ~3 to 4 Ga ago (2, 4) and the formation of a global subduction network at ~2 Ga ago related to present-day plate tectonics might offer an explanation to the diachronous stabilization of cratons between ~2.9 Ga ago and as young as ~1.6 Ga ago (40). It has been argued that plate margins associated with cratons may have actually been comparable in length to the present-day plate network (44), thus implying relative uniformity of plate tectonic processes since ~2 Ga ago, consistent with our findings.

The ultimate linking up of the world’s plates into the global subduction network by ~2 Ga ago would have irreversibly accelerated mantle cooling. Global subduction cools the mantle both directly through slab cooling (45) and indirectly by promoting more efficient large-wavelength mantle convection (46). The expected acceleration of subduction-driven planetary cooling at ~2 Ga ago can thus explain both why secular mantle cooling is insignificant until Proterozoic time and the onset of low-temperature metamorphism at this time (4749).


Seismic line

RFs from temporary short-period deployments are usually much noisier compared to broadband RFs due to their variable recording conditions and narrow frequency response (50). However, the advantage of using short-period RFs in this study is to reach an unprecedented lateral resolution (~0.5 km) over a long distance (~300 km) that is critical in imaging large-scale whole-crust architecture. The stacking of a large number of RFs with similar ray paths can greatly enhance coherent crustal RF signals (51). To compute the RFs, 16 teleseismic events with magnitude greater than 5.6 from the 1-month deployment were carefully selected (table S1). Teleseismic events with epicentral distances between 30° and 90° are often used in RF studies. In this study, to increase the number of usable events during the short deployment, a minimum epicentral distance of 17° is used. Park and Levin (52) suggested that reliable RFs can be obtained for such short epicentral distances.

Three-component recordings of distant earthquakes were band-pass–filtered (two-pole Butterworth filter) between 0.05 and 5 Hz and rotated into the z-radial-tangential (Z-R-T) ray–based coordinate system. The first round of visual inspection was performed to remove the low-quality records (lack of a clear P arrival on the Z component). The responses of our short-period sensors have a cutoff frequency at ~2 Hz, but for many earthquakes, there is still sufficient signal below 2 Hz for the purpose of crustal RF analysis. To take advantage of this low-frequency signal, we broaden the signal frequency range by deconvolving the instrument response. This enables us to reach down to nearly 10 s (figs. S1).

Radial RFs were computed using a time-domain iterative deconvolution method (23), with the Gaussian parameter of 4.0. The second round of visual inspection was then performed to remove low-quality RFs (poor consistency compared with surrounding stations). As a result, a total of 16 events contributed 7960 radial RFs to our study. A 2D CCP RF stack code based on Dueker and Sheehan (26) was then applied to migrate the time-domain RFs.

Depth migration of RFs

To construct a depth-domain crustal conversion image, the RF records were stacked and converted from time to depth using the CCP approach (22) to produce a depth section of Ps converted phases along the profile. The CCP stacking process consists of two steps. First, the energy at each point on RFs after the direct P is assumed to be generated by Ps converters along the S-wave ray path. The amplitude is back-projected from the station to the conversion point using a reference velocity model. Second, the crustal volume is divided into bins of designated size and then all amplitudes in the same bin are averaged to produce a structural image.

The velocity model used for the time to depth migration was constructed on the basis of a wide-angle reflection result from the same profile (17). The Vp/Vs ratio was set to 2.0 in the sedimentary layer where Vp < 5.8 km s−1 and 1.73 elsewhere. Figure S4 shows the CCP stacks computed using different Vp/Vs ratios (1.70, 1.73, and 1.80), which results in a Moho depth shift of up to 5 km (fig. S3). The dipping Moho feature, however, remains unchanged. The bin size was 2 km wide and 0.5 km thick, and because our code is 2D only, we mapped off-plane RFs within 15 km distance into the profile. Neighboring bins overlap by 3 km along the profile to generate a smooth 2D image. Considering the existence of the dipping interface at a distance of 20 to 60 km, we only use the RFs from the southwest in the first half of the profile and all the RFs in the second half (Fig. 2C). We also provide CCP profiles of all RFs for comparison (fig. S4).

Analysis of north-dipping Moho interface

A continuous north-dipping Ps conversion can be traced from ~5 to 8 s between a distance of 20 and 60 km for RFs from the southwest (fig. S2, A and B). A disrupted Moho at first, and then its flat continuation, can be traced at ~6 s for the distances of 60 to 300 km. In contrast, RFs from the southeast (fig. S2, C and D) exhibit a weak Moho present between a distance of 20 to 60 km, and the Moho between a distance of 60 to 300 km is consistent with RFs from the southwest. This discrepancy in the Moho signal between 20 and 60 km is attributed to the westward off-plane dipping of the Moho, which we test further below. The synthetic waveforms were constructed using the Raysum code (53), and theoretical RFs were calculated by the time domain deconvolution method using a ray parameter of 0.065 s km−1 (Fig. 3) (23). To balance between computational efficiency and BAZ coverage, theoretical RFs between 0° and 360° BAZ with 10° spacing were calculated. The amplitude of Ps conversions clearly varies as a function of BAZ (Fig. 3A). In particular, the amplitude of the Ps conversions with a BAZ between 50° and 190° is weak, which is consistent with our observations. The result indicates that if the Moho dips to the northwest at an angle near 30°, the RF Moho signal perpendicular to the dipping direction (e.g., BAZ = 233; Fig. 3B) is strong but weak from the up-dip direction (e.g., BAZ = 128; Fig. 3B).

In addition, a dipping Moho interface may cause notable azimuthal variations in the moveout and amplitude of the converted phase (24, 53, 54). In fig. S5, both radial and tangential RFs are computed using the Raysum code and are then compared with the real data. Because of the very short deployment time, we do not have a large BAZ coverage; nonetheless, systematic variations in both components can still be matched for the dipping Moho interface. For example, in the up-dip direction (BAZ, ~120°), both radial and tangential RFs have the smallest amplitude at ~6 s where the Moho signal arrives. On the tangential RFs, the only event from the northeast quadrant (BAZ, 51.4°) yields a negative Moho signal around 6 s (blue), while on the opposite side of the dipping direction, events from ~130° to 253° BAZ show that the Moho phase has a positive sign (red). The moveout of the Moho signal on both radial and tangential RFs also largely agrees. Dipping intracrustal structures and seismic anisotropy may further complicate tangential RFs (54), but a detailed investigation of their effects is beyond the capacity of our limited data.

Analysis of shallow structures

Between a distance of 0 to 50 km, three large continuous positive phases are present (fig. S6). We performed a 1D synthetic test to suggest that they are the conversions and multiples from shallow structures. We assume that the positive phases near 0.8 and 2 s a distance of 0 and 50 km were conversions from interfaces at depths of 3 km (L1) and 11 km (L2), and synthetic RFs with and without multiples were constructed using the Raysum code (fig. S6B). The multiple of L1 and the conversion of L2 combine to form a larger peak near 2 s, and the multiple of L2 combines with the Moho conversion near 5.5 s. When L1 and L2 became shallower to the north, the multiple of L2 became shallower and is separated from the Moho conversion between 50 and 80 km distance, where the direct P conversions of L1 and L2 and the multiple of L1 are combined to form a complex peak as the first arrival (fig. S6C). Therefore, the positive inclining phases at 2 to 5 s between a distance of 50 to 80 km are multiples from L2 that interfere with the Moho signal at 0 to 40 km. Note that there is a continuous near-horizontal positive phase at 3.5 s between 0 and 60 km, which seems to be connected with the inclining phases between 50 and 80 km and does not represent multiples of L2.

Paleogeographic reconstruction

Siberia, Laurentia, and Baltica are reconstructed according to Evans and Mitchell (33), which includes an intraplate rotation within Siberia discussed by those authors. This reconstruction is also supported by identification of coeval large igneous province events (55). North China is configured closely with Siberia on the basis of both geological (7) and paleomagnetic (35) correlations. Antarctica (Mawson Continent) is restored to Australia to match the basement geology of the Curnamona and Mount Isa regions (56). Updated compared to its placement by Evans and Mitchell (33), Australia is configured on the basis of a new paleomagnetic pole from 1790 Ma ago (57). At present, the positions of Amazonia, West Africa, and India are largely conjectural pending firm paleomagnetic testing (58). Additional uncertainties related to the paleomagnetic reconstruction of Nuna are addressed in the Supplementary Materials.


Supplementary material for this article is available at

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


Acknowledgments: The manuscript benefitted from discussions with M. Brown, L. Chen, J. H. Guo, P. F. Hoffman, X. F. Liang, J. M. Wang, W. J., Xiao, T. Y. Zheng, and Coffice 442 (IGGCAS). Y. D. Liu and B. Xia assisted in compiling Fig. 1. P. Cawood and an anonymous reviewer provided constructive comments that improved the manuscript. This is a contribution to IGCP 648 and contribution 1512 from the ARC Centre of Excellence for Core to Crust Fluid Systems ( Funding: This study was supported by the NSFC (grant nos. 41888101 and 91855207) and the Program of the Chinese Academy of Sciences (grant no. Y201715) to B.W. Author contributions: B.W. initiated the idea and designed the project. B.W., X.Y., and X.T. deployed the seismic array. X.Y. and X.T. processed the seismic data. H.Y. worked with modeling methods. R.N.M. and U.K. did the paleogeographic reconstruction. All authors interpreted the data. B.W. and R.N.M. shaped the scope of the study and wrote the manuscript with input from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The computed receiver function waveforms used to compute the CCP model are stored at the Center for Open Science website:, doi:10.17605/OSF.IO/4SPHU. All other data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.

Stay Connected to Science Advances

Navigate This Article