The PSCS, which has now been completely removed from the surface, is inferred to have subducted during the opening of the present-day South China Sea (SCS). This has been argued based both on a gap in tectonic reconstructions prior to this opening (e.g., Holloway, 1982) and the observation of fast slab-like anomalies below this region in tomographic models (e.g., Hall and Spakman, 2015). However, the details of this subduction remain uncertain. Traditionally, the PSCS has been envisioned as subducting southward below Borneo (e.g., Hall, 2002; Seton et al., 2012; Zahirovic et al., 2014). Based on slab-unfolding methods, Wu and Suppe (2018) instead argue that the PSCS subducted both northward and southward during the opening of the SCS, due to both the position of the subducted PSCS slab in the mantle and its subhorizontal orientation.

Unlike the PSCS, the PSP is still present at the Earth’s surface (Figure 1). But its Cenozoic history remains uncertain due to its isolation from global plate circuits by surrounding subduction zones (see discussion in Wu et al., 2016). Due to this, vast differences exist in the proposed PSP reconstructions. A common view is that the PSP was much larger in the past, and began its current subduction along the Eurasian margin before 30 Ma (e.g., Seno and Maruyama, 1984; Xu et al., 2014; Zahirovic et al., 2014). On the other hand, based on the slab anomalies underlying the PSP in the present day, Wu et al. (2016) argue for a much smaller historic PSP, which begins subducting along the Eurasian margin at ∼16 Ma. In this reconstruction, the spatial gap produced by the smaller PSP is made up by the now fully-subducted “East Asian Sea” (EAS) slabs, which have been inferred based on the imaged slabs to have subducted southward below the PSP. Therefore, under this viewpoint, the smaller PSP overrides the southward-subducting EAS slabs as it migrates northward towards its present day position. For a more detailed review of the reasoning behind these differing reconstructions, we refer the reader to the discussion of PSCS reconstructions in Lin et al. (2020) and the discussion of the PSP and EAS slabs in Wu et al. (2016) and Zahirovic et al. (2014).

The study of Lin et al. (2020) investigated the mantle flow predicted by differing tectonic reconstructions in East Asia. Reconstructions were assimilated as surface velocity boundary conditions in the global mantle convection code TERRA (Bunge and Baumgardner, 1995). This allows geological information about past plate motion to enter the flow, as downwellings develop at locations of plate convergence (e.g., Bunge et al., 2002). Mantle convection was modelled by solving the conservation equations for mass, momentum, and energy in the truncated anelastic liquid approximation (Jarvis and McKenzie, 1980a). Lithospheric thickness variations were generated self-consistently by solving the conservation equations with tectonic velocities and a surface temperature T _S = 300 K as boundary conditions. The computational domain was discretised on a regular grid based on the icosahedron, with ∼80 million grid points and a minimum resolution of ∼25 km. The mantle was heated from below and from within with a CMB temperature T _CMB = 4,200 K and a radiogenic heating rate of 6 × 10^–12 W⋅kg⁻¹. The viscosity was assumed to be Newtonian and depends on temperature (T) and depth (d) as η d , T = η 0 A d exp V ∗ d R E − R CMB − E ∗ T − T S T CMB − T S , (1) where V ^∗ = 3.976 and E ^∗ = 4.610 are non-dimensional constants controlling the degree of depth and temperature dependence respectively. An increasing viscosity with depth is imposed by the factor exp[V ^∗ d/(R _E − R _CMB )], in which the depth is normalised by the mantle thickness R _E − R _CMB . The reference viscosity is η ₀ = 10²² Pa⋅s and is multiplied with a radial pre-factor A(d) in order to impose a low-viscosity asthenosphere. The resulting maximum, minimum, and average radial viscosity profiles are reported in Supplementary Figure S1, together with the total radial pre-factor η ₀ A(d). This model rheology leads to two-sided subduction, and one-sided subduction was not forced by imposing either dipping weak plate boundaries or a prescribed dipping slab thermal structure. Slab material therefore sinks vertically unless carried passively by mantle flow.

Four different tectonic reconstructions (Models 1a, 1b, 2a & 2b) were implemented, which are designed to test different end-members of the proposed reconstructions of the PSP and the PSCS introduced above. Models labelled “a” (Models 1a & 2a) implement a southward PSCS subduction, while models labelled “b” (Models 1b & 2b) implement a double-sided PSCS subduction. Models labelled “1” (Models 1a & 1b) implement a large PSP in the style of Zahirovic et al. (2014), while models labelled with a “2” (Models 2a & 2b) implement a simplified form of the slab-unfolded model of Wu et al. (2016), featuring a small PSP that overrides a single “vanished ocean” (VO) plate. In these reconstructions, the small PSP overrides the VO plate and begins its subduction along the Eurasian margin at ∼12 Ma. The naming scheme for these models is summarised in Table 1. Figure 3 shows the reconstructed plate motions at 30 Ma, highlighting their key differences. The primary aim of Lin et al. (2020) was to test the present-day mantle heterogeneity predicted under each reconstruction against seismic tomography. Model 2b led to a present-day mantle state which had the highest correlation with seismic tomography, and so is our preferred reconstruction.

Prior to the assimilation of these reconstructions at 30 Ma (Models 1a & 1b) and 45 Ma (Models 2a & 2b), the reconstruction of Matthews et al. (2016) is used as a surface boundary condition on velocities since 410 Ma. A reference model implements the reconstruction of Matthews et al. (2016) for the entire 410–0 Ma timespan. All models therefore implement the subduction of a large Izanagi plate during the Mesozoic and the subduction of the Izanagi-Pacific ridge at ∼55 Ma. Running these models since the Paleozoic means that an Earth-like convective planform is allowed to develop, generating vigorous plumes in the Pacific domain. Thus, a realistic hot upper mantle arises in this region during the Mesozoic. Along with the implemented Izanagi subduction and downwelling slabs in East Asia, this provides a realistic scenario in which to test our asthenospheric flow hypothesis.

