The Dynamic Topography of Eastern China Since the Latest Jurassic Period

Some changes in the topography of eastern China since Late Jurassic times cannot be well explained by lithospheric deformation. Here we analyze global mantle flow models to investigate how mantle‐driven long‐wavelength topography may have contributed to shaping the surface topography of eastern China. Paleodrainage directions suggest that a southward tilted topography once existed in eastern north China in the latest Jurassic Period, which is different from that at present day (southeastward tilting). Our model dynamic topography reveals a southward tilting topography between 160 and 150 Ma, followed by southeastward tilting and rapid subsidence, which is compatible with paleodrainage directions and tectonic subsidence of the Ordos Basin. The Cretaceous anomalous subsidence of the Songliao and North Yellow Sea basins, as well as the Cenozoic anomalous subsidence of the East China Sea Shelf Basin, can also be explained by dynamic topography. An apatite fission track study in the Taihang Mountains reveals four stages of evolution: Late Jurassic fast unroofing, Cretaceous slow unroofing, early Cenozoic fast unroofing, and late Cenozoic slow unroofing. We propose that mantle flow influenced this surface unroofing because the model predicts Late Jurassic dynamic uplift, Cretaceous dynamic subsidence, early Cenozoic dynamic uplift, and late Cenozoic dynamic subsidence. Apatite fission track data from northern south China are also in reasonable agreement with predicted dynamic topography between 80 and 30 Ma. The spatial and temporal agreement between geological observations and model dynamic topography indicates that mantle flow has had a significant influence in shaping the surface topography of eastern China.


Introduction
The topography of eastern China is laterally variable, consisting of plains mostly below 500-m altitude (e.g., the Northeast Plain, the North China Plain, and the Middle-Lower Yangtze Plain) interspersed between mountains and foothills. Buried rift basins underlie most of these plains, reflecting past continental extension and diffuse deformation. The tectonic evolution of eastern China is complex. Subsequent to the final closure of the Paleo-Asian Ocean along the Solonker suture during late Permian times (Windley et al., 2007;Xiao et al., 2003), the Late Triassic amalgamation of the South China Craton (SCC) to the North China Craton (NCC) constitutes the final assembly of the major part of eastern China (S. Zhang et al., 1996). Significant intraplate deformation has reshaped the topography of eastern China since its assembly.
The NCC was stable for about 1.7 billion years since its final cratonization at 1.85 Ga (Zhao et al., 2005) until early Mesozoic deformation and magmatism, during which a continuous sequence of marine and terrigenous strata was deposited. The destruction of the eastern NCC during Late Mesozoic to Cenozoic times sets it apart from other cratons (Zhu et al., 2012). In contrast, the western NCC remained stable until the present day, and deformation was concentrated on its margin. Seismic tomography reveals a lithospheric keel exceeding 200 km in thickness (Chen et al., 2009). The destruction of the eastern NCC is marked by an intense deformation since Jurassic times (S. , in particular the extensive Jurassic folds and thrusts of the "Yanshanian Movement" (Wong, 1927). Recent studies indicated that this tectonic movement affected the whole of eastern China (S. Li et al., 2011), whereas it only affected the sedimentary cover in south China. In contrast, the metamorphic basement of the NCC is also deformed . Since Early Cretaceous times, compression transitioned to episodic extension, forming rifts or grabens in eastern China (Figure 1). This structural deformation is thought to have shaped the present-day complex topography of eastern China, as evidenced by the extensive distribution of folded mountains and fault-bounded basins. These local topography contrasts have been exhaustively studied. For example, Wu et al. (1999) mapped the existence of three-level planation surfaces in north China. The apatite fission track (AFT) method has also been used to analyze the uplift history of mountains in eastern China (e.g., Cao et al., 2015;Xu et al., 2016). In contrast, few studies have considered the influence of mantle flow on the whole region. To our knowledge, only Y. Wang, Huang, et al. (2016) computed some two-dimensional thermochemical convection models to explain the difference in heat flux and residual topography between the eastern and western NCC regions.
Generally, the difference of topography caused by isostatic equilibrium can be related to contrasts in the thickness and density of the crust and lithospheric mantle. The influence of mantle flow beneath the lithosphere on surface topography (Pekeris, 1935) is called dynamic topography. Given the limited constraints on structure and density of the continental lithosphere, dynamic topography is usually defined as the topography originating from sources beneath the upper thermal boundary layer of mantle convection (Flament et al., 2013). Anomalous topography changes that cannot be explained by isostatic effects or lithospheric deformation have been explained by dynamic topography. Examples include the development of elevated passive continental margins Müller et al., 2017), the Late Cretaceous emergence of Sundaland during sea level rise , and the asymmetric topography of the South Atlantic Ocean (Flament et al., 2014). Eastern China has been influenced by two adjacent ocean basins since Late Jurassic times: the Pacific domain and Tethyan domain (S. Li et al., 2011) and the fast convergence rate of the two domains have inevitably had a significant influence on mantle flow and dynamic topography.

Tectonics
Here we aim to study the influence of dynamic topography on the evolution of the topography of eastern China. We first compute the dynamic topography evolution of eastern China by coupling global plate reconstructions and mantle flow models. We then compare the predicted temporal and spatial evolution of dynamic topography to available geological observations in this region, including the Jurassic drainage direction of eastern NCC, the unroofing history of mountains, and the Mesozoic-Cenozoic subsidence history of major basins.

