Three-dimensional coseismic displacements and slip distribution of the 2021 Mw 7.4 maduo earthquake: Synergy of SAR, InSAR and optical images

On 21 May 2021, an earthquake of moment magnitude (Mw) 7.4 occurred near Maduo county in Qinghai Province China. This is the first major earthquake that occurred in the interior of the Bayan Har block in the Tibetan Plateau over the past 70 decades. Focusing on this event, we conducted a study on three-dimensional (3D) coseismic displacement reconstruction and its tectonic implication. We acquired both Synthetic Aperture Radar (SAR) and optical imaging satellite imagery, including SAR images from Sentinel-1 and Advanced Land Observation Satellite-2 (ALOS-2) as well as the multi-spectrum images from the Sentinel-2 (S2) satellite. We applied the Interferometric SAR (InSAR) and pixel-offset tracking (POT) techniques to coseismic SAR data pairs and reconstructed two dimensional displacements. With the constructed displacement fields in multiple viewing directions, we resolved the 3D coseismic displacements (north-south, east-west, and up components) by integrating the regional strain model with variance components estimation (SM-VCE). We recommend using the standard deviation value at each grid cell, which can be calculated in the resampling, as the initial weight. Based on the resolved 3D coseismic displacements, we further estimated the dip angles for the two segments (F4 and F5) in the east of the rupture zone. The associated moment magnitude is about Mw7.4, which corresponds to the released energy of ∼ 1.74 × 10 20 Nm.


