Mapping the b-values along the Longmenshan fault zone before and after the 12 May 2008, Wenchuan, China, MS 8.0 earthquake

Theb-value in the Gutenberg-Richter frequencymagnitude distribution, which is assumed to be related to stress heterogeneity or asperities, was mapped along the Longmenshan fault zone which accommodated the 12 May 2008, Wenchuan, MS 8.0 earthquake. Spatial distributions ofb-value before and after the Wenchuan earthquake, respectively, were compared with the slip distribution of the mainshock. It is shown that the mainshock rupture nucleated near to, but not within, the high-stress (low b-value) asperity in the south part of the Longmenshan fault, propagating north-eastward to the relatively low stress (high b-value) region. Due to the significant difference between the rupture process results from different sources, the comparison between slip distribution and pre-seismic b-value distribution leads to only conclusion of the rule-of-thumb. The temporal change ofb-value before the mainshock shows a weak trend of decreasing, being hard to be used as an indicator of the approaching of the mainshock. Distribution of b-values for the aftershocks relates the termination of the mainshock rupture to the harder patch along the Longmenshan fault to the north.


Introduction
The unexpected occurrence of the Wenchuan M S 8.0 earthquake on 12 May 2008, in Sichuan Province of southwest China, provides a unique opportunity to test some of the concepts in use in earthquake physics, among which one of the important concepts is the relation between slip distribution and the distribution of locked segments along an earthquake Correspondence to: Z. L. Wu (wuzl@cea-igp.ac.cn) fault, mapped by different measures such as the distribution of microseismicity.
The concept of asperities, strong spots, as a segment of an earthquake fault that resists faulting more than its surrounding, has been used in models for earthquake rupture since long time ago (e.g.Wyss and Brune, 1967;Lay et al., 1982;Aki, 1984).Recent results in rupture dynamics show that, as a simplified concept, asperities and barriers along an earthquake fault generally reflect the complex fault geometry as well as the heterogeneities of the Earth media (Zhang and Chen, 2006a, b).As a conceptual model of earthquake rupture, stress is accumulated within an asperity, but less so along the fault surface around it.When the asperity fails, a mainshock rupture propagates to some distance into the lessstressed fault segments (Wyss et al., 2000;Zuñiga and Wyss, 2001).Recent results also show that some earthquakes exhibit the properties of "persistent asperities" (Igarashi et al., 2003;Yamanaka and Kikuchi, 2004), probably associated with the fractal nature of an earthquake fault.
One of the approaches to mapping asperities is to use the frequency-magnitude distribution of small earthquakes, or the scaling relation of small earthquake population.Previous researches show that high ambient stress can perturb the b-value in the Gutenberg-Richter-type frequency-magnitude distribution, from its normal value b=1 to lower values, and that low stress perturbs b-value to anomalously higher values (Scholz, 1968;Wyss, 1973;Urbancic et al., 1992;Lahaie and Grasso, 1999).From laboratory experiments, pore pressure records, and analysis of earthquakes related to mining, it was proposed that locally high stress can perturb the normal b=1 to low values of 0.5 (Scholz, 1968;Wyss, 1973;Urbancic et al., 1992).Creeping segments of the fault are found to be characterized by high b-values up to b≈1.5 (Amelung and King, 1997).This relation, the relation between the variation of b-value and local stress heterogeneity, provides a feasible means in mapping the asperities, and has been applied to several cases (Amelung and King 1997;Wiemer and Wyss, 1997;Wyss et al., 2000).In recent years, with the enhancement of the quality of earthquake catalogues and the implementation of the inversion of detailed rupture process, test of such a concept has become more and more practical (Wyss, 2001;Schorlemmer et al., 2004a, b;Schorlemmer and Wiemer, 2005;Chen et al., 2006;Nakaya, 2006;Sobiesiak et al., 2007;Ghosh et al., 2008).
In this study we try to investigate the distribution of bvalue along the Longmenshan fault zone which is responsible for the great Wenchuan earthquake.Previous investigation using b-value mapping shows that fault segments with different sliding behaviors exist along the central-southern portion of the Longmenshan fault zone (Yi et al., 2006).Being a unique great thrust event in continental regions, the large extent of the earthquake fault, about 300 km long propagating from SW to NE http://earthquake.usgs.gov/eqcenter/eqinthenews/2008/us2008ryan/finite fault.phpmakes it possible to use local earthquake catalogues even if with poor completeness and poor location accuracy as compared to those regions such as San Andreas Fault.After the earthquake, rupture process was quickly available (for example, http://www.geol.tsukuba.ac.jp/ ∼ nisimura/20080512; http://earthquake.usgs.gov/eqcenter/eqinthenews/2008/us2008ryan/finite fault.php;and http://www.tectonics.caltech.edu/sliphistory/2008 e sichuan/e sichuan.html,providing detailed information which can be compared with the distribution of asperities.Moreover, deployment of mobile seismic stations in the meizoseismal region provided aftershock catalogue with pretty good quality, making it possible to compare the b-value distribution before and after the main rupture.

