Present‐Day Crustal Deformation of Continental China Derived From GPS and Its Tectonic Implications

We process rigorously GPS data observed during the past 25 years from continental China to derive site secular velocities. Analysis of the velocity solution leads to the following results. (a) The deformation field inside the Tibetan plateau and Tien Shan is predominantly continuous, and large deformation gradients only exist perpendicular to the Indo‐Eurasian relative plate motion and are associated with a few large strike‐slip faults. (b) Lateral extrusions occur on both the east and west sides of the plateau. The westward extrusion peaks at ~6 mm/yr in the Pamir‐Hindu Kush region. A bell‐shaped eastward extrusion involves most of the plateau at a maximum rate of ~20 mm/yr between the Jiali and Ganzi‐Yushu faults, and the pattern is consistent with gravitational flow in southern and southeastern Tibet where the crust shows widespread dilatation at 10–20 nanostrain/yr. (c) The southeast borderland of Tibet rotates clockwise around the eastern Himalaya syntaxis, with sinistral and dextral shear motions along faults at the outer and inner flanks of the rotation terrane. The result suggests gravitational flow accomplished through rotation and translation of smaller subblocks in the upper crust. (d) Outside of the Tibetan plateau and Tien Shan, deformation field is block‐like. However, unnegligible internal deformation on the order of a couple of nanostrain/yr is found for all blocks. The North China block, under a unique tectonic loading environment, deforms and rotates at rates significantly higher than its northern and southern neighboring blocks, attesting its higher seismicity rate and earthquake hazard potential than its neighbors.


