Normal Faulting in the 2020 Mw 6.2 Yutian Event: Implications for Ongoing E–W Thinning in Northern Tibet

: Extensional earthquakes in the Tibetan Plateau play an important role in the plateau’s orogenic evolution and cause heavy seismic hazard, yet their mechanisms remain poorly known, in particular in harsh northern Tibet. On 25 June 2020, a Mw 6.2 earthquake struck Yutian, Xinjiang, o ﬀ ering us a rare chance to gain insights into its mechanism and implications in the Tibetan extension. We used both descending and ascending Sentinel-1 images to generate coseismic deformation associated with this event, which indicates a typical extensional mechanism with a maximum subsidence displacement of 25 cm and minor uplift. The causative fault constrained with interferometric synthetic aperture radar (InSAR) data based on a ﬁnite fault model suggests that the fault plane has a strike of 186.4 ◦ and westward dip of 64.8 ◦ , and the main rupture is concentrated at a depth of 3.6–10.8 km with a peak slip of 0.85 m. Our source model indicates that the 2020 Yutian event ruptured an unknown high-angle blind normal fault with N–S striking. The total released geodetic moment yields 2.69 × 10 18 N · m, equivalent to Mw 6.23. We used dense interseismic global positioning system (GPS) measurements to reveal an approximate 7 mm / yr extensional motion in the Yutian region, but it still does not seem large enough to support high local seismicity for normal events within 12 years, i.e., Mw 7.1 in 2008, Mw 6.2 in 2012, and this event in 2020. Combined with Coulomb stress change modeling, we speculate that the seismicity in Yutian is related to the lower lithospheric dynamics.