Plate Tectonic Motions Since the Latest Jurassic Period
Plate tectonic reconstructions are essential to understand geological processes in a regional context. Plate velocities and subduction zone locations from global tectonic reconstructions have recently been used as time-dependent boundary conditions of geodynamic models , helping understand the links between surface processes and the deep Earth.
Ophiolite complexes and other geologic records show the NCC and SCC assembled at the end of the Middle Triassic Period, which represents the final assembly of the major block of eastern China (Zhang et al., 1996). The tectonic evolution of eastern China since 160 Ma has been significantly influenced by the Pacific and Tethyan tectonic domains, as well as by the subduction of Mongol-Okhotsk Ocean seafloor to the north until its final closure at 150 Ma. The lack of preserved oceanic crust from the Neo-Tethyan  Zahirovic et al. (2016). Mid-ocean ridges are shown as black solid lines, subduction zones are shown as toothed black lines, and absolute plate velocities are shown as black arrows. White stars A-D show the locations used for sampling age of the subducting oceanic crust. Cross-sections NW1, NW2, and SN are anchored to the North China Craton. Continents of Archean age are in pink, Proterozoic in yellow, Phanerozoic in cyan, and inferred continental material in blue. Ocean floor is colored by age according to the color scale.

10.1029/2017TC004830
Tectonics Ocean and the Izanagi plate makes reconstructions uncertain. In this study, we implement two alternative kinematic models Zahirovic et al., 2016) as boundary conditions of mantle flow models. The reconstruction by Zahirovic et al. (2016) is based on that by Müller, Seton, et al. (2016), with some refinements to the eastern Tethyan domain from 160 Ma to present. In Müller, Seton, et al. (2016), the East Java, West Sulawesi, and eastern Borneo are placed on the New Guinea margin prior to breakup at~160 Ma, following Zahirovic et al. (2014). However, due to the poor preservation of seafloor spreading histories, Zahirovic et al. (2014) proposed an alternative scenario, where East Java, West Sulawesi, and eastern Borneo originate from the Argo Abyssal Plain on the NW Shelf of Australia, which is adopted in Zahirovic et al. (2016). In addition, the Sepik and Philippine Arcs are attached to New Guinea before moving northward in Zahirovic et al. (2016), and the north margin was an active margin character, which is consistent with the age and geochemistry of the Central Ophiolite Belt on New Guinea (Permana, 1998) and also the magmatic and paleomagnetic history of the Philippine Arc (Balmater et al., 2015;Deng et al., 2015). Another major change in Zahirovic et al. (2016) is that the peak of (ultra) high-pressure metamorphics in the Luk-Ulo suture zone on East Java at~115 Ma is reinterpreted to be the onset of subduction of the Woyla Sea backarc, rather than the previously interpreted collision of the East Java, West Sulawesi, and the Woyla intraoceanic arc in Müller, Seton, et al. (2016). In both reconstructions, relative plate motions are calculated based on preserved seafloor magnetic anomaly identifications and their associated rotations and finally form a hierarchical chain with all plate motions linked back to Africa. Africa is linked to a hotspot reference frame  back to 70 Ma and to a true polar wander-corrected paleomagnetic reference frame  for earlier times. A 10°longitudinal shift is applied to the true polar wander-corrected reference frame at 100 Ma, and the models transition linearly from the hotspot reference frame to the paleomagnetic reference frame between 100 and 70 Ma.
The past existence of the Izanagi plate is suggested by the M-sequence Japanese magnetic lineation set to the east of the Japan subduction zone, which records a westward younging Jurassic-Cretaceous spreading system (Woods & Davies, 1982). The westward subduction of the Izanagi plate was accompanied by the growth of the Pacific Plate to the east ( Figure 2). With the assumption of no major change in spreading rate, direction, and accretion from M0 (120 Ma) to the cessation of spreading along the Izanagi-Pacific ridge, the ridge is inferred to have intersected the east China continental margin at about 60-55 Ma in the reconstructions Zahirovic et al., 2016; Figure 3), which agrees with some onshore geological constraints from Japan and Korea (Seton et al., 2015). To the south, the now completely vanished Tethys Ocean basin existed between Gondwana and Laurasia before the opening of the Indian Ocean. The evidence for this ocean basin is primarily preserved in the terranes and ophiolite complexes along southern Eurasia and the Mediterranean Sea. Tethys seafloor moved northward and transferred the Gondwanaderived continental terranes and microcontinents toward Eurasia. The breakup between India, Africa, and Antarctica occurred during the Late Jurassic Period, and the separation of India and Madagascar at about 87 Ma constitutes the final breakup of Gondwanan continental blocks . The main phase of continent-continent collision between the Indian and Eurasian plates may have occurred from 47 Ma and finished by about 40 Ma .

