Exhumation History of the Greater Khingan Mountains (NE China) Since the Late Mesozoic: Implications for the Tectonic Regime Change of Northeast Asia

The Northeast Asia region showcases the NNE (north-northeast) trending Greater Khingan Mountain (GKM) Range, which exhibits a conspicuous gradient zone of gravity anomaly and a rapid decline in crustal thickness from west to east [1]. The GKMs in their present form are predominantly the result of tectonic tensions transmitted from the subduction of both the Mongol-Okhotsk Ocean located in the northern part and the Pacific plate situated in the eastern portion since the Early Mesozoic [2-7]. The Mongol-Okhotsk Ocean is presumed to have been blocked with double-sided subduction, indicating its eventual closure during either the Jurassic or Cretaceous period [8-13]. The Mesozoic period was witness to the Paleo-Pacific plate’s complex geodynamic regime that incorporated processes involving slab roll-back, initiating the present-day Pacific slab’s subduction, and the formation of a colossal mantle wedge in the Cenozoic <55 Ma [14]. The transition time of these two dynamic systems in Northeast China and the corresponding superficial mechanism has not been clearly defined.

The GKMs record far-field responses from plate-boundary processes generated by the evolution of the Mongol-Okhotsk Ocean and the Pacific Ocean throughout the Phanerozoic. Due to its intracontinental setting, the GKMs serve as an ideal zone for evaluating the far-field effects of intraplate tectonism (Figure 1(a)) [15, 16]. The action of the two dynamical systems in the deeper parts of the GKMs region causes tectonic uplift of the shallow mountain ranges, and the low-temperature thermochronometer is a good record of this mountain-building process [17-19]. Researchers have conducted a large number of low-temperature Thermochronology studies on GKMs (Figure 1(b)). However, previous models based on apatite thermochronology studies have led to diverse views on the cooling history and tectonics of the region. While one view suggests that the exhumation of the Greater Khingan range began only at the start of the Cenozoic epoch [20], another maintains that it started already in the late Mesozoic and continued during the Cenozoic [18]. Moreover, Wu et al. (2016) suggested that major cooling events occurred during the Cretaceous (120 to 90 Ma) and at ca. 60 Ma [21]. However, Pang et al. (2020) held that the GKMs underwent two episodes of rapid exhumation, which occurred at ca. 100–60 Ma and ca. 50–0 Ma, respectively, with a regional reheating event at ca. 60–50 Ma [22]. Scholars have also debated the driving forces behind the exhumation of the GKMs. While some attribute it to mantle upwelling [23], others maintain that the phenomenon resulted from back-arc extensional processes during the combined subduction of the Mongol-Okhotsk Ocean and Pacific Oceans [22, 24, 25] or solely as a byproduct of Pacific Ocean subduction [18, 21, 22]. The development of the GKMs is closely associated with the formation of the adjacent basins [26, 27]. Consequently, in-depth research on the GKMs and their neighboring basins has been conducted [18, 20-22]. Nevertheless, the basin-mountain relationship remains incompletely understood.

Therefore, in order to clarify the transition time of the two major dynamical systems, the Mongol-Okhotsk Ocean and the Paleo-Pacific Ocean, in Northeast China and the corresponding records of thermal processes in the shallow part of the Earth’s crust, we conducted thermal dating studies (including apatite fission-track [AFT] dating and zircon U-Pb dating) of sedimentary rocks in the central GKMs and deduced the thermal evolution history of the region based on the measured thermochronology data through two-dimensional (HeFTy) and three-dimensional (Pecube) thermal history simulations [28-32]. In addition, published low-temperature thermochronological data from the GKMs and peripheral basins were analyzed to synthesize the geological information and then infer the correlation between the mountain range and the basin system and to further explore the structural consistency between the deep and shallow crust, which is controlled by geodynamic processes.

Microcontinental blocks that formed part of the Paleozoic-early Mesozoic Paleo-Asian orogenic belt are exposed in northeast China [4, 33-35]. These blocks underwent amalgamation with the northern periphery of the North China craton during the Late Permian-Middle Triassic as part of the initial phase of the closure of the western Paleo-Asian Ocean [36-40]. From the early Mesozoic onward, subduction processes were documented in the area, which was intimately linked to the cessation of the Mongol-Okhotsk Ocean and Paleo-Pacific plate. The subduction of the Paleo-Pacific plate has remarkably influenced the entire Northeast China region since the early Mesozoic, most notably during the Early Cretaceous [3, 41, 42]. The tectonic and geological development of NE China in the Cenozoic was largely regulated by the circum-Pacific tectonic realm [43-45].

Situated at the heart of the Northeast Asian Xingmeng Orogenic Belt, the Greater Khingan constitutes a pivotal tectonic entity, extending southward at the periphery of the Siberian Plate. Dominated by a principal tectonic line running predominantly northeast and north-northeast, the region presents a distinctive basin-mountain topography, featuring the Hailar-Erlian Basin and the Songliao Basin flanking either side (Figure 1) [46]. According to the division of tectonic units, the GKMs are mainly composed of the Erguna Massif, Xing'an Massif, and Songnen Massif from north to south. Large-scale magmatic events erupted in the GKMs region during the Late Paleozoic-Mesozoic, resulting in the distinctive feature of a large distribution of regional igneous rocks [40, 47]. In addition, two large deep strike-slip faults, Derbugan and Nenjiang-Balihan, have been developed on the north-western and eastern sides of the GKMs, respectively (Figure 1(a)), and a large number of NE-NNE, EW, and NW directed faults have been developed between these two faults, constituting a spatial pattern of “lattice-like faults” [18, 48]. The investigation was conducted within the central zone of the GKMs (Figure 2), where the observable sedimentary layers encompass Ordovician detrital rocks, Mesozoic volcanic rocks, and conglomerates. The predominant orientation of these geological samples is generally in the NE direction (Figure 2(a) and 2(b)).