Data used for the mapping of b-values before the Wenchuan earthquake
The region considered in this investigation is taken as 29.0 • ∼34.0 • N and 101.5 • ∼107.5 • E, including the Longmenshan fault zone as well as its surrounding regions.Note that in Chinese "shan" means mountain, this fault system is correlated with topography, as can be seen from the map. Figure 1 shows the epicenter map of the earthquakes with M L ≥2.4 in this region for the period from 1977 to the day before the Wenchuan M S 8.0 earthquake.Earthquake catalogue used is the Monthly Earthquake Catalogue provided by the China Earthquake Networks Center (CENC), compiled based on the local catalogues from regional/local seismic networks, with magnitude unified as M L .Homogeneity of the catalogue as a function of time, space and magnitude firstly need to be confirmed by the analysis of the catalogue, because of the need of as many completely-recorded events as possible for statistical resolution power (Wyss and Stefansson, 2006).Figure 2a plots the cumulative number of earthquakes with M L ≥2.4 in the catalogue, which shows a constant slop with time since 1977.This change is to much extent due to the upgrading of local seismic networks around 1976.This result about temporal homogeneity is in consistency with that of Yi et al. (2006).Figure 2b shows the magnitude-time distribution of earthquakes in this region.Figure 2c shows the frequency-magnitude distribution for earthquakes since 1977.From the fitting of Gutenberg-Richter's law, the completeness magnitude is determined as M L 2.4.Note that the analysis in Fig. 2a and c, i.e., taking of earthquakes above a cutoff magnitude to determine the period with better completeness, and taking of earth-quakes after a time to determine the completeness magnitude through Gutenberg-Richter distribution, have some trade-off with each other.As a matter of fact, selection of the starting time for analysis and the cutoff magnitude of the catalogue is a process of iteration and optimization by considering different factors.In the local catalogues, due to the limitation of location accuracy, depth cannot be determined accurately with acceptable reliability, although there have been some results of accurate location of part of the earthquakes using double difference methods (e.g.Zhu et al., 2005).Also considering that the Longmenshan fault zone is a thrust fault zone with at least three parallel faults striking along the SW-NE direction, in this study, we do not consider the depth distribution and just take all the shallow earthquakes for the analysis.

Method for evaluating the b-value
There have been several approaches to the calculation of bvalue, with different advantages for special needs in different works (Aki, 1965;Hamilton, 1967;Shi and Bolt, 1982;Marzocchi and Sandri, 2003;Amorèse, 2007;Sandri and Marzocchi, 2007).In this article, we use the maximum likelihood method which has shown to be a robust and unbiased estimation in most cases (Bender, 1983;Wiemer and Wyss, 1997).For an earthquake catalogue with mean magnitude M mean and cutoff magnitude M min , the maximum likelihood gives b = log 10 (e)/(M mean − M min ) (1) In the calculation one also needs to take into account the effect of the discrete binning of magnitudes on the value of M min .According to Utsu (1965), if the binning is 0.1 magnitude unit, An estimation of the standard deviation of the b-value can be obtained using the formulae of Aki (1965): We map the distribution of b-value in the studied region using a grid of 0.1 • ×0.1 • , with a circular window of radius 30 km, for the time duration of 31 years from 1977 to the day before the Wenchuan earthquake.Previous works used 5 km×5 km, r=40 km during the period 1983∼1999 ( Öncel and Wyss, 2000); 5 km×5 km, r=40 km during the period 1988∼1998 (Zuñiga and Wyss, 2001); 0.01 • ×0.01 • , r=20 km during the period 1980∼1999 (Wyss and Matsumura, 2002); 1 km×1 km, r=30 km during the period 1990∼2003 (Nakaya, 2006); or 0.025 • ×0.025 • , r=35 km during the period 1999∼2001 (Ghosh et al., 2008).Comparing to previous works, our grid is coarser, the time duration is longer, but the size of the sliding window is comparable.In the paper of Yi et al. (2006), the grid was taken as 0.1 • ×0.1 • , with sampling radius 20 km, time duration from 1977 to 2004, and cutoff magnitude M L 2.0.To highlight the Longmenshan fault zone, we also consider the polygon with dashed lines in Fig. 1, and further divide the polygon into 19 bins, similar to the mapping of Jiang and Wu (2006).Based on the results of Wyss et al. (2000), Zuñiga and Wyss (2001), and Chen et al. (2006), in the color bar representing b-value, the upper end is truncated at 1.5 and the lower end is truncated at 0.5.Following Bender (1983), if there are less than 25 events in the sliding window, the grid is discarded and not shown in the mapping due to its unacceptable uncertainty.

