Characterized source model of the 2013 Lushan earthquake (Mw 6.6) by the empirical Green’s function method

Ground motions near the source area of the mainshock of the 2013 Lushan earthquake (Mw 6.6) in Sichuan Province in China were reproduced using the characterized source model and the empirical Green’s function method (EGFM). The best-fit characterized source model consisted of one strong motion generation area (SMGA) and a background area. The synthesized ground motions of the characterized source model were in fairly good agreement with the observed ground motions in the frequency range from 0.5 to 30.0 Hz at ten strong motion stations. For the 2013 Lushan earthquake (Mw 6.6), both the relationships between the SMGA and the seismic moment, and those between the flat amplitude of the acceleration source spectrum in the short period and the seismic moment almost followed the empirical scaling relationships of inner fault parameters developed for crustal earthquakes. The reasons for the largest peak ground acceleration (PGA) (> 1 g) in the strong-motion observation history of China recorded at the 51BXD strong motion station were investigated from the source and site effects. We found that the directivity effect did not contribute to the largest record by comparing the effect of different positions of the rupture starting point on the synthesized ground motions. The nonlinear effect of shallow layers was negligible, as indicated by the similarity of the earthquake H/V spectral ratios between the mainshock and EGF events. A large shear-wave velocity contrast might not exist in the shallow layers as the station was situated on the slope of a small rock hill. Finally, we agreed with previous studies that the hanging-wall effect and topographic effect might be the reasons for generating the largest record at Station 51BXD.


Introduction
The 2013 Lushan earthquake (M w 6.6) struck Lushan County, Ya'an City, Sichuan Province, China at 08:02 on April 20, 2013 (Beijing time, 00:02 on April 20, 2013, UTC). The mainshock caused considerable damage to buildings; 196 people died and 21 people were reported missing. The strong ground motions during the mainshock of the 2013 Lushan earthquake were recorded by as many as 114 strong-motion stations including some near-fault stations, as shown in Fig. 1. The largest PGA was 1005 cm/s 2 calculated from the uncorrected (raw) acceleration record in the east-west (EW) component at Station 51BXD approximately 16 km from the epicenter, which was the first time that 1 g was exceeded in the strong motion observation history of China. It is important to clarify the source and site effects to interpret this large value.
To date, the source rupture process of this earthquake has been illuminated by many kinematic source models developed with the waveform inversion method using teleseismic data (Liu et al. 2013), joint inversion of strong-motion data and teleseismic data (Hao et al. 2013;Zhang et al. 2014), and joint inversion of strongmotion data, teleseismic data and GPS records (Li et al. 2017). Despite the different data sets and methods, these kinematic source models have much in common in that the mainshock ruptured a blind thrust fault in the southern part of the Longmen Shan fault belt and that large slip areas surrounded the hypocenter. The synthesized ground motions using these inverted source models agreed well with the long-period (e.g., > 1 s) observed ground motions, whereas they could not match the observed ground motions at near-source stations in the period range of engineering interests (e.g., 0.1-10.0 s). Then, the characterized source model consisting of several asperities and a background area in the seismogenic zone was applied to synthesize the ground motions. This model involves outer and inner fault parameters. The outer fault parameters are defined to outline the overall picture of the target earthquake, such as the entire source area and seismic moment, while the inner fault parameters are defined by asperities with large slips and stress drops. Strong ground motions of engineering interest are mostly generated from asperities. Therefore, the asperities are called strong motion generation areas (SMGAs).
The outer and inner fault parameters follow specific scaling relationships for crustal earthquakes (Somerville et al. 1999;Dan et al. 2001;Irikura and Miyake 2001,  . The large rectangle denotes the projection of the kinematic source model (Li et al. 2017), of which the thick line part denotes the top edge of the fault. The focal mechanisms are shown with beach balls. The star denotes the epicenter of the mainshock. Triangles and solid red circles represent the strong motion stations during the mainshock and the EGF event, respectively. Gray lines show fault traces (Deng et al. 2003) in the study area 2011; Miyake et al. 2003;Miyakoshi et al. 2020). However, the empirical scaling relationships of inner fault parameters for crustal earthquakes in China have not been adequately investigated.
In this study, we elaborate the process to construct the characterized source model of the mainshock with the empirical Green's function method (EGFM) (Irikura 1986) and then examine whether the inner fault parameters of the mainshock agree with the empirical scaling relationships thus far developed. Finally, we quantitatively discuss the source and site effects on the largest record (PGA > 1 g) at Station 51BXD.