Model Setup and Initial Condition
We solve the equations for incompressible mantle convection with the finite-element method in a spherical domain using CitcomS (Zhong et al., 2008). Surface plate velocities extracted from the above plate reconstructions are implemented in our mantle convection models as model boundary conditions. The age of the ocean floor from the reconstructions and the tectonothermal age of the continental lithosphere (Artemieva, 2006) are used to build the thermal structure of the lithosphere assuming a half-space cooling model. The subducting slabs are assimilated in the models to a maximum depth of 350 km (45°dip angle), below which convection is dynamic. The average model resolution is 50 × 50 × 15 km at the surface, 28 × 28 × 27 km at the core-mantle boundary (CMB), and 40 × 40 × 100 km in the midmantle. The viscosity depends on depth and temperature following: where η r is a depth-dependent prefactor defined with respect to the reference viscosity (η 0 ) for two layers (Table 1): 1 for mantle above 660 km (upper mantle) and 100 for mantle below 660 km (lower mantle). E η is the activation energy (E ηUM for the upper mantle and E ηLM for the lower mantle), T is the temperature, T b is the ambient temperature of the mantle, T η is a temperature offset, and R is the universal gas constant. We do not consider the activation volume in these Boussinesq (incompressible) models. The values used for activation energy (Table 1) are smaller than suggested by experimental values due to numerical limitations. The activation energy and the temperature offset are adjusted to obtain viscosity variations over 3 orders of magnitude across the range of temperatures considered. These lateral viscosity contrasts, which are lower than expected to occur within the solid Earth (e.g., Stadler et al., 2010), can be computed with a resolution at which the computation of time-dependent mantle flow simulations is possible.
The Rayleigh number is defined as where α 0 is the coefficient of thermal expansivity, ρ 0 the density, g 0 the acceleration of gravity, ΔT the temperature difference between surface and CMB, h M is the depth of the mantle, κ 0 the thermal diffusivity, η 0 the viscosity, and the subscript "0" indicates reference values (Table 1). Thermal expansivity is assumed to be constant. The CMB temperature (3100 K) is lower than estimated (> 3300 K, Lay et al., 2008). A relatively low CMB temperature is acceptable for incompressible models in which the adiabatic gradient is not considered.
Our calculations start at 230 Ma to make sure the models reach a dynamic equilibrium (which takes~50 Ma; Flament et al., 2014) before 160 Ma. The initial condition consists of a global temperature field (T b = 1685 K) with cold slabs inserted to 1,400-km depth, with a dip of 45°above 425 km and 90°below ( Figure S1 in the supporting information). In addition, slabs are initially twice as thick in the lower mantle than in the upper mantle, to account for the greater viscosity of the lower mantle. An isothermal boundary condition is applied both on the surface (T = 273 K) and at the CMB (T = 3100 K). In order to obtain long wavelength dynamic topography, the initial condition includes a basal layer to suppress the formation of mantle plumes because their temperature is too large in incompressible convection models and isolate the contribution of subduction zones to longwavelength dynamic topography. Three factors contribute to the absence of plumes in the setup: the dense layer is 4.2% denser than the ambient mantle (corresponding to buoyancy ratio B = δρ ch /(ρ 0 α 0 ΔT) = 0.5, whereρ ch is the density contrast), its initial thickness is taken equal to the thickness of the basal thermal boundary layer (laterally uniform and 113-km thick corresponding to 2% of the volume of the mantle, which is appropriate for Large Low Shear Velocity Provinces, Burke et al., 2008;Hernlund & Houser, 2008), and the model run time (230 Myr) is relatively short. Plumes occur when the ratio of compositional to thermal buoyancy is smaller (e.g., Zhong & Hager, 2003), for instance, if the thermal boundary layer is thicker than the basal dense layer (e.g., Hassan et al., 2015).
We consider two mantle flow models: Model 1 uses the reconstruction of Zahirovic et al. (2016) as timedependent boundary conditions, whereas Model 2 uses the reconstruction of Müller, Seton, et al. (2016). All other parameters are identical between the two models (Table 1).

Computation of Dynamic Topography
The dynamic topography is obtained by scaling the total normal stress σ rr on an air-loaded top surface, using no-slip boundary conditions. We ignore buoyancy and lateral viscosity variations above 250-km depth, which is the base of Archean mantle lithosphere in the initial condition of the models (Flament et al., 2014), to isolate the topography arising from mantle flow.

10.1029/2017TC004830
Tectonics The time-dependent dynamic topography at the surface is computed in 10-Myr intervals following where Δρ is the density difference between air (ρ a = 0 kg/m 3 ) and the shallow mantle ρ sm = 3,340 kg/m 3 , because surface dynamic topography is most sensitive to shallow sources of buoyancy (e.g., Hager & Clayton, 1989). Using the density of the shallow mantle increases the amplitude of dynamic topography, compared to using the reference density of the mantle ρ 0 = 4,000 kg/m 3 (Table 1).

Low-Temperature Thermochronology
AFT thermochronology has been widely used to determine the uplift and denudation rates of orogenic belts, to unravel the thermal history of sedimentary basins, to identify sedimentary provenance, and to date volcanic rocks (e.g., Donelick et al., 2005;Fitzgerald et al., 1995;Flowers et al., 2012). AFT thermochronology is based on the accumulation of radiation damage in uranium-bearing minerals due to the spontaneous nuclear fission of 238 U (Fleischer & Price, 1963). The damage (fission track) is gradually annealed under elevated temperatures and can even be reset to zero (Crowley, 1985). The annealing and growing processes of AFT can be used to reconstruct cooling trajectories and to quantitatively constrain exhumation histories.
In this study, we use recently obtained data from the Taihang Mountains in the central NCC (Cao et al., 2015) and in the northern SCC (Wang et al., 2015). The test of the apatite was completed at the China University of Geosciences (Beijing), using standard external detector method and the ζ age calibration approach (Hurford & Green, 1983). Ages were calculated using the Zeta calibration method (Hurford & Green, 1983). The cooling histories were modeled using HeFTy (Ketcham, 2005). We directly compare the burial and unroofing history revealed by the AFT study with temporal patterns of dynamic uplift and subsidence (section 4.2), to analyze the influence of mantle flow on surface topography.

Modeling Tectonic Subsidence
Subsidence analysis tracks the subsidence and uplift history of a sedimentary basin since deposition (McKenzie, 1978). The tectonic subsidence curve is the idealized subsidence history of a basin that would have existed if only water, and no sediment, filled the subsiding accommodation space. It provides a way of normalizing subsidence in different basins that have undergone very different sedimentation histories. Extensional events are characterized by significant lithosphere thinning, rapid tectonic subsidence, positive thermal anomaly, and active magmatism. Once rifting ceases, the thermal anomaly decays due to thermal reequilibration, as the mantle beneath the rift cools which results in postrift subsidence that slows according to a half-space cooling model (Jarvis & McKenzie, 1980;White, 1993). The total subsidence is directly related to the amount of thinning during the rifting, unless dynamic topographic change postdating the onset of rifting plays a significant role.
Eastern China has undergone several episodes of rifting since the Cretaceous Period, which have affected its surface topography. Most basins in eastern China are significantly controlled by faulting. Well data and some other data (e.g., seismic profiles) constrain the evolution histories of these basins (e.g., prerift stage, synrift stage, and postrift stage). We forward modeled the water-loaded tectonic subsidence history of major basins using the approach of White (1993), assuming alternative stretching factors for synrift stages. The stretching factor that best matches the onset of the synrift and postrift stages was used to calculate the inferred tectonic subsidence history. We defined anomalous subsidence for a given well as the difference between subsidence history derived from observed data and that derived from the theoretical tectonic subsidence (e.g., Xie et al., 2006). The anomalous subsidence is unrelated to the thermal cooling of the lithosphere during the postrift stages. This anomalous subsidence can be compared with predicted dynamic topography. Observed subsidence greater than predicted by a model suggests dynamic subsidence, whereas observed subsidence smaller than predicted by a model suggests dynamic uplift.

