Earth and Planetary Science Letters ‘Two go together’: Near-simultaneous moment release of two asperities during the 2016 M w 6.6 Muji, China earthquake

On 25 November 2016, a M w 6.6 earthquake ruptured the Muji fault in western Xinjiang, China. We investigate the earthquake rupture independently using geodetic observations from Interferometric Synthetic Aperture Radar (InSAR) and regional seismic recordings. To constrain the fault geometry and slip distribution, we test different combinations of fault dip and slip direction to reproduce InSAR observations. Both InSAR observations and optimal distributed slip model suggest buried rupture of two asperities separated by a gap of greater than 5 km. Additional seismic gaps exist at the end of both asperities that failed in the 2016 earthquake. To reveal the dynamic history of asperity failure, we inverted regional seismic waveforms for multiple centroid moment tensors and construct a moment rate function. The results show a small centroid time gap of 2.6 s between the two sub-events. Considering the > 5 km gap between the two asperities and short time interval, we propose that the two asperities failed near-simultaneously, rather than in a cascading rupture propagation style. The second sub-event locates ∼ 39 km to the east of the epicenter and the centroid time is at 10.7 s. It leads to an estimate of average velocity of 3.7 km/s as an upper bound, consistent with upper crust shear wave velocity in this region. We interpret that the rupture front is propagating at sub-shear wave velocities, but that the second sub-event has a reduced or asymmetric rupture time, leading to the apparent near-simultaneous moment release of the two asperities.