Methodology
The empirical Green's function method (EGFM) was first proposed by Hartzell (1978) and then revised by Kanamori (1979) and Irikura (1983Irikura ( , 1986. The ground motions of a large event are expressed as a superposition of the ground motions of a small event in Eqs. (1)-(3) (Irikura 1986): Where U(t) represents the ground motions of a large event, and u(t) is the ground motions of a small event regarded as the empirical Green's functions (EGFs). The subscripts i and j are indices increasing along the strike and dip directions, respectively. The terms r, r ij , and r 0 are the respective distances from the station to the hypocenter of the small event, from the station to the (i,j) subfault, and from the station to the rupture starting point on the fault plane of the large event. In addition, ξ ij is the distance between the rupture starting point and the (i,j) subfault. V s and V r are the S-wave velocities in the source region and rupture velocity on the fault plane, respectively. F(t) is the correction function of slip velocity between large and small events which was modified by Irikura et al. (1997). The parameters, N and C, are the ratios of fault sizes (length and width) and stress drops between the large and small events, respectively. Furthermore, n' is the division number to shift the artificial periodicity to a frequency higher than that of interests, e Euler's number, and Γ the rise time of the large event. (1) (2) The formulation of the EGFM in Eq.
(1) implies that the synthesized ground motions consider the contribution of the Green's functions from the far-field term only. Because the empirical Green's functions contain the contributions of near-, intermediate-and far-field terms, the synthesized ground motions naturally contain the contributions from all terms. According to Nozu (2006), the Fourier spectral ratio of total waves corresponding to the near-, intermediate-and far-field terms to the far-field S-wave term approaches 1 when the normalized frequency (= 2πfr/V s , r is the hypocentral distance) is larger than 8. This means that the contributions of the near-and intermediate-field terms can be neglected only if f > 8V s /2πr. As V s in the source region is 3.50 km/s (Han et al. 2014), and the shortest hypocentral distance is approximately 20 km, the far-field S-wave term in the synthesized ground motions is comparable with that of the total terms without any corrections at frequencies larger than 0.2 Hz. We adopted the method proposed by Miyake et al. (1999) to estimate the fault size of the EGF event, N and C by fitting the theoretical source spectral ratio to the observed one. The observed source spectral ratio is estimated from the ratio of the observed Fourier spectrum (vector summation of two horizontal components) removing the path effect between large and small events, while the theoretical source spectral ratio is expressed as the following equation: where M o and m o are the seismic moments of large and small events, respectively, and f ca and f cm are the corner frequencies of large and small events, respectively.
Once the best-fit theoretical source spectral ratio is determined by searching the minimum of weighted least squares (Miyake et al. 2003), the corner frequencies f ca and f cm are directly obtained according to Eq. (4). The proportional parameters, N and C, are estimated from Eq. (5). The fault size r a and stress drop ∆σ a of the EGF event are estimated by applying Eq. (6) (Eshelby 1957;Brune 1970Brune , 1971): Except for the above parameters (i.e., N, C, r a ), we determine the position indices (i, j) of the rupture starting point within the SMGAs, V r and Γ, by applying the adaptive simulated annealing (ASA) algorithm, which employs a more efficient sampling of the parameter space than the conventional simulated annealing algorithm. We gave a brief introduction to the ASA algorithm by referring to Satoh (2006). Let x(t) be the vector composed of parameters x i (t) (i.e., C, N, i, j, V r and Γ). Here, x i (t) is generated from the generating function by use of random variable ∆x i as expressed in Eq. (7). ∆x i is generated from the generating temperature T i,gen and uniform random value u as expressed in Eq. (8): where T i,gen is calculated from the temperature reduction function as where k is an annealing factor for each generating temperature, c is a constant defining the relative width of the generating distributions, and α is a factor associated with the speed of temperature decrease. The values of the parameters-T i,gen (0), c, and α are set to 1.0, 1, and 0.6, respectively, based on several numerical experiments.
Once new values for x(t) are generated, the difference of the misfit functions between new and previous is computed as where E(x) is the misfit function defined by the residual between synthesized and observed records for acceleration envelopes and displacement waveforms (Miyake et al. 1999).
The new values x(t) are accepted if ∆E < 0. If ∆E > 0, x(t) is accepted if the uniform random value u is less than P, which is expressed as the following equation: where the acceptance temperature is calculated as Eq. (9).

