Strong motion PGA prediction for southwestern China from small earthquake records

For regions without enough strong ground motion records, a seismology-based method is adopted to predict motion PGA (peak ground acceleration) values on rock sites with parameters from small earthquake data, recorded by regional broadband digital monitoring networks. Sichuan and Yunnan regions in southwestern China are selected for this case study. Five regional parameters of source spectrum and attenuation are acquired from a joint inversion by the micro-genetic algorithm. PGAs are predicted for earthquakes with moment magnitude (Mw) 5.0, 6.0, and 7.0 respectively and a series of distances. The result is compared with limited regional strong motion data in the corresponding interval Mw± 0.5. Most of the results ideally pass through the data clusters, except the case ofMw7.0 in the Sichuan region, which shows an obvious slow attenuation due to a lack of observed data from larger earthquakes (Mw ≥ 7.0). For further application, the parameters are adopted in strong motion synthesis at two near-fault stations during the great Wenchuan Earthquake M8.0 in 2008.


Introduction
Ground motion prediction equation (GMPE) is a vital field in the research of engineering seismology.A great number of research results have been published, and most of them are empirical, such as those for the western United States by the Next Generation Attenuation (NGA) project (Power et al., 2008;Bozorgnia et al., 2012;Boore et al., 2014) and those for Japan (Si and Midorikawa, 2000;Kanno, 2006).Empirical GMPEs are developed, mainly based on plenty of strong ground motion records.For the regions without enough observed data, like the mainland of China and many others, it is difficult to build GMPEs this way.Obviously, there are more records of small earthquakes from regional earthquake monitoring networks than those of strong motion in general.The main goal of this paper is to see if these data could be a basis to estimate the regional parameters in a seismology-based attenuation relationship (Hanks and McGuire, 1981;Boore, 1983;Atkinson, 1984;Boore and Atkinson, 1987).If they are really helpful, the bottleneck restriction of insufficient strong motion records will be broken completely.Five regional parameters, stress drop σ , Q 0 and η in a quality factor, and R 1 and R 2 in a geometric attenuation factor, are acquired from a joint inversion by the micro-genetic algorithm.The objective function is defined as the summation of square differences between envelopes of calculated and observed Fourier velocity amplitude spectra.The result is compared with the limited observed PGA (peak ground acceleration) data that have never been used in the inversion, to show the feasibility of this method.The regional parameters are further adopted in synthesis of strong motions at two near-fault rock stations with a hybrid source model for magnitude 8.0.The result is compared with the records during the great Wenchuan Earthquake in 2008.

Data for inversion of the five parameters
Small earthquake data recorded from January 2001 to December 2007 by Sichuan and Yunnan digital earthquake net- works are adopted for the inversion.The period is from the beginning of operations of both networks to the end of the year before the Wenchuan Earthquake.There are 29 digital seismological stations in Sichuan and 26 stations in Yunnan; their locations are shown as triangles in Fig. 1.
A total of 147 records from 82 small earthquakes in Sichuan and 863 records from 154 small earthquakes in Yunnan, with M w 3.5 ∼ 4.5 and focal depth ≤ 30 km, are downloaded from the China Earthquake Networks Center (2011) and processed using a conventional procedure.The epicenters are shown as circles in Fig. 1.
The distributions of focal depth and hypocentral distance with M w are shown in Figs. 2 and 3, respectively.The hypocentral distances distribute more or less uniformly from 50 to 300 km in the regions.
The time window of the velocity time history, with time interval 0.02 s, starts from S-wave arrives until that 80 % energy is included.