Evolution of Mantle Temperature
We consider mantle 10% colder than ambient temperature to be subducted slabs. A cross section of the temperature field shows the subduction processes of the (proto-) Pacific plate and the subduction-induced negative dynamic topography behind the trench (Figures 4 and S2). During Late Mesozoic times, the subducting

10.1029/2017TC004830
Tectonics slab in the model is the Izanagi plate, which continued subducting until the beginning of Cenozoic times (Figures 4 and S2). At about 60-55 Ma the Izanagi-Pacific spreading mid-oceanic ridge intersected the trench as shown by the very thin oceanic lithosphere at 60 Ma (Figures 4 and S2), forming a slab window beneath East Asia. The slab window beneath can still be recognized at 30 Ma. The continuous northward movement of the Philippine Sea Plate led to the formation of double subducting configuration since approximately 30 Ma, which lasted until present day (Figures 4 and S2).

Evolution of Dynamic Topography
The dynamic topography at 160 Ma in eastern China is characterized by a broad negative dynamic low, as low as À2,000 m in the south and À600 m in the north (Figures 5 and S3). The larger amplitude in the south is due to the combined effects from two subducting oceans, the Pacific and Tethys Oceans. The dynamic topography is E-W trending, southward dipping at 160 Ma, and is about À600 to À800 m in the NCC, À800 to À1,200 m in the SCC, then slowly rotates anticlockwise to become trending NE, dipping southeastward by  Figure 2). Gray dashed contours denote mantle 10% colder than ambient. Magenta lines on top of the section are air-loaded dynamic topography (black lines show mean dynamic topography, which is by definition equal to 0), with vertical scale exaggerated by a factor of 333. The numbers above the color scale denote nondimensional temperatures, and the number below the color scale denote dimensional temperatures. The black solid lines at 660-km depth denote the upper-lower mantle boundary; the red solid lines at 250-km depth denote that buoyancy is ignored above.

10.1029/2017TC004830
Tectonics 120 Ma (Figures 5 and S3). The dynamic subsidence continues to 80 Ma (Figures 5 and S3). Eastern China underwent rapid dynamic subsidence between 160 and 80 Ma, by about 500 m in eastern NCC and northeastern SCC (to about À1,300 m). The dynamic southeastward tilting reflects the increasing impact of the subduction of the Izanagi plate. Since 80 Ma, the dynamic topography has been characterized by two NE-SW trending dynamic topography lows, the western one covering most of the NCC and SCC and the eastern one immediately behind the subduction zone, covering the Yellow Sea and the East China Sea (Figures 5 and S3). The eastern dynamic low migrates eastward due to the retreat of the subduction zone since 40 Ma (Figure 1), accompanied by depocenter migration and tectonic deformation . The variation in amplitude along subduction zones reflects the dependence of slab buoyancy on the age of the subducting lithosphere (Flament et al., 2014). The amplitude is very large between 160 and 90 Ma because old lithosphere was subducting during that period, and the amplitude becomes smaller afterward until about 60 Ma when the Izanagi-Pacific mid-oceanic spreading ridge intersected the convergent margin along eastern China, at which time the dynamic topography reached its least negative value (Figures 4  and S2). Dynamic topography then became more negative again until the present day.

Comparison to Seismic Tomography
In this section we qualitatively compare the present-day thermal structure predicted by the mantle flow models to two tomography models. Positive seismic velocity anomalies imaged by seismic tomography reveal subducting slabs beneath East Asia extending from oceanic trenches into the upper mantle (Huang & Zhao, 2006;X. Liu et al., 2017). The resolution of P wave tomography models is generally better beneath continents than beneath oceans. In contrast, S wave tomography models image equally the mantle beneath oceans and continents, despite their generally lower resolution (Grand, 2002;Widiyantoro et al., 1998). We carried out a qualitative comparison of predicted present-day mantle temperature and tomography models by overlying the contours of mantle 10% colder than ambient on two tomography models, the P wave model (C. Li et al., 2008) and the S wave model of Grand (2002; Figure 6). The modeled cold anomalies predicted for both flow models match the first-order shape and depth of positive seismic velocity anomalies in the tomography models. The subducted Izanagi plate at longitude 110°E and between 1,000-and 2,000-km depth in section NW1 is imaged by both tomography models. Model 1 shows a better agreement with tomography than Model 2. There are two subducting slabs along section NW2, corresponding to the Pacific plate in the east and the Philippine Sea Plate in the west. Both mantle flow models predict a similar mantle structure, which is in reasonable agreement with tomographic images. Figure 6. P wave (C. Li et al., 2008) and S wave (Grand, 2002) tomography models along sections NW1 and NW2 (see Figure 2 for location). Magenta contours denote mantle 10% colder than surrounding temperature from Model 1, green from Model 2. Paleogeography is an important constraint for geodynamic setting and tectonic evolution. Provenance analysis of terrigenous sediments can be effectively used to assess paleogeographic reconstructions (Dickinson, 1988). Successive sequences of the Middle-Late Mesozoic Period are well preserved and exposed in the Luxi Uplift (located in inner eastern NCC), which are surrounded by Cenozoic rift basins. Through the analysis of eHf(t) values, combined with zircon dating and sandstone detrital modes, Xu et al. (2013) proposed that the Late Jurassic sediments in the Luxi Uplift were mainly derived from the mixed provenance of primarily from the northern NCC and secondarily from the southern Central Asian Orogenic Belt, implying that a provenance region of elevation higher than the NCC once developed along the northern margin of the NCC during Early to Late Jurassic times. Xu et al. (2013) reconstructed a southward flowing Late Jurassic drainage network that shed eroded materials into the inner NCC. A north-south dynamic topography profile crossing the southern Central Asian Orogenic Belt and Luxi Uplift shows the southward tilting of the topography during latest Jurassic times (Figure 7), which is compatible with the paleodrainage flowing direction inferred from the geological record. Due to the poor preservation of Cretaceous strata in this area (with only thin strata deposited in little rift basins), the evolution of this drainage is less constrained.