Relation between mainshock rupture and pre-seismic b-value distribution
Figure 3 shows the spatial distribution of b-value for earthquakes from 1977 to the day before the Wenchuan mainshock.In the same figure, the epicenter of the mainshock is also plotted as a comparison.It can be seen that the nucleation point of the mainshock locates to the north-east of the higher stress (or exactly speaking, lower b-value) region.The rupture propagated further north-eastward and penetrated another higher stress (or exactly, low b-value) region and propagated further north-eastward.To further investigate the re-lation between the mainshock rupture and the pre-seismic bvalue distribution, Fig. 4a highlights the Longmenshan fault zone as already shown by the polygon with dashed lines in Fig. 1.The completeness of earthquake catalogue in this region is shown in Fig. 2d and e.The whole polygon is divided into 19 bins.The b-value for the events in each bin is calculated, as shown in Fig. 4a, together with the slip distribution results from different sources.In the figure, we adapted the 2-dimensional slip distribution profile into a 1-D distribution of total slip and take the 70% and 50% maximum, respectively, as the first-order measure of the slip distribution or the "effective length".Hinted by previous results of temporal b-value variation before earthquakes (Imoto, 1991;Kebede and Kulhánek, 1994;Sahu and Saikia, 1994;Enescu and Ito, 2001;Cao and Gao, 2002;Ziv et al., 2003;Nuannin et al., 2005), we also try to investigate whether there is any change of b-value before the Wenchuan earthquake. Figure 4b shows the temporal variation of b-value within the polygon.A weak trend of decrease can be observed before the mainshock, albeit hard to be used as an indication of the approaching of the great earthquake, indicating why this earthquake is to much extent "unexpected".This pattern is further confirmed by considering the b-value distribution within different time periods.Figure 4c displays the results for the periods 1977∼2008 (31 years before the Wenchuan earthquake), 1993∼2008 (15 years before the earthquake), and 1998∼2008 (10 years before the quake).From the figure it can be observed that basically there is no significant variation of b-value distribution, except that near to the epicenter of the Wenchuan earthquake there was a weak trend of pre-seismic b-value decrease, or, stress increase.