Methodology
Fourier amplitude spectrum of ground motion from a point source of a given earthquake can be presented as (Boore, 1983;Atkinson and Boore, 2000;Boore, 2003) where C is a constant, as Eq. ( 2); M 0 is the corresponding seismic moment; f is the frequency; R is the distance; S( q ) is the source spectrum; G( q ) is the geometric attenuation term; D( q ) is the medium energy dissipation term; A( q ) is the near-surface transfer function with V 30 620 m s −1 for generic rock (Boore and Joyner, 1997); P ( q ) is the high-frequency filter with f max = 5 Hz from regional strong ground motion records; I ( q ) is the motion type conversion function for displacement, velocity and acceleration.The constant, proportion factor, is given by where R θφ is the radiation pattern of the shear excitation, with the assumption that the energy was equally partitioned into two horizontal components, R θφ = 0.6; F is the free surface effect, F = 2.0; V is the vectorial partitioning of shear-wave energy into two components of equal amplitude,  and Hanks, 1980); R 0 is a reference distance, R 0 = 1 km (Boore, 2003); ρ s is the density in the vicinity of the source, ρ s = 2.8 g cm −3 ; β s is the shear-wave velocity in the vicinity of the source, β s = 3.5 km s −1 (Wu et al., 2001;Xie et al., 2012).
In order to eliminate the dependence of the motion amplitude on source size, the following source spectrum model (Tao et al., 2008) is adopted rather than the well-known ω 2 spectrum (Brune, 1970) or the two corner frequency spectrum (Atkinson, 1993).
where f c is the corner frequency, estimated by Eq. ( 4) (Boore, 2003); parameter a and b could be taken as a where σ is the dynamic stress drop.
The geometric attenuation term G(R) is commonly adopted as in Eq. ( 5), for the distance ranges predominated by the shear wave, the surface wave, and by the wave guide effect in between these two (Atkinson and Mereu, 1992; Atkinson and Boore, 1995).
where R 1 and R 2 are regional geometric attenuation parameters, and are concerned with the regional structure of the earth's crust and dynamic characteristics of the crust medium.
The energy dissipation along the wave propagation path is very complicated; the medium energy dissipation term D(R, f ) can be taken as (Boore, 1983) where Q(f ) is the quality factor, and can be described as (Beresnev, 1998) where Q 0 and η are the region-related coefficient and the exponent, respectively.The above-mentioned five regional parameters, σ , R 1 , R 2 , Q 0 , and η, could not be measured directly.In general, σ , Q 0 , and η obtain values from inversion of regional far field records; R 1 and R 2 are estimated from the regional earth crust structure.In practice, there are several different values for the former three parameters from different authors with different data, even different forms of the formula in some regions, but no values are published in other regions.It is also not easy to estimate the values of the latter two in some regions with complex topography, such as in Sichuan and Yunnan.Therefore, the authors of this paper suggested a joint inversion procedure by the micro-genetic algorithm, µGA  (Carroll, 2009), to see if the records of small earthquakes are helpful.In nature, it is a search scheme.The five values are searched from the pre-chosen ranges that are large enough to cover all possible values (Atkinson and Boore, 1995;Qiao et al., 2006;Mao et al., 2005;Zhang et al., 2007;Qin and Kan, 1986;Su et al., 2006;Liu, 2005;Hu et al., 2003), by minimizing the objective function defined in Eq. ( 8).
where FA 0 is the ith amplitude of enveloped Fourier amplitude spectrum of the kth record; FA j (i, k) is the corresponding amplitude by Eq. ( 1) and the j th generation parameter values; 4096 is the total amplitude number of the Fourier spectrum from the time step 0.02 s and the corresponding frequency spacing 0.0122 Hz; N is the total number of the records.
The five parameter values of the solution, after having searched 2000 generations, are listed in Table 1.
For each given pair of magnitude and distance, the Fourier amplitude spectrum by Eq. ( 1) with the parameter values in Table 1, combined with a random phase spectrum, is transformed into a time domain.The time history is then enveloped for the amplitude nonstationarity.Then, the enveloped time history is transformed back into frequency domain, and scaled to the amplitudes directly from Eq. (1).The scaled spectrum is transformed back to the time domain again, and PGA for this magnitude distance is read from this final time history.Fifty acceleration time histories are generated and the average PGA is the stable expected PGA.

Strong ground motion prediction in Sichuan and Yunnan regions
It is a kind of check, to compare the result with regional limited strong ground motion data in Sect.2.2.These data are never adopted in the inversion.The attenuation curves for M w 5.0, 6.0, and 7.0, with data points in the intervals M w ± 0.5, are shown in Fig. 6.In the cases of M w 5.0 and 6.0, the result curves pass through the corresponding strong motion data, and the data points are distributed evenly on   both sides of the curve, which means our results present the mean ground motion levels well.In the case of M w 7.0, our result for the Sichuan region is obviously higher than the observed data, since we do not have data showing magnitudes of more than M w 6.5 there.
To get a quantitative measure, the residual between the observed value and the corresponding predicted value is calculated by where PGA is the observed value; PGA' is the corresponding predicted value.Distributions of ε for the cases M w 5.0, 6.0,The mean value and standard deviation for all distances but the three magnitude intervals are listed in Table 2.The statistic features are stable as a whole and unusual in just a few intervals.It means that there are really not enough observed data in the interval.The values of standard deviation are comparable with those of empirical prediction in the regions with rich strong motion data.The facts show that it is possible to build an attenuation relation by small earthquake data recorded by regional broadband earthquake networks, for regions without enough strong ground motion records.

Application of the regional parameters to predict near-fault strong motion
Near-fault ground motion of a large earthquake is complicated and governed by the source mechanism; therefore, it is not possible to be predicted by GMPEs, but possible by  Thirty hybrid source models are built (Liu and Tao, 2013); 10 of them are shown in Fig. 8.The motion from each subsource is generated in the same way as mentioned above.The motion from the all sub-sources on rupture plane are superposed by the time histories from every sub-source in the time domain, taking into account the time difference between the rupture in the sub-source and the initial rupture and the time difference between wave arrival times at the ground point from the sub-source and from the rupture start point.
For each station, 30 acceleration time histories are synthesized from these 30 models; the mean values of PGAs are listed in Table 3.The PGAs vary around the mean values with some randomness, and those values are close to the observed PGAs.Response spectra of the 30 motions at each station are shown as gray lines in Fig. 9, in which the black dotted line shows the mean spectrum at each station, and the black solid line shows the spectrum closest to the mean from source model 4. The acceleration time histories from this source model are chosen as samples for the two stations, shown as the top couple in Fig. 10.The two couples below are E-W and N-S components of the observed time histories at the stations.The amplitudes are similar, but the durations are different obviously.For Maoxian station, the second and the third sub-shocks are not synthesized.
The response spectra of motions at the two stations from model 4 are shown with those of the observed time histories in Fig. 11, as the black solid line and dotted lines, respectively.The spectra are close to each other at Maoxian station, but they are different between a period of 0.3 and 0.4 s and longer than 0.7 s at Pixian station.

