The Effect of the Wenchuan and Lushan Earthquakes on the Size Distribution of Earthquakes along the Longmenshan Fault

: Changes in the stress state of faults and their surroundings is a highly plausible mechanism explaining earthquake interaction. These stress changes can impact the seismicity rate and the size distribution of earthquakes. However, the effect of large earthquakes on the earthquake size distribution along the Longmenshan fault has not been quantiﬁed. We evaluated the levels of the b value for the stable state before and after the large earthquakes on 12 May 2008 (Wenchuan, M S 8.0) and 20 April 2013 (Lushan, M S 7.0) along the Longmenshan fault. We found that after the mainshocks, the size distribution of the subsequent earthquakes shifted toward relatively larger events in the Wenchuan aftershock zone ( b value decreased from 1.21 to 0.84), and generally remained invariable in the Lushan aftershock zone ( b value remained at 0.76). The time required for the b value to return to stable states after both mainshocks was entirely consistent with the time needed by the aftershock depth images to stop visibly changing. The result of the temporal variation of b values shows decreasing trends for the b value before both large earthquakes. Our results are available for assessing the potential seismic risk of the Longmenshan fault as a reference. in can


Introduction
Following the Wenchuan M S 8.0 earthquake on 12 May 2008, the Longmenshan fault zone was struck by the 20 April 2013 M S 7.0 Lushan earthquake. The Longmenshan fault zone is composed of several almost parallel thrust faults, forming a boundary fault between the Sichuan Basin and Tibetan Plateau, and controls the seismicity of the Longmenshan region ( Figure 1). The epicenters of the Wenchuan and Lushan earthquakes were approximately 87 km apart, and the focal mechanism of both events showed a thrust rupture [1].
According to the characteristics of the Wenchuan and Lushan earthquakes, whether the M S 7.0 Lushan event was a strong aftershock of the M S 8.0 Wenchuan earthquake or a new and independent event has been a topic of debate. For example, some researchers suggest that the two large earthquakes were independent events. The reasons are as follows: (1) there is no overlapping area between the Wenchuan and Lushan earthquake rupture zones [2]; (2) the Wenchuan and Lushan earthquakes were generated in different faults in the Longmenshan fault zone [3]; (3) the rupture processes of the Wenchuan and Lushan earthquakes were different, and the aftershock zones of the two events were nearly 45 km apart [4]. Alternatively, some scientists propose that the Wenchuan and Lushan earthquakes were a mainshock-aftershock sequence and note that the Lushan event struck in an area where Coulomb stress was increased due to the Wenchuan earthquake [5,6]. The controversy over the relationship between the Wenchuan and Lushan earthquakes highlights the complexity of earthquake interaction in the Longmenshan fault zone. It is widely accepted that earthquake interactions can be understood by identifying changes in static and dynamic stress states around faults [7][8][9]. The most observable effect of this stress change is a significant increase in seismicity rate, which is generally considered an aftershock phenomenon [10][11][12]. Statistically, aftershock activity is classically described by n(t) = K/(t + c)p, where n(t) is the number of aftershocks after time t and K, c, and p are constants that describe the aftershock productivity [13]. Ogata [14,15] described aftershock activity as a multigenerational branching process and proposed the epidemictype aftershock sequence model, which is a stochastic point process model of self-exciting point processes.
Changes in stress can impact the seismicity rate and the frequency size distribution, which is alternatively known as the frequency-magnitude distribution (FMD) [16] or Gutenberg-Richter (G-R) law [17] and is expressed as logN = a−b M, where N is the number of events in a given time period with magnitude greater than M, a describes the seismicity of a given seismogenic volume, and b is the slope of the FMD. Previous studies showed that b values fall within the range of 1.02 ± 0.03 on a large scale for a long time [18,19]. For regions on a smaller scale, the b values show a broad range of spatial and temporal variations. For example, the b value ranged from 0.5 to 2.5 in the Andaman-Sumatra region and California [20]. Interpretation of the variation of b values is based on several factors, including stress state [21,22], focal depth [23], faulting style [24,25], fluid pressure [26], and so on.
The earthquake size distribution generally follows a power law, with a slope of b values, which characterizes the relative occurrence of large and small events. A low b value indicates a larger proportion of large earthquakes and vice versa. Zhao et al. [27] (2008) compared the spatial footprint of b values before and after the Wenchuan earthquake in the Longmenshan fault zone, and the results showed that the b values tend to change from lower in the southern region to higher in the northeastern region. The temporal change in b values before the Wenchuan M S 8.0 earthquake showed a decreasing trend [28][29][30], and the Lushan M S 7.0 earthquake showed similar temporal change trends in b values [31]. These studies focused on the variation trend of b values before and after the two large events along the Longmenshan fault. However, the fundamental effect of the Wenchuan and Lushan earthquakes on the size distribution of earthquakes along the Longmenshan fault has not been quantified, which limits our understanding of how the apparent stress changes in the region affect the size distribution of earthquakes.
In this study, we evaluated the spatiotemporal evolution of the b values along the Longmenshan fault in the past nearly 20 years. Moreover, we estimated the levels of the b value for the stable state before and after the Wenchuan and Lushan earthquakes and quantified the effects of the two large earthquakes on the size distribution of subsequent events at different times. In addition, the spatial evolution process of the deep seismogenic environment in the Wenchuan and Lushan aftershock zones in two and three dimensions was illustrated via spatial scanning and data fitting, which can be used to analyze the aftershock activity of the two large earthquakes.