The asthenospheric build-up and subsequent influx proposed in Section 2 arises under all considered reconstructions. Transects of the reference model showing the full history of this influx from 60 to 0 Ma can be found in Supplementary Figures S2–S8. Hot asthenospheric material is driven westward through the Izanagi-Pacific slab window from ∼50–30 Ma, subsequently flowing laterally below East Asia during the late Cenozoic. Figure 4 shows four snapshots from the evolution of our preferred Model 2b, which implements a double-sided subduction of the PSCS and a smaller PSP. While this influx occurs under all considered tectonic reconstructions, the differences between the plate motion models described above affect how much hot material flows below East Asia, which parts of the asthenosphere are occupied by slabs, and the induced late Cenozoic dynamic topography. For further details of these models we refer the reader to Lin et al. (2020).

Intraplate basaltic magmatism has been widespread in East Asia since the latest Paleogene (Figure 1), spanning from the well-known eruptions in Indochina and Hainan (e.g., Yan et al., 2008), throughout the South China Block (SCB) (e.g., Gong and Chen, 2014), and extending far northward to the Changbai mountains on the Sino-Korean border (e.g., Liu et al., 2015). These basalts are predominantly of Ocean-Island-Basalt (OIB) type (e.g., SCS and Indochina, Yan Q. et al. (2018) and references therein; SCB, Ho et al. (2003); Tan-Lu fault zone, Chung, (1999); Changbai Mountains, Liu et al. (2015); inter-regional, Kimura et al. (2018)), have moderate to high mantle potential temperatures (∼1,320°C–1,440°C, see Kimura et al., 2018), and coincide spatially with seismically slow anomalies (Ball et al., 2021). Together these observations suggest that a source of hot, upwelling asthenospheric material is required in order to explain this recent phase of eruptions. We note that this doesn’t necessarily preclude a contribution from slab dehydration (e.g., Richard and Iwamori, 2010), at least for some of this magmatism.

In order to investigate a potential link between the late Cenozoic volcanism in East Asia and our asthenospheric flow hypothesis, we compare the timing and location of intraplate volcanic eruptions to the distribution of hot asthenospheric material from 30 to 0 Ma predicted by the models of Lin et al. (2020). We trace the distribution of temperature anomalies in the asthenosphere by defining the model lithosphere by the 1,300°C isotherm, and the model asthenosphere as the 200 km thick channel below this. In the case that a subducting slab is encountered (i.e., T < 1,300°C through the asthenosphere) a default lithosphere thickness of 100 km is taken, which produces the sharp jumps near subduction zones seen in Figure 5. To avoid the influence of “drip” artefacts left over in the asthenosphere from the subduction of the Izanagi plate (see discussion of model limitations in Section 5.4), we take the mean anomalous temperature in the model asthenosphere, as this best reflects the full lateral extent of the influx of hot material. We make use of the global database of intraplate volcanic samples compiled by Ball et al. (2021). While this is primarily a database of Neogene volcanism, many Paleogene samples are also included, meaning it is suitable for the 30–0 Ma timespan considered here. Additional eruptions in the SCB, Hainan, and South Korea are taken from Gong and Chen (2014) and Arai et al. (2001). The post-rift basaltic eruptions in the PRMB have been included from the map compiled by Sun et al. (2010).

Figure 5 shows the mean anomalous temperature in the asthenosphere at 6 Ma intervals from 30 to 0 Ma, overlaid with intraplate eruptions of age within ±1.5 Ma. Given that the uncertainties in dating are often on the order of ∼1 Ma although in some rare cases, were as high as ∼7–10 Ma, see database of Ball et al. (2021), and the slow pace of mantle flow, this approximation is not likely to result in any spurious connections between the modelling results and mapped eruptions. Maps for the full time-series are reported in Supplementary Figures S9–S19. We briefly note that the anomalous temperature values are systematically over-estimated in the models of Lin et al. (2020), and so we simply consider the distribution of anomalously hot material through time, rather than interpreting its values.

The predicted lateral extent of hot material beneath East Asia from 30 to 0 Ma shows a good match to the distribution of intraplate volcanism during this time (Figure 5). The volcanism in NE China is matched well by all models. Moving further south, Models 2a & 2b, which implement a small PSP, predict a hot, positively buoyant asthenosphere below SE China, Taiwan, and Penghu from 24 to 0 Ma, which is consistent with the OIB-type magmatism observed here (e.g., Chung et al., 1995; Ho et al., 2003). In contrast, Models 1a & 1b, which feature a large PSP, predict predominantly downwelling, cold material in the asthenosphere below these regions, which is inconsistent with this form of magmatism. Some of the westernmost eruptions, such as the Hannuoba and Taihang basalts (HN and TH in Figure 1), are only reached a few Myrs after the onset of magmatism (cf. 24–12 Ma snapshots of Figure 5). Others, such as the eruptions in the Indochina peninsula, are not reached in any of the considered models, although the flow is typically ≲1,000 km from these eruptions. The lateral distribution of hot asthenospheric material at the present-day time-step in Figure 5 shows that this hypothesis is also able to account for the seismically slow anomalies observed in the upper mantle throughout East Asia. Here as well the modelled asthenospheric material does not match the full westward extent of the slow anomalies observed in tomographic models (e.g., Schaeffer and Lebedev, 2013) by ∼1,000 km.