Introduction
Continental China is located on the eastern Eurasian plate and is one of the most tectonically active regions in the world due to widely distributed plate boundary deformation Yin, 2010;Zhang et al., 2004). The northward indentation of the Indian plate into the Asia continent caused widespread deformation, including mountain building of the Himalaya and Tien Shan orogenies and crustal shortening and uplifting of the Tibetan plateau, accompanied with large-scale eastward extrusion of the plateau materials (Molnar & Tapponnier, 1975;. Along the east margin of the continent the Pacific and Philippine Sea plates subduct underneath the Eurasian plate, causing reactivation of the North China craton (Zhu et al., 2011). Due to this process, the east Eurasian lithosphere broke into tectonic blocks of various sizes, and active faults developed at the block boundaries, whose ruptures during earthquakes have been constant threats to human society and civilizations (Gu, 1989; Figure 1). Understanding the crustal deformation processes and mechanisms is therefore important not only for Earth science research but also for mitigation of geohazard, particularly the earthquake hazard.
Crustal deformation research has been revolutionized during the past three decades, owing in particular to the utilization of space geodesy techniques such as GPS and InSAR for crustal deformation monitoring (e.g., Feigl et al., 1993;Massonnet et al., 1993). Early campaign GPS surveys in China and its neighboring regions conducted in the 1990s resulted in the documentation of deformation fields in various regions, such as across the Himalaya (e.g., Banerjee & Bürgmann, 2002;Bettinelli et al., 2006;Bilham et al., 1997;Jouanne et al., 1999;Paul et al., 2001), North China (e.g., Shen et al., 2000), the Altyn Tagh fault (e.g., Bendick et al., 2000;He et al., 2013;Shen et al., 2001), and the Tibetan plateau (e.g., Chen et al., 2000;Chen et al., 2004;Wang et al., 2001). Establishment of the Crustal Movement Observation Network of China (CMONOC) enabled large-scale monitoring of deformation field over the entire continental China (Niu et al., 2005;Wang et al., 2003). The data from the CMONOC network have been periodically analyzed, providing updated results of crustal motion velocity field for deformation research at various scales (e.g., Gan et al., 2007;Liang et al., 2013;Rui & Stamps, 2019;Shen et al., 2005;Wang et al., 2003;Wang et al., 2017;Zhang et al., 2004;Zheng et al., 2017). Although these velocity results have been widely used and provided crucial constraints for tectonic deformation and other studies, the usefulness of these solutions, however, has sometimes been limited due to relatively short time spans, uneven spatial coverage, and/or contamination by coseismic and postseismic deformation of large earthquakes.
In this study, we include data which are more than that have been used in all of the previous studies and produce a secular velocity solution in continental China by using a rigorous data processing algorithm and carefully modeling the coseismic and postseismic deformation caused by large earthquakes. Using the velocity data set combined with other data sets covering surrounding regions from previous studies, we derive secular strain rates, block motion rates, and major fault slip rates and quantify the patterns of deformation. The velocity solution with a high spatial resolution and consistency unmasks certain deformation patterns, which were previously not clear or known. Such patterns include westward extrusion of the Tibetan plateau, rotation of the southeast borderland of the Tibetan plateau, internal deformation within the blocks of eastern China, and strain distributions along faults such as the Zhangjiakou-Bohai seismic zone (ZBSZ) in North China. These results offer new insights for tectonic interpretations and seismic hazard assessments.

GPS Data
In this study, we gather and process GPS data from multiple sources (listed in Table 1 and shown in Figure  S1), which are described below. surveyed once every 2 to 3 years from 1999 to 2007. This network was expanded in Phase II in 2009, and 230 continuous and~1,000 "regional network" sites were added. The "basic network" sites in Phase I were either upgraded to continuous sites or downgraded to "regional network" sites in Phase II (the "basic network" and "regional network" sites are collectively referred as campaign sites hereafter). Thẽ 2,000 CMONOC I and II campaign sites were surveyed entirely once every 2 years from 2009 to 2015. During each campaign, 4 days of data were collected per site. These data were analyzed in some previous studies such as Wang et al. (2017) and Zheng et al. (2017). We add more data observed at the campaign sites of this network, located in North China, Yunnan, the eastern Tien Shan, and areas affected by large earthquakes (e.g., the 2001 Kokoxili and 2008 Wenchuan earthquakes). Inclusion of these data has enhanced the capability of detecting transient deformation caused by earthquakes and/or other sources. 2. Data from densified regional campaign GPS networks. These networks are usually located in areas with active seismic activity, such as the "North China/Capital Circle" (Shen et al., 2000), Shanxi rift zone, Sichuan, and Yunnan regions. These sites were usually built with deep-root monuments and compulsory antenna alignment devices, which are similar to that of the CMONOC campaign sites, and were observed regularly. These networks have at least five epochs of campaign surveys, and the observation time spans range between 5 and 17 years. In addition, we also gather data from several campaign surveys of the early years in the 1990s. These include nearly 100 sites measured as geodetic control sites in 1991-1994 and surveyed repeatedly by multiple organizations to study crustal deformation in China. The data used by Bilham et al. (1997), Bendick et al. (2000), Wang et al. (2001), Shen et al. (2001), Chen et al. (2004), and He et al. (2013) are also incorporated. More than 4,000 sites are included in our data processing, although not all of the sites have sufficient observations to derive the motion velocity with sufficient precision (criteria for site selection are specified below). In the early surveys, the occupation time at each site was short, and the number of satellites observed was also limited. Processing all the synchronous data together increase the phase double-difference combinations and strengthen the overall solutions and tie to the IGS global network. 3. Data from regional continuous GPS sites. Most of these sites are located in the Sichuan basin, Nepal Himalaya (Ader et al., 2012), and central Asia regions, and have long enough (several years of) observation histories. Raw data for the latter two regions are obtained from the UNAVCO data archives. We also use data from 29 sites distributed widely within the Tibetan plateau, and the data are provided by the Institute of Tibetan Plateau Research, Chinese Academy of Science and China Meteorological Administration.

Data Processing
We adopt a unified strategy to process/reprocess all raw GPS data using the GAMIT software (Herring et al., 2010a) to ensure the quality and homogeneity of the solutions. The major geophysical models and parameters used in the processing are listed in Table S1 in the supporting information (Boehm et al., 2006(Boehm et al., , 2007Lagler et al., 2013;Lyard et al., 2006;Saastamoinen, 1973). Instead of using the IGS precise orbits or daily solutions of global tracking sites from the IGS Analysis Centers, we incorporate~100 globally distributed IGS sites in our data processing to estimate the Earth orientation (polar motion and ut1), satellite orbits, and site position parameters together. This approach has the advantage of minimizing potential biases, such as the ones introduced by changes of geophysical models and parameters employed in IGS orbit production (Shen et al., 2011). For the earlier data of the 1990s, including both the global and regional data in a single solution can improve the accuracy of satellite orbits, particularly for East Asia, which in turn improves the accuracy of site positioning.
The GAMIT software adopts the double-difference approach to model the carrier-phase data, which limits the number of sites to be processed in a single solution with some computation power. For the data observed between 1991 and 1994, with a relatively smaller number of sites observed in each day, all the campaign data are processed together with the global IGS data to obtain a daily solution directly. For the data observed after 1994, the data are processed in more than one group divided by the regional site locations, with a few continuous sites included in each of the campaign solutions. The daily solutions are loosely constrained, and the parameters on Earth orientation, satellite orbits, and site positions are output along with their variance/covariance matrices. The solutions of the same day are subsequently combined using the GLOBK software (Herring et al., 2010b).

Site Displacement Modeling
We use the QOCA software (http://gipsy.jpl.nasa.gov/~qoca) to model GPS site displacements, which involve all the campaign sites and~100 continuous sites (~60 sites in China and~40 sites globally distributed) selected according to the observation history and spatial distribution. QOCA employs the Kalman sequential filter technique, whose advantage is that it can adequately accommodate temporal correlation of errors in the data time series and incorporate full variance/covariance information with the data input.
Here, we aggregate all the daily solutions into a single solution in which site position time series are modeled as a function of the initial position x 0 at time t 0 , secular velocity v, coseismic (also instrumental if needed) step D i , and postseismic displacement A i , in the form of where the first sum includes earthquake, instrumental, and other steps at a time t i , and n c is the total number of steps; the second sum includes the postseismic displacements of earthquakes in logarithmic form, t i and T i are earthquake occurrence time and decay time constant, respectively, and n p is the total number of earthquakes. H(t − t i ) is a Heaviside function. The decay time constant T i is determined to be 100 days by trial and error for all the large earthquakes. A few early postseismic data epochs may be discarded for the near-field sites if they appear to be inconsistent with the assigned logarithmic function.
For the campaign sites affected by one or more earthquakes, the observations may not be enough to constrain all the parameters associated with the earthquakes. We therefore use a priori offsets as constraints for coseismic displacements and a fraction of the a priori offsets (20% for the 2001 Kokoxili and 2008 Wenchuan earthquakes and 50% for the 2004 Sumatra and 2011 Tohoku-Oki earthquakes, respectively) for the postseismic displacements. All of the constraints are assigned with uncertainties which are in the form of ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0:6c i ð Þ 2 þ 0:3c j À Á 2 þ 0:3c k ð Þ 2 q for the ith component, with c i , c j , and c k being the a priori values for any combination of the east, north, and up components. No a priori constraints are assigned for the coseismic and postseismic parameters of the continuous sites. More than 20 earthquakes are considered for various deformation effects in our analysis ( Figure 2). Details about assignments of a priori coseismic displacements for sites and identifications of sites affected by earthquakes are described in the supporting information (Hsu et al., 2009;Jiang et al., 2014;Shen et al., 2009;Sun et al., 2013;Wan et al., 2008;Wang et al., 2006Wang et al., , 2011.
Furthermore, we enlarge the errors for the daily solutions by a factor of two to account for some of the external errors that are not well calibrated, such as systematic errors between different antenna types used and possible antenna setup errors, etc. Also, we introduce random walk type perturbations for site positions as 0.2, 0.2 and 1.0 mm/ ffiffiffiffi ffi yr p for the east, north, and up components respectively, to accommodate timedependent systematic errors, especially for the vertical component. Finally, we select 17 IGS sites with relatively uniform global distribution to align their velocities to that of ITRF2008 by transformation, to tie the solution to the ITRF2008 reference frame.
As mentioned above, only a fraction (~100) of the continuous sites are included in the QOCA solution, and the data from these sites are also decimated according to the epochs of campaign observations. We therefore produce another solution for all the continuous sites using all the daily data of 1999-2016. We use~60 IGS sites around the globe to tie their loosely constrained daily solutions to the ITRF2008 reference frame by the seven-parameter transformation of site positions, and then model the position time series using equation ((1)) with addition of annual and semiannual variation terms (in sinusoidal form). Comparing the two sets of velocity solutions of~100 common sites, we find that the results agree within 0.1 mm/yr in general with no systematic deviation. The velocity solution of the continuous sites derived from the position time series is then merged with the QOCA velocity solution to produce the final solution.

Crustal Motion Solution
The final solution of crustal motion is obtained after several iteration runs, each time the site position time series are inspected for outliers, and the outliers spotted are either fixed with corrected metadata or discarded if the problems are not fixable. The known problems include poor data quality (such as observations made using the GeoTracer receiver/antenna units in 1999 and 2000), disturbance due to anthropogenic activities (such as land subsidence caused by groundwater and oil pumping and coal mining in the North China plain), and short data time span (<3 years for sites located in eastern China, and <2 years for sites located in western China). In addition, we find that 28 CMONOC I campaign sites located in the eastern Tien Shan region show unusual displacements before 2008, and the problem may be due to systematic biases introduced by using different receiver/antenna units between surveys, particularly some substandard antennas. We then remove the troublesome epochs of data in derivation of the velocities for the affected sites.
Some of the CMONOC II continuous sites were upgraded from the CMONOC I campaign sites using the existing or rebuilt monuments, and we estimate the velocities of the two phases before and after the upgrade independently and tie them together later. We also tie the velocities of any two sites which are located within 1.5 km distance, since no discernable tectonic deformation is expected between the site pairs at such a close distance. The 89 site pairs are listed in Tables S3, and only one site for each pair is retained in the final solution. Our final solution contains 2,403 site velocities (including 38 sites located outside of the study region for the transformation of reference frames in future application), and the data are provided in Table S4.
For the convenience of tectonic interpretation, we transform the velocity solution from ITRF2008 to a stable Eurasia reference frame by solving the rotation Euler pole of 10 sites located in northern Europe and Siberia. Our estimated Euler vector for the Eurasian plate (−0.087, −0.514, and 0.741 mas/yr) agrees well with that of the ITRF2008 plate motion model produced by Altamimi et al. (2012) (−0.083, −0.534, and 0.750 mas/yr), and the slight difference in rotation rate corresponds to less than 1.0 mm/yr velocity difference at the Earth's surface in our study region and has no effect on the deformation analysis. The postfit residuals of the velocities are less than 1.0 mm/yr for the 10 fiducial sites, which demonstrates that the region occupied by these sites is indeed tectonically stable and our result is with high accuracy. We further incorporate several velocity solutions from previous studies (Kreemer et al., 2014  provide better coverage of the tectonic deformation field in the surrounding continental regions of China, especially the Himalaya and western Tien Shan regions. We transform these velocity solutions to coincide with our solution by rigid body rotation, minimizing postfit residuals at the common sites. The combined velocity field with respect to the Eurasian plate is presented in Figure 3, and the additional data are listed in Table S5.

Strain Rates
To visualize the crustal deformation field in continental China, we interpolate the horizontal GPS velocity field to derive a continuum strain rate field using the method developed by Shen et al. (1996Shen et al. ( , 2015. For a discrete data set such as GPS velocities, smoothing is required for data interpolation. In this method the degree of smoothing varies according to the in situ data strength, and is regulated with optimal selection of the characteristic distance σ. As shown in Figure S3a, the characteristic distance σ is closely related to the spatial density of GPS sites (Figures S1 and 3) and reflects the spatial scale over which strain rates are averaged. In this study, σ is smallest (~40 km) along the ZBSZ and the western Tien Shan, and is less than 100 km along the Shanxi rift zone, the Qilian Shan-Qinling region, the southeast borderland of the Tibetan plateau (i.e., the southwestern Sichuan and Yunnan regions), and most of the Himalaya. GPS site spacing is indeed densest in these areas ( Figure S1). On the other hand, σ is close to 150 km for most of South China, the Tarim basin, and the Alashan region, and greater than 250 km in NE China. The solutions for crustal deformation for these regions are not compromised much, however, due to very little tectonic Figure 2. Earthquakes detected and modeled in this study. The beachballs are earthquake focal mechanisms from the Global GMT catalog (https://www.globalcmt.org/CMTsearch.html). Colored lines denote the spatial regions affected by the quakes whose focal mechanisms share the same color. The gray area in North China denotes the subsidence region. deformation therein. Characteristic distance σ is close to 200 km in the Naqu region and 150 km for the rest of northern Tibet, where tectonic deformation is active, but the strain rate estimates have to be averaged over large areas due to sparse GPS site distribution. For other regions, σ ranges from 100 to 150 km, reflecting the moderate degree of data smoothing therein.
The horizontal strain rate field is shown in Figure 4a, in which the maximum shear strain rate τ max is plotted at the background, and the two horizontal principal strain rates τ 1 and τ 2 are displayed as conjugate pairs of vectors. Maximum convergence appears across the Himalaya, with about 60-80 nanostrain/yr normal convergence across most segments of the Himalaya and up to 120 nanostrain/yr across the thrust fault system of the eastern Himalaya. About 50-60 nanostrain/yr N-S convergence is found across the central Tien Shan. Little E-W deformation is detected within the two orogenies, suggesting pure thrust motion across the Himalaya Main Thrust (HMT) and the north and south frontal faults bounding the Tien Shan ranges. A maximum shear strain rate of 40-60 nanostrain/yr is found along the Xianshuihe-Xiaojiang fault system, which delineates the north and east boundaries of the crustal materials undergoing large-scale clockwise rotation around the eastern Himalaya syntaxis (EHS). This left-lateral shear zone cuts through the Red River fault and deflects toward the southwest . A significant shear strain rate of~30 nanostrain/yr is also found in the Lijiang and Longling fault zones, which were hit by multiple strong earthquakes during the last century (Wen et al., 2008).
The strain rates in southern Tibet are characterized by an almost equal amount of N-S shortening and E-W extension of~20 nanostrain/yr. The central Tibet region shows more WNW-ESE extension than NNE-SSW shortening at most places from west to east. The deformation is broadly distributed, but in western Tibet much of that is attributed to large-scale interpolation and needs to be cautiously interpreted ( Figure S3a). The strain rates in the Qaidam basin and the Qilian Shan region demonstrate NE-SW shortening and NW-SE extension, with the shortening rates significantly exceeding the extension rates. This deformation pattern suggests contraction across a series of thrust faults and basin-ranges, such as the North Qilian and

Journal of Geophysical Research: Solid Earth
Danghe-Nanshan mountains, the faults in between, and the basin bounding faults located at the north and south rims of the Qaidam basin . About 25 nanostrain/yr of shear strain is also detected across the Haiyuan fault, revealing its sinistral slip component of deformation. Sinistral strain on the order of 20 nanostrain/yr is shown along the Altyn Tagh and East Kunlun faults respectively, but this is likely an underestimate of the local deformation due to the lack of data spanning most of the fault length ( Figure S3a). Patches of noticeable strains appear in the North China plain, with the largest of~15 nanostrain/yr sinistral shear along the ZBSZ near the Bohai Bay. This high shear strain rate may also partially reflect the postseismic deformation of the 1978 Mw7.8 Tangshan earthquake (Figure 1). For other regions, such as in NE and South China and in the interior of the Ordos and Tarim basins, deformation rates are relatively low, on the order of a few nanostrain/yr or less.
The result of dilatation rate τ dil (τ dil = τ 1 +τ 2 ) shows strong areal contraction, up to 80 nanostrain/yr and 40 nanostrain/yr, along the Himalaya and Tien Shan orogenies, respectively (Figure 4b), which is expected due to crustal shortening across the thrust faults bounding these mountains. Most of the central Tibet experiences little areal change. Dilatation of 10-20 nanostrain/yr is found along the central and western Bangong-Nujiang suture zone and in the Naqu region. Also an expansion belt with dilatation rate of 10-20 nanostrain/yr is placed between the Jiali and Ganzi-Yushu-Xianshuihe faults, which extends south and southwestward, passing through the Lijiang fault and reaching the Longling and Dayingjiang faults near the Yunnan and Myanmar border ( Figure 4b). We also plot the second invariant of the horizontal strain rates Figure S3b which is the measure of full strain rate field. A comparison shows that the pattern of τ 2inv is very similar to that of τ max , which is not surprising as the deformation field is overall dominated by shear strain. Figure 4c shows rotation rates in the continental region of China which are quite interesting. It is understood that parts of the rotation pattern are associated with simple shear associated with strike-slip faulting, such as the sinistral slips along the Ganzi-Yushu, Xianshuihe-Xiaojiang, East Kunlun, Altyn Tagh, and Haiyuan faults, and dextral slips along the Longriba and Jiali faults. In central Tibet, the clockwise and counterclockwise rotation regions border along the Ganzi-Yushu fault and its NW extension. Large-scale clockwise rotation of up to 90 nanorad/yr (5.2°/Myr) is also found around the EHS relative to the Eurasian plate. Maximum counterclockwise rotation of~50 nanorad/yr (2.9°/Myr) occurs along the Xianshuhe-Xiaojiang fault system, which is in sharp contrast to the largest clockwise rotation at its west. More rotation results for the Tibetan borderland and blocks located outside of the Tibetan plateau will be presented below.

Localized and Broadly Distributed Deformation
In the western part of continental China, the deformation field is rather complex. We make a series of velocity profiles to examine the regional deformation fields. As the major driving mechanism of the deformation field is the collision between the Indian and Eurasian plates, we rotate the velocity field to be referenced to the Indian plate using the sites on the Indian subcontinent, and decompose the velocities into two components: along (or tangential to) and normal to (or perpendicular to) the relative plate motion. We further group the velocities into 1.5°bins along small circles centered at the India-Eurasia Euler pole of relative motion (Figure 5), and plot them in profiles as shown in Figure 6.
In the westernmost part of the network, profiles A-A′, B-B′, and C-C′ show crustal motion across the Hindu-Kush and western Tien Shan regions. Little convergence is found across the HMT. Instead,~20 mm/yr shortening is broadly distributed within the~1,300 km range, spanning the Hindu-Kush Thrust (HKT), Main Pamir Thrust (MPT), and Tien Shan ranges. The convergence rate changes from west to east. They arẽ 10,~7, and~3 mm/yr across the HKT, MPT, and Tien Shan ranges along the A-A′ profile, and~5,~8, and~8 mm/yr across the same structures along the C-C′ profile, respectively. The HKT and MPT undergo sinistral and dextral motions, respectively, accommodating~6 mm/yr westward extrusion of the Pamir plateau crust in between.
Within the interior of the Tibetan plateau, N-S convergence evolves remarkably from west to east. It is almost entirely concentrated across the HMT (~15 mm/yr) from profiles D-D′ to F-F′ (76-80°E), with little shortening within the plateau north of the mountain range. From 80°E (profile G-G′) eastward, convergence increases gradually and becomes more broadly distributed, and is accommodated across a cluster of conjugate faults in central and northern Tibet. The overall convergence is close to~40 mm/yr across the eastern 10.1029/2019JB018774

Journal of Geophysical Research: Solid Earth
part of the plateau, with~20 mm/yr convergence from the HMT to the Jiali fault along profile N-N′, and more gentle convergence north of the Jiali fault. The Tarim basin is near rigid, and minor shortening of 1-2 mm/yr is found within the basin. North of the basin, convergence is distributed across the Tien Shan and multiple fold and thrust fault systems north of the Tien Shan. The convergence across the Tien Shan decreases from~10 mm/yr at 75°E to~5 mm/yr at 85°E, and minor shortening of 1-2 mm/yr is found within the Jungger basin. About 2-3 mm/yr convergence takes place obliquely across the Altai fault located northeast of the Jungger basin.
The velocity component normal to the India-Eurasia relative plate motion direction is mainly associated with the E-W extension in the interior of the Tibetan plateau. The westward extrusion rate peaks at~6 mm/yr, which is lower than the eastward extrusion rate of~20 mm/yr, and the pivot boundary between the westward and eastward extrusion is located at about 82°E. In the central and eastern parts of the plateau, sinistral and dextral shears appear in the regions north and south of the Ganzi-Yushu fault, respectively, to accommodate the eastward extrusion of the plateau. These shear motions are smooth for most part of the region except across several major faults. About 4 mm/yr sinistral shear is found across the Ganzi-Yushu fault. The East Kunlun fault slips left-laterally with a rate of~12 mm/yr at 94°E, which diminishes gradually eastward to~10 mm/yr at 98°E and down to~5 mm/yr at 102°E ( Figure S4). Transpressional deformation is broadly distributed across the Danghe-Nanshan and Haiyuan faults, with the sinistral and shortening rates of~7 and~5 mm/yr across the Danghe-Nanshan fault and~7 and~6 mm/yr across the Haiyuan fault, respectively ( Figure S5). Minimal deformation is detected across the Karakash fault, and no significant velocity offset is found close to the fault ( Figure S6). The Altyn Tagh fault demonstrates sinistral shear from its west end to east end, with slip rates of~11 mm/yr at 83°E,~9 mm/yr at 89°E, and~5 mm/yr at 93°E, respectively ( Figure S6). At the southwestern edge of the plateau the Karakorum fault slips right laterally, and the slip rate for the southern fault segment is~3 mm/yr, while the northern segment shows no obvious slip ( Figure S7). About 3 mm/yr sinistral shear is detected across the Fuyun fault and a broad range (~500 km) north of the fault in the Altai mountains.
One of the most intriguing features of crustal deformation in continental China is the clockwise rotation of southeastern Tibet around the EHS. This rotation, as some previous studies envisioned, is relative to the deeper part of the Earth interior, whether it is the asthenosphere attached to the stable part of the Eurasian plate

Journal of Geophysical Research: Solid Earth
or eastern China (e.g., Wang et al., 2013). To visualize the localized rotation pattern, we convert the velocity field to be referenced to South China (Figure 7a). Further, we divide the rotating terrane of southeastern Tibet into eight fan-shaped sections and plot velocity profiles of sites located within the sections in Figure 7b. The horizontal velocities of sites in the fan-shaped sections are decomposed into radial (centered at the rotation pole) and tangential components, and the tangential components for sites located in the rotating terrane are expected to show a linear trend if the terrane material is undergoing a rigid body rotation with respect to South China.
The first order feature captured from the velocity tangential components of the fan-shaped sections is the strong sinistral shear along the Xianshuihe-Anninghe-Zemuhe-Xiaojiang fault system (shown in left panel of Figure 7b). The sinistral shear across the southeast section of the Xianshuihe fault is about 9-10 mm/yr and is localized within a relatively narrow range, suggesting a single strand of fault for this section of the fault system. In contrast, broad sinistral shears of~10 and~9 mm/yr are distributed within 100-150 km range across the Anninghe and Zemuhe sections of the fault system, respectively, suggesting that the Daliangshan fault, located to the east of the Anninghe and Zemuhe faults, is also experiencing sinistral shear. A sinistral shear zone of~100 km width is detected spanning the Xiaojiang fault with an overall rate of~8 mm/yr. This shear movement extends southwestward across the Red River fault zone and spreads out into broader areas with a reduced rate of~5 mm/yr. The width of the deformation zone is broadened to~200 km southwest spanning the Menglian and Mengxing faults.

Journal of Geophysical Research: Solid Earth
Within the rotating terrane, the sites located closer to the rotation center show a zone of gentle dextral shear (profiles A-A′ to F-F′), corresponding to the deformation field across the Jiali, Nujiang, and Lancongjiang fault zones. The sites located between the dextral and sinistral shear zones, however, share a common tangential motion with respect to South China within each cross section, suggesting translation motion than rotation. Because the translation motion direction changes gradually along the rotation, the large-scale translation deems to be realized by rotation and translation of a cluster of microblocks separated by active faults within the rotating terrane.
Compared with the tangential components, velocities of the radial components are scattered, especially for the ones across the southern (E-E′ to H-H′) profiles (right panel of Figure 7b). In the eastern part of the rotating terrane, about 3-4 mm/yr shorting is found across the Anninghe, Zemuhe, and Daliangshan fault system (C-C′ and D-D′ profiles). In the southwest of the region, about 5 mm/yr extension is detected across a group of faults such as the Lijiang, Dayingjiang, Longling, and Wanding faults (G-G′ to H-H′ profiles). The E-E′ and F-F′ profiles obliquely traverse the Red River fault, and no obvious deformation is identified on the fault.

Block Motions
From the strain rate map, it is evident that the deformation pattern in continental China is predominantly block-like outside of the Tibetan plateau and Tien Shan regions. To quantify the block motion rates, we delineate seven blocks outside of the Tibetan plateau and Tien Shan regions: the South China, North China, NE China (which is part of the Amurian plate; Heki et al., 1999), Ordos, Alashan, Tarim, and Jungger blocks, and compute their translation and rotation rates with respect to the Eurasian plate and their internal uniform strain rates. The results are plotted in Figure 8 and documented in Table 2. From these results we find the following: (1) The northward motion of the Tibetan plateau pushes the Tarim block moving northward at a rate of~13.8 mm/yr and the Jungger block moving north-northeastward at a rate of~6.4 mm/yr, respectively, and the differential motion is absorbed by crustal shortening across the Tien Shan ranges and clockwise rotation of the Tarim block. The rotation rate of the Tarim block is~11.0 nanorad/yr, to accommodate the eastward increase of N-S shortening across the Tibetan plateau and the westward increase of N-S shortening across the Tien Shan ranges (Figures 4a and 5a). (2) The eastern part of continental China escapes east-southeastward, with the rate increasing gradually from north to south: 1.7 mm/yr for NE China,~5.5 mm/yr for North China, and~7.5 mm/yr for South China, respectively. These three blocks also rotate counterclockwise at rates of~0.8,~1.5, and~0.2 nanorad/yr, respectively.
(3) The Ordos block moves together with the North China block, but with a rotation of~3.0 nanorad/yr counterclockwise. This rate is about twice as much as that of the North China block (~1.5 nanorad/yr), which possibly resulted from the Tibetan plateau's eastward push at its southwest border. (4) The Alashan block rotates clockwise at a rate of~1.8 nanorad/yr, which reflects a transition from the clockwise rotation (~11.0 nanorad/yr) of the Tarim block at its west to the counterclockwise rotation (~3.0 nanorad/yr) of the Ordos block at its east. Table 2 lists the mean strain rate estimates and their uncertainties for the blocks. The result shows internal strain rates at the level of a couple of nanostrain/yr for all of the blocks. Most of the strain rate estimates are much greater than their uncertainties, suggesting significance of the results. However, because errors of GPS velocities are usually strongly correlated within a block, the strain rate uncertainties derived from least squares regression are therefore grossly underestimated when such correlations are not accounted for in the estimation process. We assume an average error correlation of 60% between the sites (the value is obtained by inspecting the variance/covariance of the velocity solution) and use the F test to evaluate significance of the result, following Shen et al. (1996) and Wang et al. (2003). The F test confirms significance of the internal strain rate estimates within all the blocks at the 99% confidence level.
Crustal deformation in eastern China is much smaller than in western China, and is mostly broad and subtle. We plot six velocity profiles across the blocks in eastern China to investigate the deformation activities (Figure 9). Before projection we convert the site velocities to local reference frames, which are anchored to the Ordos (A-A′), South China (B-B′), South China (C-C′), NE China (D-D′), North China (E-E′), and South China (F-F′) blocks, respectively. The A-A′ profile strikes NNE across the Ordos block, and no obvious offsets are detected across the Hetao rift zone and the Qinling Northern Frontal fault at the northern and southern block boundaries. The B-B′ and C-C′ profiles strike NNE from South China to NE China and cross the western and eastern parts of North China and NE China, respectively. The deformation field along the western profile (B-B′) shows~2 mm/yr sharp sinistral shear across the ZBSZ, and another~2 mm/yr gentle sinistral shear across the rest of the North China plain. The deformation field along the eastern profile (C-C′) is quite smooth, and~4 mm/yr sinistral shear is broadly distributed throughout the eastern part of the North China plain and the southern part of the NE China block. Also along this profile about 1-2 mm/yr

Journal of Geophysical Research: Solid Earth
convergence is found across~400 km distance in the southern part of the North China plain. The data are more scattered in the western part of the North China block than in the eastern part, possibly due to disturbance from hydrological deformation within the North China plain (e.g., Feng et al., 2018).
The D-D′ profile strikes ESE across the NE China block, along which~1 mm/yr dextral shear is found across 200 km distance range spanning the Yilan-Yitong fault. The E-E′ profile strikes ESE from the Alashan block to the North China block, along which the Yinchuan rift zone and the Liupanshan fault experience about 2-3 mm/yr dextral shear. About 4 mm/yr extension is also found spread out over~500 km distance across the Yinchuan rift zone, Liupashan fault, and eastern part of the Alashan block. A broad dextral shear zone is found spanning the Shanxi rift zone and Taihang mountains, but no obvious extension is detected. The Tanlu fault looks quiet, with no noticeable deformation across. The F-F′ profile strikes ESE across the South China block, from the Songpan-Ganzi region to the Taiwan Strait, and the component parallel to the profile shows steady convergence across the Songpan-Ganzi region and the Longmen Shan fault system, with a total amount of~5 mm/yr shortening within~400 km range. This distributed deformation must be due to convergence across multiple strands of faults in the region. No obvious normal deformation is found for the rest of the profile from the Sichuan basin to the Taiwan Strait, indicating a rather stable and coherent interior of the South China block. The component normal to the profile shows widespread dextral shear of up to 8 mm/yr within the Songpan-Ganzi region NW of the Longmen Shan fault, within~300 km distance centered around the Longriba fault. Dextral shear of only~1 mm/yr across the Longmen Shan fault is detected.

GPS Secular Velocity Solution
In this study we have processed GPS data observed for the past two and a half decades in continental China and its vicinity in a rigorous way, to produce a crustal secular velocity solution. Compared with numerous velocity solutions published and widely used in the past (e.g., Gan et al., 2007;Liang et al., 2013;Rui & Stamps, 2019;Wang et al., 2017;Zhang et al., 2004;Zheng et al., 2017), our solution has the following major differences and improvements.

Data collection and spatial coverage
Our solution incorporates GPS data from the CMONOC I and II networks, with a temporal span of 1999 to 2016. The longer time span of data enables more accurate velocity estimates. The solution also includes more than 500 sites from other projects, and velocity solutions of most of these sites are published for the first time, such as~50 sites in southern Tibet,~60 sites in the Yunnan region, and~40 sites in the vicinity of the Shanxi rift zone, etc. (Table 1 and Figure S1). Inclusion of these sites enables better spatial resolution, particularly for the above regions with complex tectonic deformation.

Treatment for earthquake induced displacements
As shown in Figure 2, large earthquakes in China and its vicinity constantly affected the deformation field, and most part of continental China experienced postseismic deformation of at least one event for any given period of time during the past 25 years. How to accommodate such deformation fields becomes a serious issue. Instead of abandoning all the data affected by earthquakes or ignoring the effects (e.g., Liang et al., 2013;Rui & Stamps, 2019;Wang et al., 2017;Zheng et al., 2017), we model the coseismic and postseismic displacements caused by the large earthquakes. Our approach not only helps improve the interseismic velocity solution but also provides the coseismic and postseismic displacement solutions for further studies of these events and the deformation processes.

Solution quality
One way to assess the quality of velocity solution is through evaluation of the uncertainties. However, uncertainties are usually determined based on matching data to a model, which is methodology dependent and without accounting of epistemic errors (e.g., seasonal variation and antenna setup error). For a discrete data set representing a spatially continuum function such as the interseismic velocity field that are not affected by shallow fault creep or other abrupt spatial variation, smoothness of the velocity solution is crucial. We use the bootstrapping approach to measure smoothness, which is similar to that used by Zheng et al. (2017) and Kreemer et al. (2018). In this approach the velocity estimate at each site is obtained by interpolating the velocity data set without using the datum of the site , and then compared with the datum of the site to yield the residual. We conduct the statistical analysis for the residuals of all the sites. For comparison, we perform the same analysis for the velocity solutions of Liang et al. (2013), Wang et al. (2017), Zheng et al. (2017), and Rui and Stamps (2019). Figure S8 plots the velocity residual histograms and their fitting curves to a normal distribution. It shows that the bootstrapping residuals are more concentrated for our solution, with its histogram RMS significantly smaller than that of the other four solutions. Therefore, in the sense of smoothness test of the velocity field our solution is better than the other four previous solutions.

Comparison With Previous Deformation Results
A great deal of research has been done for crustal deformation and tectonic processes in East Asia using GPS data. As aforementioned, these studies include analyses of tectonic block motion and quantifications of fault slip rates and strain rates, etc. It would be a daunting task to compare the deformation result derived from this study with all of previous studies. Instead, we will make comparisons with only a few of the most recent studies, focusing particularly on the differences. We will not compare results with the ones dealing with block motions in Tibet (e.g., Loveless & Meade, 2011;Thatcher, 2007;Wang et al., 2017;Zhang et al., 2013), since we do not invoke a kinematic block motion model to represent the deformation field in Tibet, and such a comparison would not be quite objective.

Strain rate estimates
Strain rates are derived by interpolating the GPS velocity field. The strain rate estimates are useful to identify unknown deformation sources and can be used to infer seismic hazard potentials (e.g., Shen et al., 2007;Wang et al., 2015). Allmendinger

Journal of Geophysical Research: Solid Earth
Stamps (2019) derived strain rates in continental China, to which our strain rate pattern shows overall similarity. However, our result ( Figure 4) in general offers a more robust solution with high strain rates correlated with active faults and less fluctuation due to artifacts at low strain rate zones. It also demonstrates clear dilatation and rotation patterns, particularly the expansion belt located in southern and southeastern Tibet, and the contrasting rotations of clockwise rotation around the EHS and counterclockwise rotation around the Xianshuihe-Xiaojiang fault system, respectively (the latter was also revealed by Kreemer et al., 2014, based on the analysis of an early version of this velocity solution). Significant differences are also spotted in various places. For example, our result does not show large strain rate along the southernmost section of the Yilan-Yitong fault north of the Bohai Bay as seen in Rui and Stamps's (2019) result. Theirs could be the result of GPS velocity outliers associated with anthropogenic/hydrologic activities in this region. Also, high contraction rates of~80 nanostrain/yr were found for the Pamir and Tien Shan orogenies in Zheng et al.'s (2017) and Rui and Stamps's (2019) solutions, respectively, in contrast with our result of~50 nanostrain/yr. Another significant difference is the shear strain rate across the central East Kunlun fault, which is determined to be 20-30 nanostrain/yr in our solution versus 50-60 nanostrain/yr in Zheng et al.'s and Rui and Stamps's results. The anomalous difference is possibly due to their velocity solutions affected by the postseismic deformation of the 2001 Kokoxili earthquake (He et al., 2018;Ryder et al., 2011). Zhang et al. (2018) proposed 6.0 ± 1.3 mm/yr sinistral shear over the North China plain, which is significantly different from our estimate of 3.5 ± 0.3 mm/yr. Much of the difference, we infer, comes from the use of different reference frames. Zhang et al. used a GPS velocity data set referenced to the stable Eurasian plate. As we document in this study, however, eastern China rotates counter-clockwise with respect to the Eurasian plate, at rates of~0.2,~1.5, and~0.8 nanorad/yr for the South China, North China, and NE China blocks, respectively ( Figure 8 and Table 2). A part of the sinistral shear shown in Zhang et al.'s estimate has to come from the relative rigid body rotation. Considering the 1.5 nanorad/yr rotation over~1,000 km length of the North China block, the differential sinistral velocity across the block would be~1.5 mm/yr, which explains most of the 2.5 mm/yr discrepancy between Zhang et al.'s result and ours. Reduction of the shear strain rate in the region means reduced estimates of seismic moment accumulation rate and earthquake potential.

Karakorum fault slip rate
The Karakorum fault is arguably the most tectonically active fault in southwestern Tibet (Figures 5 and S7a), and there has been a long debate concerning its deformation rate. High fault slip rate of~10 mm/yr was estimated geologically by England and Molnar (1997) and geodetically by Banerjee and Bürgmann (2002), while other geological and geodetic studies showed low slip rates of about 0-6 mm/yr along segments of the fault (e.g., Brown et al., 2002;Jade et al., 2004;Kundu et al., 2014;Murphy et al., 2000;Robinson, 2009;Wang & Wright, 2012;Zheng et al., 2017). Recent GPS studies of Kundu et al. (2014) and Zheng et al. (2017) reported a decreasing rate from northwest to southeast, which is contradictory to our result. Our estimated fault slip rates are not significant from zero throughout the entire fault except across the southern central segment at 80°E, with a dextral slip rate of 3 ± 1 mm/yr ( Figure S7b). In the studies of Kundu et al. (2014) and Zheng et al. (2017), the rate estimates might be contaminated by the HMT motion at the far field, especially for the northwestern segment of the fault.

Present-Day Crustal Deformation of Continental China 4.3.1. Crustal deformation pattern
The present-day crustal deformation of continental China has been studied extensively and debated for decades. Numerous studies focused on using GPS velocity data to quantify the deformation field, such as to characterize the deformation by a block motion model (e.g., Meade, 2007;Thatcher, 2007;Wang et al., 2003;Wang et al., 2017). Some recent studies, using updated and more densified GPS velocity data, however, found that block motion modeling has its limitation to describe the deformation field, particularly within the Tibetan plateau (e.g., Wang et al., 2017;Zhang et al., 2013). Using the GPS velocity solution of this study, we carefully examine the deformation patterns inside and outside of the Tibetan plateau. As shown in Figures 5  and 6, the deformation field within the plateau is mostly continuous, which is particularly true for the velocity component along the relative plate motion direction. Except for across the HMT, this velocity component for sites located within the Tibetan plateau and Tien Shan mountains shows rather smoothly distributed shortening without sharp gradients (Figure 6a, similar to the result shown by Zhang et al., 2004). Such a deformation pattern cannot be explained by elastic deformation across a limited number of faults which are locked within the upper crust, since the deformation is broadly distributed, and not concentrated in the vicinity of a few large faults. The velocity component normal to the relative plate motion shows a combination of continuous deformation and offsets across major faults within the plateau, and contributions from the two are about equally significant. If a block motion model is invoked to describe crustal deformation in the plateau, as the spacing of faults and block size decrease, it becomes impossible to extract information about individual faults from the associated deformation field. If the block sizes and numbers follow a fractal distribution (Bird, 2003), the fractal dimension has to be high, and with the current network density within the region, the GPS data, unfortunately, are not adequate to address the issue of fractal distribution. Nevertheless, the highly continuous deformation pattern within the Tibetan plateau seems to lend support to a model with widespread viscous deformation in lower crust which in return manifests the upper crust deformation (e.g., Flesch et al., 2001;Royden et al., 1997).

Two-way extrusion of Tibetan plateau
Eastward extrusion of the Tibetan plateau has been recognized and investigated for the past several decades (e.g., Avouac & Tapponnier, 1993;England & Molnar, 2005;Gan et al., 2007;Molnar & Tapponnier, 1975;Wang et al., 2017;Zhang et al., 2004;Zheng et al., 2017). The westward extrusion of the plateau, although it has been proposed based on geological investigations and geodetic data (e.g., Banerjee & Bürgmann, 2002;Garzione et al., 2003;Tapponnier et al., 1981;Zhang et al., 2004), has not been reliably documented. In this study, we determine that the lateral extrusion patterns separate at~82°E longitude (Figure 5b). From this longitude westward the extrusion starts to gain momentum, which spans a region between the HMT and the MPT, and reaches its peak at~6 mm/yr for the Hindu-Kush region. The eastward extrusion takes the Altyn Tagh and East Kunlun faults as the north boundary, and the sinistral slip of the Altyn Tagh fault also initiates at~82°E longitude and extends eastward. The rate of extrusion increases eastward to form a bell-shaped flow, which reaches a maximum of~20 mm/yr between the Jiali and Ganzi-Yushu faults (Figures 5b and 6b). This belt is also accompanied with 10-20 nanostrain/yr extension along its path, which is consistent with a model of the crustal material flowing out of the plateau under gravitational push (Cook & Royden, 2008;Flesch et al., 2001;Royden et al., 1997). The widespread areal expansion within the plateau also suggests that the plateau may have ceased its growth and entered the phase of elevation reduction, as inferred by Ge et al. (2015).

Southeastern Tibet deformation
Deformation of southeastern Tibet is characterized by the large-scale clockwise rotation around the EHS. The rotating terrane is flanked by the Xianshuihe-Xiaojiang fault system as its northeast, east, and southeast outer boundaries and by the Jiali, Nujiang, Lancangjiang, Dayingjiang, and Longling faults as its inner boundary, respectively (Figures 4c and 7b). These two boundary zones are underlain by two bands of high electric conductivity materials in the lower crust and upper mantle (Bai et al., 2010), implying that the high strain field observed at the surface extends downward through the upper mantle. The rotating terrane is then suggested to be driven by gravitational flow, whose two flanking zones in the lower crust and upper mantle are weakened by differential motion across them (e.g., Gueydan et al., 2014;Niemeijer et al., 2010). Viscosities within the shear zones are reduced with respect to the outside, allowing for the lateral flow. The internal deformation pattern of the rotating terrane reveals that the rotation is accomplished, at least in the upper crust, by a combination of translation and rotation of small blocks within the rotating terrane, whose primary motion directions are parallel to the tangential direction of the large-scale rotation. Internal shear does exist, but mainly along faults parallel or subparallel to the rotation direction. Deformation along the NW trending faults such as the Red River, Chuxiong-Jianshui, and Wuliangshan faults seems to contribute relatively minor to regional deformation (Figure 7b).

Tectonic deformation within blocks
Our result demonstrates that blocks outside of the Tibetan plateau and Tien Shan undergo small but not negligible internal deformation (Figure 8 and Table 2). Mean strain rates of the blocks are usually on the order of a couple of nanostrain/yr. The strain rate patterns agree with regional tectonic loading environments. For example, the principal compressional strain rates within the Tarim and Jungger basins are oriented N-S, reflecting the long-range tectonic loading due to the India-Eurasia collision. The Jungger basin is surrounded by thrust nappe structures, and the relatively high internal strain rate is caused mainly by the compressional deformation near the basin margins. For the North and NE China blocks such strain rates are approximately E-W compression and N-S extension, suggesting that the tectonic origin could be convergence between the Pacific-Philippine Sea and Eurasian plates and eastward extrusion of the Tibetan plateau. If so part of the strains has to be elastic and would be released during and after large subduction-zone earthquakes along the Japan and Ryukyu trenches. Such scenario indeed took place during and after the 2011 great Tohoku-Oki earthquake, whose coseismic and postseismic displacements in the North and NE China regions demonstrated centimeter level eastward motion and E-W extension ( Figure S2e). The remaining part of the E-W shortening must be accommodated by permanent deformation.
The eastern China blocks undergo minor counter-clockwise rotation, with the North China block rotating at a significantly greater rate than the South China and NE China (Amurian) blocks. The internal strain rate of the North China block is also greater than that of the South China and NE China blocks. These results support the inference of Liu and Yang (2005), that sandwiched between the expanding Tibetan plateau and the stable Amurian plate including NE China, the North China block experiences significantly larger differential stresses than the South China block, and the deformation is further enhanced by faulting and thermal thinning of the lithosphere. That also explains why North China experiences more active seismicity than the north and south neighboring blocks.

North China Deformation
The North China plain consists of a cluster of NNE trending active faults such as the Tanlu, Jixian-Tangshan, and Xingtai-Langfang faults (Figure 9a), which slip rightlaterally and cause earthquakes periodically (Zhang et al., 2018). Driving by the E-W compressional boundary forces, the elongated subblocks between these faults rotate counterclockwise, resulting in E-W shortening and N-S extension for the region and dextral slip on these NNE trending faults. As the boundary fault between the North China and NE China blocks the ZBSZ undergoes sinistral shear motion. It, however, also accommodates additional shortening and extension at the northern end of the North China subblocks rotating with respect to the stable Yan Shan terrane, resulting in a broadened and fragmented fault zone. This revelation explains how the strain rate pattern shown in Figure 8 is associated with the regional tectonics in North China. Also as the ZBSZ transects the NNE trending Tanlu and Yilan-Yitong faults in the Bohai Bay at its eastern end, the~2 mm/yr sinistral motion across the fault zone becomes much broader and spread over~500 km for the eastern section. This is in clear contrast with the shear motion within a~100 km range for the western section of the fault zone inland (Figure 9b). This result explains much diffused seismicity in the Bohai Bay area (Teng et al., 2014) and attests different seismic hazard patterns between the eastern and western sections of the ZBSZ.