Spatio-temporal variations in the b-values of aftershocks
Assessment of aftershock probabilities can be conducted using different approaches, such as the analysis of aftershock sequence (Reasenberg and Jones, 1990;Ogata and Zhuang, 2006) and the calculation of CFS change (Stein, 1999;Toda and Stein, 2000).Mapping of b-values can provide some useful clues to the cause of the aftershocks and useful for estimation of aftershock probability.Figure 5 shows the distribution of aftershocks with magnitude larger than 3.0 for the period 12 May 2008 14:28:00 to 7 September 2008 10:00:00 local time.It can be seen that aftershocks overlap with the earthquake rupture zone, showing the characteristics of predominant single-side propagation.Figure 6 shows the magnitude-index picture (Ogata et al., 1991) of the aftershock sequence.From the figure it can be seen that in the first three days, as indicated by the portion to the left of the vertical red dashed line in Fig. 6a, there was a significant missing of small events, probably due to the poor monitoring capabilities during the first three days after the mainshock when the mobile seismic stations had not been deployed in the meizoseismal region.Due to the temporal change of completeness magnitude, "dynamic" calculation of b-value (Ogata et al., 1991) is used.In the "dynamic" calculation, the whole aftershock sequence is divided into K segments, with the cutoff  magnitude in the k-th segment being M cut k , number of earthquakes in this segment above M cut k being N k , and the average magnitude for these earthquakes being M mean k .Then the "dynamic" b-value can be calculated via In our calculation we divide the segment for each 200 events along the index axis, so N k is less than 200 for each segment.Due to the natural ending of the catalogue the last segment contains 181 events.
Figure 7 shows the distribution of b-values for the period from 12 May 2008 to 7 September 2008, together with slip distribution.There is little evidence showing the termination of the main rupture.However, seen from Fig. 7, there is a region to the north with higher stress or lower b-value.It  1 are not included in the calculation of b-value.Therefore the maximum M L here is the maximum M L of earthquakes in the b-value calculation.In the transformation of magnitudes, M S 6.0 corresponds to M L 6.3.Note that based on Eq. ( 3), number of samples determines the standard derivation of the b-value.might be reasonable to connect this low b-value region with a "harder" patch on the fault which stopped the rupture propagation.
Comparing the distribution of b-values before and after the mainshock gives little evidence for the cause of aftershocks due to lack of seismicity to the north part of the earthquake fault before the mainshock.Comparing the b-values before and after each strong aftershock with magnitude larger than M S 6.0, as shown in Table 1, also fails to give sufficient information for the 1st, the 2nd, and the 8th strong aftershock.Meanwhile the 4th strong aftershock seems not associated with significant change of b-value (0.82 versus 0.85).The number of strong aftershocks being still small, it is hard to draw any definite conclusion by the changing of b-values before and after these strong aftershocks.However, considering the change of b-values together with the location of these events, it may be seen that the strong aftershocks near the front of the rupture (No. 6 to the north and No. 3 to the south) are accompanied by the decrease of stress or increase of b-value, while those within the rupture area (No. 5 and No. 7) are accompanied by the increase of stress or decrease of b-value.This pattern is to some extent understandable in physics if the termination of the rupture near to the crack tip and the existence of barriers along the rupture zone are considered.1. Table 2 shows the statistics for each bin in details.

Discussion and conclusions
One of the well-accepted theoretical consideration related to b-value variation is that (Wiemer and Wyss, 1997;Wyss, et al., 2000;Wyss, 2001) an earthquake fault consists of locked segments that resist faulting (asperities) and unlocked parts characterized by creep; In the creeping segments stresses are largely released so they can not build, whereas in asperities stresses are concentrated.Most mainshock energy radiates from asperities, while the creeping segments may participate in a main rupture by co-seismic slip; A rupture may involve only one asperity and stop as major earthquake, or it may reach one or several neighboring asperities, and continue to generate a great earthquake.Wiemer and Wyss (1997) suggested that high stress is the most likely cause for the low b-value in asperities and that the frequency-magnitude distribution may be used as a stress indicator that can identify the location of patches on the fault under high stress.The study in this paper, therefore, is two-folded.On one hand, we try to use the case of the Wenchuan earthquake to investigate whether the above hypothesis holds for this great thrust event.On the other hand, we hope that the mapping of b-values in this region for different periods may be able to reveal some of the features related to the seismogenesis and aftershock sequence of this great earthquake.Our result shows that the main rupture initiated near, but not exactly within, a higher stress region to the south of the fault (as indicated by the lower b-value), propagating north-eastward to the region of lower stress (as indicated by higher b-value), and stopped at a higher stress or harder region (as indirectly indicated by the lower b-value of aftershocks).
It has to be cautioned that, in addition to stresses, the b-value can also be affected by geological conditions (Hatzidimitriou et al., 1985;Wang, 1988;Tsapanos, 1990;Ogata et al., 1991).Since the geological setting in the study area is quite complicated (Burchfiel et al., 2008), numerous factors can also influence the spatial and temporal variations in b-value.At present the inversion results of rupture process are still inconsistent with each other.In Fig. 4a and Fig. 7 we plotted a simplified version of slip distribution results from different agencies.From the figure it can be seen that even if for this coarse/first-order slip distribution, the three results are not consistent with each other in second order details.The only pattern which can be confirmed is that the rupture propagated within the relatively low stress (high bvalue) region.Whether or not the high stress region to the south of the earthquake rupture blocked the south-westward propagation of the rupture is still a question in need of further investigation.