There are many Mesozoic volcanic rock formations exposed, among which the Middle Jurassic Tamulangou Formation (J₂tm), Upper Jurassic Manketouebo Formation (J₃mk), Manitu Formation (J₃mn), and Lower Cretaceous Baiyingaolao Formation (K₁b) are volcanic sedimentary formations (Figure 2(a) and 2(b)). The Murui Formation (J₃mr) is a set of conglomerates that represent the bottom of the Upper Jurassic and are in angular unconformable contact with the overlying Manketouebo Formation (J₃mk), Baiyingaolao Formation (K₁b), and the underlying Luohe Formation (O₃lh), respectively (Figure 2(a)–2(b)). Conglomerates in the different regions cover different lithostratigraphic units. Most of the conglomerates are unconformable within the Tamulangou Formation and Luohe Formation, with the overlying strata being the Manketouebo Formation (Figure 2(c)). It is worth noting that the Late Mesozoic stratigraphy in Northeast Asia is widely time-penetrated, and the age of formation of the above strata varies by about 20–10 Ma in different regions [47, 49, 50]. Although this highly controversial issue has been identified by the geological community, it is currently inconclusive, and therefore, when discussing the age of formation of a particular stratum, we usually refer to the age of formation at the location of the build-up section.

In the Jurassic–Early Cretaceous, the middle section of the GKMs exhibited a well-established network of NE trending faults which were instrumental in shaping the primary structural framework of the region (Figure 2(d)). Among them, the Chaihe thrust fault, a large-scale and deeply penetrating fault, is notable for its exposed length of approximately 15 km. The fault plane’s NW-inclined dip angle, ranging between 45° and 50°, resulted in an overthrust of the Upper Ordovician Luohe Formation onto the Middle Jurassic Tamulangou Formation (Figure 3(a) and 3(b)). The thrust fault’s structural analysis provides supporting evidence for the existence of a thrust nappe, with an actual nappe distance exceeding 4.2 km [50].

The Murui Formation is characterized by the variegated dimensions of conglomerate material—initially fine, evolving toward coarseness, subsequently reversing the pattern, and transitioning back to fine grain size (Figure 3(d) and 3(e)). Provenance analysis of sandstone embedded in this formation indicates that the detrital sediments are predominantly derived from the Chaihe fold-thrust belt [46]. Given the observations made on the geological characteristics, the Upper Jurassic Murui Formation can be classified as a syntectonic deposit originating from the Chaihe fold-thrust belt [46, 50].

The Luohe Formation of the Upper Ordovician and the Tamulangou Formation of the Middle Jurassic exhibit a faulted contact (Figure 2(d)). Additionally, the Murui Formation displays an unconformable relationship with the underlying Baiyingaolao Formation (Figure 2(e)). These stratigraphic contact relationships, in conjunction with the sedimentation times, provide valuable timing constraints that can be utilized in inversion modeling.

Trilobite fossil assemblages are present in the Luohe Formation’s northeastern segment and can be attributed to the Caradoc Series [51]. While the stratigraphic of Luohe Formation’s lithological association, paleontological attributes, and stratigraphic compared to its western extent (e.g., Alatan, Inner Mongolia, Heli, Sumu, and Hambudunzhao regions), indicate Luohe Formation depositional age in the Late Ordovician [52]. Zircon U-Pb dating conducted on detrital zircons extracted from the northern sector of the study area reveals the youngest zircon age peak age to be around 460 Ma, and the minimum age limit for the formation of the Luohe Formation is the Late Ordovician [53]. In light of the examination of plant fossil specimens and zircon U-Pb dating, it is reasonable to infer that the emergence of the Luohe Formation transpired in the Late Ordovician era, encompassing a time frame with an estimated duration of 460–440 Ma.

In the examined region, the Tamulangou Formation, representing the most recent volcanic sedimentary stratum underlying the Murui Formation, is estimated to have formed approximately 147 Ma [47]. Conversely, the Manketouebo Formation above the earliest volcanic sedimentary stratum is believed to have undergone diagenesis around 146 Ma [50]. This evidence suggests that the emergence of the Chaihe thrust fault and subsequent syntectonic deposition of the Murui Formation occurred between roughly 147 and 146 Ma.

The middle section of the GKMs was chosen as the primary location for sampling. A total of five detrital rock samples were retrieved for AFT thermochronology testing (Figure 2). Specifically, one sample (D05-6) was collected from the Upper Ordovician Luohe Formation, situated in the vicinity of the hanging wall associated with the fault plane of the Chaihe thrust fault (Figure 2(a)). In order to investigate the exhumation history of the fault structure, the HeFTy inversion modeling method was employed. The remaining samples (D08-6, A01, D00-1, and P1-11) were obtained from the Upper Jurassic Murui Formation. To evaluate the thermal evolution of the GKMs, two-dimensional thermal history inversion modeling (HeFTy) and three-dimensional numerical simulation tests of the three-dimensional geothermal field in the shallow crust (Pecube) were carried out on the basis of the collected data.

In order to elucidate the inverse modeling process, U-Pb zircon dating was executed on volcanic rocks (Figure 3(c)) present both above and below the unconformity demarcating the Baiyingaolao and Murui Formations. The stratigraphically significant Sample P1-13A originating from the Baiyingaolao Formation along with Sample A01 derived from the Murui Formation were subjected to dating. This enabled the establishment of vital boundary limits for subsequent inversion modeling concerning thermal history. The sampling details have been meticulously documented in Table 1.