Conclusions
An approach to predict strong ground motion PGA in southwestern China is presented in this paper, by means of a seismology-based attenuation relationship with five regional parameters: σ , Q 0 , η, R 1 , and R 2 .The values of those parameters are estimated from a joint inversion of small earthquake data recorded by regional broadband digital earthquake networks, rather than values from some available references.The objective function of the inversion by µGA is constructed by minimizing the differences between envelopes of Fourier spectra of the small earthquake records and those calculated for the corresponding magnitude and distance of the quakes.The PGA results for M w 5.0 and 6.0 pass through clusters of the limited strong motion data observed in the regions as a whole.The result for Sichuan M w 7.0 is higher than the observed data, since there are no strong motion data in the range of M w > 6.5.From this kind of blind test, since in the process of inversion, strong ground motion data are never adopted, the feasibility of this approach is shown for regions without enough strong ground motion data.Furthermore, the regional five parameters are adopted to predict near-fault strong ground motion of a large earthquake M8.0 by motion synthesis, though for the case study at two rock stations in Sichuan, the motion is recorded during the Wenchuan Earthquake.For each station, 30 acceleration time histories are synthesized from 30 finite fault hybrid source models.The results show that PGA can be predicted very well by the mean and that the average spectra are good.The amplitudes of time histories are also similar to the records, but their envelopes are not in accordance with those observed.

Figure 1 .
Figure 1.Seismological stations and small earthquakes in Sichuan and Yunnan regions.
There are 69 strong ground motion stations in Sichuan and 36 stations in Yunnan.A total of 1234 records from 118 Sichuan earthquakes and 78 records from 27 Yunnan earthquakes, with M w ≥ 4.5 and focal depth ≤ 30 km, are adopted.The processed records are directly downloaded from the China Strong Motion Networks Center.The strong ground motion data of Sichuan are all from aftershocks of the Wenchuan Earthquake.The distributions of M w focal depth and M w hypocentral distance are shown in Figs. 4 and 5.The focal depth is from 8 to 20 km in Sichuan, and most focal depths of big earthquakes in Yunnan are less than 15 km.Hypocentral distances in Sichuan are from 50 to 300 km, and most distances in Yunnan are less than 120 km; some of them are less than 20 km.

Figure 2 .
Figure 2. M w focal depth distribution of small earthquake records in (a) Sichuan region and (b) Yunnan region.

Figure 3 .
Figure 3. M w hypocentral distance distribution of small earthquake records in (a) Sichuan region and (b) Yunnan region.

Figure 4 .
Figure 4. M w focal depth distribution of strong ground motion records in (a) the Sichuan region and (b) the Yunnan region.

Figure 5 .
Figure 5. M w hypocentral distance distribution of strong ground motion records in (a) the Sichuan region and (b) the Yunnan region.

Figure 6 .
Figure 6.Comparison of the result with strong ground motion data in the Sichuan region and the Yunnan region.

Figure 7 .
Figure 7. Residual distribution in the Sichuan region and the Yunnan region.

Figure 8 .
Figure 8.Ten hybrid source models in the case of Wenchuan Earthquake, redrawn from Liu and Tao (2013).
and 7.0 are shown in Fig.7.The two figures, at the bottom in Fig.7, show the ε distribution for all magnitude range, within bins equally spaced with respect to log 10 R.In the cases of M w 5.0 and 6.0, ε values are evenly distributed around the mean of 0.0 which means the ratios of the observed values and the predicted values are around 1. Most residuals of M w 7.0 in the Sichuan region are lower than 0.0 because of the same reason as mentioned above.No perceptible trend indicates that the path terms could represent the data trends reasonably, which is the same as the conclusion inBoore et al. (2014).

Table 2 .
Mean value and standard deviation of residuals.Maoxian station (MXT, 31.7 • N, 103.9 • E) and Pixian station (PXZ, 30.9 • N, 103.7 • E), shown as squares in Fig.1, in the areas with intensity IX and VIII, respectively, during the M8.0 Wenchuan Earthquake, are selected as examples, since strong ground motion was recorded there dur- M w = 5.0 M w = 6.0 M w = 7.0

Table 3 .
Synthesized PGA and observed PGA at two rock stations (gal).