Determination of parameters required for the EGFM
The characterized source model in this study is composed of a SMGA consisting of N × N subfaults along the strike and dip directions, and a background area in the seismogenic zone on the fault plane. Because the background area mostly contributes to long-period strong ground motions (Kurahashi and Irikura 2010), we estimated the strong ground motions only from the contribution of the SMGA neglecting the contributions from the background area. The number of SMGAs generally depends on the complexity of the source rupture process. For large magnitude earthquakes, multiple SMGAs seem to be distributed randomly (e. g. Kamae and Irikura 1998;Irikura 2010, 2013). For the 2013 Lushan earthquake, several source-rupture-process studies based on waveform inversion analyses were common in that the large slip area was concentrated around the hypocenter. Therefore, it is reasonable to suppose that one SMGA located in the large slip area of the kinematic source models is suitable for synthesizing the ground motions.  Table 1. The difference in the dip and rake angles between the mainshock and the aftershock are influential on the radiation pattern of ground motions from the source over a longer period range but less influential on those over a shorter period range, as the radiation characteristics over a shorter period are naturally smoothed due to diffraction and scattering along propagation paths (Kamae and Irikura 1992). Moreover, the site effects are more influential on the mainshock ground motions in Table 1 Information on the mainshock (Li et al. 2017) and the EGF event a refers to Han et al. (2014) b refers to Yi et al. (2016) Figure 1 also shows the projection of the kinematic source model with a large rectangle. The kinematic source model developed by the waveform inversion method (Li et al. 2017) is adopted for constructing the characterized source model, as the hypocenter is close to the position determined by highresolution relocation (Han et al. 2014). The ground motions of the M w 4.3 aftershock at the stations, where the records of both the mainshock and the EGF event were available are regarded as EGFs. We synthesized ground motions at seven strong motion stations (e.g., 51BXD, 51BXM, 51BXY, 51YAM, 51QLY, 51YAD, 51PJD) surrounding the source region by applying the EGFM. The 51BXZ near-source station is excluded, as a strong nonlinear effect is observed (see details in the "Discussion" section). Second, we estimated the inner fault parameters related to the SMGA from the source spectral ratios of the mainshock ground motions to the EGF ones. We employed the S-wave portion with a length of 40.96 s of the mainshock and EGF event to calculate the observed Fourier spectra. We followed the method of Miyake et al. (1999) to adopt the hypocenter distance for the mainshock and EGF events to remove the path effect (Q s (f) = 155f 0.6804 , Tao et al. 2016) from the Fourier spectra to obtain the observed source spectra. The average source spectral ratios and standard deviations for the target stations are shown with thick gray and dashed curves in Fig. 2. Three parameters, i.e., seismic moment ratio (M o /m o ) between the mainshock and the EGF event, corner frequencies of the mainshock (f cm ) and the EGF event (f ca ), are required to estimate the theoretical source spectral ratio as expressed in Eq. (4). The ratio M o /m o is known from the moment magnitudes of the mainshock and the EGF event, while f cm and f ca are determined when the theoretical spectral ratio best fits the observed source spectral ratio in a specific frequency range. We set the lower frequency limit for fitting the observed source spectral ratio to be 0.5 Hz considering the low signal-to-noise ratio of observed ground motions and the upper frequency limit to be 30.0 Hz, which is below the cutoff frequency of the anti-aliasing filter of the seismograph. We determined the best fitting of the theoretical source spectral ratio (red curves in Fig. 2) to the average observed spectral ratio (thick curves in Fig. 2) in the frequency range of 0.5 to 30.0 Hz when f ca = 1.90 Hz and f cm = 0.17 Hz. Subsequently, we applied Eqs. (5) and (6) to calculate N = 11, C = 2.12, r a = 0.68 km, and ∆σ a = 4.85 MPa. For small or moderate earthquakes, the fault length can be identical to the fault width; thus, the fault length or width of the EGF event is √ πr a =1.21 km. It should be noted that the ground motions from the SMGA are approximately expressed as omega-squared spectra with a corner frequency f cSMGA (Boatwright 1988;Miyake et al. 2003). Then, f cm and M o in Eq. (5) can be expressed with f cSMGA and M oSMGA (seismic moment of the SMGA), respectively (Miyake et al. 2003;Somei et al. 2020). Therefore, N and C can be regarded as the ratio of the fault size between the SMGA and the EGF event, and the ratio of the stress drop between the SMGA and the EGF event, respectively.
Third, we applied the ASA algorithm to determine the optimal values of the parameters required by the EGFM, such as i, j, V r , Γ as well as N and C determined by fitting the observed source spectral ratio which are not sufficiently accuracy considering the small signal-to-noise ratio of the observed ground motions in the low frequency range. We fixed the rupture starting point at the hypocenter (Li et al. 2017) and kept the size of the fault length/width of the EGF event constant. Under the same acceptance temperature, the iteration of the ASA algorithm was carried out 5 times, and the search of indices (i, j) of the rupture starting point was carried out 10 times for each iteration, whereas the final value of k was set to 100. Therefore, the optimal parameters were determined when the misfit function reached the global minimum among the 5000 iterations.  Fig. 2 Fitting of the average observed source spectral ratio (thick curve). Gray curves show the observed source spectral ratios at the strong motion stations. Dashed curves show the standard deviation above and beneath the average observed source spectral ratios. The red curve shows the best-fitting theoretical source spectral ratio in the frequency range of 0.50 ~ 30.0 Hz range for each parameter in the second column. The search ranges for V r and rise time-τ a of the EGF event (Γ = Nτ a for the mainshock) are approximately assigned, while the upper and lower limits for C and N are evaluated by the fitting of the theoretical source spectral ratio to the observed ones as shown in dashed curves lying one standard deviation above and below the average observed source spectral ratio. In accordance with the fitting frequency range of the observed source spectral ratio, we applied a bandpass filter of 0.5-30.0 Hz beforehand to the waveforms with a time window of 30 s, which contains the S-wave portions at the seven target stations. Because step-like pulses were observed on the baseline of the records in the EW component at Station 51QLY for both the EGF event and the mainshock, we abandoned synthesizing the ground motions in the EW component to eliminate the special process of baseline correction.

