Deformation of the Crust and Upper Mantle beneath the North China Craton and Its Adjacent Areas Constrained by Rayleigh Wave Phase Velocity and Azimuthal Anisotropy

: The North China Craton (NCC) has experienced strong tectonic deformation and lithospheric thinning since the Cenozoic. To better constrain the geodynamic processes and mechanisms of the lithospheric deformation, we used a linear damped least squares method to invert simultaneously Rayleigh wave phase velocity and azimuthal anisotropy at periods of 10–80 s with teleseismic data recorded by 388 permanent stations in the NCC and its adjacent areas. The results reveal that the anomalies of Rayleigh wave phase velocity and azimuthal anisotropy are in good agreement with the tectonic domains in the study area. Low-phase velocities appear in the rift grabens and sedimentary basins at short periods. A rotation pattern of the fast axis direction of the Rayleigh wave together with a distinct low-velocity anomaly occurs around the Datong volcano. A NW–SE trending azimuthal anisotropy and a low-velocity anomaly at periods of 60–80 s are observed subparallel to the Zhangbo fault zone. The whole lithosphere domain of the Ordos block shows a high-phase velocity and counterclockwise rotated fast axis. The northeastern margin of the Tibetan plateau is dominated by a low-velocity and coherent NW–SE fast axis direction. We infer that the subduction of the Paleo-Paciﬁc plate and eastward material escape of the Tibetan plateau mainly contribute to the deformation of the crust and upper mantle in the NCC. velocities are estimated in the southwestern Zhangbo fault