Comparing Dynamic Topography With AFT Thermochronology Results
In order to compare the burial and unroofing history inferred from two recent AFT studies (Cao et al., 2015;Wang et al., 2015) with the temporally and spatially dependent dynamic topography, we plot them together over the same time interval in Figure 8. We first analyze the result of the Taihang Mountains located in the central NCC then the Chencai area in northern SCC.

.1. Taihang Mountains
The Taihang Mountains are a north-south trending mountain range located in the central NCC, belonging to Trans-North China Orogen according to a tectonic subdivision of the Precambrian basement (Zhao et al., 2005). In an earlier study, we carried out AFT sampling, analysis, and modeling along the middle and south parts of the Taihang mountains (Figure 1; Cao et al., 2015). The modeled burial depth-time paths from two samples can be divided into four stages (Figures 8a and 8b): (1) rapid unroofing before 140 Ma (the onset time of the unroofing is unknown because AFT thermochronology has no constraint to temperature higher than 110°C); (2) a slow-to-no unroofing stage occurred until 65 Ma; (3) another rapid unroofing event was recorded until about 50 Ma followed by relative slow unroofing until about 35 Ma; and (4) the paths of the two samples are different after 35 Ma, one recording rapid unroofing until the present day ( Figure 8a) and the other showing almost no unroofing (Figure 8b). The first sample (thm188) was collected on the east flank of the Taihang Mountains and the second (thm206) on the west flank. The east flank sample might be influenced by the postrift subsidence of the Bohai Bay Basin since the early Miocene (Cao et al., 2015), which has little influence on the west flank sample that is located further away from the basin. As a consequence the AFT paths of the west flank sample better show the influence of mantle flow, which are selected to compare with dynamic topography.
Similarly, the predicted vertical motion history of the Taihang Mountains can also be divided into four stages, showing (1)  The dynamic uplift predicted for stages 1 and 3 corresponds to the rapid unroofing revealed by AFT data, and the dynamic subsidence predicted for stages 2 and 4 correspond to slow-to-no unroofing in the AFT data. The elevation of the samples (from mountains) should still be higher than the surrounding area during dynamic subsidence periods, so that little sedimentation and hence little change in AFT burial depths are expected during stages 2 and 4 (Figures 8a and 8b). The transition from dynamic uplift to dynamic subsidence predicted by Model 1 (Model 2) is 20 Myr (10 Myr) earlier than cessation of the first rapid unroofing record by AFT data from stages 1 to 2 and 5 Myr earlier for the transition from stages 2 to 3. The possible reasons for this offset are discussed in section 4.4. Nevertheless, the predicted evolution of dynamic topography is compatible with the AFT results.
The magnitude of unroofing history is several times larger than that of dynamic topography. Note that the change of burial depth is not a perfect proxy for the amplitude of subsidence or uplift on the surface (Flowers et al., 2012), because crustal unloading during unroofing enhances the uplift history. In the absence of direct relationship between dynamic topography and unroofing, the correspondence of timing of rapid unroofing with dynamic uplift and slow-to-no unroofing with dynamic subsidence is encouraging.

Chencai Area
We selected two samples 433 and 438 in the Chencai area (Wang et al., 2015), northern SCC. The thermal history model results showed similar cooling histories, revealing fast unroofing during 80-55 Ma, slow unroofing between 55 and 30 Ma, renewed rapid unroofing between 30 and 20 Ma, followed by slow-to-no unroofing afterwards. The dynamic topography shows rapid dynamic uplift between 80 and 50 Ma followed by slower dynamic uplift between 50 and 40 Ma, which is compatible with unroofing history, then slow dynamic subsidence from 40 Ma to the present day. The dynamic subsidence during 40-30 Ma and 20-0 Ma is compatible with the slow unroofing history, while it is at odds with rapid unroofing between 30 and 20 Ma (Figures 8c and 8d). The rapid unroofing might be related to the thermal subsidence of the East China Sea Shelf Basin (C. Li, Zhou, et al., 2009) or could have been caused by structural deformation that is not considered in our models.

Comparing Dynamic Topography With Tectonic Subsidence
The wide emplacement of rift basins in eastern China reflects episodic extension since Cretaceous times, which has contributed to the evolution of regional surface topography. We investigated the influence of mantle flow on the basins by comparing anomalous subsidence with dynamic topography evolution. We compiled available tectonic subsidence histories from some basins in eastern China, including the Ordos Basin, the Songliao Basin, the North Yellow Sea Basin, and the East China Sea Shelf Basin. The Ordos Basin experienced no structural deformation, while the other three basins underwent rifting during Meso-Cenozoic times (Figure 9).