Stratigraphic data sets hold crucial information about past mantle flow. Mapping conformable and unconformable contacts on inter-regional scales can yield a proxy for the large-scale uplift and subsidence associated with dynamic topography, which may in turn be linked to convective motion in the Earth’s mantle (Friedrich et al., 2018). This method of “hiatus mapping” has been carried out on the scale of single (Vibe et al., 2018; Carena et al., 2019) and multiple (Friedrich et al., 2018; Hayek et al., 2020, 2021) continents, and has provided important information on the influence of mantle flow on large-scale topographic changes since the Mesozoic. However, likely due to the potential tectonic influence of nearby subduction zones, East Asia has not been the focus of any hiatus mapping studies to date—although an inter-regional Cretaceous–Paleocene sedimentary hiatus in Southeast Asia has been noted and linked to mantle flow (Clements et al., 2011). In our preceding comparison, it was found that the proposed influx of hot asthenospheric material approximately matches the distribution of late Cenozoic intraplate volcanism in East Asia. Due to its positive buoyancy, this influx is also expected to induce a coeval dynamic uplift. We may therefore use the hiatus mapping concept to trace this uplift and compare it to the model predictions of Lin et al. (2020), allowing us to test our hypothesis against the sedimentary record over a large region.

In order to gauge large-scale topographic changes during the late Cenozoic, we make use of generalised stratigraphic columns for a group of sedimentary basins (both onshore and offshore) in eastern China. We map columns for the basins surrounding Hainan Island (Shi et al., 2011); the Pearl River Mouth Basin (PRMB, Shi et al., 2008); the Taiwan region (including both exposed strata on Taiwan Island and the rift basins in the Taiwan Strait (TWS), Lin et al., 2003; Shi et al., 2008; Huang et al., 2012); the East China Sea Shelf Basin (ECSSB, Shi et al., 2008; Wang et al., 2019); the Subei Basin (Zhou et al., 2019); the South Yellow Sea Basin (SYSB, Yang et al., 2020; Zhang R. et al., 2020); the Hefei Basin (HF, Jiaodong et al., 2012); the North Yellow Sea Basin (NYSB, Wang et al., 2017); the Bohai Bay Basin (BBB, Tang et al., 2019); and the Songliao Basin (SLB, Song et al., 2018). A map of these basins can be found in Figure 1.

Stratigraphic columns are mapped at the resolution of geological series (ten to a few tens of Myrs, see the definition of series as a unit of chronostratigraphy, e.g., Cohen et al., 2013; Ogg et al., 2016), given that stratigraphic columns with stages resolution were available only in a few offshore basins in eastern China. We follow the approach and terminology of Friedrich (2019), Carena et al. (2019) and Hayek et al. (2020), Hayek et al. (2021), in which, for example, “Base of Miocene” refers to the conformable/unconformable surface below this series (i.e., Oligocene sedimentation/hiatus), and a series is marked as hiatus if and only if there are no preserved sediments of this series (i.e., missing stages within a series are not marked as hiatus). In some cases, stratigraphic columns were available for individual sub-basins and so these too have been mapped where available. Stratigraphic columns are mapped as single points, approximately central in a given basin or sub-basin. While not portraying the true scale of a given basin, this approach tends to avoid issues with intraplate deformation as stratigraphic data are reconstructed back in time to their appropriate palaeolocations. In some cases, however, we still run into issues of this kind. For example, the generalised stratigraphic columns for Taiwan Island represent exposed strata which were deposited on the continental margin during the early Cenozoic. These sediments have been significantly deformed and moved from the location where they were deposited (Suppe, 1980) during the Taiwan Orogeny since ∼6.5 Ma—which is not accounted for when these points are reconstructed back in time. Due to this, the reconstructed point data will not reflect the true palaeolocations of deposition, and so this must be kept in mind when interpreting the resulting maps. Given the size of the Songliao Basin, and the lack of nearby basins, a representation of its stratigraphic column as a single point is particularly misleading (for other larger basins, stratigraphic columns for each sub-basin were available, or nearby basins filled out the surrounding space). We therefore make the choice to instead plot the generalised stratigraphic column of the Songliao Basin as a group of points—one for each sub-unit of the basin (sags and uplifts). This better represents the true area of the data used in NE China. Given that the Songliao Basin has not undergone significant deformation since the Paleocene (e.g., Song et al., 2018), the locations of the points relative to each other will have remained relatively constant throughout the late Cenozoic. This means that the reconstruction, which keeps these points fixed relative to each other, will result in a reasonable representation of their palaeolocations.

As can be seen in Figure 4, the modelled influx of hot asthenospheric material is present below far eastern China from the mid/late Eocene and flows below this region (i.e., SLB, BBB, S/SYSB, HF, ECSSB, TWS, TNB, see Figure 1) during the Oligocene. This coincides temporally with an inter-regional late Eocene–Oligocene sedimentary hiatus recorded in these basins, the exact timing of which depends somewhat on location. To investigate the potential link between this large-scale hiatus and the inflowing Pacific material, we compare our Base of Miocene hiatus map to the modelled change in dynamic topography during the Oligocene. A comparison between our Base of Pliocene hiatus map and the predicted change in dynamic topography during the Miocene can instead be found in Supplementary Figure S20, along with an associated discussion in Supplementary Text S4.

