Stress drop assessment of the August 8, 2017, Jiuzhaigou earthquake sequence and its tectonic implications *

By using a broadband Lg attenuation model developed for the Tibetan Plateau, we isolate source terms by removing attenuation and site effects from the observed Lg-wave displacement spectra of the M 7.0 earthquake that occurred on August 8, 2017, in Jiuzhaigou, China, and its aftershock sequence. Thus, the source parameters, including the scalar seismic moment, corner frequency and stress drop, of these events can be further estimated. The estimated stress drops vary from 47.1 kPa to 7149.6 kPa, with a median value of 59.4 kPa and most values falling between 50 kPa and 75 kPa. The estimated stress drops show significant spatial variations. Lower stress drops were mainly found close to the mainshock and on the seismogenic fault plane with large coseismic slip. In contrast, the highest stress drop was 7.1 MPa for the mainshock, and relatively large stress drops were also found for aftershocks away from the major seismogenic fault and at depths deeper than the zone with large coseismic slip. By using a statistical method, we found self-similarity among some of the aftershocks with a nearly constant stress drop. In contrast, the stress drop increased with the seismic moment for other aftershocks. The amount of stress released during earthquakes is a fundamental characteristic of the earthquake rupture process. As such, the stress drop represents a key parameter for improving our understanding of earthquake source physics.