Introduction
The North China Craton (NCC) is the oldest Archean craton in China where the earliest crustal rock can be dated back to 3.8 Ga [1]. In the late Paleoproterozoic (~1.85 Ga), the North China Craton was amalgamated by the collision between the eastern NCC (ENCC) and western NCC (WNCC) along the Trans-North China orogen (TNCO) until its final cratonization [2] (Figure 1a). The NCC, in the eastern Eurasian plate, is composed of several tectonic blocks including horsts and grabens shown in Figure 1b. In addition, the NCC is currently surrounded by orogenic belts formed in different geological episodes. To the north, the central Asian orogenic belt (CAOB) is a product of the closure of the Paleo-Asian ocean in the Paleozoic-Mesozoic [3]. To the southwest, the Qilian orogen (QLO) has been assembled by Paleozoic subduction-accretion [4]. To the south, the Qinling-Dabie orogenic belt (QDO) has been formed by the collision between the Yangtze plate and the NCC in the Triassic [2].  [5] and Zhu et al. [6]). The thick blue rectangle indicates the study area. The bold black lines mark the boundaries of the major tectonic units. The cyan, blue and red circles are the large historic earthquakes (M ≥ 6, [7,8]). The black arrows represent the movements of the different blocks or plates. The magenta arrows show the GPS velocity [9]. The Paleo-Pacific plate was subducted beneath the Eurasian plate in a NNW direction (Figure 1a), and then turned gradually to a westward direction from~50-42 Ma describing a counterclockwise rotation [6]. Previous studies have shown that this subduction of the Paleo-Pacific plate below the Eurasian plate transformed the dynamic deformation of the NCC from N-S shortening to almost E-W intraplate extension, resulting in a series of largescale extensional structures, detachment faults and rift grabens in the NCC [10][11][12][13]. The stagnant slab dehydration of the Paleo-Pacific plate in the mantle transition zone led to an unsteady asthenospheric flow, partial melting of the overlying lithosphere and magmatic activity [12,[14][15][16]. Especially in the eastern NCC, the lithospheric thickness has been thinned by at least 100 km, and the original cratonic structures and properties have been significantly modified due to the lithospheric partial melting mentioned above [12,13]. The uplift and outward growth of the Tibetan plateau have been caused by the collision between the Indian and the Eurasian plates since~50 Ma [17,18]. The Qinling-Dabie orogen between the Sichuan basin and the Ordos block was proposed as acting as an important channel for the eastward material escape of the Tibetan plateau [17,19]. This extrusion of the Tibetan plateau even affected the deformation of the crust and mantle in the Trans-North China orogen [20][21][22]. Therefore, the NCC suffers the westward subduction of the Paleo-Pacific plate in the east and eastward extrusion of the Tibetan plateau in the southwest, leading to significant deformation of the crust and upper mantle structure [10,11,23,24]. In general, cratons are generally considered to be underlain by a thick, cold and refractory mantle lithosphere [13][14][15]. Lithospheric deformation is likely to cause prominent changes in velocity and anisotropy of the crust and mantle. Thus, the velocity structure and anisotropy can effectively constrain the deformation of the crust and mantle structure in turn, which help us discuss the tectonic evolution of the NCC and its adjacent areas.
Recently, many studies revealed the velocity structure and seismic anisotropy beneath the crust and upper mantle of the NCC [23,[25][26][27]. However, these results on understanding the lithospheric thinning and deformation of the NCC are still controversial. For example, the low-velocity anomalies of the asthenosphere revealed by Xu et al. [16], related to the lithospheric thinning, are segmented by irregular high-velocity anomalies in the eastern NCC, which are obviously inconsistent with the results obtained by full-waveform tomography [23], indicating that the eastern NCC is dominated by extremely low velocities. The upwelling of hot material beneath the Datong volcano is indicated by low-velocity anomalies in previous results [20,21,25,27], but its origin and mechanism are still under debate. Huang et al. [27] proposed that the asthenospheric upwelling came from the northeastern margin of the Ordos block. Yao et al. [20] confirmed this hypothesis and attributed it to lithospheric stretching. However, the results of S-wave [21,25] suggested that the hot materials originated from the asthenosphere in the vertical direction rather than the northeastern margin of the Ordos block. Chen et al. [25] further speculated that the upwelling of the asthenosphere beneath the Datong volcano would flow into the eastern NCC. Seismic anisotropy is often used to investigate the deformation of the crust and upper mantle. The S-wave splitting studies provide useful information about anisotropy in the crust and upper mantle [24,26]. However, the depth of the anisotropic structure is usually difficult to obtain by S-wave splitting. Moreover, studies on lithospheric deformation are still discussed due to the different interpretations of the NCC. Based on the fast axis directions of the local S-wave and SKS/SKKS phases, Gao et al. [28] suggested that the deformation of the crust and upper mantle was not a simple coupled pattern. The fast axis direction of the SKS/SKKS phases correlates well with the absolute plate motion (APM) direction, implying the APM-driven model of the upper mantle anisotropy in the NCC [29]. This model disagrees with the hypotheses inferred from stress-aligned microcracks, frozen fossil anisotropy and asthenospheric flow [22,25,26]. Compared with S-wave splitting, the surface wave azimuthal anisotropy estimated from different periods has a better vertical resolution and can provide a depth independent anisotropy to constrain the deformation of the crust and upper mantle.
Using teleseismic data recorded by the permanent stations in the NCC and its adjacent areas, the goal of this study is to further improve constraints on the deformation of the crust and upper mantle. We estimate simultaneously Rayleigh wave phase velocity and azimuthal anisotropy at periods of 10-80 s. By comparing the previous results, we focus on investigating the tectonic deformation and geodynamic mechanisms forming the NCC and its adjacent areas.

Data Acquisition and Processing
The seismological observation system in China has achieved rapid development [30] and the permanent seismic stations cover densely the NCC and its adjacent areas (Figure 1b In this study, the classical two-station method is used to measure the Rayleigh wave phase velocity. This method is not affected by errors such as the epicenter, the earthquake origin time and the initial phase, but only requires that the event and two stations are roughly on the same great circle path [25]. In order to measure the high-quality and reliable dispersion curves, we require that: (1) the magnitude of the event is not less than 5.5 and the focal depth is not greater than 50 km, which may ensure the signal-to-noise ratio of the Rayleigh wave; (2) the epicenter distance is between 20 • and 120 • to reduce the influence of the near-source effect and high-order surface wave interference; (3) the interstation distance should be larger than 3 wavelengths to decrease the errors in the dispersion measurement; (4) calculating the cross-correlation function with the seismograms of the same event recorded by the selected station pair after down sampling and removing the means, trends and instrument response. According to the above processing, a total of 1347 earthquake events were selected in this study ( Figure 2). We found that the distribution of these earthquakes has good epicentral distance and azimuth coverage, but a considerable number of earthquakes are concentrated in areas such as the western Paleo-Pacific subduction zone, resulting in multiple great circle paths passing through a certain station pair. Therefore, we adopt the time-frequency phase weighted stacking method [31,32] to stack the waveforms of the cross-correlation function with the same great circle path for the station pair, which can improve the signal-to-noise ratio of the cross-correlation functions and measurement accuracy of Rayleigh wave phase velocity dispersion curve. Then, the automated frequencytime analysis method [33] is used to estimate the Rayleigh wave phase velocity dispersion curves in the period range of 10-80 s for all station pairs. The final dispersion curves are selected by the criteria of signal-to-noise ratios greater than 20 and phase velocities less than twice the standard deviation at each period.

Surface Wave Tomography Inversion
We use the tomographic method of Barmin et al. [34] to invert Rayleigh wave phase velocity and azimuthal anisotropy on a 0.5 • × 0.5 • spatial grid beneath the NCC and its adjacent areas. This method is based on the damped least squares inversion to minimize a penalty function composed of a linear combination of data misfit, model smoothness and the perturbation to a reference model. For the weakly anisotropic medium, the relationship between Rayleigh wave phase velocity c(T) with a period of T and the azimuthal angle θ can be expressed as [35]: Here, c 0 is the isotropic phase velocity of the reference model, and c 1 , c 2 and c 3 , c 4 are 2θ and 4θ anisotropic components, respectively. Generally, 4θ can often be ignored in weakly anisotropic media, and then the strength of anisotropy, A, and the azimuth angle of the fast wave, θ, can be obtained by the following two formulas: During the inversion, the two parameters specifying the weight of smoothing and the width of the smoothing area control the damping factor, which is determined by a number of tests using different combinations of parameters by considering data misfit, model resolution and model smoothness [36]. Meanwhile, we refer to the processing steps of Guo et al. [37] to perform the surface wave tomography inversion. We first invert Rayleigh wave phase velocity and azimuth anisotropy at each period. We then remove the dispersion curves with large misfits, update the data set and then perform the inversion again to obtain the final model of Rayleigh wave phase velocity and azimuth anisotropy in this study.

Robustness Analysis
After the quality control of dispersion curves in Section 2.1, we plot Rayleigh wave phase velocity dispersion curves as a function of period in Figure 3a. The fewest dispersion curve numbers at a 10 s period are still as high as 8760. At a period of 55 s, the ray number is higher than 45,000. A considerable number of dispersion curves provide quite dense and crossed ray paths and ensure the high resolution of the final inversion. Here, we also take two station pairs as an example to check the phase velocity dispersion curves ( Figure 3b). The paths of the station pairs WUL-LZH and JJ-BJT span different geological blocks in the study area ( Figure 1b). The station pair WUL-LZH is from the south of the Qilian orogen, and across the Ordos block and the Trans-North China orogen. The path between stations JJ and BJT lies mostly within the Bohai Bay basin. The Rayleigh wave phase velocities at short (<15 s) and long (>40 s) periods along the path of JJ-BJT are obviously lower than that along the path of WUL-LZH because of the sediments and thin lithosphere in the Bohai Bay basin of the eastern NCC. In contrast, in the period range of 15-40 s, the phase velocities are much higher in the Bohai Bay basin, indicating the colder and stronger lower crust and uppermost mantle between stations JJ and BJT. Histograms of data misfits to measurements from the final tomography at different periods are presented in Figure 4. The standard deviation displayed in the upper right corner of each panel increases gradually from 1.38 s to 2.43 s. This rising standard deviation is mainly due to the increase of random error and decrease of the signal-to-noise ratio in Rayleigh wave phase velocity measurements. However, the Gaussian distributions of data misfits are nearly concentrated between ±4 s, which implies that the final models obtained by this inversion fit the observed data well.
In principle, the resolution and reliability of tomographic results depend on the dense and crossed paths of Rayleigh wave phase velocity measurements. The study area is dominated by a good path coverage ( Figure 5), although the rate of data rejection at each period is more than 30% based on the strict selection criteria in Section 2.1. The resolution maps are simultaneously determined during tomography ( Figure 6). As described by Barmin et al. [34], the resolution is estimated from the resolution matrix and defined by a complicated function of path density. In Figure 6, we find that the better resolution gradually shrinks from the edge to the central part of the study area due to the path coverage. Overall, most of the study area is covered by a robust resolution of~1 • -2 • (~100-200 km).     (Figures 7 and S1). This low-velocity anomaly extends slightly northwest (Figure 7c,d), which is also observed by the corresponding low S-wave velocity [20,27]. The low-velocity anomalies beneath the northeastern margin of the Tibetan plateau are obvious but not continuous. At long periods of 60-80 s (Figure 7e,f), the Ordos block and the northeastern part of the Alax block show significant high-velocity anomalies. The NW-SE trending low-phase velocities are estimated in the southwestern Zhangbo fault zone (Figures 7e,f and S2). Similar low-velocity anomalies were revealed by Chen et al. [25] and inferred to be related to the upwelling asthenosphere below the Datong volcano.

Azimuthal Anisotropy
We denote the azimuthal anisotropy of the Rayleigh wave phase velocity with a black error bar, where the length and orientation stand for the anisotropic magnitude and direction of the fast axis, respectively. The statistical direction of the fast axis in each tectonic block, e.g., the Alax block, the Qilian orogen and Songpan-Ganzi terrane, the south Ordos block, the north Ordos block, the Lüliang and Taihang mountain, the Bohai Bay basin, and the Datong volcano region, are represented by the red rose diagram at each period ( Figure 7).
Two prominent rotations of azimuthal anisotropy are observed at periods of 15-40 s in the study area (Figure 7a-d). One is centered on the Ordos block, where its surrounding fast axis direction rotated counterclockwise. Another is present under the Datong volcano region (Figures 7c,d and S3). Bounded by the latitude of 38 • , the azimuthal anisotropy of Rayleigh wave phase velocity in the north and south side of the Ordos block is obviously different. The direction of the fast axis of the Ordos block is NE-SW in the north, while it turns to NW-SE in the south. At periods of 60-80 s (Figure 7e,f), the fast axis in the Ordos block seems to be parallel to the direction of the SKS/SKKS splitting [24,38] and the APM [39] (Figure 9). In the northeastern margin of the Tibetan plateau, the fast direction is NW-SE and parallel to the strike of the Qilian orogen (Figure 7), which is consistent with the azimuthal fast direction of the Ps and SKS/SKKS phases splitting [24,40]. The direction of the fast axis in the Trans-North China orogen is NE-SW at periods of 10-30 s and parallel to the strike of the Lüliang and Taihang mountains. At periods of 40-80 s, the fast axis of the southern Trans-North China orogen aligns along a NWW-SEE direction that is different from the direction of the northern Trans-North China orogen where the fast axis trends almost N-S at 40 s period and NNW-SSE at 60-80 s periods. The direction of the fast axis is NW-SE in most of the Bohai Bay basin region. Especially, the magnitude of anisotropy in the southwestern Zhangbo fault zone is very strong at periods of 60-80 s (Figure 7e,f).

Discussion
The Rayleigh wave phase velocity is most sensitive to the S-wave velocity structure of the Earth's medium in the depth range of about 1/3 to 1/2 the wavelength, although it is also affected by the temperature and mineral composition of the medium. The sensitivity kernels of Rayleigh wave phase velocity are shown in Figure 8. According to the sensitivity kernel curves, the short-period Rayleigh waves have a narrow sensitive depth range and a high resolution while long-period Rayleigh waves shows the opposite pattern. The depth of the sensitive kernel overlaps with the increase of the period. The Rayleigh wave phase velocity of 80 s is mainly sensitive to the S-wave velocity structure at 100-150 km. Therefore, maps of Rayleigh wave phase velocity and azimuthal anisotropy (Figure 7) in this study can reliably characterize the structure and deformation within the depth range of 150 km.
The cause of azimuth anisotropy of the Rayleigh wave phase velocity is relatively complex. It is generally believed that the shape preferred orientation and the lattice preferred orientation mainly contribute to the seismic anisotropy in the crust and lithospheric mantle [41,42]. The shape preferred orientation is considered to be caused by geological structures such as oriented microcracks, faults, rifts and orogenic belts [41]. The lattice preferred orientation is produced by the alignment of anisotropic minerals due to structural deformation [42]. Under the regional tectonic stress field, a series of faults and orogenic belts will be formed on the surface of the Earth, resulting in the alignment of microcracks in the superficial crust. Rayleigh wave phase velocity azimuthal anisotropy of the upper crust is controlled by the surface tectonic evolution so that the direction of the fast axis is parallel to the strike of faults or orogenic belts [41,43]. In the mid-lower crust and lithospheric mantle, the Rayleigh wave phase velocity azimuthal anisotropy is the "fossil" anisotropy frozen by the preferentially aligned anisotropic minerals during the recent tectonic movement [44,45]. Quartz, mica, amphibole and other anisotropic minerals are widely distributed in the mid-lower crust. Olivine and pyroxene are the main mineral constituents of the upper mantle. These minerals are anisotropic. Especially olivine has the strongest ability of alignment, which is considered to be the main cause of upper mantle anisotropy [24,42]. The Rayleigh wave phase velocity azimuthal anisotropy of the asthenospheric mantle may be related to the flow of mantle material and is usually parallel to the flow direction [46,47]. Therefore, measurements of Rayleigh wave phase velocity azimuthal anisotropy are helpful to study the regional tectonic evolution and interaction among adjacent tectonic blocks and possible perturbations of the mantle flow.

Lithospheric Thinning and Volcanism of the NCC
Previous studies have shown that the present lithospheric thickness of the NCC gradually decreases from west to east and the lithospheric mantle gradually varies from an old cratonic type to modified and young types [48]. The composition and nature of the Archean mantle lithosphere beneath the eastern NCC has transformed since the Mesozoic, leading to a significant craton rejuvenation [12,48]. This evolution of the mantle lithosphere is thought to cause the geological deformation responses in the crustal and upper mantle, such as the extensional basins and faults in the upper crust, the lithospheric deformation and the intracontinental volcanism [12,13,48].
Most of the eastern NCC is occupied by the Bohai bay basin. The Taihang mountain and the Sulu orogen serve as the boundaries of the Bohai bay basin along its western and eastern edges, respectively (Figure 1). The Bohai bay basin is an extensional basin with several secondary depocenters. The extension of the Bohai bay basin reached a climax between 42 and 32.8 Ma, with NNW-SSE-directed extensional rates of 0.2-0.4 cm/yr [6]. In this study, we observe a low-velocity anomaly at a period of 15 s below the Bohai bay basin that reflects the low-velocity structure of the sedimentary basin. The 15 s fast axis under the Bohai bay basin trends NNW-SSE and is roughly consistent with the direction of the present GPS velocity [9] (Figure 1a, magenta arrows) and the past extensional velocity of the basin [6].
The Trans-North China orogen and western NCC are interpreted to represent the old cratonic stable lithosphere since the Precambrian based on the thick crust [49,50] and lithosphere thickness [51] and low heat flow [52]. However, the recent high-resolution fullwaveform images [23], the detailed lithospheric thickness by S-wave receiver functions [15] and the geochemical study on Cenozoic basalts [53] suggest that the lithosphere was thinned in the eastern NCC, and this thinning may also affect the Trans-North China orogen. The rejuvenated craton is generally characterized by a shallow lithosphere-asthenosphere boundary (LAB) and a low-velocity below the lithosphere-asthenosphere boundary depth. The Rayleigh wave phase velocities of 60-80 s correspond to the S-wave velocity structures at depths of 70-150 km (Figure 8). Referring to the lithosphere-asthenosphere boundary depth (Figure 9b), Figure 7e,f demonstrates the structure of the mid-lower lithosphere and the uppermost asthenosphere beneath the NCC. The distinct high-velocity anomalies beneath the Ordos block and the northeastern part of the Alax block indicate that the thick lithospheric root and stable craton are still preserved in the western NCC. In contrast, the low-velocity anomalies on both sides of the North-South Gravity Line beneath the northern Trans-North China orogen support the assumption that the lithosphere would be modified locally beneath the Trans-North China orogen [15,23,53]. Our results provide a new understanding that this lithospheric thinning may occur only in the northern Trans-North China orogen.
The widespread Early Cretaceous magmatism associated with the Paleo-Pacific plate subduction in the NCC, especially in the Bohai bay basin and its surrounding regions, implies that the thinning of the eastern NCC lithosphere mostly took place during this period [11][12][13]. Beneath the Datong volcano, the prominent low-velocity zone indicates the upwelling of hot materials ( Figure S1). The asthenospheric upwelling probably originates from the mantle transition zone beneath the Datong volcano based on the prominent low P-and S-wave velocity anomalies from the surface down to approximately 660 km [23,54]. Another important phenomenon is that the upwelling asthenosphere would flow horizontally after encountering the upper lithospheric layer, flowing southeastward into the eastern NCC (Figures 7e,f and S2). The hot materials continue to rise by thermally eroding the overlying lithosphere and extending in the lower crust (Figure 7c,d), driving surface extension and volcanism ( Figure 10). This upwelling and its horizontal flow probably causes the observed azimuthal anisotropy [25,55]. The low-velocity anomalies beneath the Datong volcano at periods of 30-40 s (Figure 7c,d) shows the intrusion northwestward in the uppermost mantle with weaker azimuthal anisotropy due to the small range expansion in the strong lithosphere. However, the rotating fast axis around the Datong volcano (Figures 7c,d and S3) agrees well with the results obtained from Ps and SKS/SKKS splitting [22,24], affecting the whole lithosphere, which may represent stress-induced anisotropy caused by the mantle upwelling [22,55]. The low-velocity anomalies of 60-80 s along the southwest side of the Zhangbo fault zone (Figures 7e,f and S2), a noticeable seismic zone with strong azimuthal anisotropy (Figure 7e,f), suggest that the large-scale flow of the upwelling materials enhanced the observed azimuthal anisotropy. Previous studies speculated that the Zhangbo fault zone plays a role as the channel for the upwelling asthenosphere [25,56]. Our results provide new insights into the range of the low-velocity anomaly and the corresponding azimuthal anisotropy, supporting that the upwelling mate-rials of the Datong volcano flow into the low-velocity asthenosphere of the eastern NCC along the southwestern Zhangbo fault zone (Figures 7e,f and S2). Combined with the previous hypotheses [12][13][14][15], that the thinning of the cratonic lithosphere in the eastern NCC resulted from slab rollback related to the subduction of the Paleo-Pacific plate, we further infer that the asthenospheric flow from the mantle upwelling of the Datong volcano also contributed to this thinning ( Figure 10).  [22,26,40] and at period 40 s vs. the SKS/SKKS phases splitting [24,38] (a) Distribution of Moho depth (after Li et al. [50] and Xu et al. [22,40]). The black arrows indicate the absolute plate motion (APM) direction of the GSRMv2.1 model [39]. The abbreviations of geological units are defined in the caption of Figure 1. (b) Distribution of lithosphere-asthenosphere boundary depth (after Chen et al. [15]).

Deformation Patterns of the Northeastern Margin of the Tibetan Plateau and the Western NCC
The continental collision between the Indian and Eurasian plates has resulted in the uplift and outward growth of the Tibetan plateau for~50 Ma [17,18]. Counteracted by the western NCC, the northeastern margin of the Tibetan plateau has a steep topography and strong lithospheric deformation. Thus, it is essential to propose a suitable model for understanding the mechanism of the growth and deformation in the northeastern margin of the Tibetan plateau. The crustal channel flow [5,57] and vertical coherent deformation of the entire lithosphere [58,59] have been debated recently to predict the very different seismic anisotropy and velocity structures in the northeastern margin of the Tibetan plateau [24,40,60,61].
According to the continuous low-velocity layer in the mid-lower crust associated with the observed seismic anisotropy, the crustal channel-flow model may explain the lateral crustal deformation along the southeastern margin of the Tibetan plateau [62][63][64]. In the Songpan-Ganzi terrane, the crustal thickening and deformation were also considered to be driven by the crustal channel flow based on the existence of regional low-velocity zones [60,61]. However, the low-velocity anomalies of Rayleigh wave phase velocity at periods of 30-40 s (Figure 7c,d), sampling the mid-lower crust (Figures 8 and 9), are discontinuous, coinciding with the distribution of S-wave velocity revealed from ambient noise tomography [65]. The discontinuous intracrustal low-velocity anomalies are most likely caused by partial melting and may facilitate the shortening between the Tibetan plateau and the western NCC [65][66][67]. In particular, the coherent crustal and upper mantle anisotropy ( Figure 9) estimated from the Ps and SKS/SKKS phases splitting contradicts the crust-mantle decoupling due to crustal channel flow [24,40]. The statistical rose diagrams of the northeastern margin of the Tibetan plateau at periods of 15-80 s (Figure 7) show a NW-SE trending azimuthal anisotropy, which agrees well with the S-wave splitting in the crust and upper mantle (Figure 9). Our results in this study conform to the prediction of a vertical coherent deformation of the entire lithosphere. Moreover, the prominent lowvelocity anomalies at 20-40 s periods are not interrupted along the collision belt between the northeastern margin of the Tibetan plateau and the western NCC, and invade partially the interior of the Alax block. The encroachment of the Tibetan plateau depends on the strength of the surrounding blocks. This may explain the observed low-phase velocity anomalies in this study. Therefore, we propose a scenario where the vertical coherent lithospheric deformation in the northeastern margin of the Tibetan plateau is produced by the northeastward expansion of this domain upon the western NCC.

Rotation of the Ordos Block
The Ordos block located at the western NCC is a stable block but its peripheral areas with high-strain rate and seismicity have been active since the Cenozoic [19,68,69]. The Ordos block is surrounded by four extensional rift systems, e.g., the Weihe graben to the southeast, the Fenhe graben to the east, the Hetao graben to the north and the Yinchuan graben to the northwest, and the Liupan fold and thrust belt to the southwest (Figure 1). Therefore, the Ordos block and its adjacent regions play a significant role in the present tectonic deformation. A kinematic counterclockwise rotation of the Ordos block was proposed to accommodate tectonic driving forces from its surrounding regions [19,24,68].
The recent field investigations [70], GPS observations [68,69] and high-resolution remote sensing imagery [71] reported that a series of right-lateral strike-slip faults, with a rate of about 0.8-2.6 mm/yr, occur along the extensional grabens surrounding the Ordos block. In addition, the earthquake focal mechanism solutions in these grabens support the pattern of superficial crustal deformation mentioned above [8]. The rotated direction of Paleo-Pacific plate subduction (Figure 1a) might be coupled to the crustal deformation pattern that provide the extensional-transtensional stress environment for promoting the counterclockwise rotation of the Ordos block [6,10]. This counterclockwise rotation pattern was also observed by the fast axis directions of the Ps and SKS/SKKS phases splitting ( Figure 9) in the mid-lower crust and upper mantle [22,24,40]. It is concluded that the fast axis directions from S-wave splitting are produced by the extrusion of the lithospheric material of the Tibetan plateau [22,24,40]. The azimuthal anisotropy of the Rayleigh wave phase velocity around the Ordos block in the crustal and mantle lithosphere, corresponding to the periods of 15-40 s (Figure 7a-d), is quite consistent with the proposed model of rotation. The low Rayleigh wave phase velocity (Figure 7a,b) as well as the S-wave velocity [27] anomalies along the periphery of the Ordos block are the manifestation of the weak crustal structure that is more likely to produce the transtensional shear deformation inside the crust of the extensional grabens. Our results seem to support the previous hypothesis that the eastward extrusion of the Tibetan plateau and the subduction of the Paleo-Pacific plate caused the shortening rate as high as 3.0 mm/yr in the Liupan thrustfold belt and extensional rate of 1.3-1.7 mm/yr across the Weihe and Fenhe grabens [69], respectively, and jointly contributed to the rotation of the Ordos block [20][21][22].

Conclusions
In this study, Rayleigh wave phase velocity and azimuthal anisotropy at periods of 10-80 s was estimated from the teleseismic data recorded by 388 permanent stations in the NCC and its adjacent areas. Our results reveal the detailed crustal and upper mantle structures of the NCC and its adjacent areas, which provides a new reference for understanding the lithospheric deformation. The prominent low-velocity anomalies with strong NW-SE azimuthal anisotropy from the Datong volcano to the eastern NCC indicate that the lithospheric thinning of the eastern NCC may be also caused by the asthenospheric upwelling of the Datong volcano, except for the subduction of the Paleo-Pacific plate. Our observed coherent crustal and upper-mantle azimuthal anisotropy beneath the northeastern margin of the Tibetan plateau conforms to the prediction of a vertical coherent deformation of the entire lithosphere. Combined with the previous geological and geophysical studies, the eastward extrusion of the Tibetan plateau and the subduction of the Paleo-Pacific plate may contribute to the counterclockwise rotation of the Ordos block.  Figure S3: Refined Rayleigh wave phase velocity and azimuthal anisotropy at 30 and 40 s periods indicating the rotating fast axis around the Datong volcano. The abbreviations of geological blocks are defined in the caption of Figure 1. The black error bar denotes the magnitude and fast wave direction of azimuthal anisotropy. The average Rayleigh wave phase velocity (Ave.) at each period is shown in the corresponding panel. Figure S4: Refined Rayleigh wave phase velocity and azimuthal anisotropy at 30 period indicating the discontinuous low velocities in the northeastern margin of the Tibetan plateau. The abbreviations of geological blocks are defined in the caption of Figure 1. The black error bar denotes the magnitude and fast wave direction of azimuthal anisotropy. The average Rayleigh wave phase velocity (Ave.) at 30 s period is shown in the figure.
Author Contributions: X.X. and D.Z. conceived and designed the paper; X.X. inversed the results and wrote the manuscript; X.H. and X.C. participated in the data processing; D.Z. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.