Introduction
Earthquake hazard assessment requires an understanding of potential fault rupture length, an important factor controlling the earthquake size. How far ruptures can propagate during earthquakes is limited by intrinsic physical properties of the fault, such as geometrical complexities, stress variations and frictional conditions. Surface rupture observations and numerical models suggest that fault step-overs of ∼5 km are sufficient to stop earthquake rupture propagation (e.g., Wesnousky, 2006Wesnousky, , 2008Biasi and Wesnousky, 2016). Typically, the rupture of asperities on the same fault is assumed to follow a 'cascading model' in which adjacent fault patches sequentially slip seismically as the rupture propagates in its main direction (e.g., Yu et al., 2010). Recent studies highlight a rupture style involving multiple faults during a large magnitude event in the same local tectonic system, via dynamic and/or static * Corresponding author.
E-mail addresses: l .bie @liv.ac .uk (L. Bie), s .hicks @soton .ac .uk (S. Hicks). stress transfer (Hamling et al., 2017;Hicks and Rietbrock, 2015;Hill et al., 2015;Nissen et al., 2016). Nissen et al. (2016) found a time gap of only 19 s between the initial M w 7.0 rupture, and near-instantaneous triggered M w 6.8 rupture (referred to as an aftershock) on different faults over distances greater than 50 km. In this scenario, a physical gap between fault segments/asperities may not completely stop ruptures.
It is common that rupture models of large earthquakes (magnitude greater than 7) show multiple patches of slip along a single fault segment, or along adjacent fault segments. Examples of such ruptures include the 2008 M w 7.9 Wenchuan earthquake (e.g., Fielding et al., 2013), 2010 M w 8.8 Maule earthquake (e.g., Lin et al., 2013) and the 2011 M w 7.1 Van earthquake (Elliott et al., 2013;Zahradník and Sokos, 2014). However, rupture behavior involving near-simultaneous failure of multiple asperities is less commonly reported for earthquakes with moderate magnitudes (i.e. M w <7).
A recent example is the 2015 M w 6.5 Lefkada earthquake in western Greece. By solving for multiple point sources for this rupture, Sokos et al. (2016)  as two patches of slip in the kinematic slip model derived independently from inversion of geodetic observations by Bie et al. (2017). This type of rupture effectively results in a higher seismic moment release over a smaller amount of time. The simultaneous rupture model has important implications for seismic hazard models, earthquake early warning systems, and earthquake forecasting. A simultaneous rupture on two asperities would also be difficult to discern during routine earthquake monitoring due to the similar faulting mechanisms, similar size and limited gap in space and time between two sub-events. Whilst basic information about rupture source complexity (especially timing) can be gained from teleseismic data (e.g., Vallée et al., 2011), more detailed source imaging requires seismic data at regional and local distances (Zahradník and Sokos, 2014).
On 25th November, 2016 (14:24:30 UTC), a M w 6.6 earthquake struck the westernmost part of China, close to the border with Tajikistan and Kyrgyzstan. This event was preceded by a M w 5.2 foreshock that occurred about five minutes earlier (14:18:59 UTC), as reported by the USGS National Earthquake Information Center (NEIC). Two M w 5.1 aftershocks were recorded within one day following the main-shock. Field observations after the earthquake reported limited surface rupture (Chen et al., 2016), suggesting a buried rupture on the Muji fault.
The Muji fault is part of the complex Pamir fault system that accommodates the India-Asia plate collision (Burtman and Molnar, 1993). To the west, the SW-dipping Muji fault likely branches from the Main Pamir Thrust (Robinson et al., 2004), which accounts for as much as 300 km of shortening due to N-S Pamir-Tianshan collision (e.g., Strecker et al., 1995;Arrowsmith and Strecker, 1999). The eastern end of the Muji fault bends to SSE and connects with the extensional Kongur Shan normal fault system (Fig. 1). Based on the analysis of fluvial terrace offsets, Chevalier et al. (2011) found a minimum dextral slip rate of 4.5 ± 0.2 mm/yr for the Muji fault.
As shown by instrumental seismicity catalogues, the regional area of the Pamir mountain range is abundant in moderate to large earthquakes (Fig. 1). Most crustal seismicity occurs in the upper 15 km of the crust (Schurr et al., 2014). The most recent earthquake was the M w 6.4 Sary-Tash, Kyrgyzstan earthquake, which occurred near the southern boundary of the Alai Valley, and probably ruptured the southern basin-bounding thrust fault. The largest event in the last century to have occurred within 100 km of the Muji fault is the 1974 M w ∼7.3 Markansu Valley earthquake (Storchak et al., 2013), ∼20 km to the NW of the Muji earthquake.
Various focal mechanisms were proposed for the 1974 event. Ni (1978) suggested a thrust faulting mechanism, while Jackson et al. (1979) argued for rupture on a right-lateral strike-slip fault, and an aftershock sequence with more than one type of faulting. This discrepancy likely indicates a structural complexity of regional tectonics as mentioned by Jackson et al. (1979), where thrust and strike slip faults jointly accommodate strain due to the India-Asia plate collision. In the Muji basin, according to the ISC-GEM catalogue (Storchak et al., 2013), an M w ∼5.6 earthquake occurred in 1957, ∼20 km SE of the 2017 earthquake.
The 2016 earthquake was the first large instrumentally recorded earthquake on the Muji fault. Aftershocks following the Muji earthquake occurred both east and west of the main-shock epicenter, and mainly to the south of Muji fault trace (Fig. 1). This pattern suggests that the Muji fault dips to the south. In this study, to investigate the spatial distribution of slip during the Muji earthquake, we use surface displacement measurements by interferometric synthetic aperture radar (InSAR) from multiple satellites (Sentinel-1A/B and ALOS-2). Furthermore, based on the segmented static rupture behavior revealed by InSAR observation and inversion, source inversions based on regional seismic waveform data are explored independently to reveal the dynamic history of asperity failure process during the earthquake. We find that the 2016 Muji earthquake ruptured two asperities >5 km apart, almost simultaneously. We combine our geodetic and seismic source inversions to infer what might lead to the near simultaneous failure of two asperities during earthquakes.

InSAR data and processing
The epicentral area of the Muji earthquake was well covered by SAR observations (Fig. 2). We obtained two tracks of Sentinel-1 data in TOPS (Terrain Observation by Progressive Scans) mode with different viewing geometries from European Space Agency. Furthermore, a pair of ascending SAR acquisitions in strip-map mode from ALOS-2 was ordered from Japan Aerospace Exploration Agency (JAXA). Both Sentinel-1 and ALOS-2 data were processed using the open-source GMTSAR software developed and maintained at the Scripps Institution of Oceanography (Sandwell et al., 2011). To remove the topographic contribution to the interferograms, we simulated synthetic fringes using the three arc-second (∼90 m) Shuttle Radar Topography Model (SRTM) digital elevation model (Farr et al., 2007) and subtracted it subsequently from original interferograms. We use two Sentinel interferograms in addition to the ascending ALOS-2 interferogram for subsequent inversions. To reduce the computational effort, we down-sampled each interferogram using a quadtree algorithm (Jónsson et al., 2002). Before down-sampling the interferogram, unwrapping errors to the north of Muji fault trace were manually removed. For the purpose of computational efficiency, we only invert for slip using observations in the area between 38.7 • N-39.8 • N and 73.4 • E-74.8 • E. Excluding further areas where no data is available, only 4808 points are left. The satellite azimuthal and incidence angles are averaged to the same resolution of LOS displacement.