The change in dynamic topography is calculated in a Lagrangian sense. That is, we calculate this change in the reference frame of the Eurasian plate as it moves with respect to the underlying time-dependent dynamic topography signal. The modelled asthenospheric flow predicts a broad N-S trending band of uplift during the Oligocene (Figure 6), which matches the distribution of sedimentary hiatus during this time to first-order. The BBB and NYSB (highlighted as blue stars in Figure 6) both stand out as outliers amongst the sedimentary hiatus in eastern China. It is worth noting, however, that rifting continued into the late Oligocene in these basins (Qi and Yang, 2010; Wang R. et al., 2020), and so rapid syn-rift tectonic subsidence is likely to largely overwrite any dynamic topography signal in the stratigraphic record. In addition to this, while syn-rift sedimentation occurred in these basins, the late Oligocene strata are clearly eroded, and have been completely removed in some shallower parts of the basins. Indeed, this aspect of the development of the BBB has been investigated further using Apatite Fission Track Thermochronolgy (AFT), which found a rapid cooling event from 27 to 16 Ma, removing an estimated 1.5–1.8 km of sediment (Tang et al., 2019). This uplift event has been noted too in the NYSB (Wang R. et al., 2020), and while it has been attributed to tectonic compression, the stratigraphic cross-sections of the basin are also consistent with a regional uplift characteristic of dynamic topography. This uplift is thought to have continued until ∼14 Ma in the NYSB (Wang R. et al., 2020). And so, while a full sedimentary hiatus was not induced in these basins, they remain consistent with the inter-regional picture of uplift and erosion during the Oligocene and early Miocene.

Additional misfits can be found in some offshore regions, such as the eastern ECSSB and the northern margin of the SCS, in which the predicted dynamic uplift covers regions of Oligocene sedimentation. However, the early Cenozoic rifting and offshore setting of these regions means that it is somewhat less likely that a dynamic uplift will lead to a series-long sedimentary hiatus. Indeed, while a series-long hiatus has not been recorded, there nevertheless exists evidence of late Cenozoic uplift and erosion on the northern margin of the SCS (Section 5.2.3). Intense erosion has been noted too in the Taiwan region (e.g., Fuh, 2000; Shi et al., 2008), but the onset of sedimentation during the late Oligocene means that this event is not highlighted when mapping is carried out at the resolution of geological series (Figure 6).

Similarly to our comparisons with intraplate volcanism, Figure 6 shows that the downwelling induced by the large PSP of Models 1a & 1b leads to an inconsistency with the geological record. These models predict a large-scale dynamic subsidence in the ECSSB during the Oligocene, which is inconsistent with the observed sedimentary hiatus. In Models 2a & 2b, instead, the influx of hot asthenospheric material is allowed to predominate the dynamic topography signal over a larger region, leading to a consistency with the sedimentary hiatus in the ECSSB.

As can be seen in Figure 6, the SCB was excluded from our mapping. Cenozoic deposits are largely absent from this region (Figure 1), meaning that the hiatus mapping concept is not useful in gauging large-scale topographic changes during this time. Cenozoic sedimentation occurred mainly in small intermontane basins, which are not likely to be sensitive to the changes in relative sea level induced by dynamic topography. We therefore take a different approach in the SCB, which is discussed in the following section.

The SCB has garnered significant attention from thermochronological studies (e.g., Yan et al., 2009; Wang et al., 2015; Li and Zou, 2017; Su et al., 2017; Tao et al., 2019; Wang Y. et al., 2020; Qiu et al., 2020), likely due to the regional history of exhumation inferred from the widespread present-day exposure of Mesozoic and Paleozoic rocks. That this exhumation has continued into the Cenozoic is revealed by the exposure of late Cretaceous granites formed during the flat-slab subduction of the Izanagi plate (Charvet et al., 1994), indicating at least kilometre-scale erosion has occurred during the Cenozoic (Yan Y. et al., 2018). While thermochronological studies have uncovered multiple phases of heating and cooling since the beginning of the Mesozoic (Tao et al., 2019; Wang Y. et al., 2020), we focus here on the late Cenozoic cooling history of the region. This allows us to investigate the dynamic uplift expected under our hypothesis over a larger region, by filling a spatial gap in the SCB (Figure 6).

A significant rapid cooling event during the Oligocene and Miocene has been found consistently in numerous studies (see those compiled by Qiu et al., 2020), with samples covering a large (∼1,000,000 km²) region. To test whether this regional-scale erosion of the SCB during the late Cenozoic is linked to the dynamic uplift expected under our hypothesis, we compare the dynamic topography predictions of each model to the AFT-derived cooling histories published by Wang Y. et al. (2020), Yan et al. (2009), and Li and Zou, (2017). A total of 21 samples are considered, of which 9 have been reported in Figure 7 to highlight key characteristics in different regions. Clouds of “acceptable” (95% confidence limit) fits have been included to indicate the uncertainty of the AFT cooling histories. The full set of samples is included in Supplementary Figures S21–S41.

The AFT cooling histories generally show a distinct increase in cooling rate during the early Miocene (see, e.g., samples 20GN-58 and 01JF-181 in Figure 7) or late Miocene (samples 02QSH-1 and 01JF-214). In northwesterly regions (samples 02QSH-1, 01JF-68 and 01QT-03), the model predictions of dynamic topography are broadly similar, and we find a good match between the observed cooling histories and the increase in dynamic topography induced by the inflowing asthenospheric material from ∼20 Ma. In easterly samples (NA06-5, 01JF-214, 01JF-217), the closer proximity to the Ryukyu Trench leads to significant differences in the model predictions based on their implementations of the PSP history. In line with our previous comparisons, we find that the large PSP implemented in Models 1a & 1b leads to a dynamic subsidence which is inconsistent with the observed cooling histories. In contrast, the smaller PSP implemented in Models 2a & 2b allows the influx of hot asthenospheric material to dominate the Neogene dynamic topography in this region, leading to a better prediction of the increased cooling rate during this time.