. Anomalous Vertical Motion of the Ordos Basin
The Ordos Basin lies in the western NCC. The basin is underlain by a thick (>200 km) lithospheric mantle root (Chen et al., 2009) and has remained structurally stable and "cratonic" since the amalgamation of the NCC at 1.85 Ga (Zhao et al., 2005). This basin has experienced episodic uplift and subsidence, which resulted in the deposition of sediments of several thousands to more than ten thousands of meters, including early Paleozoic carbonate platforms, late Paleozoic marine to continental strata, Mesozoic lacustrine inland strata, and thick Cenozoic loess. The structure within the basin is rather simple, with little to no deformation. However, the basin is surrounded by structural belts with different structural styles and orientations (Zhang et al., 2007).
The subsidence history is compiled from the southwestern Ordos Basin (Xie & Heller, 2009); however, the exact location of the backstripped well is unknown to us. Dynamic topography is expected to evolve similarly across this area, given its long-wavelength spatial nature. We thus extracted the dynamic topography from a point near the central southwestern Ordos Basin to compare with tectonic subsidence. We plotted the time-dependent dynamic topography curve together with the subsidence curve over the same time interval (Figure 9a). In order to compare the trends and magnitude between the curves through time, all of them are shifted to 0 m at 80 Ma (a time that best shows the similarity or difference between the curves). The subsidence curve shows that the basin underwent little subsidence between 160 and 148 Ma, followed by an increase in subsidence rate until 100 Ma, at which time the subsidence rate declined slightly until 70 Ma. The basin then experienced long-term erosion until Pliocene times, since when thick loess deposited, so the subsidence curve only lasts to 70 Ma. The basin underwent 300 m of subsidence between 148 and 70 Ma. The dynamic topography predicted by Model 1 is in good agreement with the subsidence curve in both time and magnitude for most of the time period, except for time periods 165-160 Ma and 80-70 Ma. Model 2 is similar to Model 1, except that it overestimates the dynamic subsidence by about 250 m.

. Anomalous Postrift Subsidence of the Songliao Basin
The Songliao Basin is located in Northeast China. The basement of the Songliao Basin is composed of several microplates initially assembled during the early Paleozoic and then welded together with the NCC by the Solonker Suture Zone at the end of the Permian Period . We compiled one water-loaded subsidence curve calculated from three borehole successions for the basin (P. Wang, Mattern, et al., 2016). The subsidence curve indicates that the synrift stage lasted between 150 and 103 Ma, with fault-bounded grabens filled by volcanogenic successions and sediments. The following postrift stage from 103 to 54 Ma was characterized by an abnormally high subsidence rate, accompanied by flood basalt eruption. Cenozoic strata are almost absent from the basin that mostly contains Quaternary sediments, so the subsidence curve stops at 55 Ma. As documented in section 2.5, we calculated the amount of anomalous subsidence by subtracting the theoretical thermal subsidence derived from forward models from that derived from the well data. We assumed continental crust initially 35-km thick, similar to the surrounding area with little deformation since 160 Ma revealed by receiver function imaging (Zheng et al., 2007). We varied stretching factors between 1.2 and 1.8 during the rifting phase between 150 and 103 Ma in the model, and a factor of 1.31 best fitting the observed data was used to calculate anomalous subsidence.
We compared the anomalous subsidence curve to the predicted dynamic topography (Figure 9b), with all curves shifted to 0 m at 54 Ma, as in the Ordos Basin. A phase of rapid anomalous subsidence by more than 400 m is identified between 103 and 80 Ma. After a rapid transient change of topography between 80 and 79 Ma, about 20 m of anomalous uplift occurred until 72 Ma, then 50 m of anomalous subsidence until 64 Ma, followed by anomalous uplift lasting until about 54 Ma. Both Models 1 and 2 agree temporally with the rapid anomalous subsidence from 103 to 80 Ma, the cessation in rapid anomalous subsidence by about 80 Ma, and the gentle uplift from 64 to 55 Ma. Model 1 better matches the magnitude than Model 2. The predicted dynamic topography is compatible with the anomalous subsidence for the Songliao Basin. However, both models are at odds with anomalous subsidence between 80 and 64 Ma for which there must be an explanation other than deep-mantle dynamic topography.

Anomalous Vertical Motion of the North Yellow Sea Basin
The North Yellow Sea Basin is located in the southeast part of the NCC. The basin first formed during Late Jurassic times. Seismic profiles indicate that the basin experienced intense rifting during Paleocene-Oligocene times and thermal subsidence in the Miocene-Neogene times. The tectonic nature of the basin during Jurassic-Cretaceous times is debated, with some workers proposing it was a rift basin with active faults (Cao et al., 2007), while others proposed based on seismic profiles that these faults had no influence on thickness of strata (Chen et al., 2008). Here we compare three subsidence curves from three sags of the basin (W. Li, Zeng, et al., 2009) to dynamic topography (Figure 9c). We extract representative dynamic topography at the center of the basin, because long-wavelength model dynamic topography does not vary much across this relatively small basin. Both observed data and dynamic topography are shifted together at 59 Ma, to focus on the prerift phase (rifting occurred between 59 and 23 Ma). Well data reveal that the basin experienced a phase of tectonic subsidence between 160 and 100 Ma followed by uplift at 100-59 Ma. The predicted evolution of dynamic topography is compatible with anomalous subsidence of the North Yellow Sea Basin for the period 160-59 Ma. The basin experienced rifting between 59 and 23 Ma, followed by structural inversion and then subsidence until present day.