Kinematic slip modeling
In the kinematic inversion, we modeled line-of-sight (LOS) displacements jointly from the Sentinel-1 and ALOS-2 observations as induced by slip on rectangular faults buried in a uniform elastic half-space (Okada, 1985). The shear modulus is 33 GPa based on an upper crust P-wave velocity of 6.0 km/s, S-wave velocity of 3.52 km/s and density of 2.72 g/cm 3 as described in Crust1.0 in this area (Laske et al., 2013), which is similar to a regional 1-D velocity model for the Pamir region (Sippl et al., 2013a). The geomorphological lineament along the Muji mountain front, in addition to the change of sign in the LOS displacements across the fault, especially in the descending data, provide good constraints on the fault strike, which was fixed to 106.4 • in subsequent inversions. Firstly, we attempt to resolve the possible trade-off between fault dip and slip direction (rake). When the optimal dip angle was found, we further derive a distributed slip model with variable rake.
The trade-off between fault dip angle and rake were first investigated by varying dip at 4 • intervals between 64 • and 88 • , and rake ±8 • from −176 • . This generates 42 combinations of fault dip and rake. We constructed the fault plane for each dip angle before subdividing it into patches with dimension of 2 km by 2 km. Green's functions were computed assuming uniform slip on every patch. The non-negative least squares method was then applied to derive distributed slip at the fixed rake and parameters accounting for planar residual orbital errors. The root-mean-square (RMS) misfit in each run was computed and taken as indicator of how well model prediction fits the observation. After getting the minimum RMS, we further refined the searching grid at 2 • intervals between 68 • and 76 • for the dip, and between −178 • and −174 • for the rake. Finally, we took the dip angle that generates minimum RMS from the last step as optimal, and jointly inverted for slip in dextral and reverse directions. Direction and magnitude of slip on each patch was then solved through vector calculation.
To evaluate the robustness of our preferred distributed slip model, we performed checkerboard tests and also conducted error estimation. In the checkerboard tests, synthetic slip distributions with various patch sizes were used as input to forward-model LOS displacements at the InSAR data locations, and the synthetic displacements were subsequently inverted for distributed slip using the same regularization as we applied in the inversion of real data.
We further conducted statistical error estimation to assess the potential effect brought by sources such as the atmospheric noise. A variance-covariance matrix was constructed from the residual phase for each interferogram (Funning et al., 2007), and used to generate a weighting matrix and 100 sets of synthetic spatiallycorrelated noise (e.g., Cervelli et al., 2001;Bie et al., 2014). The synthetic noise was added to the data and weighted least square algorithm was then applied to obtain distributed slip for the 100 sets of perturbed data sets. Finally, one-sigma errors were derived for each patch of the distributed slip.