Synthesized ground motions using the characterized source model
We obtained the best-fit SMGA of the characterized source model when the residuals between synthesized and observed records for the misfit function consisting of acceleration envelopes and displacement waveforms (Miyake et al. 1999) reached the global minimum throughout all iterations. The optimal values of the parameters, such as C, N, i, j, V r and Γ, are listed in the third column in Table 2. The SMGA is 118.6 km 2 (10.89 × 10.89 km), the stress drop is 12.9 MPa, and the seismic moment of SMGA is 6.85 × 10 18 N·m. The location of the SMGA, as shown in Fig. 3, nearly coincides with the large slip area of the inverted kinematic source model (Li et al. 2017), which suggests that the SMGA can be regarded as equivalent to an asperity. Figure 4 shows the comparisons of waveforms between synthesized and observed acceleration, velocity and displacement in the NS and EW components at seven target stations. Figure 5 shows the comparisons of the acceleration Fourier spectra between synthesized and observed ground motions.
The Fourier spectra are smoothed using a logarithmic window function with a band width parameter of b = 50 (Konno and Omachi 1998). Overall, not only the waveforms but also the Fourier spectra of the synthesized ground motions are in fairly good agreement with the counterparts of the observed ground motions at the target stations except for some stations, where synthesized acceleration waveforms are overestimated as nonlinearity of surface geology is not included in the EGFM. We also applied the same set of parameters of the best-fit SMGA to synthesize the ground motions at three nontarget stations, i.e., 51PJW, 51DJZ, and 51CDZ, where the observed ground motions are not considered in the above inversion analysis. The good consistency between the synthesized and observed ground motions is confirmed (Additional file 1: Fig. S1 and Additional file 2: Fig.  S2), which further suggests that the characterized source model consisting of one SMGA is suitable for synthesizing the ground motions of the mainshock. Subsequently, we examined the effectiveness of the inner fault parameters, such as the SMGA and acceleration level, in the best-fit source model. Here, the SMGA is considered equivalent to the combined area of an asperity. Figure 6 shows the empirical scaling relationship between the combined area of asperities with large slips estimated from waveform inversions for long-period (> 1.0 s) motion and the seismic moment (M o ) proposed by Somerville et al. (1999) and Miyake (2001, 2011) in the left panel, and the empirical scaling relationship between the flat amplitude of the acceleration source   (Wang et al. 2017) and the 2017 Jiuzhaigou earthquake (M w 6.5) (Wu et al. 2020), as shown with square, diamond and solid circle, respectively, in Fig. 6. It should be noted that Kurahashi and Irikura (2010) constructed the characterized source model only in the southern part of the fault model, and the seismic moment in the southern part was assumed to be half of the seismic moment for the entire fault. The relationships between the SMGAs and M o for these crustal earthquakes in China agree well with the empirical scaling relationship (Somerville et al. 1999;Miyake 2001, 2011;Miyake et al. 2003;Miyakoshi et al. 2020). The A 0 values for the 2014 Ludian earthquake and the 2008 Wenchuan earthquake are slightly larger than those estimated from the empirical scaling relation, which is acceptable, as the values fall in the range indicating factors of 1/2 and 2 for the average (shown in dotted lines). Therefore, the relationships between A 0 and M o for these crustal earthquakes also follow an empirical scaling relationship (Dan et al. 2001

Discussion
The EGFM is applicable for synthesizing the ground motions on the premise of linear superposition of the EGF, while the nonlinear effects on the ground motions are excluded for applying the EGFM. To examine the influence of nonlinearity on the surface geology of observed ground motions, we compare the earthquake H/V spectral ratio of the mainshock with that of the EGF event at seven target stations in Fig. 7. The acceleration Fourier spectra in three components are calculated from the S-wave portions with a time length of 20.48 s, and then the H/V spectral ratios are evaluated following the method of Kawase et al. (2011). The degree of nonlinearity (DNL) (Noguchi and Sasatani 2008) at each station is evaluated to represent the extent of the nonlinear effect. A large DNL suggests a strong nonlinear Mo(N.m) Somerville, et al., 1999Miyake, et al., 2003Miyakoshi, et al., 2019 This study Kurahashi and Irikura, 2010Wang, et al., 2017Wu, et al., 2020 5.0 5 .5 6.0 6 .5 7.0 7 .5 8. Mo(N.m) Dan et al., 2001 This study Kurahashi and Irikura, 2010Wang et al., 2017Wu et al., 2020 Fig. 6 Empirical scaling relationships of inner fault parameters between the combined area of asperities and the seismic moment (left), and between the flat amplitude of the acceleration source spectrum in the short period and the seismic moment (right). The thick lines represent the regression equations derived by Somerville et al. (1999) in the left panel and by Dan et al. (2001) in the right panel. The dashed lines represent extrapolation of regression equations. The dotted lines denote the factors of 1/2 and 0.5 for the average effect, and vice versa. The similarity of the earthquake H/V spectral ratios between the mainshock and the EGF event and the small values of DNLs indicate that nonlinear effects are sufficiently weak at some target stations.
Although the observed PGA of the mainshock is quite large at Station 51BXD, the smallest DNL value suggests that the nonlinear behavior of the site is negligible there. This is consistent with the fact that the site condition is rock, and the field survey also reported that the station was situated on the slope of a small rock hill (Wen and Ren 2014). The moderate values of DNL at Stations 51BXY and 51QLY suggest that the nonlinear behaviors of the sites are moderate. The synthesized acceleration ground motions are slightly larger than the observed ground motions at these stations due to the linear premise of the EGFM. An exception is the 51BXZ station, where the largest value of DNL suggests that nonlinear behavior is quite strong, which should be excluded when synthesizing the ground motions of the mainshock. The synthesis through the EGFM taking the nonlinear effect into account could be accomplished by introducing certain parameters (Nozu and Morikawa 2003), which is beyond the scope of this study. Figure 8 compares the observed PGAs at target stations with those predicted by the ground motion prediction equations (GMPEs) (Zhao et al. 2016). The site conditions of the target stations can be categorized  Black lines represent the waveforms of the observed ground motions. Red, blue, and purple lines represent the waveforms of the synthesized ground motions for the best-fit model, Model A and Model B, respectively. Acceleration Fourier spectra of the observed ground motions and the synthesized ground motions of three models are shown in the right panel. The Fourier spectra are smoothed by use of a logarithmic window function (Konno and Omachi 1998) forward direction of rupture propagation (Model A). Moreover, the synthesized waveforms of Model A are in better agreement with the observed waveforms than those of the best-fit characterized source model. However, if we compare the synthesized ground motions with the observed ones at other stations, e.g., 51YAM, as shown in Fig. 11, the synthesized ground motions of Model A are in poorer agreement with the observed ones than those of the best-fit SMGA model. In terms of the fitting of the synthesized ground motions at all target stations, the best-fit SMGA source model is preferred among the three models. Furthermore, the amplitude of the acceleration Fourier spectrum of Model B is quite larger at periods of 1.0 ~ 2.0 s than that of the other two models at Station 51YAM. This indicates that the directivity effect is noticeable at Station 51YAM, which is in the forward direction of rupture propagation (Model B). In contrast, the directivity effect is not significant at Station 51BXD even though the station is in the forward direction of rupture propagation (Model A). This is because the SMGA dips toward the northwest and the distances between each subfault of the SMGA and the 51BXD station are not as variant when the rupture propagates toward the 51BXD station. We concluded that the directivity effect could not be responsible for the largest PGA in the EW component at Station 51BXD station for the best-fit SMGA source model of which the rupture starting point is close to the center of the SMGA. In addition, the synthesized ground motions for the best-fit SMGA of the characterized source model agree with the observed ground motions in the EW component at Station 51BXD during the mainshock, as shown in Fig. 4 and Fig. 5. The largest PGA can be attributed to the large peak of the Fourier spectrum at approximately 0.1 s which was observed during the mainshock and the EGF event, as shown in Fig. 12. On the other hand, if strong site amplification contributes to the largest record, it must be caused by the large shear-wave velocity contrast between the soft and hard soils in the shallow layers. However, soft soil might not exist, as the station is Observed and synthesized acceleration Fourier spectra in the EW component at Station 51BXD. The black, red, and thin lines denote the acceleration Fourier spectra for the observed and synthesized ground motions of the best-fit characterized source model, and the observed ground motions of the EGF event, respectively. The Fourier spectra are smoothed by use of a logarithmic window function (Konno and Omachi 1998) situated on the slope of a small rock hill. It is appropriate to attribute the reason to the hanging-wall effect and topographic effect which have been investigated in the previous studies (Dai and Li 2013;Xie et al. 2014;Wen and Ren 2014).