A key consideration to be made when interpreting this widespread cooling effect is that it occurs outside of the well-constrained Apatite Partial Annealing Zone (PAZ, 110°C–60°C). Generally, one cannot confidently interpret cooling histories below the lower bound of 60°C, due to the lack of fission track annealing below this temperature. However, we note that the samples considered often leave the PAZ at ∼20–10 Ma, and given their present-day exposure at surface temperatures (∼20°C), their Neogene cooling histories remain well constrained. For example, samples 20GN-58 and 01JF-181 (Figure 7) cool through the lower bound of the PAZ at ∼10 Ma, and so we may infer that ∼40°C of cooling has occurred since this time. Taking a lower bound on the geothermal gradient of ∼20°C/km, we may then estimate ∼2 km of erosion has occurred with a rate of ∼0.2 km/Ma since ∼10 Ma. We therefore find the recent regional cooling to be a robust effect in the derived cooling histories.

We note additionally that the predicted uplift is generally significantly lower than the eroded sediment thickness inferred from AFT cooling histories (Figure 7). For instance, the 2 km exhumation of samples 20GN-58 and 01JF-181 since 10 Ma, as estimated above, corresponds to only ∼250–500 m of predicted dynamic uplift. A key reason for this is that isostatic rebound has not been accounted for in these comparisons, which provides an additional component of uplift as overlying rock is eroded. Part of this disparity may also arise due to the modelled asthenospheric influx not flowing far west enough in significant amounts, in keeping with our comparison to the magmatic record (Section 4.1). Nevertheless, we find an encouraging match between the recent regional cooling and the modelled increase in dynamic topography during the Neogene.

Subsidence analysis may be used to make quantitative estimations of the history of dynamic topography affecting rifted sedimentary basins. Traditionally, this is achieved by estimating strain rates during rifting, for instance by inverting the backstripped subsidence (White, 1994), and using these in conjunction with lithospheric stretching models (McKenzie, 1978; Jarvis and McKenzie, 1980b) in order to forward model the post-rift tectonic subsidence. This forward modelled tectonic subsidence may then be subtracted from the subsidence gleaned from backstripping (Sclater and Christie, 1980), in order to find the anomalous (i.e., dynamic) uplift or subsidence which has affected a given basin during its post-rift phase. A refined method of this kind has recently been proposed by Zhao (2021), which sought to remove several of the underlying assumptions of earlier methods. Zhao (2021) used fault growth rates derived from seismic cross sections to allocate the total strain of a basin, allowing spatially and temporally varying strain rates to be calculated. The use of strain rates which are independent of the backstripped subsidence means that one separates the tectonic and dynamic components of subsidence from the outset. This allows for the full (i.e., including the syn-rift phase) history of anomalous subsidence to be derived, and the tectonic subsidence due to secondary rifting phases to be accounted for. This method was applied to the Qiongdongnan Basin (QDNB), which, given its location on the northern margin of the SCS (Figure 1), provides us with an additional observation with which to test our asthenospheric flow hypothesis. Zhao (2021) found that the QDNB was anomalously uplifted by ∼1 km beginning at 28.4 Ma, which decreased significantly after 11.6 Ma. Our intraplate volcanism comparisons in Figure 5 show that the influx of hot material reaches the QDNB during the Oligocene, and may therefore provide an explanation for this uplift. We therefore compare the anomalous uplift since 45 Ma with the modelled change in dynamic topography since this time. The resulting comparisons, and a map of the relevant borehole locations, is shown in Figure 8.

To first order, we find a good match between the modelled timing of dynamic uplift and subsidence affecting the QDNB and the anomalous uplift at boreholes P1–P8. The modelled uplift tends to have an earlier onset than the results of Zhao (2021) in Models 2a & 2b, and has a particularly high amplitude in Model 2b. Given the double-sided PSCS subduction in Model 2b, the northern subduction zone leads to additional downflow in the QDNB region, pulling more asthenospheric material towards the basin from the east. The northern subduction zone additionally leads to a larger initial dynamic subsidence in the QDNB, and so the rebound from this initial state following the subduction of the PSCS provides an additional component of uplift. The inflowing hot asthenospheric material comes to predominate the dynamic topography signal in this basin in all considered models during ∼30–25 Ma (i.e., the dynamic topography signal becomes net positive). This timing aligns with the onset of uplift found by Zhao (2021), although the modelled amplitude of the uplift tends to fall short by ∼500 m (excluding P1 and P8, in which we find an excellent match to the reconstructed amplitude). As noticed in Section 4.1, the modelled influx of asthenospheric material does not flow far west enough. This is particularly apparent in the regions of Hainan and Indochina (Figure 5), which are adjacent to the QDNB. A more vigorous influx of hot material in the geodynamic models would presumably lead to a better match to the observed amplitudes.

Another difference we observe between the results of Zhao (2021) and the model predictions of Lin et al. (2020) is the spatial variability of the dynamic uplift signal: the models predict a maximum of ∼250 m of variation of the peak of uplift between the borehole locations, while subsidence analysis predicts it to be ∼500 m. This may have resulted either from a lack of small-scale structure in the modelled dynamic topography, or from underlying assumptions in the analysis of Zhao (2021). At short enough wavelengths, the viscous stresses which would otherwise induce dynamic topography are instead opposed by the strength of the lithosphere. This filtering effect is modulated by the flexural rigidity of the lithosphere and is strongly wavelength-dependent (see, e.g., Turcotte and Schubert, 2002). To account for this effect, the dynamic topography predicted by the models of Lin et al. (2020) has been expressed as a spherical harmonic expansion truncated at l _max = 40, which filters out short-scale (i.e., ≲500 km) variability. While this is a reasonable first-order approximation of lithospheric behaviour, it might filter out too much of the mantle signal, particularly in regions of young oceanic lithosphere where the flexural rigidity might be low. The analysis of Zhao (2021), on the other hand, assumes local isostasy in both the backstripping routine and the lithospheric stretching model used to forward model the tectonic subsidence. It is therefore assumed that both the sediment load and crustal thinning during rifting vary smoothly with respect to the elastic thickness of the lithosphere. If present, short-scale variability in these features leads to additional spatial variability in the derived subsidence histories (instead of being compensated via regional isostasy). The gradient in dynamic topography of ∼10 m/km between boreholes P1 and P3 at ∼25 Ma is particularly high, suggesting that the assumption of local isostasy might not be adequate.