Multi-source moment tensor inversion
Whilst InSAR data provides excellent spatial coverage of ground deformation, we also analyzed the Muji earthquake using seismic waveforms in order to constrain the temporal-spatial distribution of slip. This method is free of any a priori constraints on rupture velocity or rupture direction.
We obtained seismic waveforms from regional broadband seismic stations located within 450 km epicentral distance of the Muji earthquake ( Fig. 1). Most stations come from the Kyrgyz Digital Network. The mean signal-to-noise ratio over all stations in the frequency range 0.016-0.08 Hz is ∼1800 allowing a detailed waveform modeling approach. Unfortunately, we had to exclude the station KSH of the China National Seismic Network, located ∼180 km to the east of the epicenter, from the inversion since waveforms were clipped. We also excluded waveforms that had long period disturbances related to tilt (Zahradník and Plesinger, 2005) or non-linear behavior, which can occur at seismic broadband stations located close to a large rupture, and if present, can bias the inversion. From these checks, we focus on records from five stations (each with three components) being used for the moment tensor inversions (Fig. 1). The waveforms were corrected for instrument response and integrated to displacement.
We used iterative deconvolution (ID) to invert for deviatoric moment tensors . In this approach, each asperity is approximated by a point-source, that allows us to directly compare with the InSAR co-seismic slip model. We used the 1-D P-wave and S-wave velocity model for the Pamir region of Sippl et al. (2013a). This is a suitable model because it is based on arrivals at stations located just to the west of the epicenter, which is where most stations used for the moment tensor inversion are located. All inversions were carried out with a high-pass frequency corner of 0.016 Hz (43 s); the low-pass corner was chosen based on whether we were inverting for a single source or a multiple-source rupture. We assessed the effect of data errors and poorly constrained Green's functions by jackknifing stations and individual components from the inversions. For multiple-point source inversions, we prescribe a fixed length moment-rate function (triangle) for all sources, which was found to produce the best overall variance reduction (VR).
We also used the non-negative least squares (NNLS) method of ISOLA (Zahradník and Sokos, 2014) to 1) semi-independently verify the result from ID, and 2) to calculate a resulting source-time function. For NNLS, we prescribe the double-couple source mechanism given by the result from ID. At each given trial point source position, we compute the moment rate from a prescribed set of equidistantly shifted isosceles triangles with a given duration. We then invert for the NNLS weights of each triangle.

Description of the displacement
The ground displacement caused by the Muji earthquake is well mapped by the Sentinel-1 and ALOS-2 satellites. Figs. 2a, d and g show Sentinel-1 ascending, descending and ALOS-2 ascending interferograms, respectively. From the two ascending interferograms, which are sensitive to vertical movement in this case, we can clearly see two areas with similar displacement patterns (marked by dashed black boxes in Fig. 2a). The displacement size and amplitude for the two areas are different. The western one is smaller in peak LOS displacement, but larger in area, compared to the eastern area. This similarity in displacement pattern but disagreement in amplitude and area affected indicates a consistent mechanism but differences in slip depth and magnitude. Slip depth of the western asperity must be deeper to produce a longer wavelength of displacement. This interpretation can also be inferred from the Sentinel-1 descending interferogram (Fig. 2d), which mainly maps horizontal displacement here. In comparison to the eastern side, the area of significant displacement on the western side is further away from the fault trace (Fig. 2d).

Inversion results
From the first step of our inversion procedure, we found that a model with dip angle of 70 • and rake of −176 • produces the lowest RMS misfit of 1.522 cm (Fig. S2). The dip angle is within several degrees of the NEIC reported focal mechanism. Right-lateral slip dominates the slip model, consistent with previous field studies on the kinematics of the Muji fault (e.g., Chevalier et al., 2011). Distributed slip with variable rake was obtained by weighted least square inversion. As expected, two major asperities were recovered (Fig. S3). The distributed slip model supports our inference from analysis of interferograms regarding the kinematic slip pattern. The western asperity is longer, wider and deeper than the eastern one. For the western asperity, the maximum slip is ∼0.9 m at a depth of ∼8.6 km, while for the eastern asperity, the maximum slip is ∼1.31 m and shallower at 4.76 km of depth. Taking the 40 cm contour line of slip as the boundary of each asperity, a ∼6 km long gap lies in between, consistent with distributed slip shown by Feng et al. (2017). The horizontal distance between maximum slip patches of the two asperities is ∼25 km. Both asperities are dominated by dextral slip. The total moment of the two major slip patches is 8.27 × 10 18 Nm, equivalent to a moment magnitude of 6.54, consistent with the seismic moment magnitude from NEIC and GCMT. Most of the rupture is buried, with only a limited number of patches with slip on them reaching the surface. This finding agrees with preliminary results from field reconnaissance of surface rupture after the earthquake (Chen et al., 2016). Our inversion result shows deep slip below 20 km depth along the fault, dominated by normal slip. The checkerboard test shows that slip down to a depth of ∼10 km is well recovered (Fig. S4). Considering further the relatively large error below a down-dip depth of 20 km (Fig. S3b), the deeper slip patch with normal slip may not be well-constrained. However, given its location below the major co-seismic ruptures, it could partially relate to after-slip, but this is not within the scope of this study.
The second column of Fig. 2 shows the predicted LOS displacements from the distributed slip model with variable rake, and the third column shows the residual. The predictions fit the data very well. The residual crossing fault trace was probably caused by our assumption of a planar fault near the surface. Considering the variation in strike along the fault trace may further reduce the residual, but success was not guaranteed as shown by Wang et al. (2017). The residuals near the eastern end of fault were likely caused by post-seismic processes or their combinations, which requires further investigation.