The process of mineral separation was performed at the Fission Track Testing and Analysis Laboratory, Northwest Institute of Eco-Environment and Resource, CAS. The apatite grains were meticulously extracted by utilizing conventional magnetic separation techniques. Subsequently, the extracted apatite grains were engraved in epoxy resin, polished to a high degree, and then immersed in a 5 N HNO₃ solution at 20°C for a period of time. This action was followed by a 20-second etching process to amplify the naturally induced fission tracks. The age of Durango (31.4 ± 0.5 Ma) and Fish Canyon Tuff (27.9 ± 0.7 Ma) apatite grains served as age-standardized. The samples, age standards, and IRMM540R dosimeter glass underwent radiation within the Oregon State University thermal reactor, with a nominal neutron flux of 9 × 10¹⁵ n•cm⁻². Upon radiation, the mica detector was subjected to etching in a 40% HF (hydrofluoric acid) solution at a temperature of 20°C to magnify the fission tracks caused by neutrons in the apatite crystals. Fission track analysis was executed in the Fission Track Testing and Analysis Laboratory, Northwest Institute of Eco-Environment and Resource, CAS. The ζ calibration method [54] was used to determine the age of fission tracks; the individual ζ value in this study was 257.56 ± 17.67.

Zircon sorting was carried out at the Hebei Regional Geological and Mineral Investigation Institute. Uniformly crushed whole-rock samples of 80–100 mesh size were subjected to an electromagnetic method for the separation of heavy minerals. Subsequently, zircon crystals were selected with the aid of a binocular lens following elutriation. The selection process involved the identification of zircons with a well-formed crystal shape and free from inclusions or cracks.

The zircon crystals were embedded on a resin surface and prepared for imaging via transmitted light, reflected light, and cathodoluminescence. The imaging work was performed at Beijing Zirnian Pilot Technology Co., Ltd., and the resulting images were used for U-Pb dating using laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS). The U-Pb dating work was carried out at the Key Laboratory of Mineral Resources Evaluation in Northeast Asia, Ministry of Natural Resources of China, utilizing an Agilent 7900 ICP-MS instrument and a Coherent GeoLasPro 193 nm ArF excimer laser. The synthetic silicate glass reference material NITS610 was employed for instrument state adjustment, while the international standard zircon 91,500 (1062Ma) [55] was utilized for the isotopic composition external standard. The Plesovice zircon, with a monitored age of 337 Ma [56], was used as the monitoring zircon standard. The resulting data were processed using ICPMSDataCal [57] and corrected via the method detailed in Reference 58 for common lead isotope ratios. Subsequently, the zircon U-Pb harmonic diagram showing a relatively young age (<1000 Ma) was calculated and drawn using Isoplot [59]. An age error of 1σ was taken into consideration, and the ²⁰⁶Pb/²³⁸U age value was adopted. The obtained zircon U-Pb ages for the samples are presented in online supplementary Table S1.

AFT dating results are shown in Table 2. Samples D08-6 and P1-11 failed to sort out enough apatite particles to be used for testing, resulting in dispersed single-particle ages, and their age results are indicative only, while the remaining samples obtained sufficient numbers of apatite particles for reliable age results (Table 2). We used the χ² test P (χ²) to analyze the age dispersion of apatite grains: (1) P (χ²) > 5% indicates that the track age distribution is more concentrated [8] and (2) P (χ²) < 5% indicates that the track age distribution is more dispersed [60]. Three of the five samples failed the χ² test (the χ² test with P (χ²) < 5%). The samples AFT ages were scattered, with the center age distributed between 61.8 ± 6.4 Ma and 260 ± 36 Ma and the average length of the track length between 11.62 ± 1.09 µm and 12.85 ± 2.01 µm. Based on the correlation between the peak ages of the sample models, the peak ages were divided into three groups, namely, P1 (between 41.9 and 16.2 Ma), P2 (between 159 and 108 Ma), and P3 (378 Ma; Table 2). We used Density Plotter 7.0 [61] to draw the length-frequency distribution of fission tracks (Figure 4), and all track lengths have been corrected by the C-axis.

Before interpreting the AFT age of the clasts, it is necessary to evaluate the possibility of resetting or partial resetting caused by the burial depth or thermal disturbance of the isotherm after sample deposition [62]. The A01, D00-1, and P1-11 grain ages for the Murui Formation samples (Upper Jurassic) are scattered, while the D08-6 grain ages are relatively concentrated (Figure 4). Among them, 68 apatite grains in D00-1 are used for AFT dating (Figure 4(c)), and only five apatite grains have AFT ages older than the sedimentary age of the stratum, indicating that the samples are almost completely thermally reset and the thermal history evolution information of the samples from the sedimentary stratum is retained. Although the number of grains used for dating sample A01 is small (Figure 4(a)), it also shows a similar age trend with D00-1, and the sample is almost completely thermally reset. However, unlike the above two discussed samples (D00-1 and A01), the AFT age of most particles in sample P1-11 is older than the sedimentary age of the stratum (Figure 4(g)), indicating that the sample has not undergone thermal reset and mainly records the ancient thermal history information of the provenance area. The grains from sample D08-6 are few (Figure 4(i)), and the grain age almost does not overlap with the stratum age, indicating that it may be composed of older grains, so it has not been reset by heating processes. The D05-6 age distribution for the Upper Ordovician samples is concentrated (Figure 4(e)), and 85 apatite grain ages are all smaller than the lower limit age of the formation, indicating that the samples have experienced thermal reset and retained the thermal history of the sedimentary formation. It is worth noting that the track length of the sample presents an obvious bimodal distribution (Figure 4(f)), which also directly indicates that the sample has experienced a complex thermal history in the stratum, resulting in partial annealing of the AFT length [63].