Fig. 1 .
Fig. 1.Epicenter map of the Longmenshan fault zone for the period from 1977 to the day before the Wenchuan M S 8.0 earthquake.The region under study is displayed in the indexing figure to the upper right.Blue dots show the locations of events with M L ≥2.4.Red circle indicates the epicenter of the M S 8.0 earthquake on 12 May 2008.Gray lines show active faults.The polygon with yellow dashed lines indicates the "highlighted" region of the Longmenshan fault zone, see Fig.4.

Fig. 2 .
Fig. 2. (a) Temporal distribution of the cumulative number of earthquakes with M L ≥2.4 from 1970 to the day before the Wenchuan M S 8.0 earthquake, in the region 29.0 • ∼34.0 • N, 101.5 • ∼107.5 • E. Vertical dashed line represents the year 1977.(b) Temporal distribution of earthquakes in the studied region since 1970.Vertical dashed line corresponds to the year 1977.(c) Frequency-magnitude distribution for the studied region since 1977, shown by cumulative distribution.Cutoff magnitude is determined in the figure as M C =2.4.Black bold line shows the fitting of Gutenberg-Richter's law.In all the figures M L is used as the unified magnitude.(d) Temporal distribution of the cumulative number of earthquakes with M L ≥2.4 from 1970 to the day before the Wenchuan M S 8.0 earthquake within the polygon in Fig. 1.The vertical dashed line represents the year 1977.(e) Frequency-magnitude distribution for earthquakes from 1977 to the day before the Wenchuan M S 8.0 earthquake, within the polygon in Fig. 1.The cutoff magnitude of completeness is M L 2.4.

Fig. 3 .
Fig. 3. Mapping of b-value in and around the Longmenshan fault zone for the period from 1977 to the day before the Wenchuan M S 8.0 earthquake.Black circle shows the epicenter of the M S 8.0 earthquake on 12 May 2008.

Fig. 4 .
Fig. 4. (a) Spatial distribution of b-value along Longmenshan fault zone.Here the fault zone is depicted by the polygon shown in Fig. 1.The whole fault zone is divided into 19 rectangles.The estimate of b-value uses earthquakes occurring in each bin during the period from 1977 to the day before the Wenchuan earthquake.Slip distribution (USGS: http://earthquake.usgs.gov/eqcenter/eqinthenews/2008/us2008ryan/finitefault.php;CALTECH: http://www.tectonics.caltech.edu/sliphistory/2008 e sichuan/e sichuan.html;TSUKUBA:http://www.geol.tsukuba.ac.jp/ ∼ nisimura/20080512) is also mapped onto the plot.See text for details.(b) Temporal variation of b-value for the events in the polygon.Sliding step is 3 months, and sliding window is one year before the time of the sliding.Error bars show the standard deviation of b-value.Vertical dashed line indicates the M S 8.0 mainshock.

Fig. 4 .
Fig. 4. (c) Spatial distribution of b-value along the Longmenshan fault zone using three different time windows preceding the occurrence of the Wenchuan M S 8.0 earthquake.

Fig. 5 .
Fig. 5. Distribution of aftershocks with M L ≥3.0 for the period from 12 May 2008 14:28:00 to 7 September 2008 10:00:00 local time.Data is from the Fast Reporting Catalogue provided by Earthquake Administration of Sichuan Province.For captions see Figs. 1 and 2. Polygon circulated by dashed lines is identical to that in Fig. 1.Aftershocks in the polygon are subjected to b-value analysis and highlighted in the figure.

Fig. 6 .
Fig. 6.(a) Magnitude versus index in the catalogue of the aftershock sequence.Vertical dashed line in red shows the time three days after the mainshock.(b) Density distribution.Black solid line shows the "dynamic" cutoff magnitude in the aftershock sequence.See text for details.

Fig. 7 .
Fig. 7. Mapping of b-value of the aftershock sequence from 12 May 2008 to 7 September 2008, comparing with slip distribution and the distribution of b-value before the mainshock.For captions see Fig. 4a.Numbered yellow dots show the location and temporal order of aftershocks with M S ≥6.0, as listed in Table1.Table2shows the statistics for each bin in details.
S Depth/km Pre-event b-value Post-event b-value

Table 2 .
Statistics of Fig. 7 in details.Pre-mainshock Post-mainshock No. of bins b-value No. of events Maximum M L b-value No. of events Maximum M * Note: * The 8 strong aftershocks listed in Table