Data and Postulates
The earthquake catalog we used here was documented by the regional seismic network and then verified by the China Earthquake Networks Center (CENC) along the Longmenshan fault during the period from 1 January 2000 to 1 January 2019. It is a relatively complete catalog containing the Wenchuan-Lushan earthquake sequence. Figure 2 shows the magnitude-time distribution of earthquakes in the Wenchuan source region and Lushan source region. The locations of the earthquakes in this catalog were corrected for accuracy. Each event includes the time, location, magnitude, and depth. Homogeneity of the catalog was iterated and optimized for subsequent research and analysis. Stress change has been widely used to interpret the triggering of the mainshockmainshock and mainshock-aftershock events [7,[32][33][34][35][36]. However, the exact measurement of stress states is difficult; thus, a relationship between the stress state and b value has been proposed [16,24,37,38]. Schorlemmer et al. [24] demonstrated that the b value could be regarded as a stress indicator that depends inversely on differential stress. Therefore, changes in the stress state on faults lead to the variation in b values, which is followed by time-dependent recovery. In general, the larger the magnitude of the earthquake, the greater the stress changes and the larger the b value fluctuations. No event of magnitude larger than M S 7.0 had been reported in the historical record of the Longmenshan fault zone and the catalog we used in this paper contains more than 80,000 events and only two large earthquakes greater than magnitude 6.5, that is, the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes. Thus, these two events are mainly responsible for the apparent changes in stress state along the Longmenshan fault zone in the past 20 years.
Therefore, there are two postulates: first, without a large earthquake perturbation, the b value will remain in a stable state with a small fluctuation range for a long time; second, after the perturbation of a large earthquake, the b value may recover to another stable state. We evaluated the levels of the b value for the stable state before and after these two large earthquakes in the study region.

Methods
The main methods we used include estimation methods (MaxCurvature for the estimation of the completeness magnitude; maximum likelihood estimation (MLE) for the b value estimation), test methods (Akaike information criterion (AIC) for the variation trend of b values; the nonlinearity index (NLIndex) for the linearity assessment of frequencymagnitude distribution) and kriging interpolation to describe the spatial and temporal evolution images of aftershock focal depth.

Completeness Magnitude (M C ) and b Value Estimation
The estimation of the completeness of earthquake catalogs is essential to the computation of b values, and the lowest magnitude of all earthquakes that are reliably detected in a space-time volume is defined as the completeness magnitude (M C ) [39]. The lower the M C , the higher the detection capability. Here, we use the MaxCurvature technique, which estimates the M C by locating the magnitude that is the highest frequency of events in the FMD. Mignan [40] showed that the MaxCurvature technique underestimates the M C in cases involving gradually curved FMDs and postulated that this underestimation tendency arises from spatiotemporal heterogeneities within the earthquake monitoring network. Therefore, we used the corrected MaxCurvature method with a correction factor of +0.2 [41], and the uncertainties were determined by bootstrapping.
The least-squares method and maximum likelihood estimation are often used to calculate the b value, and the latter approach is considered more stable. In this work, the maximum likelihood estimation used to calculate the b value and its standard deviation [42,43] where M is the average magnitude of earthquakes with M ≥ M C ; M C is the cutoff magnitude. The confidence limit of the b value is expressed as follows: where N is the number of earthquake cases of the given sample.

Estimation of the Frequency-Magnitude Distribution (FMD) Extrapolation
The nonlinearity index (NLIndex) can be used to assess whether the extrapolation of a given high-magnitude FMD is likely an overestimate or underestimate of the probable rates for large events [44]. If NLIndex ≤ 1, the FMD is regarded as linear, and if NLIndex > 1, the FMD is not linear. The slope of M cut is clearly positive or negative, respectively, indicating that the FMD overestimates or underestimates large M rates.

Akaike Information Criterion
To quantify the changing trend of b values, the P test was conducted for b values in two sample windows based on the Akaike information criterion (AIC) [45]. Hypothesis 1: the b values in the two sample windows are the same; Hypothesis 2: the b values in the two sample windows are different, represented as b 1 and b 2 , respectively. The hypothesis of the difference of the AIC leads to the difference ∆AIC [46]: where N i is the number of events in the sample windows and b i is b values in the sample windows. P b represents the probability that the events in the two sample windows come from the same population and can be derived from the AIC as follows: The b value in the sample window represents a significant change when ∆AIC ≥ 2 (P b ≈ 0.05) and is highly significant when ∆AIC > 5 (P b ≈ 0.01) [47].

Kriging Interpolation
Kriging interpolation is the most commonly used geostatistical approach for spatial interpolation. With this method, a semivariogram is used to express the spatial relationship of the distance between samples. This technique depends on the spatial model between samples to predict attribute values at unsampled locations [48]. As a widely used interpolation method, kriging takes into account the distance between unknown positions and the sample locations as well as the distance between sample locations, effectively reducing the interference of clustering in samples on the accuracy of the interpolated estimates [49]. We used the kriging interpolation algorithm to produce maps incorporating anisotropy and underlying trends from irregularly spaced data. The exponential semivariance model with the smallest prediction errors was chosen over the Gaussian and spherical models for the spatial interpolation of focal depth data.

Completeness Magnitude (M C ) and Linearity Assessment of Frequency-Magnitude Distribution (FMD)
As shown in Figure 3, the results of the corrected MaxCurvature method show that the M C of the earthquake catalog used in this work is M C = 1.5. This result is consistent with the results of previous studies on the Wenchuan earthquake zone [30,50]. Fang et al. [51] described in detail the aftershock performance and analysis of the Lushan earthquake based on the combined data from permanent and temporary seismic stations. They concluded that the minimum complete magnitude was M = 1.0. To unify the consistency of the M C of the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes in the Longmenshan fault zone, we selected events with magnitudes of M ≥ M C = 1.5.
We performed a linearity check on FMD, and the results are shown in Figure 4. The NLIndex (red) is shown for different cutoff magnitudes (upper inset) and the NLIndex ≤ 1 for all cut off magnitudes; thus, the linear FMD is accepted as the best M C .

Time-Space Analysis of b Values
Earthquake frequency will increase immediately within a short time after a large event and may exceed the recording capacity of the seismic network. Before establishing the time-space series of b values with aftershocks, we should eliminate the events documented in the early catalog, which is somewhat heterogeneous and incomplete in small events [41]. In this work, the exclusion period depends on the magnitude of completeness over time. Therefore, we first removed the events documented in the initial catalog within two months after the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes, a period for which the data are highly incomplete. Then, we calculated the spatiotemporal distributions of b values before and after two large events that occurred from 2000-2019 along the Longmenshan fault by selecting events with M ≥ M C = 1.5 and using a time window and spatial grid to calculate the b values. In this computation, the window lengths were set to at least 500 events in the Wenchuan aftershock zone and 200 events in the Lushan aftershock zone. Each window was moved forward by one event at a time. Figure 5a,b display time series of the b value in source regions, and the overall change trend conforms to our postulation that b values will undergo relatively significant changes in a period of time before and after a large earthquake. Specifically, b values show a decreasing trend before the occurrence of both large earthquakes in both zones. To ensure that this trend is statistically significant, we quantitatively assessed the temporal variation in b values using the P parameter test and selected three windows before the M S 8.0 event (W 1 , W 2 , and W 3 ) and the M S 7.0 event (L 1 , L 2 , and L 3 ). Window selection was based on the significance of changes in b values (Figure 5a,b). The results are shown in Table 1. The b value in the sample window has significantly changed when ∆AIC ≥ 2 (P b ≈ 0.05). Table 1 shows that the b values decreased before both large earthquakes with statistically significant variations.
After the Wenchuan M S 8.0 event, the b values in the Wenchuan aftershock zone experienced a period of dramatic fluctuation (indicated by the pink shading lasting not more than one year) before gradually stabilizing within a small fluctuation range (Figure 5a), which was similar to the range of b values in the third period (Figure 5c, b = 0.84). After the Lushan M S 7.0 event, the b values in the Lushan aftershock zones increased rapidly and then slowly dropped to a stable state (Figure 5b), which was similar to the FMD in the first period (Figure 5d, b = 0.76). As shown in Figure 5b (red shading), the b value required less than ten months to return to a stable state.
The  (Figure 5a). However, earthquakes that were greater than M S 4 struck the entire Wenchuan aftershock zone after the M S 8.0 event (Figure 1). In addition, there were no strong aftershocks above magnitude 6.5 along the faults, and only six events greater than M S 4 occurred within two months after the mainshock. These phenomena indicate that without the continuous perturbation of strong aftershock, the b value gradually stabilized to the state mainly determined by the background earthquakes, which tended to shift to larger events following the M S 8.0 event in the Wenchuan aftershock zone.
After the M S 8.0 Wenchuan earthquake, the Lushan aftershock zone also experienced a "seismic quiescence" of approximately two years, and only one event greater than M S 4 occurred before the M S 7.0 Lushan earthquake during the period when b values were significantly decreasing (Figure 2b). Moreover, only a few events greater than M S 4 occurred within two months after the mainshock; the subsequent events basically returned to the magnitude of the background earthquakes before the mainshock, which shows that the b values in the Lushan aftershock zone eventually returned to the background level (Figure 5b).
To analyze the spatial footprints of the changes in the b values, we divided the two study regions into 0.1 • × 0.1 • grids, sampled the 300 events nearest to each grid node within radius of 30 km, and re-estimated the M C in each node. For this purpose, we used a bootstrap approach to sample the events 1000 times randomly. The spatial footprints of the changes in the b values are consistent with the FMDs (Figure 5c,d).    The b value in the southern part of the Wenchuan source region and Lushan source region was lower than that in the northern part of the Wenchuan aftershock before the M S 8.0 Wenchuan event (Figure 5e,h). This pattern is consistent with the characteristics of the Longmenshan fault, which is a strike-slip fault in the north and a thrust fault in the southern section. It is generally considered that the b value is inversely proportional to stress, and the b values of different types of faults are as follows: b (normal) > 1, b (strike-slip)~1, and b (thrust) < 1 [24,25]. The conditions changed markedly after the M S 8.0 Wenchuan event; the b values decreased in the Wenchuan source region and Lushan source region (Figure 5f,i). This finding illustrates the effect of the Wenchuan earthquakes on stress change along the Longmenshan fault. Figure 5g,j illustrate the stable state of the b value in the Longmenshan fault zone.

Evolution of Images of Aftershock Activity Depicted by Focal Depth
Spatial scanning was performed using events with depth data in the catalog, i.e., at a step size of 0.1 • for both longitude and latitude. For all earthquakes in each 0.1 • × 0.1 • grid point, the average depth was used as the depth value of the grid point, and then kriging interpolation was applied to all the grids. We counted at least ten events in each grid node in the Wenchuan aftershock zone as samples (as well as five events for the Lushan aftershock zone) to prevent the average depth of grid points from being affected by too few events. The contour lines of depth distribution at different periods after the mainshock were obtained and superimposed with regional faults (Figures 6 and 7).
To illustrate the evolution of the deep seismogenic environment in the Wenchuan aftershock zone. Figure 6 shows the spatial evolution of the aftershock depth at one day, one week, one month, six months, one year, and three years after the mainshock. The analysis shows that the focal depth of the Wenchuan aftershock zone spread along the direction of the fault. With Mianyang as the boundary, the aftershock zone can be divided into the southern section and the northern section. The depth distribution is generally deep in the southeast and shallow in the northwest, and the average aftershock focal depth is 10-15 km.
A comparison of the aftershock activity images depicted by the focal depth information shown in Figure 6d,e revealed that the image formed one year after the mainshock did not show a significantly different pattern in the following two years. This finding indicates that after the mainshock, the aftershock frequency tends to be stable one year later, which means that the active aftershock period of the Wenchuan M S 8.0 was less than one year. Figure 7 shows the spatial evolution of aftershock depths at one day, one week, five months, ten months, one year and three years after the mainshock in the Lushan aftershock zone. The analysis shows that the focal depths of the Lushan aftershock are distributed around the fault, with the fault as the boundary, with a trend of deep in the southeast and shallow in the northwest. Moreover, a comparison of the regional aftershock activity images depicted by the focal depth information in Figure 7d,e revealed that the pattern of the image formed ten months after the mainshock presented limited changes in the following one year, which indicates that the aftershock frequency tended to be stable ten months after the Lushan mainshock. Therefore, the active aftershock period of the Lushan M S 7.0 earthquake was less than ten months.

Discussion
Earthquake interaction can change the stress state of faults, which is reflected in both earthquake activity rate and earthquake size distribution. The relationship between stress state and b value can be used to quantify the effect of large earthquakes that can significantly change the stress states and, therefore, the earthquake size distribution. In contrast to previous studies on the Wenchuan and Lushan earthquakes, our paper focuses on a quantitative analysis and discussion of the effect of the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes on the size distribution of the earthquakes along the Longmenshan fault at different times.
Interpretation of the b value and its variability according to physical mechanisms has received considerable attention and discussion. In most cases, the observation of spatial and temporal b value variability can be caused by several factors: (i) Process of estimation: homogeneity of catalog and method of calculation can affect the results. All the data in this work are from the Longmenshan fault and its surroundings, and each event includes the time, location and magnitude depth. The maximum likelihood estimation and least-squares regression method [30,[41][42][43][52][53][54][55] are used to estimate the b value and its uncertainty, but the latter is excessively affected by the largest earthquake magnitude. Marzocchi et al. [56] measured the bias on b values caused by the magnitude binning and catalog incompleteness when the b value is estimated by the maximum likelihood estimation and provided guidance to reduce the likelihood of being misled by b value variation. (ii) Stress conditions: the b value and its variation represent stress buildup and release. The differential stress is inversely dependent on the b value has been observed in laboratory experiments [57,58] as well as in the field [59,60]. The stress acting on a fault may control the variation in the b value in space and time. Parsons et al. [61] calculated the regional Coulomb stress changes on major faults surrounding the rupture resulting from the Wenchuan M S 8.0 and showed that significant stress increased in the Lushan aftershock zone. Other studies obtained similar results [5,6,62]. The spatial variation in the b values throughout the Lushan aftershock zone decreased after the M S 8.0 mainshock (Figure 5i), and the same effect occurred in the southern part of the Wenchuan aftershock zone after the M S 7.0 event (Figure 5g). These findings suggest that the b value is negatively correlated with stress and indicate the effect of the earthquake interaction along the Longmenshan fault zone. (iii) Crustal tectonics: the variation in b value can be interpreted according to the tectonic characteristics, i.e., rock heterogeneity [63], focal depth [23], pore pressure [26], and fault types [24,25]. Previous studies have shown that the b value in different types of faults is b (normal) > 1, b (strike-slip)~1, and b (thrust) < 1 [26,64]. As shown in the spatial footprints in Figure 5, the b value of the southern part of the Longmenshan fault zone is lower than that in the northern part. This pattern is consistent with the tectonic characteristics of the Longmenshan fault, which is a strike-slip fault in the north and a thrust fault in the southern part [65,66].
Stress changes seem to be a key factor that affects the b value and its variation. Excepting the approach of estimation, all other factors are secondary because they are directly or indirectly affected by the stress [67]. Therefore, the observed falls in the b values shown in Figure 5a,b was interpreted as changes in the related stress conditions, which could be precursors to large earthquakes. However, these temporal variations may occur over a timescale ranging from months to years, and the timeliness and effectiveness of this variability as an indicator are difficult to guarantee. Additionally, there is usually an insufficient number of events to accurately calculate the b value before large earthquakes. Gulia and Wiemer [41] pointed out that the period following a moderate earthquake is rich in such data, with thousands of events occurring within a short period. These events may allow real-time monitoring of the evolution of b values. The authors claim that the probability of a larger earthquake following a moderate earthquake increases by several orders of magnitude if the b value remains the same or drops significantly rather than increases. However, Brodsky [68] suggested that the observed pattern revealing the changes in b values is a statistical effect rather than deterministic and that researchers need more cases to test this claim.
In general, thousands of aftershocks occur in the period following a large earthquake. Based on these abundant data, there are two typically common operational aftershock forecasting models used in aftershock hazard assessment, namely the short-term earthquake probability (STEP) model [69] and epidemic-type aftershock sequence model [14,15]. Gulia et al. [16] reported that these models forecast a high probability for a repeat of the mainshock rupture and thus substantially overestimate the aftershock hazard. This paradox can be resolved by taking into account the stress changes and their effect on the earthquake size distribution.
Our results showed that the time series of the b value in the Longmenshan fault zone after the mainshocks exhibited a period of significant fluctuation before returning to the stable state in both the Wenchuan aftershock zone (one year) and the Lushan aftershock zone (ten months). Figures 6 and 7 show that the time required for the b values to return to a stable state after both mainshocks was entirely consistent with the time required for the aftershock depth images to cease changing visibly. The spatial footprints of the changes in the b value reveal that the southern part of the Longmenshan fault zone is lower than the northern part. This finding demonstrates that the Wenchuan and Lushan events did not change the pattern of higher stress in the southern part of the Longmenshan fault zone than in the northern part. However, the most obvious change is that after the mainshock, the size distribution of the subsequent earthquakes in the Wenchuan source region shifts toward relatively larger events (lower b values).
The Longmenshan fault zone began to develop in the Late Triassic, and severe tectonic deformation occurred during the Indo-China and Himalayan movements, forming a combination of thrust and strike-slip displacement [66,70]. Previous studies did not comprehensively quantify the stable state for the Longmenshan fault zone before and after the two large events in a long time series [5,6,27,29,30,55,71]. The b value is a measurable indicator of earthquake size distribution within a specified region and period of time and is dependent on stress. With our present results, we reported the temporal and spatial variation in b values before and after two big earthquakes and fitted the source depth in time and space to quantify the stress changes and their effect on the earthquake size distribution in the Longmenshan fault zone.

Conclusions
Based on the tectonic characteristics and potential seismicity surrounding the aftershock zones of the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes, we studied the spatial and temporal variation of b values in two source regions from 2000 to 2019. In addition, the spatial evolution process of the deep seismogenic environment in the Wenchuan and Lushan aftershock zones was drawn by spatial scanning and depth data fitting.
The results depict the decreasing trends of b values before the two large earthquakes in the study region. Additionally, the b value in the Wenchuan aftershock zone took approximately one year to enter a new stable state (b values ranging from 1.21 to 0.84), while the b value in the Lushan aftershock zone took approximately ten months to return to its original stable state (b = 0.76). Moreover, the major aftershock active periods of the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes were less than one year and ten months, respectively, which are consistent with the time required for the b value to return to a stable state. The spatial footprints of the changes in the b values results reveal that the Wenchuan M S 8.0 and Lushan M S 7.0 events did not change the pattern of high b values in the north and low b values in the south along the Longmenshan fault zone.
We quantified the effect of the Wenchuan M S 8.0 and Lushan M S 7.0 earthquakes on the size distribution of earthquakes along the Longmenshan fault. Future studies can focus on how to quantify the effect of large earthquake size distribution across different tectonic regimes and apply the findings in potential seismic risk assessment.