We have summarized the distribution frequency of apatite grains age for five detrital samples (Figure 4(j)), which show three peaks: G1 (36.8 ± 3.3 Ma), G2 (112.3 ± 5.3 Ma), and G3 (277 ± 29 Ma). About 70% of the apatite grains have ages concentrated at 112.3 ± 5.3 Ma, which is consistent with our peak age separation results for five samples (Table 2).

The annealing of the track will reduce the track density, resulting in lower dating results [64]. In addition to temperature, chemical composition, and mineral anisotropy, the main factors affecting apatite track annealing are the D_par value, the thermal history information of the source region retained in the sample is more likely to be erased, and conversely, the larger the D_par value, the easier it is for the sample to preserve the cooling information of the source region [63-66]. Therefore, discussing the D_par value can reflect the annealing resistance of the sample and thus better determine the true thermal history experienced by the rock sample. Generally speaking, the smaller the D_par, the faster the annealing rate [19, 65]. The D_par values of the five samples show a good correlation with apatite grain age and track length (Figure 5).

The D_par value for sample D00-1, from the Upper Jurassic Murui Formation, is low. This indicates that the AFT was easier to anneal, resulting in the young sample age and short track length (Figure 5). The A01 sample has a relatively large D_par value and was more resistant to annealing. The dating results also show an aging trend, and the track length is also relatively dispersed (Figure 5). The distribution of D_par values for samples D05-6, P1-11, and D08-6 are in the middle state, where the grain age distribution of P1-11 is generally older, and the distribution range of D05-6 sample track length at each interval is very average.

A total of 103 zircon grains were collected from volcanic rock sample P1-13A and detrital rock sample A01. Most of the magmatic zircon crystals are subangular, slender cylindrical shapes, with sizes between 50 and 120 µm and obvious oscillatory ring structures (Figure 6(b)). Detrital zircon crystals are characterized by a high degree of idiomorphic and large particle size, indicating short transportation distance and rapid deposition (Figure 6(e)).

From sample P1-13A, from the Cretaceous Baiyingaolao Formation, 20 zircon crystals were selected for U-Pb dating. The age distribution is relatively concentrated, ranging from 140 to 100 Ma, with a concordant age of 115 ± 1.24 Ma (MSWD (mean standard weighted deviation) = 0.25; Figure 6(a)).

Sample A01, from the Late Jurassic Murui Formation, contains detrital zircon crystals. In this sample, 83 zircon grains were selected for U-Pb dating. This gave an age peak at 140 Ma (Figure 6(c)), while the weighted average age of the youngest three zircon grains is 121.72 ± 2 Ma (MSWD = 0.056; Figure 6(d)).

The zircon dating results for samples P1-13A and A01 constrain the unconformity occurrence time between the Baiyingaolao Formation and Murui Formation to 122–115 Ma, providing an important constraint for subsequent thermal history modeling.

By inputting geological constraints and random search boxes, HeFTy inverse modeling can be used to predict the thermal history of rocks and the AFT age for each inversion path and verify the inversion results [31, 67].

Thermal history inversion modeling of detrital samples is not a simple matter. It is necessary to ensure that the samples used for testing have enough apatite particles and track length to provide robust results. Second, there should be sufficient geological information to provide precise constraints for any inversion model. To allow for thermal history inversion modeling, we selected two suitable samples (Figure 2)—sample D05-6 (from the Upper Ordovician Luohe Formation) and sample D00-1 (from the Upper Jurassic Murui Formation). For each sample, several path tests with different constraints are set for comparison (Figures 7 and 8). The Monte Carlo search method was used for the model and random substitution spacing was utilized. Kuipers Statistics was used for the GOF (good fitting) method, and 2,000,000 test paths were set for the model. See Table 3 for the constraints and explanations in the modeling process. The modeling information for samples D05-6 and D00-1 is given in online supplementary Table S2.

Sample D05-6 was taken from the area near the Chaihe, a thrust fault plane. Two models were utilized to check whether the sample records the rock lifting process from thrusting or was subject to thermal disturbance from the fault. Only one stratum sedimentation time limit constraint and an extensive random search box were inputted into the PATHA1 model—allowing HeFTy to freely explore possible thermal history paths (Figure 7(a)). Subsequently, the PATHA2 model was set up, which includes two geological constraint frames, namely, a stratigraphic sedimentation time limit constraint and a fault unconformity constraint (Figure 7(b)). Of course, it is not ruled out that the thermal disturbance caused by faults may not be lower than 120℃, but this is not recorded by AFT.

The D00-1 sample is taken from the Murui Formation in the east section (Figure 2(b)). Time constraints for the unconformity between this formation and the Baiyingaolao Formation are defined by zircon U-Pb dating (Figure 6). Two paths are set for D00-1 thermal history inversion. The constraints from stratigraphic sedimentation time limits and unconformity information were set in PATHB1 (Figure 8(a)). In fact, five AFT ages in the D00-1 sample are earlier than the sedimentary age of the stratum. To exclude the influence of this factor on the real thermal history of the sample, another group of models PATHB2 (Figure 8(b)) is set up, which adds an extensive search box to explore the thermal history before deposition.

The inversion simulation results for sample D05-6 indicate that both models show good consistency with regard to the thermal history of the rocks after the folding and thrusting event occurred, with a rapid cooling period recorded too. The two models show that the thermal history prior to the Early Cretaceous was divergent (Figure 7), and such a large time span may require more address constraints to limit the thermal history experienced by the rock. In the absence of clear fold-thrust geological constraints (PathA1), 413 good-fit paths were generated by inversion modeling, and the GOF reached 0.98. Inversion results for the best-fit path show that the sample has rapidly passed through part of the annealing zone (~83.5℃) since ~44 Ma (Figure 7(a)). The PATHA2 model, with two geological constraints, shows that the sample may have been affected by thermal disturbance from the fault, which makes the thermal history inversion insensitive to the temperature change during the thrusting process (Figure 7(b)). The best-fit path only experienced a slight cooling period at 147 Ma, and 151 good-fit paths were generated from the inversion models. Starting from ~35 Ma, it quickly passed through the partially annealed band at ~82℃.