Introduction
On 21 May 2021, an earthquake of moment magnitude (Mw) 7.4 occurred near Maduo county in Qinghai Province China (Figure 1). The determined epicenter is located at 34.59°N, 98.34°E at a depth of~17 km according to the China Seismological Network (CSN, https://ceic.ac. cn/). This Mw 7.4 earthquake occurred inside the Bayan Har block (Figure 1), which is one of the most active regions with strong earthquakes (Mw ≥ 7.0) in the Tibetan plateau (Jia et al., 2021). Most of these large earthquake events occurred along the boundaries of the Bayan Har block, e.g., the 1997 Mw7.5 Manyi earthquake (Wang et al., 2007), the 2001 Mw 7.8 Kokoxili earthquake (Vallée et al., 2008), the 2008 Mw7.9 Wenchuan earthquake (Diao et al., 2010), and the OPEN ACCESS EDITED BY Ziyadin Cakir, Istanbul Technical University, Türkiye 2008 Mw7.1 Yutian earthquake (Song et al., 2019). Whereas the Mw 7. 4 Maduo earthquake is the first major earthquake to occur in the interior of the Bayan Har block since the 1947 Dari earthquake .
The overall length of the coseismic rupture zone is~148 km (Yuan et al., 2022). It is located to the southeast of the big bend of the Kunlun fault , which is seismically active and outlines the northern boundary of Bayan Har Block (See Figure 1). In the past year, there have been several geophysical studies (e.g. He et al., 2021;Zhao et al., 2021;Jin and Fialko, 2021;Chen et al., 2022;Tong et al., 2022) and geological studies (e.g., Ren et al., 2022;Yuan et al., 2022) to determine and investigate the characteristics of the seismogenic fault of the Mw7.4 Maduo earthquake. The geological studies were mostly focused on mapping fault geometry and coseismic surface ruptures. For example, Yuan et al. (2022a) surveyed the horizontal and vertical displacements in the field and via Unmanned Aerial Vehicle (UAV) imageries. They suggested that coseismic vertical displacements were generally small compared with the strike-slip displacements, except at a~3 km zone to the west of the Yematan bridge (see Figure 1). According to the complex vertical displacement of surface ruptures, they also found a contrasting sense of fault dipping directions, especially in the west of the epicenter. Such complex geometry conditions are also supported by the study with relocated post-seismic sequences . Another feature worth noting is a~5 km secondary fault rupture strand parallel to the main rupture strand towards the east section of the coseismic rupture . This secondary fault merges towards the northwest on the main earthquake fault at around 99 E°through a blind fault segment (Yuan et al., 2022). According to the Global Positioning System (GPS) derived three-dimensional (3D) coseismic displacements, Wang et al. (2022) revealed a nearly vertical fault dipping (overall~87°) to the north with two main slip areas above 15 km in depth. However, they have omitted the aforementioned secondary branch given the limited spatial resolution of GPS.  reconstructed 3D displacements using imagery acquired by Synthetic Aperture Radar (SAR) satellites and suggested that the optimal solution of dip angle was 90°for both the main strand and secondary strand. Jiang et al. (2021) conducted a joint investigation with displacement resolved from both GPS and interferometric SAR (InSAR) and suggested the spatial diversity of the dip angle as well as an obvious segmental characteristic of the seismogenic fault. The previous studies mainly agree that the Mw 7.4 Maduo earthquake was dominated by the left lateral slip, but there may be some puzzles about other kinematic parameters of the seismogenic fault. Nevertheless, further studies of the Maduo earthquake, a rare large earthquake (Mw > 7) that occurred in the interior of a major tectonic block within the Tibetan plateau, would provide information to better understand the deformation style and earthquake potential of the Tibetan plateau.

FIGURE 1
Tectonic setting of the 2021 Maduo earthquake. The red star denotes the epicenter location of the 2021 Maduo earthquake. Blue arrows in subplot (B) represent the coseismic horizontal GPS displacements . The pink line in (A) represents the surface rupture traces. In subplot (A), the black lines are active faults from Deng et al. (2003) and the colored points represent the focal depth from aftershocks of the Maduo earthquake . The red boxes in subplot (B) are the coverages of ALOS-2 SAR data. The yellow dashed boxes represent the coverages of ascending and descending Sentinel-1 acquisitions. The purple dash boxes are denoted as coverages of Sentinel-2 data. The white box in subplot (C) is the region as shown in subplot (B).

Frontiers in Earth Science
frontiersin.org 02 In this study, we focused on the 3D coseismic displacement (north-south, east-west, and up components) reconstruction and its tectonic implication. The experiment employed both Sentinel-1 A/B (S1) and Advanced Land Observation Satellite-2 (ALOS-2) acquisitions and applied the Synthetic Aperture Radar Interferometry (InSAR) techniques, which are sensitive in estimating 3D earthquake location and magnitude (Lohman and Simons, 2005). Besides that, coseismic signals were retrieved via subpixel correlation from pre-and post-earthquake image pairs that were acquired by SAR and optical satellite sensors, respectively. The Sentinel-2 (S2) multi-band images are acquired in a view geometry differing from the side-looking SAR. It is therefore helpful in improving the constraint in 3D inversion (Zhang et al., 2017;Bacques et al., 2020;Lu et al., 2021). A strain model based on variance component estimation (SM-VCE) (Liu et al., 2019;Hu et al., 2021) was adopted in the reconstruction of 3D displacements. We also discovered the optimal weight scheme of the initial weight in SM-VCE to achieve robust and efficient 3D displacement estimations. Afterward, we further investigated the geometry and kinematic parameters of the seismogenic fault of the 2021 maduo earthquake based on the Okada model (Okada 1985) in an elastic half-space with the resolved 3D displacement field.

Data preparation
We acquired the S1 SAR imagery in both ascending and descending orbits, ALOS-2 SAR imagery in two adjacent descending orbits, and multi-band optical imagery acquired by the S2 satellite.
The acquisition dates of SAR coseismic pairs are listed in Table 1. The C-band S1 SAR imagery was acquired in Terrain Observation by Progressive Scans (TOPS) mode. The width of a standard swath is about 250 km. The L-band ALOS-2 SAR data has a longer wavelength of 23.6 cm, which is beneficial to resist decorrelation and measure displacements with a large gradient. The ALOS-2 ScanSAR mode imagery has a swath width of 350 km. The footprints of SAR and optical imagery are shown in Figure 1.
S2 is a high-resolution multispectral satellite with a standard granule size of 100 × 100 km 2 . Considering that the quality of optical images relies on cloud coverage, we selected coseismic pairs acquired on 12 October 2020 and 17 October 2021 based on the cloud amount criteria (<3%). Ortho-rectified Level-1C (L1C) products of S2 in two adjacent tracks were obtained to fully cover the coseismic rupture zone (Figure 1).
Additionally, GPS coseismic measurements published by Wang et al. (2022) were also used in the rest of the studies for validating the reconstructed 3D displacement field. They processed GPS time series ranging from 18 May 2021 to 24 May 2021 and resolved coseismic displacement at 58 continuous sites, which are located within a~300 km radius of the epicenter. The GPS-derived coseismic displacement only has minor impacts from postseismic deformation that are no more than 3 mm at the nearfield sites . As shown in Figure 1, GPS site QHAJ shows a large horizontal displacement, which may be deteriorated by localized displacements .
The Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) (Rodriguez et al., 2005) was used in the differential interferogram generation. We applied the 90 m SRTM DEM products to ALOS-2 pairs and 30 m SRTM DEM products to S1 pairs. This is because the reported spatial resolution of ALOS-2 ScanSAR data (~60 m resolution; JAXA, 2012) is lower than that for the S1 dataset (~22m; ESA, 2022b).
2.2 Extraction of 1D/2D displacement fields from multi-source remote sensing imagery 2.2.1 Sentinel-1 A/B pairs processing We processed the S1 coseismic pairs in both ascending and descending orbits. Both differential InSAR and pixel offset tracking methods were applied to reconstruct the onedimension (1D) displacement in Line-of-sight (LOS) direction as well as two-dimensional (2D) displacements in range/azimuth directions, respectively. In InSAR processes, the single look complex (SLC) SAR images were multi-looked by a factor of 10 in range and two in azimuth. The pixels with coherence (<0.6) were masked prior to the phase unwrapping. The Minimum Cost Flow algorithm (Werner and Wegm, 2002) was used to obtain a reliable unwrapped phase.
To conduct the pixel-offset tracking (POT) process (Michel et al., 1999), a window of 128 × 128 was adopted to estimate offsets with an oversampling factor of two between two SAR images, and the offsets were estimated from peaks of sub-pixel cross-correlation in range and azimuth direction respectively. This setup can guarantee a precise co-registration of the coseismic pair. The estimated local offsets were then considered as the displacement signals caused by the earthquake. These 2D displacement fields include movements between satellite and ground targets (range direction), as well as those along the satellite flight direction (azimuth direction). Although the POT-derived range movement is in the same geometry as that in InSAR, POT can achieve complete displacements even if the displacement gradient is large (e.g., in the vicinity of the earthquake rupture zone).