Regional centroid moment tensor inversion results
As a first test of the usefulness of complementary seismic data, we found the best fitting single-source CMT solution of the Muji earthquake using regional broadband waveforms, filtered to a lowpass of 0.04 Hz. As a first run, we used the main-shock origin epicenter given by the USGS NEIC (39.273 • N, 73.978 • E) to search for the best-fitting depth position using a grid search. This location is similar to those reported by other agencies and uses phase arrivals from the regional stations used in this study. For the single-source inversion the best-fitting source depth was found at 14 km with a moment tensor solution that agrees well with the published results (e.g., USGS-NEIC, GEOSCOPE), and the VR for this inversion is reasonably high (0.81).
Next, keeping the depth fixed and searching along a line of trial-point-sources along the strike orientation of the WNW-ESE nodal plane (108 • ; determined from the previous test) we searched for the best centroid position in terms of latitude and longitude, keeping the focal mechanism fixed to the result found from the first run. The best centroid position was found to lie 12 km east-south-east of the USGS hypocenter, also at a bestfitting depth of 14 km. This centroid position agrees very well with the location of the greatest slip during the earthquake in the InSAR-derived co-seismic slip model as this patch produces the greatest moment release at low frequencies. This relocation of the centroid improves the fit to the waveforms (VR increases by 11% to 0.9 with respect to the first run using the USGS epicenter). The condition number of CN = 2.8 indicates that the Green's functions matrix is dominated by non-singular elements; therefore, the moment tensor inversion is reliable in terms of the source-station configuration, frequency range and crustal model used (Sokos and Zahradník, 2013). Fig. S5 shows the waveform fit from single source model.
To investigate source complexity, we increased the high frequency corner of the inversion. We found that an elementary triangle source-time function length of 18 seconds produces a peak in VR of 0.74 (Fig. S6). The lower VR value is due to the increased frequency, which now enables us to study the source rupture process in greater detail. With the evidence for source complexity obtained from InSAR inversion, we carried out a two-point source inversion using the ID method. We continued using the trial point source geometry to replicate the geometry of the fault (strike 108 • , dip of 78 • , along-strike spacing of 5 km, and along-dip spacing of 4 km) and inverted waveforms in the frequency range 0.016 Hz to 0.08 Hz. The upper frequency limit was chosen to mitigate any effect of un-modeled 2-D and 3-D structural heterogeneity on our inversions. Using this frequency range ensures that we are negating the effect of heterogeneity on the scale of <80 km, assuming a average crustal velocity of 6.2 km/s for the region (Sippl et al.,  . Heterogeneity on such scale lengths has been imaged in the region by Sippl et al. (2013b). In order to stabilize the inversion from potential biases introduced by the imperfect azimuthal coverage, we kept the double couple source mechanism fixed to the single-source low frequency result, as there is no evidence from the InSAR observations that the fault geometry changes drastically along strike. In the best fitting model, the first sub-event has a moment magnitude of 6.4 and has a centroid time of 7.9 s after origin. The second sub-event has a smaller magnitude (M w 6.1) and occurs later with a centroid time of 10.5 s after origin (2.6 s after Sub-event 1). The overall VR given by this two-point source model is 0.74 (Table 1). Fig. S6 shows the centroid time difference as a function of triangle length for the ID method. For solutions that are within 5% of the maximum VR (at 18 s), we find that time differences may range from 1.6 to 2.8 s. Furthermore, if we run a jack-knife test by removing one station at a time from the inversion, and take those solutions with VR within 5% of the best VR, we have time differences ranging from 2.2 to 2.6 s. We assessed the statistical significance of including the second sub-event via the F-test by following the approach of Hicks and Rietbrock (2015), Sokos et al. (2016) and Dreger and Woods (2002). Based on the number of waveforms used in the inversion (five stations, each with three channels), and assuming that samples are correlated over the period of the 0.08 Hz low-pass filter for the dominant part of the seismogram (∼180 s), this yields 260 degrees of freedom. The VR ratio for the two sub-events is 1.17, which suggests that model improvement by using two sub-events is statistically significant to the 90% level. This confidence level is acceptable given the weaker, and therefore lower amplitude signal from Sub-event 2, along with the InSAR evidence of two discrete asperities.
The relative separation between the two asperities given by the InSAR inversion and seismic waveform inversion is similar.
However there is a difference of ∼5-8 km in the location of the second centroid from the two methods. Therefore, we investigated whether the precise centroid locations as determined from the InSAR slip model could better constrain the seismic waveform inversion and subsequent grid search. We carried out the same waveform inversion as per our preferred solution. Instead of grid searching for the best centroid positions and times, we instead fixed the location of each sub-event to that determined by InSAR (source positions 53 and 45; see Table 1). We used the best fitting time shift at each position. The resulting VR is 0.62 for Sub-event 1 alone, and 0.70 for Sub-events 1 and 2, which a 5% reduction in VR compared to our best-fitting model. The overall time difference at the sub-events at the fixed locations is 2.9 s compared with 2.6 s with the best model. Compared to our best-fitting model, the model improvement using the locations from InSAR is not significant at the 90% confidence level. Overall, this finding suggests that our inferred relative timing between sub-events is a robust feature of the Muji rupture. As no detailed seismic velocity model for the Muji fault area is available, the absolute positions of the two sub-event locations might be biased at scales of 5-8 km, assuming that the two methods "see" the same locations in moment release maxima.
By inputting the source mechanisms and centroid times into a NNLS inversion, we calculated the moment-rate function using the same frequency range as the multiple-source ID result. Each shifted triangle has a duration of 18 s (as given by the result from ID). The NNLS inversion result is stable, with a VR similar to that of the ID inversion. The resulting source-time function (Fig. 3d) is composed of the sum of the two sub-events and confirms the relative strength and timing of the two sub-events given by the ID result. Given the wide elementary triangles, there is no discernable separation between the two events. This matter is further discussed in the next section.

Discussion and conclusions
A key component in earthquake hazard assessment is estimating the maximum possible rupture length, and therefore, the maximum earthquake magnitude (e.g., Field et al., 2009). If the Muji fault failed all the way from the epicenter of the 2016 event to the fault's eastern end where the fault strike changes drastically (Fig. 1), a M w 7.0 earthquake could be expected based on scaling relationships (Wells and Coppersmith, 1994). It is noted that, this estimation is apparently a lower limit, without taking the segment to the west of the epicenter into consideration. In the database compiled by Mohadjer et al. (2016) of central Asian quaternary faults, the total length of the Muji fault reaches ∼133 km, capable in resulting in a M w 7.3 earthquake.
Partial failure of an active fault is not uncommon. Several factors could act as barriers impeding the rupture propagation, thus limiting earthquake magnitude. Such factors include variable fault geometry along strike, variation in fault frictional properties and locking/creeping status, rupture directivity, and past rupture history. We note that, the Muji fault west of the 2016 earthquake changes strike gradually from east-south-east to north-north-east. However, the largest fault bend is approximately 26 km to the west of the 2016 earthquake epicenter. It is therefore likely that the changing fault strike did not play a key role in impeding the western propagation of rupture. More likely, the overall rupture geometry can be attributed to the unilateral rupture propagating to the east.
Overall, Fig. 4 shows that the rupture did not reach the eastern termination of the Muji fault, where the fault strike changes drastically connecting the northern end of Kongur Shan fault. At least 16 km of the Muji fault was not broken co-seismically, leaving a potential future asperity. Similar to the western part of Muji fault, some other factor(s) rather than changing fault strike can be invoked to explain why the rupture stopped in the east. A likely explanation for this is that the eastern unbroken segment exhibits a velocity-strengthening frictional behavior. Thus, aseismic slip there is favored. An alternative explanation is that the Muji fault composes of multiple segments with variable seismic cycle periods. The eastern unbroken part has not reached the critical failure stress, even with stress perturbations from the 2016 earthquake. A similar explanation may also be applicable to the >5 km gap (as given by the 0.4 m slip contour) between the two asper-ities (Fig. 4). We note that the local earthquake catalogue from 2008-2011 of Schurr et al. (2014) shows a small, discrete cluster of events in this area. The correlation of this cluster with the slip gap indicates that this part of the fault is probably creeping. However, this catalogue likely has poorly constrained locations in the Muji area owing to the lack of seismic station coverage there. Monitoring of the post-seismic deformation following the Muji fault, especially the near field of unbroken segments, combined with a better understanding of the past rupture history along the fault, will help to reveal whether these segments are creeping aseismically, or remain locked. For this purpose, we processed several postseismic interferograms spanning first three months after the mainshock (Fig. S7). Although profiles show changes in magnitude of ground displacement across the fault trace (Fig. S7b), and the displacement increases with time, caution should be taken when interpreting how much of the deformation is actually related to aseismic slip on the fault. As pointed out by Feng et al. (2017), the probable fault-related postseismic signal can be easily masked by topographic error, atmospheric noise and seasonal changes in glaciation. Detailed noise analysis of the continued observation of the post-seismic deformation may help clarify the locking status of the un-ruptured segments. This in turn will have a significant impact on the seismic hazard assessment of the Muji fault, and other similar faults.
In terms of the detailed evolution of seismic slip during the 2016 Muji earthquake, the rupture appears complex, with two clear slip patches from geodetic data, and confirmed as two separate sub-events by seismological data, which ruptured almost simultaneously in time. The combined source time function from both sub-events (Fig. 3d) shows two major moment releases with peak-to-peak gap (centroid time difference) of 2.7 s. A synthetic test inverting the synthetic waveforms as real data, yields a similar result but with a slightly shorter time gap (2.2 s) between the two sub-events, indicating that the limited station coverage seems to bias the event timing, probably due to directivity effects. Therefore, we believe that a peak-to-peak gap ranging from 2.2 s to 3.2 s is a robust observation. Furthermore, in our centroid moment tensor and moment-rate inversion, we are constrained by having to provide the same duration of elementary isosceles-triangle functions for all sub-events. We expect that the grid-search for the best triangle length is most sensitive to the sub-event with the greatest moment release (Sub-event 1; M w 6.4). Therefore, and because the integral of the moment-rate function scales with moment magnitude, it is reasonable to assume the triangle half-duration for Subevent 2 is shorter than eight seconds. This second sub-event could therefore be 'masked' within the source of the first sub-event and therefore not be discernible at teleseismic distances. The resulting source-time function is similar to the one provided by Geoscope using the SCARDEC method (Vallée et al., 2011). It implies that the pulse at ∼1 s in the SCARDEC source time function is related to nucleation of the rupture at the hypocenter. With these factors in mind, we also carried out a test for the NNLS inversion using a triangle half-duration of four seconds. The resulting momentrate function has a VR that is within 2% of that for the original NNLS inversion, suggesting that a 4 s half-duration can also fit the data. This alternative moment-rate function is also shown in Fig. S8, and in this case, shows a clearer separation between the two sub-events. Another explanation is that the overall sourcetime function of the rupture is asymmetric, with a short rise time and a gradual fall-off, especially if the rise time is similar, or less than the elementary triangle function length prescribed for the NNLS inversion (1 s). Such a regularized Yoffe source-time function has been proposed for the dynamic interpretation of kinematic slip models (e.g. Tinti et al., 2005), and cannot be ruled out for the Muji earthquake. We also paid attention to the waveform fits at two stations located parallel to the strike, and west-north-west, of the strike of the fault (DRK and BTK). These stations should see the greatest time delay in waveforms between the two sub-events since they are located in the opposite direction to the sequence of asperity failure. A shorter moment-rate function of 12 s results in a 5% higher variance reduction at station BTK compared with our best-fitting model (Table S1). This is one complementary piece of evidence that might suggest a shorter rupture time of Sub-event 2, compared with our best-fitting moment rate length of 16 s at all stations. However, this effect cannot be fully resolved due to the unknown effect of directivity, resulting from the non-ideal station coverage.
Given the short time difference (2.2-3.2 s) and distance between the sub-events (25 km between the centroid locations), we explore the possible rupture mechanisms. Firstly, we consider the rupture behavior in terms of dynamic triggering (although some contribution of static stress changes cannot be completely ruled out). We consider that shear waves are the main controlling factor behind the dynamic triggering since for a strike slip (right-lateral) rupture the S-wave amplitudes in direction of the fault are at a maximum and the P-wave contribution should be minimal. Considering the large distance (∼25 km) between the maximum slip patch of each of the asperities (from the co-seismic slip model), and assuming a rupture velocity of 3-3.5 km/s, we can not explain the observed time difference of 2.2-3.2 s (Fig. 3d), as ∼8 s gap should exist between two moment release impulses, if the failure of the eastern asperity only starts after the completion of western rupture.
An alternative explanation is that the Muji earthquake involved super-shear rupture. Bouchon and Karabulut (2008) found that aftershocks following super-shear rupture tend to appear off the main fault trace, a characteristic pattern indicating failure of secondary structures. Aftershocks of the Muji earthquake seem to follow this characteristic pattern, although this may be due to the fault dipping to SSW and location uncertainties. Considering a rupture velocity of 6.0 km/s, equivalent to P-wave velocity, it will take at least 4 s for the rupture to propagate from centroid location of Sub-event 1 to that of Sub-event 2. This time gap is larger than the time difference of ∼2.5 s we found (also maybe at the upper end); we therefore conclude that even a rupture velocity approaching super-shear speeds is also not a likely explanation for the small gap between sub-events, and we therefore rule out super-shear rupture scenario.
As pointed out above, our synthetic test indicates that Subevent 2 might have a shorter source duration as Sub-event 1.
Sub-event 2 lies ∼39 km from the epicenter; an average rupture velocity of ∼3.7 km/s is expected, given the observed time delay of 10.5 s. Such rupture speed would be consistent with a normal sub-shear rupture propagation based on upper crust S-wave velocity in the region reported by Sippl et al. (2013a). We note that, the rupture velocity estimated here is an upper bound, due to fact that the centroid position of Sub-event 2 from multiple moment tensor inversion is slightly to the east of the eastern slip patch derived by InSAR observation. Asymmetric source time functions, which cannot be ruled out as discussed above, might also influence the inferred rupture velocity.
To reconcile this apparent contradiction between the observed time centroid time delay of 2.2-3.2 s, indicating super-shear rupture, and sub-shear rupture as derived from the time delay between the origin time and the centroid time of Sub-event 2, we carried out a simple synthetic test assuming a temporally more compact triangle moment-rate function of 12 s duration for Subevent 2, with a centroid time delay of 2.5 s. The subsequent inversion using an 18 s triangle leads to an almost identical moment rate function (Fig. S9) as per our preferred inversion result. This suggests the inversion is not able to sufficiently resolve the true time function of the second sub-event. Our explanation for the rupture behavior from our observations of the Muji earthquake therefore is that the rupture front is propagating in the sub-shear regime but that the Sub-event 2 has a temporally more compact moment rate function, leading to the apparent near-simultaneous moment release of the two asperities and significant overlap of the moment rate functions due to the fixed triangle length constraint. Furthermore, we cannot rule out differences between the momentrate rise and fall time for each Sub-event.
Our study of the Muji earthquake is unique in terms of the observation that multiple asperities fail near-simultaneously. In comparison to earthquakes with similar magnitude that involve failure of multiple asperities, e.g. the 2015 M w 6.5 Lefkada earthquake in Greece (Sokos et al., 2016), the Muji earthquake shows shorter temporal gap but slightly larger spatial gap between the sub-events. To our knowledge, the Muji earthquake has smallest temporal gap between two sub-events being reported for the intermediate sized earthquakes. The short centroid time gap between the two sub-events does not favor a typically assumed cascading rupture or super-shear rupture of similar asperities. Our observation suggests that even for an intermediate sized earthquake, the heterogeneity of fault properties might have a significant impact on the moment release rate and therefore the associated possible near fault damage. Although similar observations may not have been documented before for other intermediate sized earthquakes, given the resolving capability for teleseismic inversions to image closely spaced sub-events in both time and space (e.g., Hicks and Rietbrock, 2015), it may be a feature of other ruptures. Our work implies that the potential for enlarged rupture areas due to the failure of multiple asperities should be considered when calculating seismic hazard for a particular fault zone.