Anomalous Postrift Subsidence of the East China Sea Shelf Basin
The east China continental margin was stretched during early Cenozoic times, forming the China Sea Shelf Basin. We compiled two well data sets from the west part of the basin (Figure 1; C. Li, Zhou, et al., 2009). The well data sets indicate that the synrift stage lasted between 65 and 40 Ma then the rifting ceased and transitioned to subsidence. We assumed continental crust initially 30-km thick, similar to that of the coastal area (Q. Li et al., 2013) to the west with little deformation since the Late Mesozoic. A stretching factor of 4.2 for the synrift stage (65-40 Ma) best fitting the observed data was used to model subsidence.
We plotted the dynamic topography together with the anomalous subsidence with respect to the present day over the same time interval (Figure 9d). Backstripped well h1 reveals an anomalous uplift of about 130 m between 40 and 23 Ma, then subsidence until 17 Ma, followed by a transient uplift until 12 Ma, returning to continued subsidence until the present day. The postrift subsidence history of well c5, located at the margin of a Paleogene-Oligocene sag, suggests a different history, starting with a very short uplift phase immediately followed by a rapid subsidence until 34 Ma, then a phase of anomalous uplift until 6 Ma, transitioning to subsidence until present day. The predicted dynamic topography for the two wells is similar for each model due to their close proximity. The predicted vertical motion for Model 1 reproduces the anomalous subsidence of well h1, with a smaller magnitude and simpler history. Model 2 is at odds with the anomalous subsidence of well h1 between 40 and 31 Ma but agrees reasonably with the later stage. Both models differ from well c5 data, possibly because the well is located at the margin of the rift, and recorded a nonrepresentative stratigraphic history.

Model Uncertainties and Mismatches Between Geologic Observations and Model
Dynamic Topography 4.4.1. Uncertainties of the Models We considered two different plate reconstructions because past plate motions are increasingly uncertain back in time (e.g., Torsvik et al., 2010). The difference in the predicted dynamic topography (compare Figures 5 and S3) reflects that model results are sensitive to imposed plate motions since relatively small differences in tectonic reconstructions can lead to locally significant differences in model predictions (in this case in the northwest part of the study area). The onshore geological record can provide additional constraints for older times due to the poor preservation of oceanic crust. Zahirovic et al. (2016) refined the eastern Tethys based on analysis of the volcanism, ophiolite obduction events, which improved the fit of Model 1 to these geological observations. The considered reconstructions do not account for lithospheric deformation, which may influence the surface velocities and the locations of subduction zones to some extent (e.g., Flament et al., 2014). This effect should be small compared to the absolute velocities of the plate motions. In addition, we do not consider the effect of asthenospheric flow since the presented dynamic topography is derived from sources deeper than 250 km. We note that small-scale upper mantle convection is limited in our time-dependent global mantle flow models in which lateral viscosity variations are limited to 3 orders of magnitude. In addition, considering the depth dependence of thermal expansivity (e.g., Tosi et al., 2013) or the compressibility of the mantle would reduce the effect of deep (lower mantle) sources of thermal buoyancy on surface topography. The initial condition of the models is also uncertain. We inserted slabs down to 1,400-km depth in the initial condition ( Figure S1). Flament et al. (2017) showed that inserting slabs to more than 800-km depth in the initial condition better reproduces the present-day structure of the lower mantle and that results are similar when slabs are initially inserted to more than 1,100-km depth.
Anomalous subsidence for a given well is the difference between subsidence history derived from observed data and that derived from the theoretical tectonic subsidence. The theoretical tectonic subsidence is calculated with assumptions of the initial continental crust thickness and onset and end time of synrift stages. The duration of synrift stages can be well constrained from well data and other geological records, while the crustal thickness prior to rifting is poorly constrained. The crustal thickness we used for the Songliao Basin and the East China Sea Shelf Basin are from the surrounding areas with little deformation since the Late Mesozoic times. However, the result are not very sensitive to the initial crustal thickness. Taking the Songliao Basin as an example, the theoretical subsidence between 150 and 50 Ma increases by only 40 m if the initial crustal thickness is increased by 5 km (from 35 to 40 km). Another assumption we made is that the postrift subsidence of a basin is only affected by cooling of the lithosphere and dynamic topography, which may be an oversimplification (Flament et al., 2013). Other factors, such as transient structural deformation, may contribute to anomalous subsidence, even though deformation is usually limited during postrift stage.
Because of information loss, limited precision of measurements, and some other factors, more than one timetemperature history is consistent with a given ending condition, so the cooling paths modeled using the software HeFTy (Ketcham, 2005) for a given sample consist of a set of thermal histories that are consistent with the measured data, including many good fit paths (represented by the magenta envelope in Figure 8) and a best fit path. Cao et al. (2015) assumed the time-temperature histories for all samples to be monotonic, because the samples are collected from mountains with little sedimentation since the latest Jurassic. The uncertainty of the cooling paths becomes larger for older times, which is reflected by the broadening of the "good fit" envelope back in time. Following Flowers et al. (2012), we compared the trend rather than the magnitude of the dynamic topography and cooling paths of the samples (see section 4.2.1), and we note that although absolute values are uncertain, the trends are resolved beyond uncertainty in the thermochronology and are consistent across mantle flow model cases (Figure 8).