Inverse modeling results for sample D00-1 have different thermal history paths compared with results obtained from sample D05-6. Both models (PATHB1 and PATHB2) record two rapid cooling events after the unconformity period. However, the inversion results do not show a clear cooling path during the thermal history (i.e., from deposition to the period of the unconformity). The PATHB1 model, with both stratigraphic constraints and unconformity constraints, has generated twenty-nine good-fit paths, and the best GOF inversions were 0.92 and 0.45, respectively. The best-fit path shows that the first rapid cooling event occurred at ~116–78 Ma, from 115℃ to 83℃, and the second rapid cooling event occurred at 4.8 Ma, from 83℃ to part of the annealing zone (Figure 8(a)). The addition of the widely searched Box in PathB2 has had some effects on the model, making the thermal history of the sample after the unconformity period record higher temperatures and a shorter cooling time (Figure 8(b)). PathB2 shows that the sample was cooled from ~139℃ to ~82℃ at ~121–87 Ma, with a second cooling event occurring since 1.25 Ma (Figure 8(b)).

Pecube is a software package designed to interpret thermochronological data by solving the heat transport equation in a 3D scenario, in a crustal or lithospheric block. Given the fixed and dynamic parameters in the surface process model, the tectonic and geomorphic evolution can be evaluated from a thermal history perspective. The time-temperature history, geothermal field, and exhumation rate of rock mineral particles from the deep crust to the surface can be simulated by using forward simulations and NA (Neighborhood Algorithm) inversion simulations, and the ages of different thermal chronological systems can be predicted [68-70]. NA inversion is based on the Bayesian NA [71, 72], which finds the optimal value of given model parameters, within the specified range, to minimize the Misfit function. The Misfit function is the norm of the difference between the observed and predicted values. See details in Braun [62, 68] and Braun et al. [69] [63, 68, 69].

Sample D00-1, in the eastern section of our study area, has a thermal history simulation that indicates a cooling event during the Cretaceous (Figure 8). To verify whether the whole region was cooled and exhumed during this period, we selected samples D05-6 and D08-6, which passed the χ² test, from the westernmost side of our study area, for NA inversion and forward modeling using the Pecube code. We introduced in section 2.2.2 that the Murui Formation is a syntectonic deposit of the Chaihe thrust fault, so the AFT of sample D08-6 from the Murui Formation may retain some cooling information in the marginal zone (Chaihe Thrust Fault). Two scenarios are utilized to explore the inflection time when the stripping rates changed. A 0.2° × 0.2° range is used as the simulation operation area (Figure 9). For the NA inversion, the current terrain data were used for the upper boundary of the model. The data were from the GTOPO30 digital elevation model (resolution: 30 m), which is sufficient to distinguish the landform difference for each sample. The crustal thickness selected to simulate the lower boundary of the model was 40 km [73]. The model was set so that the topographic relief was stable through time. From the previous thermal history inversion results and geological information, such as the formation time of the Chaihe fold-thrust belt, the time boundary for the model was stipulated between 150 and 70 Ma, and the basement temperature was set at 800℃ [74, 75]. See Table 4 for other thermal kinematic parameters and Table 5 for dynamic parameters. The simulation time for each model was set to 11,200 times. Finally, the best inversion results were carried forward, and the predicted values for the sample age and elevation were obtained.

The inversion results from Scenario 1 and Scenario 2 have very good constraints on the start time (T1) and end time (T2) of the exhumation event and the exhumation rate. The two-dimensional scatter distribution shows good convergence, and the probability distribution curve of the kernel function also shows obvious peaks. The lowest Misfit values for the two models are 0.83 and 0.75, respectively. The range of low Misfit values has a strong correlation with the peak values for each parameter and good convergence.

The distribution state of the kernel function probability curve, for the starting time parameter (T1) for Scenario 1, shows a single peak, with the maximum probability occurring at 147.6 Ma (Figure 10(a)). The denudation rate increased from 0.027 km/myr to 0.142 km/myr after 147.6 Ma (Figure 10(a) and 10(b)). The peak value of the kernel function for the exfoliation rate (E1) is between 0.02 and 0.04 km/myr, and the kernel function of E2 presents an obvious single peak feature.

The distribution state of the kernel function probability curve for the termination time parameter (T2), for Scenario 2, shows a single peak with the maximum probability occurring at 70.1 Ma (Figure 10(c)). After 70.1 Ma, the stripping rate decreased from 0.094 to 0.021 km/y (Figure 10(c) and 10(d)). The kernel function probability curves for the stripping rates E3 and E4 converge well, with a single peak value.

The three-dimensional forward crustal thermal model shows the temperature field and denudation direction within the crust (online supplementary Figure S1). In general, the forward modeling results are consistent with the thermal evolutionary history of the region. The fitting age, of the best-fitting samples from the two scenarios, should be less than the measured age obtained for sample D08-6 (Figure 11). This is because the detrital apatite particles in D08-6 reflect the cooling information of the source area, although there is clear evidence indicating that the source is near the thrust fault [46]. However, Pecube can only track the thermal evolutionary history of bedrock, which causes this phenomenon where the fitting age is smaller than the observation age. The scene-fitting age for sample D05-6 is slightly larger than the observation age. We propose that AFT annealing of the sample has caused the measured age to be smaller, and the bimodal distribution characteristics of the D05-6 sample track length support this conclusion (Figure 4(f)) [64]. Second, the samples were taken from the Upper Ordovician Luohe Formation. The distribution characteristics of apatite grain ages and track lengths show that the samples have recorded the thermal history since deposition in the stratum (Figure 5). Because the time boundary for the model was between 150 and 70 Ma (Table 5), this reduced sample annealing time, causing the fitting value to be larger than the measured value.