ALOS-2 pairs processing
We processed ALOS-2 coseismic pairs in two adjacent descending orbits. The bandwidth of acquired SLCs in ScanSAR mode is 28 MHz with a swath width of 350 km, which is large enough to cover the entire rupture fault with a single standard frame. ALOS-2 SAR images were multi-looked by a factor of eight in range and 30 in azimuth. We applied the standard differential interferogram generation approach (Rosen et al., 2000) to each sub-swath of the ScanSAR data to then mosaic them in the geographic coordinate. Note that we did not apply the POT approach to ALOS-2 pairs. This is because: (1) ALOS-2 SAR data is in the L band which has better capability to measure a large displacement; and (2) the ionospheric artifact can largely affect the azimuth displacement estimation in the POT method (Brcic et al., 2010;Meyer, 2011). We used the same InSAR process scheme as that applied to S1 pairs to resolve the unwrapped phase of ALOS-2 pairs. A second-order polynomial was then applied to unwrapped interferograms to model and mitigate the ionospheric artifact.

Sentinel-2 pairs processing
The L1C product of the S2 satellite is an orthoimage product that is geometrically refined. The Sentinel Application Platform (SNAP) developed by ESA was used to conduct the atmospheric correction of L1C. Afterward, we used MicMac software (Rosu et al., 2015;Galland et al., 2016;Rupnik et al., 2017) to conduct sub-pixel correlation and reconstruct the lateral displacement between the pre-earthquake and post-earthquake images of each individual band. The MicMac software is developed to compute the correlation of image pairs in the spatial domain by using a regularization algorithm (Rosu et al., 2015;Galland et al., 2016). The derived displacement fields for visible bands (B2, B3, and B4 bands) and near-fared band (B8) were then stacked to reduce the overall noise level. The outcome is a two Dimensional (2D) ground deformation field in both East-West (EW) and North-South (NS) directions, given the orthographic projection of S2 products. However, the NS component was not included in the rest of the 3D reconstruction. This is because of the existence of large residuals and orbital errors in the NS measurements.
To sum up, the SAR coseismic products are limited to range displacement and azimuth displacement in the satellite side-looking geometry. Additionally, SAR measurements are largely ambiguous with regard to NS movement (Wright et al., 2004). This is related to the sun-synchronous orbit of SAR satellites, which would cause further difficulties in resolving 3D displacement given these SAR data shared a similar side-looking geometry. Therefore, we included the displacement products derived from S2 images in this study. The S2 coseismic pairs captured the displacement in another viewing direction that can introduce an extra constraint in the estimation of the 3D coseismic field.