Introduction
On August 8, 2017, a devastating earthquake occurred in the Jiuzhaigou scenic region in Sichuan Province, China, causing 25 fatalities, 525 injuries and 6 missing people (http://www.xinhuanet.com/politics/szzsyzt/scjzgdz 170808/index.htm). As reported by the China Earthquake Network Center (CENC), the magnitude 7.0 Jiuzhaigou earthquake occurred at 21:49:24 local time (13:19:46 UTC) and was located at 103.82°E, 33.20°N with a depth of 20 km (Figure 1). Waveform inversions indicated that the seismic moment magnitude M W was approximately 6.5, and the focal mechanism solutions indicated dominantly strike-slip faulting (Han et al., 2018;Liang et al., 2018b;Yang et al., 2017;Yi et al., 2017). The results published by different authors were comparable, despite slight differences (Table 1). Based on coseismic slip models inverted from geodetic observations, the rupture process of the Jiuzhaigou earthquake was dominated by a left-lateral slip with a deficit of shallow slip (Chen et al., 2018;Ji et al., 2017;Shan et al., 2017;Sun et al., 2018;Wang et al., 2018). The epicentral region features complicated and intricate fault systems, which are mainly composed of the Kunlun, Tazang, Huya and Minjiang faults ( Figure 1a). However, the main seismogenic fault of the Jiuzhaigou earthquake is not directly linked to any of these faults, and a field survey failed to identify any surface ruptures or the seismogenic fault itself (Xu et al., 2017b). Aftershock relocation demonstrated that the aftershocks were confined within a narrow zone that extended from the epicenter of the mainshock both northwest toward the juncture of the Tazang and Minjiang faults and southeast to the northern end of the Huya fault (Fang et al., 2018;Liang et al., 2018a;Liang et al., 2018b;Yi et al., 2017). The distribution of aftershocks outlined a near-vertical and south-southeast-trending fault plane coinciding with one of the nodal planes of the focal mechanism solutions for the mainshock (Han et al., 2018;Liang et al., 2018b;Yang et al., 2017;Yi et al., 2017). The morphology and kinetics of this seismogenic fault were consistent with the south-southeast strike and sinistral nature of the Huya strike-slip fault, indicating that the Jiuzhaigou earthquake might have occurred on a young fault representing a blind extension of the Huya fault (e.g., Sun et al., 2018;Yang et al., 2017;Yi et al., 2017).
The 2017 Jiuzhaigou earthquake, with a magnitude of 7.0, was the third largest earthquake on the eastern margin of the Bayan Har block since the 2008 Wenchuan M S 8.0 and 2010 Lushan M S 7.0 earthquakes. The Bayan Har block is one of the major components involved in the eastward expansion of the Tibetan Plateau, and its eastward movement is accommodated by significant (~10-12 mm/a) sinistral slip along its northern boundary, the Kunlun fault (Zhang et al., 2004). The slip rate on the  seismic Lg-wave Q distribution at 1.0 Hz (Zhao et al., 2013), overlaid on the locations of faults and sutures (black lines) and stations (triangles), with the relevant legend. The focal mechanisms of the Jiuzhaigou earthquake (red) and nearby large earthquakes (black) were retrieved from GCMT (http://www.globalcmt.org), the smaller aftershocks (circles) are plotted according to their magnitudes and depths, respectively. Kunlun fault was found to decrease toward the east and was partially absorbed by deformation surrounding the fault tip (Kirby et al., 2007). The left-lateral motion is also transferred to the horsetail-shaped splay faults associated with the Kunlun fault, including the Tazang, Huya and Minjiang faults (Xu et al., 2017b). Subsurface crust images from P-wave receiver functions suggest that these processes may result in a massive thrust nappe structure beneath the Minjiang fault and in the uplift of this region (Liu et al., 2017). The intense crustal deformation and complicated fault system produce intense seismic activity on the eastern margin of the Bayan Har block. In particular, the Jiuzhaigou region was located in a historical earthquake gap and has been identified as an area at risk for potentially large earthquakes (Xu et al., 2017a).
∆σ Previous investigations of the 2017 Jiuzhaigou earthquake have provided information on the mainshock and have shed light on the regional tectonic setting as well. Defined as the difference between average stresses on a fault before and after an earthquake, the stress drop is one of the fundamental source-scaling parameters used to characterize an earthquake, and it provides information on the physics of the rupture processes and exhibits regional variation related to the tectonic setting (Abercrombie and Rice, 2005;Allmann and Shearer, 2009;Shearer, 2009;Shearer et al., 2006). Given the priori assumptions made in the source model (e.g., Brune, 1970), estimating the stress drop is not easy because it is closely related to the shape of the source spectrum, which can be distorted by seismic attenuation, especially at higher frequencies. To retrieve the source spectrum from the observed spectrum of a large earthquake, the empirical Green's function (EGF) method takes the spectra of nearby small earthquakes, which have sufficiently high corner frequencies and are flat over the observation band, as an approximation of Earth's response to a delta-function source and can then correct the spectrum of the target event for attenuation and other path effects (e.g., Allmann and Shearer, 2009;Shearer et al., 2006). However, the EGF method requires small earthquakes close enough to the target one and cannot remove the effects of near-source attenuation that may vary across regions. Another approach for large datasets in which an event is recorded by multiple stations is to directly solve for the attenuation term using generalized inversion methods and then to compensate the spectrum using the calculated attenuation (e.g., Oth, 2013;Oth et al., 2011). However, this approach can result in difficulties due to the trade-off between the source and attenuation terms. The eastern Tibetan Plateau is characterized by strong attenuations with significant lateral variations. Therefore, a reliable high-resolution broadband attenuation model could greatly improve the accuracy of the stress drop measurements for the Jiuzhaigou earthquake and its aftershocks.
In this paper, based on corrections using a regional Lgwave Q model covering the eastern Tibetan Plateau with a resolution of approximately 1 degree and a broad frequency band from 0.05 Hz to 10.0 Hz (Zhao et al., 2013a(Zhao et al., , 2013b, we obtained the Lg-wave source spectra for the 2017 Jiuzhaigou earthquake and its aftershock sequence. Adopting the Brune (1970) source model, we estimated the seismic moments and corner frequencies and calculated the regional average stress drop for the Jiuzhaigou area. The temporal and spatial variations in the stress drop were investigated, improving the understanding of seismic energy release during the Jiuzhaigou earthquake. Because source spectra were retrieved independently at individual frequencies without prior assumptions, such as the constant stress drop model, which leads to conclusions of self-similar source physics, we were able to investigate the self-similarity of the Jiuzhaigou earthquake sequence, which could provide hints regarding the tectonic background and physical mechanism of these earthquakes.

Data and methods
The 2017 Jiuzhaigou earthquake and its aftershock sequence generated abundant seismograms recorded by broadband digital seismic stations at regional distances. Based on the CENC catalog, 194 aftershocks occurring before 6 October 2017 with magnitudes greater than 2.0 were selected, and their origin times, epicenters and depths were obtained from the catalog. We collected seismic data from 166 CENC stations from the China National Digital Seismic Network (CNDSN) to investigate the source spectra for the 2017 Jiuzhaigou earthquake and its aftershocks.
To model the source spectra using Lg waveforms, we followed the method used by Zhao and his coworkers (Zhao and Xie, 2016;Zhao et al., 2010Zhao et al., , 2013aZhao et al., 2013b). First, the Lg waveforms were extracted using a group velocity window of 3.6-3.0 km/s, and a time window with the same length was used before the first Pwave arrival to pick the pre-event noise. Spectra for both Lg-wave and noise were calculated and sampled at 58 individual frequencies log-evenly distributed between 0.05 Hz and 10.0 Hz. Next, we calculated the signal-to- noise ratios at individual frequencies and dropped out low-quality data with a signal-to-noise ratio less than 2.0. Last, by conducting the noise correction, , where A(f) is the spectral amplitude and the subscripts S, O and N denote the true signal, the observations from the raw Lg waveforms, and the noise series, respectively.

A ki
For the observation of event k at station i, the Lg-wave spectral amplitudes can be written in the following form where S is the source term, G is the geometrical spreading factor, is the attenuation term that can be expressed by an integral along a great circle ray path from event k to station i, P is the site response term and R is the computational error and cumulative random effects in the Lg propagation between station i and event k. The geometrical spreading is , with representing the epicentral distance from event k to station i and representing a reference distance fixed at 100 km. In this study, the path attenuation was computed using a regional Lg-wave attenuation model (Zhao et al., 2013b). This Lg-wave attenuation model was derived using tomography based on the single-station and two-station Lg amplitude, and the inversion was done at individual frequencies independently without any assumptions on frequency dependence. Strong lateral variations correlated to the geological structures and tectonics were revealed in this attenuation model, which provided useful information on the deformation of the crust in the Tibetan Plateau. According to formulas from Zhao et al. (2010), the path attenuation from an event to a station can be directly calculated given an Lg-wave Q model. Therefore, with known attenuation and geometric spreading terms, we can strip the propagation effects from the raw Lg-wave spectra from individual stations. Then, by fitting spectra from multiple stations with the theoretical source model using the least-squares orthogonal factor decomposition method (LSQR, Paige and Saunders, 1982), we obtained the source spectrum S. In practice, to take site effects into account, we allowed small perturbations in the attenuation to deal with the trade-off between attenuation and site terms, and we treated the site response as unresolved error. In this way, we hope to reduce potential effects from uncertainty in site response. In addition, the majority of the stations used in this study are located off the Tibetan Plateau, and these stations are generally located in surrounding high-Q regions and have similar tectonic backgrounds. Consequently, the site response probably varies little among these stations. However, a simplified assumption regarding site effects could still limit the reliability of our measurements.

Seismic moment and source spectrum corner frequency
By removing attenuation effects from the raw spectra using the high-resolution broadband Lg-wave Q model, we obtained the Lg-wave source spectral amplitudes at 58 discrete frequencies for the 2017 Jiuzhaigou earthquake and its 166 aftershocks. For these source spectra, their high-frequency amplitudes fell off at a slope of approximately -2 on the log-log scale. Therefore, the commonly used source model (Brune, 1970) could provide a good fit to Lg-wave spectral data from the Jiuzhaigou earthquake sequence. Based on the theoretical model, the source spectrum can be written as where is the seismic moment, is the corner frequency and and are the average density and shear wave velocity in the source region (Herrmann and Kijko, 1983;Street et al., 1975). A typical value of 2.7 kg/m 3 was assigned to the density , and the shear wave velocity km/s was obtained from Bao et al. (2015). and can then be determined by minimizing the L2 norm of the residuals between the theoretical source function and the network-determined source spectral data. The bootstrap method (Efron, 1983) was used to provide the standard deviation of the two parameters. Examples of the inverted source spectral amplitudes and the best-fit source functions for selected events are shown in Figure 2, and the inverted and for all events are listed in Table 2. The seismic moment and corner frequency for the 2017 mainshock were N·m and Hz. For aftershocks with local magnitudes between 2.0 and 4.9, the seismic moments and corner frequencies ranged from N·m to N·m and from 0.95 Hz to 0.062 Hz, respectively. Based on the linear regression, the relationship between the magnitude and seismic moment is ( Figure 3).

Stress drop
For a circular fault, the stress drop is given by where r is the radius of the fault (Eshelby and Peierls, 1957) and can be expressed as where is the observed corner frequency, v is the shear wave velocity in the source region and k is a constant depending on the specific theoretical source model. Assuming the rupture velocity to be 0.9v, Madariaga (1976) performed finite differences calculations and found that k is 0.21 for S-waves and 0.32 for P-waves given an high-frequency falloff rate. Combining equations (3)-(4), the Brune-type stress drop can be estimated from the best-fit seismic moment and corner frequency using The stress drops calculated from spectral data obtained for the Jiuzhaigou earthquake sequence by using equation (5) showed large variations from 47.1 kPa for a low stress drop aftershock to 7.1 MPa for the mainshock. The median value of the stress drop was 59.4 kPa, and most values were between 50 kPa and 75 kPa. The 7.1 MPa stress drop for the mainshock was close to the 6.0 MPa median value

Figure 2
Retrieved Lg-wave source excitation spectra for the August 8, 2017, Jiuzhaigou earthquake and selected aftershocks. Crosses are directly inverted source functions from observed Lg spectra after correction for the propagation effect. Blue solid lines are best-fitting source models, and pink shaded areas are the standard deviations. The resulting M 0 and f c values are labeled with their standard errors in each plot  (Wang et al., 2017a). The difference may be due to several reasons, including (1) the Lg-wave and strong motion data cover different frequency bands; (2) the estimated stress drop is sensitive to uncertainties in the and measurements; (3) the theoretical source model may be oversimplified (Gallovič and Valentová, 2020). In this study, the source parameters, e.g., the seismic moment, corner frequency and stress drop, were measured for the entire sequence using the same method, which allows us to investigate the relative variations in these parameters among the mainshock and its aftershocks and to explore their spatial and temporal variations during the development of the sequence.

Discussion
The stress drops diminished as the seismicity decreased in intensity, and the temporal variation highlights the overall decay of the stress drop, especially the abrupt decrease since the mainshock (Figure 4). This pattern possibly indicates that the accumulated energy was relatively fully released by the mainshock (Wang et al., 2018). The map of stress drop reveals very large variations occurring over short length scales, which may indicate complicated source mechanisms for the 2017 Jiuzhaigou earthquake and its aftershocks ( Figure 5). Close to the mainshock, which had the highest stress drop of 7.1 MPa, the aftershocks were relatively sparser, and their stress drops were less than the median value. Furthermore, the less-than-median stress drops seemingly concentrated near the trace of the major seismogenic fault, which could be inferred from the major axis of the aftershock zone, while higher stress drops appeared off the trace of this fault. The coseismic displacements inverted from geodetic data suggested that the slip on the seismogenic fault was primarily concentrated at approximately 7-15 km depth (Ji et al., 2017;Wang et al., 2017b;Wang et al., 2018). On the fault plane with a dip angle of 70° obtained from focal mechanism solution (Wang et al., 2017b), onto which all events were projected, there appeared to be no clear or strong correlation between the stress drop and depth, but we observed a tendency that aftershocks with stress drops higher than the median value of 59 kPa were more likely to occur at depths of ~15 km (corresponding to a distance of 13 km along dip), just below the depth of significant slip ( Figure 6). These observations could be related to the fact that the energy accumulated within the seismogenic fault associated with the Jiuzhaigou earthquake was relatively thoroughly released by coseismic displacement, and thus, the aftershocks off the seismogenic fault and outside of the zone of large slip released larger amounts of stress than those on the fault and within the zone of large slip.
The assumption of earthquake scaling self-similarity, suggesting that the physics of earthquakes is independent of their sizes, implies that the stress drop remains constant over earthquakes with a wide range of magnitudes (Aki, 1967). It can be predicted that on a log-log plot, the shape of all source spectra will be identical with offset along a line (Allmann and Shearer, 2009;Prieto et al., 2004;Shearer, 2009). The self-similarity appears to be roughly true for most earthquakes except for very large ones. Shearer et al. (2006) found little dependence of stress drop on moment and suggested that over 60,000 earthquakes in southern California with magnitudes ranging from 1.5 to 3.1 were self-similar. Allmann and Shearer (2007) found that the stress drops of over 40,000 earthquakes remained nearly constant with moment, implying selfsimilarity across earthquakes with magnitudes between 0.5 and 3.0 in central California. At the global scale, selfsimilarity in moderate to large earthquakes has also been suggested (Allmann and Shearer, 2009). However, a strong dependence of stress drop on earthquake size was observed at small scales, which could be related to the local stress regime (Oth, 2013). Moreover, Zhao et al. (2011)  that the stress drop increases with earthquake magnitude, favoring a nonconstant stress drop model and thus rejecting the assumption of self-similarity. In this study, a constant stress drop can be calculated from linear regression between and with a fixed slope of -3, and we found that the constant average stress drop for the 2017 Jiuzhaigou earthquake sequence was 93.0 kPa (Figure 7). However, by careful investigation, these events can be further grouped into two clusters in Figure 7, with blue symbols representing events with stress drops lower than the regional average of 93.0 kPa and red symbols representing events with stress drops higher than the regional average (Figure 7). Performing linear regression on each of the clusters of events, we obtained values of 57.4 kPa and 318.9 kPa, which provided a better fit for both clusters of events. For the blue cluster, the values of stress drop remained nearly constant versus the seismic log 10 f c log 10 M 0 moment over 3 orders of magnitude, indicating that these aftershocks were self-similar in their source physics. These events could be related to the left-lateral slip motion driven by the eastward expansion of the Tibetan Plateau (Sun et al., 2018). For the red cluster, the values of stress drop were systematically higher and more scattered, especially for the mainshock, which had a significantly high stress drop. This may be related to the very large magnitudes that caused the rupture behaviors of these events to differ from those of other moderate events in terms of the aspect ratios. The two groups of events were characterized by different levels of self-similarity, which can be attributed to a nonconstant stress drop (Taylor et al., 2002). The level can be quantified by the deviation of slope from -3 when performing linear regression between and ( Kanamori and Rivera, 2004). Overall fitting produced a slope of -2.52, but separate linear regression provided better fit for both clusters of events ( Figure 8). For the red cluster, the slope was less than -3, demonstrating that the stress drop values increased with earthquake size. In contrast, the slope for the blue cluster was -2.97, which are very close to -3, suggesting that the stress drop of events in this cluster was nearly constant and that these events were self-similar. The electrical structure obtained using magnetotelluric imaging suggested that anomalies of high conductivity extend across the northern part of the Huya fault, with its top boundary at depths of ~10-20 km in the Jiuzhaigou earthquake source region (Sun et al., 2017;Zhao et al., 2012). The majority of events in the blue cluster (i.e., 86 of 120 events) were located deeper than 10 km and were very likely to be affected by the highconductivity structure. This improves our understanding of the observed low stress drop values of events in the blue cluster because the presence of elevated fluid content, as indicated by the high-conductivity anomaly, could reduce the stress applied on the faults by increasing the fluid pore pressure (e.g., Sumy et al., 2017). However, because the epicentral depths of the Jiuzhaigou earthquake sequence were obtained from a catalog and because the depth resolution of the magnetotelluric imaging might be limited, detailed investigations into the impact of the presumed viscous crustal flow on the seismogenic dynamics of local seismicity are still required. These investigations should include analysis based on high-precision relocation results and further studies on the fine-scale local electrical structure as constrained by high-resolution magnetotelluric images.

Conclusions
The stress drop is one of the fundamental parameters Each circle represents an earthquake and is filled according to its stress drop value. Higher-than-median stress drops more likely appear for aftershocks at depths deeper than ~15 km (corresponding to ~13 km along dip) and in areas where coseismic slip was less than ~ 40 cm. The red hexagon denotes the mainshock, and the light gray area on the right represents a lack of coseismic slip data coverage Assuming the scaling relationship of and a constant stress drop model, linear regression yields an average stress drop for the region near the Jiuzhaigou earthquake of 93 kPa (black line). These events can be divided into two clusters according to the average stress drop (blue circles denote events with stress drops lower than the average, and red circles denote events with stress drops higher than the average). Separate linear fitting for each cluster yields a stress drop of 57.4 kPa for the events in the blue cluster (blue line) and a stress drop of 318.9 kPa for the events in the red cluster (red line) The grouped red and blue clusters of earthquakes are same as those in Figure 7. The results of linear regression between log 10 f c and log 10 M 0 for all events, blue cluster and red cluster are denoted by black, blue and red lines, respectively. Linear regression suggests that the events in the blue cluster, with a scaling factor very close to -3 and a nearly constant stress drop, exhibit self-similarity. The events in the red cluster are characterized by a scaling factor less than -3, which means that the stress drop increases with increasing seismic moment instead of exhibiting self-similarity, a possible indicator of different source physics characterizing the dynamic rupture of an earthquake. It is also useful in investigating the process of stress field adjustment after a major earthquake through the analysis of its aftershock sequence. The stress drop of a seismic event can be estimated from seismic observations through the seismic moment and corner frequency of the source spectrum. However, effects from the attenuation along the propagation path can distort the shape of the source spectrum, especially at higher frequencies and in regions with strong lateral attenuation variations. In this study, we used a high-resolution broadband Lg-wave Q model to compensate for the effect of attenuation in Lg-wave source spectra. The corrected data were fitted with the theoretical source model to obtain source parameters for the 2017 Jiuzhaigou earthquake and its 166 aftershocks with magnitude greater than 2.0. The seismic moment and corner frequency were then determined using the bootstrap method, and the stress drop was estimated for individual events. The obtained seismic moments were linearly related to the magnitude . The Brune-type stress drop was calculated for individual events. The stress drop for the mainshock was 7.1 MPa, while the median value of stress drops from all aftershocks was less than 60 kPa. The temporal variation in the stress drop demonstrated an abrupt decay after the mainshock. The aftershocks very close to the mainshock, near the major axis of the aftershock zone, or within the zone with large coseismic slip exhibited relatively low stress drops. The above observations suggested that the stress energy accumulated by the major seismogenic fault was largely released by the mainshock. Based on the relationship of the seismic moment versus the corner frequency on the log-log scale, the aftershocks can be separated into two groups, one of which was characterized by self-similarity with a nearly constant stress drop of approximately 58 kPa and a scaling relation of . The events in the second cluster exhibited increasing stress drop with increasing seismic moment, possibly suggestive of different source physics related to dissipation of the energy on the major seismogenic fault and localized stress heterogeneities resulting from complicated tectonic activities.