The thermal history process of a given sample can be explored by leveraging various annealing models with HeFTy inversion, utilizing information such as grain age, track length, and Dpar [31, 67]. Pecube facilitates the interpretation of low-temperature thermal chronology data by transforming geologic and geomorphic scenarios into a thermal evolutionary history. Additionally, Pecube can harness data derived from HeFty or QTQt software for an increased understanding of thermal processes [68, 69].

The results of the HeFTy and Pecube 3D inversion models suggest that rapid cooling and exhumation events occurred during two separate time periods with associated exhumation rates of 0.08 and 0.11 km/y for ~121–87 Ma and ~35 Ma, respectively. Moreover, the Pecube 3D inversion model suggests that rapid exhumation occurred between 147.6 and 70.1 Ma, and if crustal heat conduction was uniform, the exhumation rates would range from 0.094 to 0.142 km/y when assuming an average geothermal gradient of 20℃/km. The assumed ground temperature gradient of 20°C/km is considered reliable based on a study of the current heat flow structure in mainland China [75].

In essence, the results of the thermal history simulations, coupled with the broader geological framework [46, 50], signify that the GKMs were subjected to episodes of cooling and exhumation during the period ranging from 147 to 70 Ma, as well as from 35 Ma to the present.

The comparison of the age versus length distributions of AFTs in the GKMs displays a “Boomerang” pattern (Figure 12) [17]. As the AFT age gets younger, the AFT length tends to decrease and reaches a minimum at about 80 Ma, after which the track length becomes longer again, which most likely represents the occurrence of a new stripping event. In addition, the low-temperature thermochronological data (apatite and zircon fission tracks) from different areas of the GKMs also show different age and trace length distributions (online supplementary Table S3), which suggests that the GKMs region has experienced a complex thermal evolutionary history. In order to reconstruct the thermal evolutionary trajectory of the GKMs from the late Jurassic onward, both zircon U-Pb dating and detrital AFT dating were employed. Subsequently, two separate intervals of rapid cooling were identified in the GKMs since the Late Mesozoic via two independent thermal model simulations.

A prominent peak in detrital apatite FT ages ranging from 159 to 108 Ma is observed in Figure 13(a), with a considerable number of FT ages being proximal to ~112 Ma (G2). Results from HeFTy thermal history inversions, in conjunction with 3D numerical simulations, indicate the occurrence of a cooling event between 147 and 70 Ma. Pang (2020) proposed [22], based on model simulations of bedrock thermal history, that the rapid uplift of the GKMs took place between 100 and 60 Ma. Furthermore, AFT thermal history simulations on samples obtained from the Mohe basin (Figure 13(a)) in the northern segment of the GKMs unveil an exhumation event that took place between 120 and 90 Ma [21]. These findings provide solid evidence that northeast Asia has undergone a complex tectonic system transformation since the Mesozoic era, which triggered geological feedback mechanisms such as volcanic activity and mountain uplift in the GKMs.

The final suturing of the Mongol-Okhotsk Ocean is believed to have occurred sometime in the Jurassic or Cretaceous, with the ocean closing in a scissor-like or segmental fashion [8-13, 76]. A series of northeast-trending faults (Figure 2(a)), represented by the Chaihe fold-thrust belt, developed in the GKMs region during the Late Yanshan period (between 147 and 146 Ma), which may be the far-field counterpart of the Mongol-Okhotsk Ocean closure. The crust was deformed due to tectonic stresses in an NW-SE orientation [46, 50]. At the same time, the northern edge of the Amur basin, near the Mongol-Okhotsk suture zone, had a compressional uplift during the Early Cretaceous [77]. At this time, the adjacent mountain uplift was significant. The gravels in the syn-tectonic sedimentary layer (Murui Formation) developed during the same period as the Chaihe fold-thrust belt are up to boulder size (Figure 3(e)), suggesting that tectonically active activities at that time led to the rapid accumulation of detrital materials. NA inversion Scenario 1 also shows that the denudation rate increased to 0.142 km/myr after 147.6Ma (Figure 10(b)), indicating that the mountain was rapidly uplifted and denuded.

After the complete suturing of the Mongol-Okhotsk Ocean, the compressional orogeny came to a close, thereby opening a new phase of postorogenic extension for the GKMs. The extension-driven formation of metamorphic core complexes near Ganzhuermiao and Xinkailing within the GKMs occurred between 130 and 120 Ma [78-81]. Additionally, the identification of highly differentiated I-type and A-type granites and bimodal volcanic rocks from the Early Cretaceous period points toward the presence of oceanic crust material recycling in this region [47, 82-87]. The central Songliao Basin also shows evidence of an extensional detachment structure on the eastern flank of the GKMs, with an Ar-Ar age of 126.7 ± 1.54 Ma [88].

In the Early Cretaceous period, noteworthy NNW-SSE extension occurred in the southern section of the Mongol-Okhotsk suture zone, as documented in sources [26, 27, 89]. Additionally, the Songliao Basin was subject to ongoing development within an extensional environment [90]. It is observed that the swift exhumation of the GKMs in tandem with the rapid subsidence of the Songliao Basin and Hailar Basin leads to the conclusion that the formation of the basin and mountain exhumation are both outcomes of regional extension.