A multitude of hypotheses have been proposed to explain the recent inter-regional phase of mafic volcanism in East Asia (Figures 1, 5) on smaller, regional scales. Of the mechanisms proposed, a hypothesised mantle plume below Hainan Island on the northern margin of the SCS is by far the most prevalent (e.g., Yan et al., 2008; Lei et al., 2009; Wang et al., 2012; Huang, 2014; Xia et al., 2016; Yan Q. et al., 2018; Zhang Y et al., 2020, etc.). The “Hainan plume” has been proposed to explain the basaltic eruptions in Leizhou, Hainan, Indochina and the SCS basin since the late Miocene (e.g., Yan et al., 2008), as well as the seismically slow anomalies observed below the SCS region in numerous tomographic models (e.g., Lei et al., 2009). Limited tomographic resolution in the lower mantle means it is not yet possible to conclusively prove the existence of a Hainan plume extending down to the Core Mantle Boundary (CMB), and while a seismically slow anomaly in the upper mantle appears to be a robust feature, similar anomalies are observed throughout East Asia, as discussed in Section 1. It is worth noting too that a plume in this region has not been resolved in global tomographic studies which have resolved plumes (e.g., French and Romanowicz, 2015; Hosseini et al., 2020), although it may be the case that tomographic resolution is not yet high enough to resolve the proposed Hainan plume. Additionally, the magmatic record in this region lacks the clear arrival of a plume head, which is evidenced at well-known plume locations by the formation of a Large Igneous Province (LIP) (Richards et al., 1989). The westward migration of a small plume into this region has been found to be a plausible geodynamic mechanism using mantle circulation models (Zhang and Li, 2018). As noted above, it remains unclear, though, whether the magmatic and tomographic evidence supports this hypothesis.

The lack of a clear plume signal is true, similarly, for the eruptions on inter-regional scales. In addition to the absence of LIPs, the eruptions don’t show the linear age progression characteristic of plume-fed hotspots, and don’t align with absolute plate motions, with contemporaneous eruptions spanning over large (∼5,000 km) distances. The pattern of these eruptions (Figure 1) is instead reminiscent of those observed in the far-field of mantle upwellings, for instance in Africa (Lees et al., 2020). An east–west younging in the onset of eruptions is expected under our hypothesis, and is indeed observed in several regions. For example, the onset of volcanism in SE China and Taiwan is earlier than in Hainan, which in turn has an earlier onset than the main phase of eruptions in Indochina. This east–west younging can also be observed in NE China, with an earlier onset of volcanism in Changbai than the more westerly regions of Halaha and Dalenuor (Liu et al., 2001). However, there are clear exceptions to this trend. For instance, the onset of volcanism in Hannuoba is earlier than the more easterly eruptions in the Jiaodong peninsula, which in turn is earlier than the eruptions in the Korean peninsula. It is worth noting, however, that the study region is interspersed with offshore areas within which the magmatic record is highly biased towards the most recent eruptions. Post-rift intraplate eruptions in, for example, the ECSSB and SYSB (e.g., Yao et al., 2020) must be mapped and dated using seismic stratigraphy and borehole data due to recent sediment cover. As such, it is less likely that such eruptions have been mapped and dated accurately in the current literature. This, more generally, leads to a bias towards onshore regions. Indeed, this bias can be clearly observed in Supplementary Figures S9–S19, with the vast majority of mapped eruptions lying within onshore regions. It may therefore be the case that a clearer east–west younging will emerge as the magmatic record in these offshore regions becomes more complete.

Another prominent line of thinking invokes the influence of the Pacific slab in sourcing the intraplate volcanism (e.g., Liu et al., 2001) and underlying seismically slow anomalies (e.g., Ma J. et al., 2019) in NE China, through deep slab dehydration (e.g., Richard and Iwamori, 2010; Kuritani et al., 2019) or the interaction of the slab with a hydrous Mantle Transition Zone (MTZ) (Yang and Faccenda, 2020). While we don’t necessarily preclude some contribution from these mechanisms, we note that they are not required under our hypothesis. Relatively high mantle potential temperatures (∼1,380°C–1,440°C) have been noted in NE China (Kimura et al., 2018), indicating that a hot mantle source is required in this region. In addition to this, geochemical analyses have found no evidence of a typical slab-derived isotopic signature in these basalts (Chen et al., 2007; Ward et al., 2021). In any case, the distribution of seismically slow anomalies and intraplate volcanism spans far into regions which are not underlain by slabs (e.g., the SCB, see Gong and Chen, 2014; He and Santosh, 2021), suggesting that these mechanisms are not sufficient on inter-regional scales.