Introduction
The Tibetan Plateau, with an average altitude of 5000 m over a lateral area of thousands of kilometers, exhibits complex large-scale intercontinental deformation and serves as a prototype for intracontinental orogenic process on the Earth (e.g., [1][2][3], Figure 1). The N-S shortening and lithospheric thickening in Tibet is well established by the India-Asia collision, but the mechanism of extensional faulting along N-S striking grabens remains speculative and has been widely disputed [4][5][6][7][8][9][10][11]. Indeed, the cause of extension in Tibet is directly related to the fundamental issue for the plateau rise and deeply affects climate and environment changes [6]. To address the question of what the mechanism of normal faulting is in the Tibetan Plateau, numerous hypotheses have been proposed, e.g., gravitational collapse of high topography [7], thermal evolution of the thickened lithosphere [8], convective removal of Beach balls represent focal mechanism solutions of events with magnitude larger than Mw 6 during 1976-2020 from the Global Centroid-Moment-Tensor (GCMT), light blue represents historical events and red the 2020 event. Yellow points are aftershocks (M ≥ 4) within one week after the main shock from the US Geological Survey (USGS). Blue and black lines indicate regional active faults from [12,13], respectively. Blue rectangles depict footprints of Sentinel-1 A/B images used in this study. The insert map shows the location of the epicenter and the study region in Tibet Plateau. The yellow star is the epicenter, and the dashed blue rectangle is our study region.
Geodetic observations throughout the plateau interior verify a significant ESE-WNW extension at a rate of 23 ± 3 mm/yr, which is accommodated half by strike-slip and half by normal faulting [2,7,14]. Different from the strike-slip concentrated on several major faults (>10 mm/yr; e.g., Altyn Tagh, Karakorum, Kunlun, and Xianshuihe faults), these normal faults with low deformation magnitude are associated with numerous narrow-width rift flank uplifts and basins [3] and are poorly constrained due to sparse global positioning system (GPS) sites and time series interferometric synthetic aperture radar (InSAR) data. Therefore, the active deformation and relative mechanism of normal faulting remain barely known. Present characteristics of normal-faulting activities are mainly delineated by focal mechanisms from earthquake catalogues, suggesting high seismicity with normal faulting in southern Tibet but a low level in northern Tibet divided by the boundary of the Jiali fault (i.e., Bangong-Nujiang suture) [15,16]. In contrast to the south, the Tibet extension has commonly been accepted as being responsible for the coupling of the Indian lower crust [17]; the mechanism of Beach balls represent focal mechanism solutions of events with magnitude larger than Mw 6 during 1976-2020 from the Global Centroid-Moment-Tensor (GCMT), light blue represents historical events and red the 2020 event. Yellow points are aftershocks (M ≥ 4) within one week after the main shock from the US Geological Survey (USGS). Blue and black lines indicate regional active faults from [12,13], respectively. Blue rectangles depict footprints of Sentinel-1 A/B images used in this study. The insert map shows the location of the epicenter and the study region in Tibet Plateau. The yellow star is the epicenter, and the dashed blue rectangle is our study region.
Geodetic observations throughout the plateau interior verify a significant ESE-WNW extension at a rate of 23 ± 3 mm/yr, which is accommodated half by strike-slip and half by normal faulting [2,7,14]. Different from the strike-slip concentrated on several major faults (>10 mm/yr; e.g., Altyn Tagh, Karakorum, Kunlun, and Xianshuihe faults), these normal faults with low deformation magnitude are associated with numerous narrow-width rift flank uplifts and basins [3] and are poorly constrained due to sparse global positioning system (GPS) sites and time series interferometric synthetic aperture radar (InSAR) data. Therefore, the active deformation and relative mechanism of normal faulting remain barely known. Present characteristics of normal-faulting activities are mainly delineated by focal mechanisms from earthquake catalogues, suggesting high seismicity with normal faulting in southern Remote Sens. 2020, 12, 3012 3 of 18 Tibet but a low level in northern Tibet divided by the boundary of the Jiali fault (i.e., Bangong-Nujiang suture) [15,16]. In contrast to the south, the Tibet extension has commonly been accepted as being responsible for the coupling of the Indian lower crust [17]; the mechanism of the northern Tibet extension far from the Indian plate is more enigmatic and impedes our understanding of the overall Tibet extension [16].
On 25 June 2020 at 21:05 UTC (13:05 local time), a Mw 6.3 earthquake (Ms 6.4 according to the China Earthquake Networks Center (CENC]) struck Yutian County, Xinjiang Province, China recorded by the US Geological Survey (USGS) earthquake catalogue. This event occurred in a junction area of the Altyn Tagh, Karakax, and Kunlun fault systems in northwestern Tibet, 164 km south of Yutian County. Within one week, five M 4.0+ aftershocks occurred along a lineament striking approximately N-S. As it is a sparsely populated region, no casualties were reported in this event. The USGS, Global Centroid-Moment-Tensor (GCMT) [18], and CENC agree that the primary focal mechanism for this event was normal faulting. Note that this junction region has suffered high-level seismicity in recent times, consisting of one strike-slip event (Mw 6.9 in 2014) and three normal-faulting events (Mw 7.1 in 2008, Mw 6.2 in 2012, and Mw 6.2 in 2020) reported from the GCMT archive ( Figure 1). As an important window to obtain insights into extension activity in northern Tibet, it has drawn the attention of scholars both at home and abroad [11,12,[19][20][21][22]. However, the harsh terrain, high altitude, and extreme weather conditions make this remote region barely accessible, and as a result only a few field investigations were conducted for the 2008 Yutian event [12]. In addition, InSAR images experience severe decorrelation for the 2008 event [11,20], and are not available for the 2012 event due to SAR satellite updating. That is, the knowledge of these causative faults remains limited in this region. Consequently, the 2020 event gains much more significant value in terms of exploring the evolution of this normal event to incorporate the ongoing northern Tibet extension. Earthquakes with good source modeling would yield much insight into their tectonic structure and dynamic mechanism [23]. In this study, we analyze both descending and ascending C-band Sentinel-1 satellite data to determine the coseismic displacements of the 2020 Yutian event. Based on the InSAR observations, we implement source modelling to determine fault geometry and slip rupture on the causative fault. In addition, we collect previous interseismic GPS measurements to simply analyze the local extensional setting. Finally, we carry out Coulomb stress change modeling to explore the seismic behavior in northern Tibet.

Tectonic Setting
It is well recognized that the present deformation of the Tibetan Plateau is accompanied by not only mega thrust-slip and strike-slip faults, but also ample interior normal faults. As shown in Figure 2, a large number of normal faults associated with north-trending rifts have been deduced from structural relief in the plateau. Note that normal faults in the northern part of the Jiali fault zone exhibit small extensional features and isolated grabens, in contrast to the southern part, resulting in much less knowledge in northern Tibet [24]. The GCMT earthquake catalogue has documented many normal slip events that occurred in Tibet in the past few decades. These numerous normal slip events are of moderate magnitude, resulting in the sparse GPS measurement in Tibet. With European Remote Sensing (ERS-1/2), Envisat, and ALOS SAR images, coseismic and postseismic deformation were derived to constrain the coseismic rupture and lithospheric rheological structure, e.g., the four normal-faulting earthquakes (Mw 5.7-6.2) that occurred in the Pumqu-Xainza rift (PXR) in 1993 and 1996 [25], the 2004 Mw 6.2, 2005 Mw 6.3, and 2008 Mw 6.7 Zhongba events [26], the 2008 Mw 5.9 and 6.4 Gaize events [11,27,28], and the 2008 Mw 6.3 Damxung event [29]. Yao et al. [30] noted that the Mw 6.2 event in Zhongba in 2004 could have been triggered by the remote Mw 9.1 and Mw 8.6 earthquakes in Sumatra in 2004 and Nias in 2005, respectively. Most of these events that occurred in southern Tibet and some near central Tibet have attracted a lot of attention. In contrast, the Yutian region in the northwestern corner of Tibet, where the recent Mw 7.1 earthquake in 2008, the Mw 6.2 earthquake in 2012, and this Mw 6.2 earthquake in 2020 occurred, has been considered as a high-seismicity exception for normal faulting. Furthermore, the 2008 event was mapped by incomplete coverage of InSAR images [11,20], and the 2012 event escaped from available SAR image, resulting from the harsh terrain [31]. To some extent, the 2020 Yutian event is the first major normal-slip event well imaged by InSAR in this region. Therefore, it adds significant value for the investigation of the Tibet extension.
active blocks, such as the western Kunlun, Qaidam-Qilian, Bayan Har, and Qiangtang blocks ( Figure  1) [12]. In the differential movement between the Altyn-Tagh and Longmu-Gozha co-fault systems, E-W-dipping extensional stress is generated. Since 2008, four events (Mw > 6) have been recorded by the GCMT catalogue, involving three normal-slip events on different faults and one strike-slip on the Altyn-Tagh fault [32]. From the 10 m resolution Sentinel-2 optical image ( Figure 2, insert map), we find an appropriate location of this Yutian event lying between the Chung Muztagh mountains in the west and Heishibei Lake in the east. The peak of the highest Chung Muztagh mountain and an unknown lakefront mountain are covered by some ice sheets. Between the two mountains, a relatively flat landscape is characterized by some alluvial fans and small gullies, and the major gully strikes in the N-S direction, which corresponds to a common north-trending rift. Note that a known normal fault in the eastern flank of the Chung Muztagh mountains with east dipping has been inferred from field geological surveys after the Mw 7.1 Yutian event in 2008 [12], but whether it is the causative fault for this 2020 event is still unclear. Figure 2. Seismicity in Tibet and local geomorphic features (after [11]). Beach balls represent focal mechanism solutions of events with magnitude larger than Mw 6 during 1976-2020 from GCMT, in which normal, strike, and thrust events are rendered in green, blue, and black, respectively. The normal Yutian event in 2020 is plotted with red. The insert map shows near field geomorphic features for this event from the Sentinel-2 optical image. Red symbols with black triangles show locations of gullies.

Coseismic Displacement Analyzed by InSAR Data
Our study benefits from complete coverage of SAR imagery over the epicenter region. Two pairs of terrain observation with progressive scans (TOPS) images from the C-band Sentinel-1 radar satellite, operated by the European Space Agency (ESA), were employed to derive the coseismic displacement associated with the Yutian event. The novel constellation design and imaging modes Figure 2. Seismicity in Tibet and local geomorphic features (after [11]). Beach balls represent focal mechanism solutions of events with magnitude larger than Mw 6 during 1976-2020 from GCMT, in which normal, strike, and thrust events are rendered in green, blue, and black, respectively. The normal Yutian event in 2020 is plotted with red. The insert map shows near field geomorphic features for this event from the Sentinel-2 optical image. Red symbols with black triangles show locations of gullies.
Tectonically, northwestern Tibet, where the 2020 Yutian earthquake occurred, is characterized by conjugate mega-strike-slip faults with minor normal faults [12]. The regional tectonics consist of the mega-strike-slip Altyn Tagh fault to the north and en echelon Longmu-Guozha co-faults spreading from the Altyn Tagh fault to the south, which originated from an interaction of several active blocks, such as the western Kunlun, Qaidam-Qilian, Bayan Har, and Qiangtang blocks ( Figure 1) [12]. In the differential movement between the Altyn-Tagh and Longmu-Gozha co-fault systems, E-W-dipping extensional stress is generated. Since 2008, four events (Mw > 6) have been recorded by the GCMT catalogue, involving three normal-slip events on different faults and one strike-slip on the Altyn-Tagh fault [32]. From the 10 m resolution Sentinel-2 optical image ( Figure 2, insert map), we find an appropriate location of this Yutian event lying between the Chung Muztagh mountains in the west and Heishibei Lake in the east. The peak of the highest Chung Muztagh mountain and an unknown lakefront mountain are covered by some ice sheets. Between the two mountains, a relatively flat landscape is characterized by some alluvial fans and small gullies, and the major gully strikes in the N-S direction, which corresponds to a common north-trending rift. Note that a known normal fault in the eastern flank of the Chung Muztagh mountains with east dipping has been inferred from field geological surveys after the Mw 7.1 Yutian event in 2008 [12], but whether it is the causative fault for this 2020 event is still unclear.

Coseismic Displacement Analyzed by InSAR Data
Our study benefits from complete coverage of SAR imagery over the epicenter region. Two pairs of terrain observation with progressive scans (TOPS) images from the C-band Sentinel-1 radar satellite, Remote Sens. 2020, 12, 3012 5 of 18 operated by the European Space Agency (ESA), were employed to derive the coseismic displacement associated with the Yutian event. The novel constellation design and imaging modes make the Sentinel-1 satellite advanced in terms of a shortened revisiting period and large-scale space coverage but with similar pixel resolution to previous generations (i.e., ERS-1/2 and Envisat satellites), which greatly strengthen its interferometric capacity to capture surface deformation signals. The analyzed SAR dataset involves one descending pair (T165D) on 17 June 2020 and 29 June 2020, and one ascending pair (T158A) on 22 June 2020 and 4 July 2020 (detailed InSAR information listed in Table 1). The two image pairs have full coverage for the epicenter region ( Figure 1). All of these InSAR data were processed by the commercial GAMMA software platform following a mature two-pass differential InSAR approach [33]. We adopted a multi-look factor of 20:4 in the range and azimuth direction for the Sentinel-1 images and removed the phase contribution of topography with 1 arc-sec resolution digital elevation data [34]. Then the interferograms were filtered, unwrapped, and geocoded within a common processing flow. The processed interferograms are shown in Figure 3a,b, which show the main deformed region featuring a pattern of two lobes along the N-S direction, and a flat valley surface between the Chung Muztagh mountains to the west and Heishibei Lake to the east. In addition, some decorrelations in the interferograms are caused by ice sheet and lake water cover. Note: A and D indicate ascending and descending orbit, respectively. σ is the standard deviation of line-of-sight (LOS) measurements in non-deforming regions of interferograms and α is the e-folding correlation length scale of 1D covariance functions of interferograms. All data errors for original and corrected interferograms are listed.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 18  To refine the quality of interferograms, we considered the atmospheric error corrected by the Generic Atmospheric Correction Online Service (GACOS) for InSAR [35]. As the major error source of InSAR, the atmospheric effect can reach up to several centimeters on an interferogram [36], and it should be taken into account in applications of small-scale deformation activities, e.g., mining subsidence, landslide, volcanic deformation, tectonic movement, and even coseismic deformation caused by a small or moderate event [35,37]. The atmospheric correction model related to topography has been successfully applied in many cases (e.g., [38,39]), but this empirical model hardly distinguishes error from actual deformation and leads to overestimation [40]. The GACOS data make up a high-resolution tropospheric delay model developed by the European Centre for Medium-Range Weather Forecast, and advanced medium-long-scale (>10 km) atmospheric error correction for InSAR [41]. This study shows relative atmospheric error ranging 2-4 cm for the whole interferogram and 1-2 cm for near field inferred from the GACOS model, and an average 1 cm phase standard deviation reduced after correction (atmospheric phase contribution shown in Figure 3c,d). We estimated the precision of the corrected interferograms with a simple 1D covariance function, which suggests a standard deviation of 0.45 cm and 0.8 cm for T165D and T158A, respectively ( Table 1).
The final interferograms shown in Figure 4 exhibit clear deformation signals in a complex tectonic region, and their displacement ranges from −22 to 10 cm (Figure 4a,b). The interferograms indicate a maximum negative line-of-sight (LOS) displacement value of 19.2 cm and 21.9 cm for T165D and T158A, respectively (negative LOS value represents increasing distance from the satellite), with a dominant vertical component. In comparison to the negative displacement, the positive displacement is very small. Figure 4c,d shows the LOS coseismic displacement profiles corresponding to the dashed line (AA') drawn in Figure 4a,b. The profiles from descending and ascending images show similar feature, major subsidence displacement accompanied by minor uplift. In addition, the profiles in Figure 4 show no surface rupture caused by this event, and as a result it is difficult to directly relate any of the previously known faults to the causative fault.   With the two interferograms in different look vectors, we can implement two-dimensional decomposition to resolve the east-west and upward-downward components by neglecting the north-south component, also called a 2.5-D surface deformation map. The LOS displacement is less sensitive to the north-south displacement, which can be neglected in the simple 2.5-D decomposition method. In addition, the 2020 Yutian event occurred on a normal fault with approximate N-S strike, which consists mostly of a vertical component. Therefore, the 2.5-D decomposition method is suitable for this case. As shown in Figure 5a, there exists maximum eastward displacement with a value of 10.7 cm and westward displacement with a value of 7 cm. In the vertical component from Figure 5b, a clear subsidence displacement with a maximum of 25 cm is observed on the left side, and some minor uplift on the right side. This 2.5-D surface deformation field suggests that the 2020 Yutian event is predominantly normal slip.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18 With the two interferograms in different look vectors, we can implement two-dimensional decomposition to resolve the east-west and upward-downward components by neglecting the north-south component, also called a 2.5-D surface deformation map. The LOS displacement is less sensitive to the north-south displacement, which can be neglected in the simple 2.5-D decomposition method. In addition, the 2020 Yutian event occurred on a normal fault with approximate N-S strike, which consists mostly of a vertical component. Therefore, the 2.5-D decomposition method is suitable for this case. As shown in Figure 5a, there exists maximum eastward displacement with a value of 10.7 cm and westward displacement with a value of 7 cm. In the vertical component from Figure 5b, a clear subsidence displacement with a maximum of 25 cm is observed on the left side, and some minor uplift on the right side. This 2.5-D surface deformation field suggests that the 2020 Yutian event is predominantly normal slip.

Source Modeling
Given that coseismic displacement is considered to be responsible for the subsurface rupture slip within a homogeneous elastic half-space [42], most recent earthquake mechanisms have been described with a finite rectangular dislocation model by geodetic observation. The analytical solution of a rectangular dislocation model involves fault geometry (longitude and latitude of fault location, length, width, depth, strike, dip) and relative kinematic parameters (slip magnitude and rake direction). As a complex geophysical inversion problem, a mature two-step inversion strategy including nonlinear uniform and linear distributed slip models was adopted to solve these parameters in previous works (e.g., [31]), and it is also followed in this study.
Before implementing the source model, we utilized a resolution-based downsampling method for the interferograms to generate the input dataset for inversion. Each Sentinel-1 interferogram yielded 10 7 -10 8 measurements in line of sight from the satellite, providing dense coseismic displacement observations. Indeed, it would cause heavy burden on computational efficiency and disturbance on unknown parameters if all of these points were used for inversion calculation. Therefore, it was necessary to resample the interferograms with the proper number points and retain the main deformation feature. Gradient-based quadtree sampling is the most widely used method for InSAR downsampling [43], but it is difficult to retain a suitable distribution and number of points when there is limited deformation magnitude caused by a small or moderate event. Conversely, a resolution-based resampling method developed by [44] can well distinguish the information in near and far fields for interferograms, and additional information of rough fault location is needed.

Source Modeling
Given that coseismic displacement is considered to be responsible for the subsurface rupture slip within a homogeneous elastic half-space [42], most recent earthquake mechanisms have been described with a finite rectangular dislocation model by geodetic observation. The analytical solution of a rectangular dislocation model involves fault geometry (longitude and latitude of fault location, length, width, depth, strike, dip) and relative kinematic parameters (slip magnitude and rake direction). As a complex geophysical inversion problem, a mature two-step inversion strategy including nonlinear uniform and linear distributed slip models was adopted to solve these parameters in previous works (e.g., [31]), and it is also followed in this study.
Before implementing the source model, we utilized a resolution-based downsampling method for the interferograms to generate the input dataset for inversion. Each Sentinel-1 interferogram yielded 10 7 -10 8 measurements in line of sight from the satellite, providing dense coseismic displacement observations. Indeed, it would cause heavy burden on computational efficiency and disturbance on unknown parameters if all of these points were used for inversion calculation. Therefore, it was necessary to resample the interferograms with the proper number points and retain the main deformation feature. Gradient-based quadtree sampling is the most widely used method for InSAR downsampling [43], but it is difficult to retain a suitable distribution and number of points when there is limited deformation magnitude caused by a small or moderate event. Conversely, a resolution-based resampling method developed by [44] can well distinguish the information in near and far fields for interferograms, and additional information of rough fault location is needed. Considering the magnitude of the 2020 Yutian earthquake, we carried out the latter method to resample the interferograms, and retrieved 680 and 713 points for T165D and T158A, respectively (retained points are shown in Figure 6).
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 18 Considering the magnitude of the 2020 Yutian earthquake, we carried out the latter method to resample the interferograms, and retrieved 680 and 713 points for T165D and T158A, respectively (retained points are shown in Figure 6). In the uniform slip model inversion, we inverted the fault geometry by assuming the subsurface rupture with a uniform slip. Differently from previous direct-search methods that only consider the observed data error, a Bayesian inverse approach developed by [45] determines optimal model parameters with the posterior probability density functions (PDFs) of source parameters as well as the prior knowledge of observations. With this Bayesian inversion method, the ranges of the initial model are not restrictive (e.g., the strike and dip angles can vary from 0° to 360° and -90° to 90°, respectively), implying that all possible fault geometry and kinematics parameters can be explored (detailed initial and optimal model solutions can be found in Table 2). After 10 6 iterations, we removed a burn-in period of 1 × 10 4 iterations and retained the results with the converged Markov chain. The retained candidate solutions are illustrated in Table 3 and Figure 7. Besides the parameter of width having a strong correlation with some other parameters, e.g., depth and X center, all other parameters demonstrated good statistical distribution, suggesting a tight constraint for model solution with the uniform slip model. In the uniform slip model inversion, we inverted the fault geometry by assuming the subsurface rupture with a uniform slip. Differently from previous direct-search methods that only consider the observed data error, a Bayesian inverse approach developed by [45] determines optimal model parameters with the posterior probability density functions (PDFs) of source parameters as well as the prior knowledge of observations. With this Bayesian inversion method, the ranges of the initial model are not restrictive (e.g., the strike and dip angles can vary from 0 • to 360 • and −90 • to 90 • , respectively), implying that all possible fault geometry and kinematics parameters can be explored (detailed initial and optimal model solutions can be found in Table 2). After 10 6 iterations, we removed a burn-in period of 1 × 10 4 iterations and retained the results with the converged Markov chain. The retained candidate solutions are illustrated in Table 3 and Figure 7. Besides the parameter of width having a strong correlation with some other parameters, e.g., depth and X center, all other parameters demonstrated good statistical distribution, suggesting a tight constraint for model solution with the uniform slip model.    The optimal uniform slip model (Table 3) demonstrates that the 2020 Yutian event ruptured a fault plane with a size of 11.9 × 7.37 km 2 , striking 186.4 • , and dipping at a high angle of 64.8 • . The main slip consists of a normal slip of 0.78 m and left lateral slip of 0.13 m, which is concentrated between a depth of 3.3 and 10.0 km and does not reach the surface. Note that the nodal plane solutions from USGS and GCMT along N-S striking are east-dipping with angles of 42 • and 50 • , respectively. These double solutions from USGS and GMCT with seismic wave give opposite fault dip trends in comparison to our source model, suggesting an inconsistency between seismological and geodetic results. Assuming a shear modulus of 33 GPa [19], the modelled geodetic moment is 2.3 × 10 18 N·m, less than the GCMT estimation of 3.29 × 10 18 N·m.
As illustrated in Figure 8, the synthetic and relative residual displacements demonstrate that the InSAR datasets fit our uniform model well. It is noted that the standard deviation of residuals of the geodetic data for T165D and T158A is 1.1 cm and 1.4 cm for the model, respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 18 between a depth of 3.3 and 10.0 km and does not reach the surface. Note that the nodal plane solutions from USGS and GCMT along N-S striking are east-dipping with angles of 42° and 50°, respectively. These double solutions from USGS and GMCT with seismic wave give opposite fault dip trends in comparison to our source model, suggesting an inconsistency between seismological and geodetic results. Assuming a shear modulus of 33 GPa [19], the modelled geodetic moment is 2.3 × 10 18 N·m, less than the GCMT estimation of 3.29 × 10 18 N·m. As illustrated in Figure 8, the synthetic and relative residual displacements demonstrate that the InSAR datasets fit our uniform model well. It is noted that the standard deviation of residuals of the geodetic data for T165D and T158A is 1.1 cm and 1.4 cm for the model, respectively.  Table 3. The black line indicates the fault line projected on the surface along its dip plane upward, and the black rectangular frame is the uniform slip plane in the subsurface.
In the next step, we estimated the variable slip distribution on a planar fault by fixing the fault geometry shown in Table 3. To avoid the fault edge effects, we extended its length and width to 24 km × 20 km, and then divided it into 1 km × 1 km patches (for a total of 480 patches). We considered the slip rake angle varying -130° to -60°. After that, we estimated the variable slip distribution with the steepest descent method (SDM) [46]. A normalized smoothing factor of 0.12 was selected as a trade-off of the slip roughness and root mean square (RMS) misfit. As shown in Figure 9, the preferred slip distribution displays a simple oval pattern with one asperity, and its maximum slip is 0.85 m. The main slip (>0.4 m), with a size of ~11 × 7 km 2 , is confined to a depth of 3.6-11.4 km. The rupture slip is dominated by normal-slip, with a little strike-slip. The distributed slip model yields a geodetic moment of 2.69 × 10 18 N·m, equivalent to Mw 6.23, larger than that inferred from the uniform slip model.  Table 3. The black line indicates the fault line projected on the surface along its dip plane upward, and the black rectangular frame is the uniform slip plane in the subsurface.
In the next step, we estimated the variable slip distribution on a planar fault by fixing the fault geometry shown in Table 3. To avoid the fault edge effects, we extended its length and width to 24 km × 20 km, and then divided it into 1 km × 1 km patches (for a total of 480 patches). We considered the slip rake angle varying −130 • to −60 • . After that, we estimated the variable slip distribution with the steepest descent method (SDM) [46]. A normalized smoothing factor of 0.12 was selected as a trade-off of the slip roughness and root mean square (RMS) misfit. As shown in Figure 9, the preferred slip distribution displays a simple oval pattern with one asperity, and its maximum slip is 0.85 m. The main slip (>0.4 m), with a size of~11 × 7 km 2 , is confined to a depth of 3.6-11.4 km. The rupture slip is dominated by normal-slip, with a little strike-slip. The distributed slip model yields a geodetic moment of 2.69 × 10 18 N·m, equivalent to Mw 6.23, larger than that inferred from the uniform slip model. With the slip distribution model, the fitting residual errors for T165D and T158A are within a standard deviation of 1 cm (Figure 10). The synthetic and relative residual displacement are improved in comparison to the uniform slip model, even though there is some residual displacement near the fault.  With the slip distribution model, the fitting residual errors for T165D and T158A are within a standard deviation of 1 cm (Figure 10). The synthetic and relative residual displacement are improved in comparison to the uniform slip model, even though there is some residual displacement near the fault. With the slip distribution model, the fitting residual errors for T165D and T158A are within a standard deviation of 1 cm (Figure 10). The synthetic and relative residual displacement are improved in comparison to the uniform slip model, even though there is some residual displacement near the fault.

Local Extension Determined by Interseismic GPS
Since the 1990s, GPS measurements have been implemented in China and its neighboring regions, providing crucial quantitative constraints for crustal motion research at various scales [47][48][49].
Given that previous GPS measurements were derived from different regional GPS networks and reference epochs by different institutions, Wang et al. [49] reprocessed the raw GPS data included in all previous studies and their own unpublished sites by using a rigorous data processing algorithm, to produce a secular velocity solution with high spatial resolution for continental China. In this study, we collected interseismic GPS observations from [49], as shown in Figure 11. The new GPS velocity solutions unmask the Tibetan Plateau featuring a western extrusion rate of 6 mm/yr and eastward extrusion rate of 20 mm/yr, divided by the pivot boundary at about 82 E • [49], close to the 2020 Yutian event.
In order to reveal both strong extensional motion and strike-slip in the Yutian area, we selected GPS sites in two profile planes, A and B, throughout the area, and simply replotted their E-W and N-S components, as shown in Figure 11. Profile A is mainly associated with the E-W extension and shows a relative motion of~7 mm/yr in the E-W direction and a few in the N-S direction. Profile B reveals relative motion of 2-3 mm/yr in the E-W direction and~8 mm/yr in the N-S direction. The N-S motion of profile B is mainly associated with strike-slip motion attributed to the left-lateral strike-slip Altyn-Tagh fault. We conclude that the Yutian region had a concentration of both strong strike-slip and extensional motions.

Local Extension Determined by Interseismic GPS
Since the 1990s, GPS measurements have been implemented in China and its neighboring regions, providing crucial quantitative constraints for crustal motion research at various scales [47][48][49]. Given that previous GPS measurements were derived from different regional GPS networks and reference epochs by different institutions, Wang et al. [49] reprocessed the raw GPS data included in all previous studies and their own unpublished sites by using a rigorous data processing algorithm, to produce a secular velocity solution with high spatial resolution for continental China. In this study, we collected interseismic GPS observations from [49], as shown in Figure 11. The new GPS velocity solutions unmask the Tibetan Plateau featuring a western extrusion rate of 6 mm/yr and eastward extrusion rate of 20 mm/yr, divided by the pivot boundary at about 82 E° [49], close to the 2020 Yutian event.
In order to reveal both strong extensional motion and strike-slip in the Yutian area, we selected GPS sites in two profile planes, A and B, throughout the area, and simply replotted their E-W and N-S components, as shown in Figure 11. Profile A is mainly associated with the E-W extension and shows a relative motion of ~7 mm/yr in the E-W direction and a few in the N-S direction. Profile B reveals relative motion of 2-3 mm/yr in the E-W direction and ~8 mm/yr in the N-S direction. The N-S motion of profile B is mainly associated with strike-slip motion attributed to the left-lateral strike-slip Altyn-Tagh fault. We conclude that the Yutian region had a concentration of both strong strike-slip and extensional motions. Figure 11. GPS velocities in Tibet (after [49]). Blue arrows represent horizontal GPS-derived interseismic velocities with respect to stable Eurasia. Green and pink rectangle zones show GPS profiles A and B, which were selected to analyze strike-slip and extension throughout the northwestern corner of Tibet. The insert map shows E-W and N-S motion components for profiles A and B. In profile A, green triangular points and gray square points indicate E-W and N-S motion components. In profile B, pink triangular points and gray square points indicate N-S and E-W motion components. Figure 11. GPS velocities in Tibet (after [49]). Blue arrows represent horizontal GPS-derived interseismic velocities with respect to stable Eurasia. Green and pink rectangle zones show GPS profiles A and B, which were selected to analyze strike-slip and extension throughout the northwestern corner of Tibet. The insert map shows E-W and N-S motion components for profiles A and B. In profile A, green triangular points and gray square points indicate E-W and N-S motion components. In profile B, pink triangular points and gray square points indicate N-S and E-W motion components.

Coulomb Stress Change Modeling
A slip that occurs on faults within the coseismic process changes the stress field in the surrounding crust, thereby promoting or inhibiting local seismicity. The Coulomb stress failure criterion has been widely employed to explore earthquake interactions since the 1990s, and has become a crucial quantitative index to interpret stress field changes and seismic hazards [50]. The general hypothesis of Remote Sens. 2020, 12, 3012

of 18
Coulomb stress triggering is that a positive value promotes failure and a negative value inhibits failure. A low Coulomb stress change of 0.1 bar was considered enough to trigger a future earthquake in a previous study [51]. As Yutian is a strong seismic region, the Mw 6.2 event in 2020 was the fourth Mw > 6 event in the region since 2008. To explore the relationship between the 2020 earthquake and the other three events, we used the Coulomb v3.4 package developed by the USGS for Coulomb stress change modeling in this study [52]. Note that there were no finite fault planes for the 2012 and 2014 events constrained by geodetic data; previous studies only employed the fault parameters derived from seismologic data for Coulomb stress change analysis [52,53]. In addition, the 2008 event was covered with incomplete InSAR data, resulting in fault parameters from previous studies having large differences [11,20,22]. Therefore, we used the fault parameters for the 2008, 2012, and 2014 events referring to [53,54] (detailed fault parameters can be found in Table 4). To some extent, these fault parameters are roughly limited to present observations. Taking the causative fault of the 2020 mainshock as a receiver fault, we calculated the static Coulomb stress change from the previous three events and all four Mw > 6.0 events. Figure 12a shows a positive Coulomb stress change region along the WNW direction between the 2008 and 2014 events. The epicenter of the Mw 6.2 earthquake in 2020 falls in the increased Coulomb stress region, suggesting that previous events occurred to promote this event. Although this event partially released the stress in the Yutian region (Figure 12b), the positive Coulomb stress change region still exists.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 18 A slip that occurs on faults within the coseismic process changes the stress field in the surrounding crust, thereby promoting or inhibiting local seismicity. The Coulomb stress failure criterion has been widely employed to explore earthquake interactions since the 1990s, and has become a crucial quantitative index to interpret stress field changes and seismic hazards [50]. The general hypothesis of Coulomb stress triggering is that a positive value promotes failure and a negative value inhibits failure. A low Coulomb stress change of 0.1 bar was considered enough to trigger a future earthquake in a previous study [51]. As Yutian is a strong seismic region, the Mw 6.2 event in 2020 was the fourth Mw > 6 event in the region since 2008. To explore the relationship between the 2020 earthquake and the other three events, we used the Coulomb v3.4 package developed by the USGS for Coulomb stress change modeling in this study [52]. Note that there were no finite fault planes for the 2012 and 2014 events constrained by geodetic data; previous studies only employed the fault parameters derived from seismologic data for Coulomb stress change analysis [52,53]. In addition, the 2008 event was covered with incomplete InSAR data, resulting in fault parameters from previous studies having large differences [11,20,22]. Therefore, we used the fault parameters for the 2008, 2012, and 2014 events referring to [53,54] (detailed fault parameters can be found in Table 4). To some extent, these fault parameters are roughly limited to present observations. This study Taking the causative fault of the 2020 mainshock as a receiver fault, we calculated the static Coulomb stress change from the previous three events and all four Mw > 6.0 events. Figure 12a shows a positive Coulomb stress change region along the WNW direction between the 2008 and 2014 events. The epicenter of the Mw 6.2 earthquake in 2020 falls in the increased Coulomb stress region, suggesting that previous events occurred to promote this event. Although this event partially released the stress in the Yutian region (Figure 12b), the positive Coulomb stress change region still exists.

Discussion
Geodetic measurements from multi-geometric-view InSAR interferograms enable us to conclude that the mechanism of the 2020 Yutian event is mainly normal-slip rupture. Different from GPS measurements, InSAR only provides displacement observations in one dimension in the LOS direction, and as a result it is hardly direct in judging the source mechanism of an earthquake by the coseismic displacement itself. For the first-and second-generation SAR satellites (e.g., ERS-1/2 and Envisat), most of the time only one geometric-view interferogram is available, such as only descending interferograms for the four normal-faulting earthquakes in 1993 and 1996 [39]. Benefiting from the design of systematic global monitoring, the third-generation Sentinel-1 SAR satellite has large-scale coverage and a short revisit interval and its ascending and descending acquisition have been greatly improved, which means we can observe coseismic deformations with different geometric views with high quality [55]. In addition, multi-geometric-view observation provides more important information than 1D regarding fault behavior and models [56]. Integrating the descending T165D and ascending T158A interferograms from Sentinel-1 SAR data shows that the 2020 Yutian earthquake was dominated by vertical deformation, involving a maximum subsidence of 25 cm, but minor uplift. According to elastic dislocation theory, the 2020 Yutian event suggests a typical extensional mechanism similar to the Mw 6.5 Norcia earthquake in Italy in 2016 [57].
Our source model indicates that the 2020 Yutian event ruptured an unknown high-angle blind normal fault with N-S striking. Our uniform slip model with a global Bayesian inverse approach indicates that this event dominated by normal slip ruptured a fault plane with a steep dip angle of 64.8 • striking 186.4 • , which is not consistent with either nodal plane from USGS or GCMT. Besides this 2020 event, there are two other instrumentally documented normal events in the Yutian region, the Mw 7.1 and Mw 6.2 events in 2008 and 2016. There are many debates about the fault geometry of the 2008 Yutian event, but most of them agree on a high dip angle between 47 • and 62 • (e.g., [11,22]). The Mw 6.2 event in 2012 was out of geodetic constraints and there is no published investigation, but its source mechanism from GCMT shows a nodal plane striking in the N-S direction with a dip angle of 55 • . These normal faults in northern Tibet with high angles are different from normal faults with low angles in the Yadong-Gulu and Lunggar rifts in southern Tibet. Kapp et al. [15] proposed that the Tibetan rifts initiated as high-angle normal and evolved at low angle in response to increasing extension and the rift basin being uplifted and eroded. Therefore, the normal fault in Yutian should be an early stage of rift formation in Tibet. In addition, the distributed slip model in Figure 9 shows the rupture slip is concentrated at a depth of 3.6-11.4 km, with no surface rupture. The main rupture depth of the 2020 event is similar to the depth of 0-15 km derived from the Mw 7.1 event in 2008 [11,22], but its released moment was not enough to reach the surface, as the 2008 event was. From the geological investigation of regional active faults by [12], it can be seen that the location of the causative fault for the 2020 event was an unknown fault between Chung Muztagh mountains and Heishibei Lake. In comparison with the known normal fault with east dipping on the east flank of Chung Muztagh mountains, the length of the causative fault close to Heishibei Lake is much shorter but west-dipping. Given the relative position and dipping trend of these two normal faults, they outline a complete graben structure in this region.
The strong extensional setting in Yutian is supported by high seismicity with normal events. Previous time series InSAR applied in the Yutian region exhibited only strong lateral-slip motion on the Altyn-Tagh and Longmu-Gozha-co faults, but few extensional components [58,59]. However, the GCMT earthquake catalogue verifies that the Yutian region has strong extensional motion by frequent normal events. We speculate that previous time series InSAR results, limited by their observation accuracy, cannot distinguish present extensional motion. Present dense interseismic GPS measurements reveal both western and eastward extrusion in Tibet, and their pivot boundary is just within the Yutian region [50]. We determined that those GPS data exhibit a total extensional motion of~7 mm/yr in the Yutian region. If all these extensions imposed one of these normal slip events, i.e., the Mw 7.1 event in 2008 (with normal slip >5 m; e.g., [11,19,22]), the Mw 6.2 event in 2012 (with normal slip of 0.49 m; refer to [53]), or the Mw 6.2 event in 2020 (with normal slip of 0.8 m), the shortest recurrence period of characteristic earthquakes is >70 years.
Note that the theoretical earthquake cycle (>70 years) is much longer than the 12 years for these normal events. Two possible explanations have been proposed to address this problem. One explanation is that there is seismic clustering behavior. If so, most of the Yutian region will lead to a long-lasting quiet spell after these events. With the Coulomb stress change modeling, the 2020 event falls into a region of positive stress change, suggesting that it would have been triggered by the previous events. However, Bie et al. [53] inferred that the 2012 earthquake could not be explained in terms of Coulomb stress change due to the 2008 Yutian earthquake, as its epicenter fell into the negative stress change zone. Assuming that these Coulomb stress change values are reliable, the present normal events in the Yutian region do not yield a clustering behavior. An alternative explanation is that there should be lower lithospheric dynamics to simultaneous affect the nucleation of these normal events. From a previous geological study [1], the large strike-slip faults in Tibet die out westward in the Yutian region, suggesting an asthenosphere convection process. Seismic imaging profiles across Tibet demonstrates that the Yutian region is located in the lithosphere-asthenosphere boundary between the Indian and Asian tectonic plates [60]. In addition, present GPS measurements exhibit widespread areal expansion within both southern and northern Tibet, suggesting that the plateau may have ceased growing and entered the phase of reduced elevation [49]. Therefore, we are more inclined to accept the latter explanation that normal faulting in northern Tibet is related to the lower lithospheric dynamics to some extent. Recently, another Mw 6.3 normal-slip earthquake (July 2020) struck Nima County, northern Tibet (Figure 2). Ongoing E-W thinning is also occurring in northern Tibet, which should be received more attention in future work.

Conclusions
Surface deformation associated with the Mw 6.2 earthquake in Yutian in 2020 is revealed by both descending and ascending Sentinel-1 interferograms, with LOS displacement ranging from −22 to 10 cm. We used the finite fault dislocation model in an elastic half-space to estimate the causative fault parameters. Our preferred model indicates that this event occurred on a N-S striking plane dipping to the east at a high angle, and the main rupture dominated by normal slip was concentrated at a depth of 3.6-11.4 km with a dip angle of 64.8 • . Combined with previous regional faults, this unknown blind normal causative fault reveals graben tectonics between the Chung Muztagh mountains and Heishibei Lake. With present interseismic GPS observations, a~7 mm/yr E-W extension is inferred in the Yutian region. We speculate that northern Tibet is presently experiencing ongoing E-W thinning, and its normal faulting may be partly contributed by deep lithospheric evolution.