Reconstruction of 3D coseismic displacements
Based on the above-generated 1D/2D displacement fields, we can reconstruct the 3D field of the Mw 7.4 Maduo earthquake. To begin the reconstruction, we first sampled the 1D/2D field into the same geographic grid. Each input coseismic displacement field was resampled into a 0.01°× 0.01°grid size using uniform averaging. The associated standard deviation of each individual grid cell during this averaging process was also calculated and used as an initial weight in the rest of the 3D field estimation.
A strain model (SM) based 3D displacement reconstruction approach (Liu et al., 2018;Liu et al., 2019; was applied here. First of all, a design matrix (A) that connects observations in various viewing directions and unknown 3D displacement at each pixel was constructed. The mathematic expression of the inversion model at a single pixel is a simple linear equation that can be written as follows: in which d 3d is the unknown 3D displacement components in the east-west (d e ), north-south (d n ), and up (d u ) direction and the unit is meter; L denotes the coseismic displacement components in different satellite viewing directions that were derived in the last subsection.
From left to right, L s1−ra is the range offset of S1 data derived from offset-tracking; L s1−azi is the azimuth offset of S1 data derived from offset-tracking; L s1−los is the LOS displacement of S1; L 41 a2−los is the LOS displacement of ALOS-2 data in orbital 41. L 42 a2−los is the LOS displacement of ALOS-2 data in orbital 42. L s2−ew is the east-west displacement of S2. Note that capital letters A and D denote ascending and descending orbits, respectively.
Given the SAR/optical remote sensing reconstructed coseismic components are the projection of 3D deformation components into a certain direction, the design matrix A can be generated from the satellite-looking vector of each component. The geometry vectors of each input displacement field were calculated according to corresponding satellite orientation parameters (e.g., heading angle or/and incidence angle). The basic expression of A can be written as , where a and b are calculated from incidence angle and satellite heading angle depending on the satellite viewing directions, as listed here: in which iα k azi is the satellite direction of kth pixel. θ k inc is the incident angle of the kth pixel.
Afterward, we can further build the functional relationship between displacement observations and 3D surface displacements based on the SM model. The detailed steps of integrating a strain Frontiers in Earth Science frontiersin.org 04 model can be found in Liu et al. (2018) and Liu et al. (2019). Because the SM method is constrained by adjacent points, it is necessary to determine an optimal window to properly determine participating pixels (Shen and Liu, 2020;Wang and Wright, 2012). Taking into consideration that a large window may result in smooth but unrealistic 3D estimation, therefore, we determined an optimal value of window size via a trade-off curve. It was generated between misfit with metric of root mean square errors (RMSE) and different window sizes. Additionally, we took the aforementioned standard deviation calculated during the uniform averaging of each individual grid cell as their initial weight. Given that the input coseismic displacements were measured from multiple sources and techniques, it is necessary to determine a realistic weight for different input observations. Hence, a variance component estimation (Liu et al., 2018) was added to the strain model, as was the so-called SM-VCE method. We conducted an iterative process to determine the relative weight among the eight input observations (Eq. 2). The correction value of each iteration was calculated from the variance component estimation. Based on our experience, the number of iterations is limited to 35 to achieve an efficient convergence.
Note that the SM-VCM method is based on the assumption that the displacements are continuous and smooth in the space. The Maduo earthquake was mainly a strike-slip event, and the two sides of the seismogenic fault were moved relatively in different directions. Thus, the assumption is invalid for pixels in the near field of the rupture zone. Therefore, it is necessary to guarantee that the deformation signals of pixels within a certain window are spatially continuous when we used these surrounding pixels to form a strain model. To solve this issue, we set up a polyline barrier based on geologically mapped rupture information (Yuan et al., 2022). For each window, we then excluded the pixel from the other side of the barrier prior to the inversion 3D inversion.

3D surface displacements
Based on the strategy described in Section 2.1, the coseismic displacement fields in the LOS direction via InSAR are demonstrated in Figure 2. The InSAR-derived range displacement

FIGURE 2
Satellite LOS displacement fields derived from ascending (A) and descending (B) of S1 data as well as track 41 (C) and track 42 (D) of ALOS-2 satellite in descending orbit.

Frontiers in Earth Science
frontiersin.org fields show a similar spatial pattern. Their maximum displacements are both within 1 m. This is mainly because S1 and ALOS-2 have similar viewing geometries. However, the C-band S1 interferograms were decorrelated in the zone to the east of the epicenter, especially in the ascending orbit ( Figure 2A). This decorrelation is related to the large displacement, as this zone is where the maximum coseismic deformation occurred . Benefiting from the longer wavelength, the L-band ALOS-2 results have achieved better coherence in the near field. The coseismic displacements derived from the POT method are shown in Figure 3. The overall POT estimated displacement in range direction has similar patterns to that estimated from InSAR. The west part of the coseismic deformation field shows larger displacements than that in the east. As mentioned above, POT can measure offsets even when the displacement gradient is large regardless of data wavelength. Therefore, the coseismic movement near the fault rupture zone has been fully resolved. However, the azimuth displacement derived via POT shows a high noise level for both ascending and descending orbits. The sun-synchronous orbit and side-looking viewing geometry of SAR satellites, thus limits their ability in measuring NS movement. Additionally, the NS movement is relatively minor for the Mw 7.4 Maduo earthquake (Yuan et al., 2022). Both of these factors contribute to the high noise level of POT azimuth products ( Figures 3B,D), which were therefore excluded in the rest of the 3D inversion.
The measured horizontal displacements from the multi-band optical remote sensing imagery are shown in Figure 4. The spatial pattern differs greatly from the range or LOS movement in Figure 2 and Figure 3, as NS displacements did not diminish towards the far field of the rupture zone ( Figure 4A). This is likely because the displacements in the range direction are a combination of NS, EW, and UP components, while Figure 4B demonstrates only the NS component. Additionally, the spatial reference to EW movement in Figure 4 is also different from that of InSAR or POT. The measured coseismic displacement in the NS and EW directions range betweeñ 2 m. This is reasonable for movement in the NS direction. However, it is too large for that in the EW direction. We considered that the EW estimation with the S2 pair was likely impacted by large residual orbital errors.
Finally, taking into consideration that the strike of the Maduo earthquake in the NW-SE direction is approximately perpendicular to the flight direction of the descending Sentinel-1 satellite, the Range (A-C) and azimuth (B-D) displacement fields of Maduo earthquake calculated from the ascending and descending S1 pairs via POT approach.

Frontiers in Earth Science
frontiersin.org azimuth displacement result from the descending S1 pairs shows large noises, so does the S2-derived NS displacement. Therefore, they were excluded in 3D inversion. Afterward, we conducted a joint estimation with the solution described in Section 2.3. Additionally, we tested the different sizes of windows (from~2 km to~14 km) in the SM-VCE inversion. As shown in Figure 5, the L curve was generated based on each tested window size and the associated global misfit. In the end, an optimal window of 2.5 km was determined and applied to solve the 3D displacement field. The reconstructed 3D coseismic displacements are shown in Figure 6. Due to the limited spatial coverage of Sentinel-2 images, the resolved 3D displacements have a smaller spatial size that is nevertheless enough to cover approximately the entire rupture zone. Obviously, the ground surface movement of the 2021 Mw 7.4 Maduo earthquake is dominated by the EW displacements, which are in a range of -2m~2 m. The NS displacement is much smaller (ranging from -0.8m-0.8 m), which also seems not continuous from east to west of the rupture zone. The UP component is in about -0.5m-0.5 m that is mainly concentrated in the near field of the rupture zone. The existence of a vertical movement field implies that this earthquake is also associated with some minor dip-slip movement, especially in the east and west sections of the rupture zone, where the strike direction starts changing.
We plotted four profiles across the coseismic zone ( Figure 6A). As shown in Figure 7, there is a fast decaying on both sides of the Jiangcuo fault, which is the suggested seismogenic fault of the Maduo earthquake . The largest EW displacement is captured at profile DD'. Additionally, also along profile DD', there are two jumps in the UP component, which may be related to the eastern two branches of the rupture fault. The displacement signal for the UP component for these four profiles shows overall different patterns. It is likely because the dip angle varies from east to west as suggested by Yuan et al. (2022) and Wang et al. (2022).

1) Validation with GPS
To validate the quality of reconstructed 3D displacement fields, we computed the difference between GPS-derived displacements  and above-derived 3D displacements. We determined pixels within a radius of 0.02°of each GPS station and averaged them to represent the SM-VCE-derived 3D displacement. We then further computed the RMSE of differences between SM-VCE and GPS measurements in three directions (Table 2). In the same table, we also listed the result from Liu et al. (2021), which reconstructed the 3D field with multisource SAR imagery. The result indicates an improvement in the accuracy of the resolved NS and EW components. It indicates the advantages of adding displacement derived from S2 imagery and therefore highlights the importance of fusion SAR and optical products in 3D displacement reconstruction.  Frontiers in Earth Science frontiersin.org  Frontiers in Earth Science frontiersin.org 08 Yuan et al. (2022) surveyed horizontal displacement at a total of 68 sites in the field. In this study, we compared reconstructed displacements with their measurements from the geological study. We adopted the same strategy used above to select pixels around each measurement site, while we took the maximum value from the selected pixel for the comparison. The result can be visualized in Figure 8. The horizontal displacements derived by the SM-VCE method are consistent with that surveyed in a field within a range of 0.5 m. We found that the SM-VCE result shows an overall underestimation in the west section, while an overestimation in the east section compared with the in-field measurements. This is likely related to the variation in the width of the coseismic rupture zone. The presented study measured overall displacement across the entire rupture zone, while the infield measurement surveys were offset at a fixed location with limited crossing fault distance (e.g., a few or tens of meters). Another impact could be the post-seismic deformation, given that the temporal resolution of remote sensing data was limited by the revisiting time of satellites. Their products therefore may contain some level of post-earthquake displacements. As suggested by He et al. (2021), the west section hosts larger postseismic deformation compared with the other part of the fault.