Tectonics
The model setup suppresses the occurrence of active mantle plumes. This assumption is reasonable considering that no plume-related magmatism has occurred in eastern China since the latest Jurassic Period. The extensive rifting of eastern China during Early Cretaceous times led to the emplacement of felsic rocks and intermediate-felsic rocks . Lithospheric thinning caused by intraplate rifting (X. Li, 2000), thermomechanical and chemical erosion (Xu, 2001), and lithospheric delamination  has been proposed to explain this crustal melting and magmatic underplating event. The rare occurrence of mafic rocks suggests that mantle flow under eastern China is largely dominated by the subduction of the surrounding oceanic plates and that active mantle upwelling has not been significant since the Jurassic Period under eastern China. 4.4.2. Mismatches Between the Amplitude, Temporal and Spatial Evolution of Geologic Observations, and Model Dynamic Topography The above comparisons reveal that there are some mismatches between predicted dynamic topography and geological observations. As an example, the transition of dynamic topography from stages 1 (uplift) to 2 (subsidence) is 20 Myr earlier than the cessation of the first rapid unroofing recorded by AFT data for the Taihang Mountains ( Figure 8b) and 5 Myr earlier for the transition from stages 2 (subsidence) to 3 (uplift). The transitions of dynamic topography possibly occur earlier than unroofing history because mantle flow is insensitive to lithospheric and surface processes. The first temporal offset could be due to the ongoing erosion of the topography created by the early phase of rapid uplift, and the second temporal offset could reflect ongoing slow-to-no erosion until an elevation difference is created. The temporal offsets might also reflect the decreasing resolution of AFT data and plate motions for older times, as documented in section 4.4.1.
There are also mismatches between dynamic topography and anomalous subsidence. The dynamic topography is of smaller amplitude than the anomalous subsidence for the Songliao Basin and East China Sea Shelf Basin. This might be related to the uncertainties of our models, particularly the plate reconstructions incorporated in the mantle flow models that ignore lithospheric deformation. In addition, the dynamic topography misses short-lived subsidence changes for the two basins, such as the fast anomalous subsidence and uplift at 80-79 Ma in the Songliao Basin. This is not unexpected because long-wavelength topography evolves over millions of years (e.g., Flament et al., 2013), implying that there might be other factors (e.g., transient structural deformation) influencing the subsidence process. Despite these mismatches and limitations, the evolution of long-wavelength dynamic topography is in reasonable first-order agreement with geological observations covering a large area.

The Different Influence of Lithospheric Deformation and Mantle Flow on Surface Topography
The topography of eastern China has been significantly influenced by episodic intraplate deformation since Jurassic times. This complex deformation makes it difficult to investigate the influence of mantle flow on surface topography. A distinct difference between the topography controlled by deformation and mantle flow is its wavelength. Subduction-driven, long-wavelength dynamic topography usually occurs at continental scale, with a wavelength ranging from thousands to tens of thousands of kilometers, while lithospheric deformation usually controls topography more locally, at the scale of a basin or a mountain. The amplitude of isostatic topography associated with lithospheric deformation is typically several times larger than that of longwavelength dynamic topography (~± 1 km).
The above comparisons show that from latest Jurassic times (160-150 Ma) until the end of Early Cretaceous times (90-80 Ma), the Ordos Basin underwent 300 m of subsidence. To the east, the Taihang Mountains experienced slow-to-no unroofing during this time period. Further east, the North Yellow Sea Basin underwent rapid subsidence by 600 m. To the north, the Songliao Basin also underwent 400 m of subsidence between 100 and 80 Ma. This simultaneous change of topography by no more than 1 km over a large region suggests an influence of mantle flow as opposed to structural deformation, and it is compatible with our model dynamic topography predicting a rapid subsidence between 160 and 80 Ma (section 3.2). Anomalous subsidence then reverted to uplift at about 90 Ma (lasting to 60 Ma with an amplitude of 200 m) in the North Yellow Sea basin, 80 Ma (fluctuating to 54 Ma with an amplitude <200 m) in the Songliao Basin, fast unroofing during 80-60 Ma in the Chencai area, and slow-to-no unroofing transitioned to rapid unroofing at 65 Ma in the Taihang Mountain. The westward younging of the transition from 10.1029/2017TC004830 Tectonics anomalous subsidence to small amplitude uplift is compatible with the subduction of younger oceanic lithosphere from east to west and thus large-scale mantle flow. The westward depocenter migration in the Songliao Basin between 80 and 64 Ma (Wang et al., 2013) could have been caused by~300 m of dynamic uplift between 80 and 60 Ma (20 m/Myr; Figure S4) in the eastern part of the basin due to the subduction of increasingly younger oceanic lithosphere (Figures 3 and 4) and the opening of a slab window. The general absence of Upper Cretaceous strata in eastern China might also reflect this uplift episode. In addition, the coherent eastward motion of the NCC and the Izanagi subduction zone (Figure 1) relative to the mantle between 80 and 60 Ma ( Figure S5) caused the NCC to override the sinking Izanagi slab and subsidence to propagate westward (Figures 5 and S4). During early Cenozoic times, most basins in the region experienced rifting, while unroofing simultaneously occurred in the Taihang Mountains and in the Chencai area. This again suggests the influence of large-scale mantle convection on surface topography. In summary, the spatial and temporal agreement between predicted dynamic topography and geological observations indicates that mantle flow played an important role in shaping the long-wavelength topography of eastern China.

Conclusion
Significant intraplate deformation has reshaped the topography of eastern China since the Late Jurassic Period, while analyses of paleodrainage, anomalous subsidence of major basins, and unroofing of mountains indicate large-scale topography change caused by mantle flow, induced by the time dependence of rapid convergence of the Pacific and Tethyan domains. The Late Jurassic drainage indicates southward tilted topography, which can be related to fast subduction of Tethys seafloor. Then the topography tilted southeastward due to the increasing influence of the Izanagi Plate, accompanied by the simultaneous subsidence of the Ordos Basin, Songliao Basin, and North Yellow Sea Basin between 160 and 90 Ma with small amplitudes, which is in reasonable agreement with model dynamic topography. Subsequently, the transition from subsidence to uplift in these basins with a westward younging trend from 90 Ma in the North Yellow Sea Basin to 65 Ma on Taihang Mountains and the westward migration of depocenter in the Songliao Basin between 80 and 64 Ma may have been driven by rebound due to the westward subduction of younger ocean floor. During early Cenozoic times, unroofing simultaneously happening in the Taihang Mountains and in the Chencai area matches the uplift of dynamic topography. During late Cenozoic times, mantle flow caused enhanced subsidence of the East China Sea Shelf Basin. Due to the relatively slow change in mantle flowdriven dynamic topography through time and the difficulty in considering asthenospheric flow and surface processes, the model dynamic topography cannot account for short-lived subsidence rate changes and underestimates the subsidence of eastern China. Nevertheless, the evolution of long-wavelength dynamic topography is compatible with regional geological data.