Analysis of detrital samples has shown the emergence of peak age, denoted as P2 at 42–16 Ma (Table 2). Thermal history modeling of HeFTy inversion further exhibits a cooling event that has continued since 35 Ma (Figure 13(a)). This outcome is consistent with the existing studies on low-temperature and thermal chronology that indicate a rapid exhumation event occurred in the GKMs since the Eocene (Figure 13(b)) [18, 20-22]. In Northeast Asia, the Zhangguangcai Mountains and the Lesser Khingan Mountains experienced Cenozoic cooling events [91, 92]. Examination of low-temperature geochronological ages reveals a decrease in age from east to west in the Songliao basin. The rate of uplift is positively correlated with the convergence rate of the Pacific plate [45, 93, 94]. The AFT thermal history of the Weijing intrusion, located near the Hailar and southern Songliao basins, also indicates exhumation in the Cenozoic during the GKMs cooling phase [45, 94, 95]. Based on these findings, we hypothesize that the tectonic activity was responsible for large-scale mountain denudation and basin sedimentation in Northeast Asia.

The convergence between the subducting Pacific plate and the Eurasian plate transpired prior to the commencement of the Cenozoic. A surge in the rate of convergence was observed during the early Cenozoic (Figure 13(c)). Noteworthy is the adaptation of the Pacific plate subduction angle to an NWW direction within the ~50–42 Ma epoch under the Eurasian plate (Figure 13(d)). This period was characterized by the occurrence of crustal cracking and faulting on the eastern periphery of Northeast Asia. The Sea of Japan ensued as a result of the processes that took place during this duration (~22–15 Ma). The region simultaneously underwent vast back-arc extension, a phenomenon that persisted until the Quaternary [96, 97].

In order to evaluate the spatiotemporal distribution of FT ages, a compendium of FT data from the GKMs and surrounding basins has been assembled (Figure 14). To ascertain the fidelity of the compiled FT ages in depicting the cooling history of the rock, exclusively those FT data that stem from bedrock or from samples that were thoroughly annealed in sedimentary strata were included in the comparative analysis. The pertinent data are itemized in online supplementary Table S4.

The age and spatial distribution of low-temperature thermochronology (FT) results exhibit distinct features as follows: First, the elder age group (>105 Ma) is predominantly distributed in the Mohe Basin, Hailar-Erlian Basin, and on the western slopes of the GKMs proximal to the Mongol-Okhotsk suture zone, while no such FT age group has been identified in the Songliao Basin (Figure 14(a)). Second, the younger age group (105, 45 Ma) is widely distributed across the GKMs and adjacent basins, with a distinct peak value (Figure 14(a) and 14(c)). Finally, the youngest FT age group (<45 Ma) is mostly present in the Mohe Basin and Songliao Basin. We propose that the differential impact of regional tectonic activity is the primary cause for the observed variations in the distribution profiles of FT results.

The Mohe Basin was in a tectonic context of contraction and deformation during the Late Jurassic-Early Cretaceous, and the low-temperature thermochronology of the sedimentary and volcanic rocks in the basin records a rapid denudation event at 144–132 Ma [98]. In addition, the Upper Jurassic Kaikukang Formation in the Desert Basin was subjected to fold deformation and uplift under compressive stress, generating a series of NEE-directed reverse faults [99]. At the same time, the extrusion deformation event was recorded in both the Hailar and Erlian basins [100-102]. Seismic reflection profiles in the Hailar Basin show that the Upper Jurassic Tamulangou Formation in the basin was cut by a series of NEE-trending reverse faults and truncated by the Lower Cretaceous strata, suggesting that uplift of the basin occurred between the Late Jurassic and Early Cretaceous [48]. Extrusive deformation in the Erlian Basin occurred after the stratigraphy of the Upper Jurassic Xing’an Formation, which caused fold deformation and uplift of the Xing'an Formation and erosion of the [48]. The Late Jurassic-Early Cretaceous extrusion and deformation events in the above basins were contemporaneous with the development of the NE-trending Chaihe fold-thrust belt in the GKMs, and combined with the distribution of low-temperature thermochronological age data for the GKMs and the surrounding basins (Figure 14), we believe that at this time, the GKMs was in the same extrusion tectonic environment as the Mohe Basin and the Hailar-Erlian Basin.

Since the Late Jurassic, the cooling ages of apatite and zircon fission-track thermochronometers from the GKMs have exhibited a progressive west-to-east younging trend, as reported by previous studies [18]. Additionally, studies on the peak U-Pb ages of magmatic zircons from the Hailar Basin, GKMs, and Songliao Basin also reveal the same trend of west-to-east younging [103-105]. Notably, the Hailar-Erlian basin represents an extensional faulted basin, which formed in the Early Cretaceous, and its geodynamic evolution was orchestrated by the same tectonic and magmatic processes that fashioned the GKMs, thereby demonstrating a similar structural history [26, 27, 89]. During this period, although the Songliao Basin had intermittent magmatic activity, the GKMs remained exempt from such processes [25, 106, 107]. Based on the results of a magnetotelluric survey conducted on the structural belt between the GKMs, Hailar Basin, and the Songliao Basin, it was observed that the eastern edge of the Hailar Basin corresponds to the Xing'an Terrain, while the contact zone between the western edge of the Songliao Basin and the easternmost end of the GKMs occurs as a thrust wedge system [108]. It is thus inferred that segments of the GKMs and Songliao Basin were affiliated with distinct tectonic domains during this period.