The asthenospheric flow hypothesis presented here unifies the inter-regional observations of intraplate OIB-type magmatism into a single picture. The approximate match observed between the timing and location of these eruptions and the models of Lin et al. (2020) shows this to be a viable new mechanism in explaining these eruptions. An open question under our hypothesis is the differing eruption volumes observed between regions—for example the low-volume eruptions observed in SE China and Taiwan compared to those observed in Indochina or NE China (Figure 1). There are many factors which may explain this disparity in future studies, such as local variations in lithospheric thickness, the influence of underlying slabs, or as noted by Lin et al. (2020), local passive return flow in the mantle—which may play a role in facilitating particularly voluminous eruptions in certain regions. Nevertheless, the correspondence we have found between the distribution of Pacific asthenospheric material and the timing and location of intraplate volcanism in East Asia provides a starting point for future investigation into the local variations observed in the erupted magmas, such as volume and geochemistry.

The comparisons presented in Section 4 provide new constraints on the uncertain Cenozoic tectonics of East Asia introduced in Section 3.1. The earlier onset of PSP subduction along the Eurasian margin at ∼30 Ma in Models 1a & 1b leads to predominantly downwelling, cold material in the asthenosphere below SE China and the ECSSB. This downwelling corresponds poorly with the timing of OIB-type magmatism (Figure 5) and induces large-scale dynamic subsidence which is inconsistent with the observations of inter-regional uplift in far eastern China (Figure 6) and the SCB (Figure 7). In contrast, the later onset of PSP subduction at ∼12 Ma in Models 2a & 2b allows a hot, positively buoyant asthenosphere to dominate over a larger region during the late Cenozoic. This provides a good fit to the observations of OIB-type magmatism and inter-regional erosion presented in Figures 5–7. The comparisons with the results of Zhao (2021) show that Models 1a & 1b provide a better fit to the onset of uplift in the QDNB (Figure 8). However, the peak of uplift predicted by Models 1a & 1b tends to come later than that found by Zhao (2021), with the latter generally lying between the peaks of Models 1a & 1b and Models 2a & 2b. As such, it is not possible to make a strong distinction between the PSP implementations based on this form of evidence. Based on Figures 5–7, we therefore find new support for the small PSP and existence of the EAS slabs proposed by Wu et al. (2016).

While the PSCS subduction style was the focus of the original study of Lin et al. (2020), it has a largely negligible effect on the asthenospheric flow considered here. This is highlighted by the minor differences between models which implement a single- (Models “a”) or double-sided (Models “b”) subduction of the PSCS in Figures 5–7. The location of the QDNB on the northern margin of the SCS means that we see more significant differences between the models based on their PSCS implementation in this region. However, the comparisons shown in Figure 8 indicate no clear preference for a particular reconstruction. Of Models 2a & 2b, the double-sided PSCS subduction (2b) provides a better match to the amplitude of uplift; whereas of Models 1a & 1b, the single-sided PSCS subduction (1a) provides a better match to the amplitude of uplift. And so, overall, it is not possible to make any clear distinctions between the implemented PSCS reconstructions based on the comparisons presented here.

Several of our comparisons in Section 4 show that the asthenospheric flow predicted by the models of Lin et al. (2020) does not flow far west enough to explain certain observations. The comparisons in Figure 5 show a ∼1,000 km shortfall in matching both the westernmost volcanic eruptions considered and the present-day distribution of seismically slow anomalies (e.g., Schaeffer and Lebedev, 2013). In addition to this, Figure 8 shows that, while the models predict well the timing of uplift and subsidence of the QDNB, the influx of hot material does not reach this basin in significant enough amounts to match the amplitude of anomalous uplift found by Zhao (2021) by ∼500 m.

One key limitation of these models which may reconcile this disparity in future studies is the implemented mantle viscosity profile. The asthenosphere features only a moderate (200x) viscosity drop with respect to the lithosphere, which is achieved with a gentle gradient from 50 to 250 km depth (Supplementary Figure S1). Due to limitations on the maximum per-gridpoint viscosity variations in the TERRA mantle convection code, this viscosity profile is a compromise between viscosity drop and channel thickness. These models therefore implement a relatively large viscosity in the asthenosphere, which is coupled with a temperature-dependent rheology without a yield stress. As a consequence, following the cessation of the subduction of a particular plate (e.g., the subduction of the Izanagi-Pacific ridge), slabs gradually break off at mid-asthenospheric depths, leaving large “drips” of cold slab material which hang from the lithosphere into the asthenosphere. This effect can be clearly observed in Figure 4, in which a large band of cold material is left hanging in the upper mantle at 45 Ma following the subduction of the Izanagi-Pacific ridge at 55 Ma. This can be observed, too, in Figure 5, in which the bands of cold material observed throughout China represent the drips left over from old subduction zones.

The band of slab material left over from the Izanagi subduction has a significant effect on the asthenospheric flow studied here. The inflowing material from the Pacific domain is forced to push through the base of this high-viscosity barrier and largely flows at depths of ∼300–500 km (Supplementary Figures S2–S8). The blockage of this material from large parts of the low viscosity asthenosphere, as well as the intervening slab material, leads to a significant slow-down in flow velocity. Consequently, the inflowing Pacific material cannot flow far west enough during the late Cenozoic. A larger viscosity drop in the asthenosphere, making it compatible with the new constraints of ∼ 1 0 18 –10¹⁹ Pa⋅s (e.g., Chen et al., 2021), would also increase asthenospheric flow velocities (assuming the same pressure gradient) and may help to reconcile this issue. These limitations are highlighted by the low ∼3 cm/year velocity of this influx, when compared to recent inferences of ∼15 cm/year pressure-driven asthenospheric flow velocities (Hartley et al., 2011; Colli et al., 2014; Parnell-Turner et al., 2014; Chen et al., 2021). Taken together, it becomes clear that the implementation of refined viscosity profiles in the next generation of exascale-enabled mantle convection codes (e.g., TerraNeo, Bauer et al., 2020) could plausibly allow the influx of hot material reach further into East Asia during the late Cenozoic.