Optimal weight scheme of SM-VCE
During the study, we also found a strong dependence on reconstructed 3D displacements and initial weight at each input pixel in SM-VCE estimation. We tested a case run using equal weight as the initial weight scheme for each grid pixel. The corresponding result is presented in Figure 9, showing that the quality is obviously compromised with the equal initial weight. Moreover, the processing time for solving 3D displacements was nearly doubled, and it encountered more difficulties in reaching convergence. Therefore, we want to emphasize the importance of determining the proper initial weight of each input pixel when using the SM-VCE method. In Section 2.2, the applied initial weight was calculated from the standard deviation in the uniform averaging, which is recommended here. While any other weighting schemes that can properly represent the quality of the input dataset would also be applicable.

Coseismic slip distribution and structure of seismogenic fault
An elastic dislocation model (Okada 1985) with a Poisson's ratio of 0.25 was used here to discover the kinematic parameters of seismogenic fault as well as the coseismic slip distribution. The steepest decent method was applied to solve the parameter estimation (Wang et al., 2013). The relocated aftershock catalog (ranging from 22 May 2021 to 30 May 2021) and fault trace of remote sensing images were employed to constrain the fault strike, and we, therefore, construct a fault model with five segments (see Figure 6A). The average strike of F1, F2, F3, F4, and F5 are 272°, 285°, 285°, 89°, and 118°, respectively. The top and bottom of the fault depth are set at a value of 0 and 20 km, respectively. The fault plane is discretized into patches with a size of 2 × 2 km. The rake angles of slip vectors are applied to vary from -45°to 45°. A quad-tree algorithm (Simons et al., 2002) was applied to sub-sample 3D displacements. The down-sampled final dataset consists of 7548 points.
The coseismic slip distribution of the Mw 7.4 Maduo earthquake constrained by the resolved 3D displacements is shown in Figure 11. Note that the dip angles of F1, F2, and F3 which are pre-defined according to the result published by Xu et al. (2021) are 75°, 70°, and 80°respectively. While the fault dip was estimated for the rest of the two segments (F4 and F5). We set up a search space for dip angles that ranges from 85°to 115°for F4 and 50°-115°for F5 with a step of 5°. By changing the dip angle from 0°to 180°, the dipping direction changes from North to South. The optimal dip was then determined by minimizing the global misfit. In the end, the resulting dip angle of F4 and F5 are 85°and 80°towards the south direction, respectively