During the Early Cretaceous epoch, the Hailar-Erlian Basin and the Mohe Basin, located south of the Mongol-Okhotsk suture zone, underwent continuous NNW-SSE extension leading to the development of normal faults, rifts, and depressions [26, 27, 89]. At the same time, a series of positive faults gradually formed a spatial tectonic framework of “lattice rupture” under the action of tensile stress in the GKMs area, and the low-temperature thermochronological data in the southern part of the GKMs analyzed that stretching and stripping occurred in the middle of the Early Cretaceous [109]. The extension likely resulted from the preceding closure of the Mongol-Okhotsk Ocean, an event possibly attributed to slab break-off and/or delamination processes [24]. Fission-track age distribution analysis reveals that from the late Jurassic to the middle of the early Cretaceous period, the same tectonic domain controlled both the GKMs and the surrounding basins to the west and north. The rapid denudation and erosion of the GKMs provided material input into the basins, ultimately resulting in the formation of a basin-mountain terrain [22].

During the Late Jurassic and Early Cretaceous, the west region of the Songliao Basin was mainly distributed with magmatic activity, porphyry Cu and Mo ore deposits, and Pb-Zn polymetallic ore deposits [26, 106, 110]. However, a shift occurred in the middle and late Early Cretaceous, leading to widespread deposits of these types across Northeast China [106, 111]. During this time, the Songliao Basin continued to evolve in an extensional environment, with rifting and depression dominating the area. Subsequently, two tectonic reversals and uplift and exhumation events occurred in the basin since the Eocene, driven by the change in subduction direction and angle of the Pacific plate [44, 90, 112]. Taking into account geological evidence and the time-space characteristics of low-temperature thermal chronology, we propose that the GKM basin-mountain system and the Songliao Basin were both affected by the same tectonic stress field post-middle and late Early Cretaceous. This ultimately led to rapid denudation of the mountains and the concurrent formation and evolution of basins that were coherently connected both in time and space.

During the Late Jurassic to the Middle-Late Early Cretaceous period, spanning from 147 to 105 Ma, the Songliao Basin and GKMs exhibited distinct geological characteristics pertaining to volcanic activity, deposit distribution, and regional structure. Based on seismic tomography data, it is evident that subduction of the Paleo-Pacific plate did not play a significant role in driving expansion in the region during this period, possibly due to the fact that the Paleo-Pacific Ocean did not extend beneath the GKMs [108]. However, it is suggested that crustal thickening resulting from the closure of the Mongol-Okhotsk Ocean may have triggered gravitational collapse and mantle convection, leading to the upwelling of mantle material and contributing to the extension in the GKMs and basin systems (Figure 15(a)). The discovery of adakite in the GKMs also supports the hypothesis that crustal melting occurred during this period, potentially facilitated by the over-thickening of the crust due to oceanic closure [7, 41]. Although the rollback of the Paleo-Pacific plate at approximately 165–100 Ma may have presented an opportunity for regional expansion [113], further investigation is required to determine the extent of its influence.

Following the Cretaceous period (<105 Ma), Northeast Asia predominantly experienced stress from the subduction of the Paleo-Pacific Ocean beneath the Eurasian continent, constituting a significant tectonic system that would have governed the tectonic evolution of basin-mountain systems located around the GKMs. For the duration of the past 35 Ma, the GKMs have undergone extension (see Figure 15(b)). Between 50 and 42 Ma, the Pacific plate underwent an alteration in subduction orientation, lasting approximately 10 Ma, from NNW to NWW. Since the Eocene period, the overall convergent rate of the Pacific and Eurasian plates has increased [114]. This progression regulated back-arc spreading and dextral strike-slip faulting in Northeast Asia.

The following conclusions have been drawn from the dating and analysis of fission traces of detrital apatite in the GKMs region of northeastern China and from the analysis of fission trace data from the GKMs as a whole and from the basins on both sides of the GKMs:

1. Utilizing detrital AFT analysis to infer thermal history inversions, it has been demonstrated that the GKMs underwent two stages of tectonic exhumation during the late Jurassic period spanning from 147 to 70 Ma, as well as a subsequent event approximately 35 Ma ago.

2. In the realm of geologic evolution, the fission-track ages of the GKMs, Hailar-Erlian Basin, and Mohe Basin are grouped in a similar manner at greater than 105 Ma. However, unlike these aforementioned locations, the Songliao Basin displays a lack of correlation. This divergence can be attributed to variations in the tectonic domains that governed the region during this time frame. Nevertheless, between the ages of 105 and 45 Ma, the AFT ages reached their pinnacle and were widely dispersed throughout the basin-mountain system, including the Songliao Basin. These findings suggest a unified tectonic regimen was operational in the region during this period.

3. The delamination and concomitant crustal thickening resulting from the Mongol-Okhotsk Ocean closure are decisive factors in the basin-mountains system of the GKMs. Regional expansion during this period was supported by asthenospheric upwelling due to mantle convection, driven by the rollback of the Paleo-Pacific Plate. From the middle to the end of the Cretaceous period (<105 Ma), the tectonic evolution of Northeast Asia was defined by the subduction of the Paleo-Pacific Ocean under Eurasia.

We thank the personnel of the the Fission Track Testing and Analysis Laboratory, Northwest Institute of Eco Environment and Resource, CAS for their help with apatite fission-track (AFT) analysis. We also give thanks to Professor An Yin of the University of California, Los Angeles, for his suggestions on revising the article. Thank you also to the three anonymous reviewers for their key comments on the revision of the manuscript. This work was supported by the National Nature Science Foundation of China (grant number 41872234) and the Opening Foundation of Key Laboratory of Mineral Resources Evaluation in Northeast Asia, Ministry of Natural Resources (grant number DBY-ZZ-18-08).

The authors declare no conflicts of interest relevant to this study.

The thermal history simulation software used in this study is available via Ketcham et al. 2009; https://doi.org/10.1016/j.epsl.2009.05.015 for, Pecube numerical simulation program by Braun, J, 2003; https://doi.org/10.1016/S0098-3004(03),00052-9.