Another limitation of these models is the unrealistic slab morphology induced by two-sided subduction. For instance, the stagnant Pacific slab beneath East Asia (e.g., Ma J. et al., 2019) is not reproduced (Supplementary Figures S2–S8), which could have a significant effect on the predicted dynamic topography signal during the late Cenozoic (Supplementary Text S4). A variety of methods have aimed to reproduce this slab geometry using mantle circulation models (MCMs) (e.g., Mao and Zhong, 2018; Ma P. et al., 2019; Peng et al., 2021). However, these approaches have tended to come at the expense of significant additional approximations, such as the suppression of mantle upwellings. As such, it is not yet apparent whether current models can reproduce stagnant slab morphologies within an Earth-like convective planform. A conceptually simple method of generating more realistic slab behaviour is the use of a temperature-dependent viscosity coupled with a pseudo-plastic rheology (Moresi and Solomatov, 1998). The use of this rheology in MCMs produces more concentrated deformation at subduction zones, reducing the presence of “drip” artefacts, and inducing flat-slab subduction following fast trench retreat (Bello et al., 2015). Furthermore, the weak zones which arise between plates would presumably lead to more realistic slab window geometries in the uppermost mantle, which could prove vital when using MCMs to study slab window asthenospheric flow more generally. And so, while the use of a non-linear rheology comes with additional computational difficulties (e.g., Fraters et al., 2019), the self-consistent generation of realistic slab morphologies makes this an attractive option for future studies of this kind.

In any case, MCMs in general are fundamentally limited in their ability to reconstruct past mantle flow, as the initial state from which this flow is derived is unknown. This is partially overcome in MCMs by assuming an initial state, and assimilating tectonic histories as boundary conditions at the surface (Bunge et al., 2002). The models of Lin et al. (2020) illustrate nicely the utility of this technique when reconstructions are assimilated over a significant > 400 Ma timescale, as mantle upwellings are reproduced in approximately the right regions at approximately the right times. A significant improvement on this can be obtained by recasting the reconstruction of past mantle flow as an inverse problem aimed at finding an optimal initial condition via the adjoint method (Bunge et al., 2003). This approach requires as input the present-day mantle state inferred from global tomographic models, which may be used in conjunction with the adjoint equations for mantle flow (Bunge et al., 2003; Ismail-Zadeh et al., 2004; Horbach et al., 2014; Ghelichkhan and Bunge, 2016, 2018) to iteratively improve upon a “first guess” of the initial state. The assimilation of tectonic reconstructions is also required in this method to ensure convergence (Vynnytska and Bunge, 2015), and to avoid the divergent nature of chaotic mantle flow when retrodicting further back in time (Colli et al., 2015). One may therefore explicitly reconstruct the history of hot regions in the mantle inferred from tomographic models. East Asia has not been considered in any adjoint studies to date (e.g., Colli et al., 2018; Ghelichkhan et al., 2021), likely due to the significant uncertainties of the inputted tectonic histories. Following further refinement of these histories, however, the history of the seismically slow anomalies observed beneath East Asia—linked here to asthenospheric flow—may be retrodicted using these methods.

Our asthenospheric flow hypothesis has significant implications for the more southerly regions of Southeast Asia which were not considered here. Following the future model improvements discussed above, it becomes apparent that the influx of hot material studied here could flow into these region during the late Cenozoic. Due to its proximity to the Sunda Trench, Southeast Asia has traditionally been expected to be dominated by dynamic subsidence during the late Cenozoic (e.g., Wheeler and White, 2002). However, growing evidence implies the influence of hot material during this time, such as the dynamic uplift and basaltic magmatism in Borneo since the Neogene (Roberts et al., 2018), the Neogene uplift (Fyhn et al., 2009) and intraplate basalts (e.g., Yan Q. et al., 2018) in Indochina, and the seismically slow anomalies observed in the uppermost mantle throughout these regions (e.g., Schaeffer and Lebedev, 2013). While Roberts et al. (2018) found a Neogene–present-day swell of dynamic uplift around Borneo, and dynamic subsidence in the adjacent SCS, we note that the dynamic influence at the surface of this hot material will depend on the depth, wavelength, and lateral distribution of slabs in these regions (Colli et al., 2016). Future convection models which implement refined viscosity profiles may therefore reconcile the disparity between convection models (which often predict dynamic subsidence) and observations (which find dynamic uplift) highlighted by Roberts et al. (2018).

While we have focussed our analysis on a single ridge subduction event during the early Cenozoic, we note that the process proposed is generic. The hypothesis rests on the idea that hot asthenospheric material can be driven from high-pressure upwellings to low-pressure subduction zones by a combination of pressure- and plate-driven flow (Höink and Lenardic, 2010; Colli et al., 2018), which may then spill over into adjacent regions following ridge subduction events and the opening of slab windows. This mechanism may therefore be more widely applicable in future studies of dynamic topography and intraplate volcanism, as similar events could have occurred throughout Earth’s history. Indeed, this process has been proposed in order to explain the intraplate volcanism (Zhou et al., 2018) and dynamic uplift (Zhou and Liu, 2019) in the western United States during the late Cenozoic, following the opening of tears in the Juan de Fuca slab and the Juan de Fuca-Pacific slab window. Furthermore, the asthenospheric build-up stage of this process may be invoked at circum-Pacific margins in the present day, and evidenced by seismically slow anomalies beneath subducting slabs in Japan (Ma J. et al., 2019; Bogiatzis et al., 2019), Cascadia (Hawley et al., 2016; Bodmer et al., 2018, 2020), and South America (Portner et al., 2017).