FIGURE 8
Comparison between the horizontal displacements based on SM-VCE method and infield measurement.
Frontiers in Earth Science frontiersin.org (see Figure 10). The correlation between the observed and modeled dataset reached~0.98.
The determined slip distribution is demonstrated in Figure 11. And the synthetic observation, simulation, and residual from the coseismic slip model inversion explains the rationality and reliability of the inversion result ( Figure 12). The inversion result shows that the slip was mainly concentrated at a depth ≥15 km and the largest slip is approximately 6 m, which is located at the eastern section of the fault at a depth of 5 km. The moment magnitude obtained through inversion is about Mw 7.4 and the released energy is~1.742 × 10 20 Nm. The slip distribution also indicates that the slip at the top of the fault reaches~3m, which is consistent with inverted 3D displacements. The maximum slip is in the east of the epicenter, e.g., the slip primarily ranged from 1.5 m to 6 m at a depth of 0 km-15 km in segment F3. Note that the western part of F3 lacks postseismicities. This zone also corresponds to the area, where the continuous coseismic surface rupture was absent (Yuan et al., 2022). We, therefore, suspect the accumulation stress could be partially released during the main shock, making it far more below the fault strength (the maximum shear strain of the fault plane sustained). Wu et al. (2022) discussed the tectonic stress evaluation on the seismogenic fault of the Maduo earthquake with numerical simulation. They suggested that the low post-  Frontiers in Earth Science frontiersin.org earthquake seismicity zone might relate to the low tectonic stress condition prior to the Maduo earthquake.
We also plotted the geodetic strain field for the 2021 Maduo earthquake ( Figure 13). As the associated strain parameters were simultaneously calculated in SM-VCE when resolved from 3D deformation. According to Figure 13, the maximum shear strain is located at a zone around 99°E. This is where that fault extends to the east and orientation changes from NS to NW. Wu et al.  Frontiers in Earth Science frontiersin.org 11 (2022) suggest that the tectonic stress accumulated the most at this orientation-changing zone (around 99°E~99.5°E) at the hypocenter depth. Although a large coseismic strain has been released here compared with the other segments along the fault (Figure 13), the stress condition of the east section of the rupture zone is questionable after the Maduo earthquake. This point is important in order to evaluate the seismic hazard potential in this area. Therefore, further study on post-seismic activities is needed in the future.

Conclusion
In this study, we focused on the 3D coseismic displacement reconstruction and its tectonic implication for Maduo Mw 7.4 earthquake. With derived displacement in range/offset directions from ALOS-2 and S1 pairs as well as NS direction from S2 pairs, we resolved the 3D coseismic displacement by using the SM-VCE method. We discovered the optimal weight scheme of the initial weight in SM-VCE to achieve robust and efficient 3D displacement estimations. We recommend using the standard deviation calculated in uniform averaging as the initial weight or other weighting schemes that can properly represent the quality of the input dataset. And adding displacement derived from S2 imagery can improve the accuracy of estimated NS and EW components. Based on the resolved 3D coseismic deformation fields, we further investigated the dip angles for the two segments (F4 and F5) to the east of the epicenter. The optimal dip angles were determined towards the south with values of 85°and 80°r espectively. The moment magnitude obtained through inversion was about Mw7.4 corresponding to released energy of 1.742 × 10 20 Nm.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.