Source parameters, path attenuation and site effects from strong-motion recordings of the Wenchuan aftershocks (2008–2013) using a non-parametric generalized inversion technique

Secondary ( S ) wave amplitude spectra from 928 strong-motion recordings were collected to determine the source spectra, path attenuation and site responses using a non-parametric generalized inversion technique. The data sets were recorded at 43 permanent and temporary strong-motion stations in 132 earthquakes of M s 3.2–6.5 from 2008 May 12 to 2013 December 31 occurring on or near the fault plane of the 2008 Wenchuan earthquake. Some source parameters were determined using the grid-searching method based on the omega-square model. The seismic moment and corner frequency vary from 2.0 × 10 14 to 1.7 × 10 18 N · m and from 0.1 to 3.1 Hz, respectively. The S -wave energy-to-moment ratio is approximately 1.32 × 10 − 5 . It shows that the moment magnitude is systematically lower than the surface wave magnitude or local magnitude measured by the China Earthquake Network Center. The seismic moment is approximately inversely proportional to the cube of the corner frequency. The stress drop values mainly range from 0.1 to 1.0 MPa, and are lognormal distributed with a logarithmic mean of 0.52 MPa, signiﬁcantly lower than the average level over global earthquake catalogues. The stress drop does not show signiﬁcant dependence on the earthquake size and hypocentre depth, which implies self-similarity for earthquakes in this study. The ε indicator was used to determine the stress drop mechanism. The low stress drop characteristic of the Wenchuan aftershocks may be interpreted by the partial stress drop mechanism, which may result from remaining locked sections on the fault plane of the main shock. Furthermore, we compared the stress drop distribution of aftershocks and slip distribution on the fault plane of the main shock. We found that aftershocks with higher stress drop occurred at areas with smaller slip in the main shock. The inverted path attenuation shows that the geometrical spreading around the seismogenic region of Wenchuan earthquake sequence is weak and signiﬁcantly dependent on frequency for hypocentre distances ranging from 30 to 150 km. The frequency-dependent S -wave quality factor was regressed to 151.2 f 1.06 at frequencies ranging from 0.1 to 20 Hz. The inverted site responses provide reliable results for most stations. The site responses are obviously different at stations in a terrain array, higher at the hilltop and lower at the hillfoot, indicating that ground motion is signiﬁcantly affected by local topography.


D ATA S E T
A total of seven issues (Issues 12-18) of uncorrected strong-motion acceleration recordings have been officially issued in China by the China Strong Motion Network Center since the China National Strong Motion Observation Network System (NSMONS) formally began operation in 2007. The analogue recordings in China before 2007 were published in Issues 1-11. Issue 12 covers recordings from the Wenchuan main shock. Issues 13 and 14 cover recordings from the Wenchuan aftershocks obtained by the permanent and temporary stations, respectively. Issues 15, 16 and 18 cover other recordings collected during 2007-2009, 2010-2011 and 2012-2013, respectively. Issue 17 covers recordings from the 2013 Lushan earthquake sequence. Strong-motion recordings in Issues 13-16 and 18 from earthquakes that occurred on or near the rupture fault of the Wenchuan earthquake are used as the data set in this study. It is composed of more than 2000 strong-motion recordings from 383 M s (M L ) 3.3-6.5 earthquakes from 2008 May 12 to September 30 recorded at 76 permanent stations of NSMONS in Gansu and Sichuan provinces (Issue 13), 2214 strong-motion recordings from 600 M s (M L ) 2.3-6.3 earthquakes from 2008 May 14 to October 10 recorded at 83 temporary stations (Issue 14, Wen et al. 2014), and 355 additional strong-motion recordings from 57 M s (M L ) 3.1-5.5 earthquakes from 2008 October 1 to 2013 December 31 recorded at 86 stations of NSMONS (Issues 15,16 and 18). Deviations between the surface wave magnitude M s and the local magnitude M L for earthquakes in China and adjacent regions measured by the China Earthquake Network Center (CENC) were ignored . In this paper, M s was used to represent the measured magnitude by CENC. The baseline correction and a Butterworth bandpass filter between 0.1 and 30.0 Hz were performed. Fig. 1 shows the hypocentre distance (R) and geometric mean of the peak ground acceleration (PGA) for the two horizontal components (east-west and north-south) of the strong-motion recordings in this data set. PGAs of these strong-motion recordings mainly vary from 2.0 to 100 cm s −2 . Hypocentre distances of most recordings from Issues 13, 15, 16 and 18 generally range from 30 to 200 km. However, hypocentre distances for many recordings from Issue 14 are less than 30 km, with the minimum approaching 1.0 km. Recordings from Issue 14 were obtained by temporary stations deployed as close to the seismogenic fault as possible (Wen et al. 2014). Very few recordings from Issues 15, 16 and 18 were obtained from earthquakes of M s > 5.0 because very few aftershocks with large magnitudes occurred in the seismogenic area of the Wenchuan earthquake during 2009-2013. We selected available recordings from this data set according to the following criteria proposed by Ren et al. (2013): (1) 30 km ≤ R ≤ 150 km; (2) 2 cm s −2 ≤ PGA ≤ 100 cm s −2 ; (3) each selected earthquake should be recorded by at least four stations, each of which should collect at least four recordings that match (1) and (2). Finally, we employed 928 strong-motion recordings from 132 earthquakes of M s 3.2-6.5 at 43 strong-motion stations. Earthquake epicentres and strong-motion stations considered in this study are shown in Fig. 2. Most stations are located in the mountains, west of the Longmenshan fault belt. In contrast to data used by Ren et al. (2013), we added more earthquakes and strong-motion stations in this study. The S waves of the two horizontal components of the strong-motion recordings were extracted, according to studies of Husid (1967) andMcCann (1979). A cosine taper was applied at the beginning and end of the S-wave window, and the length of each taper was set at 10 per cent of the total trace length (Hassani et al. 2011;Ren et al. 2013). The Fourier amplitude spectrum of the S wave was calculated and smoothed using the windowing function of Konno & Ohmachi (1998) with b = 20. The vector synthesis of the Fourier amplitude spectra from two horizontal components was used to represent the horizontal ground motion in frequency domain.

M E T H O D O L O G Y
We applied a two-step non-parametric GIT (Castro et al. 1990;Oth et al. 2008Oth et al. , 2009 to separate attenuation characteristics, source spectra and site response functions. In the first step, the dependence of the spectral amplitudes on the distance at frequency (f) can be expressed as: where O ij (f,R ij ) is the spectral amplitude observed at the jth station resulting from the ith earthquake, R ij is the hypocentre distance, M i (f) is a scale dependent on the size of the ith earthquake and A(f, R ij ) is a non-parametric function of distance and frequency accounting for the seismic attenuation (e.g. geometrical spreading, anelastic and scattering attenuation, refracted arrivals, etc.) along the path from source to site. A(f, R ij ) is not supposed to have any parametric functional form and is constrained to be a smooth function of distance with a value of 1 at reference distance R 0 .
Once A(f, R ij ) is determined, the spectral amplitudes can be corrected for the seismic attenuation effect. In the second step, the corrected spectra are divided into source spectra and site response functions: where G j (f) is the site response function at the jth station and S i (f) is the source spectrum of the ith earthquake. The trade-off between the site and source is resolved by selecting station 62WIX as a reference site, where the site responses are constrained to be 2.0 around all frequencies (Ren et al. 2013). . The inverted source spectra for four typical earthquakes representing four magnitude levels. The dark lines represent the inverted source spectra using the total recordings in this study. The grey lines represent the inverted source spectra from 100 bootstrap inversions. The name of the earthquake is composed of the date and time of this event. Eq.
(1) can be turned into a linear problem by taking the natural logarithm and expressing it as a matrix formulation: In eq. (3), the hypocentre distance ranges are divided into N bins with a 5 km width. R 1 , R 2 . . . , R N is a monotonically increasing sequence of hypocentre distance. The weighting factor ω 1 is used to constrain A(f, R 0 ) = 1 at reference distance R 0 and ω 2 is the factor determining the degree of smoothness of the solution. The reference distance was set to 30 km, which is the smallest hypocentre distance considered in this study.
We calculated the residuals between the observed data and the synthetic results from the product of the inverted source spectra, site responses and path attenuation, as shown in Fig. 3. The residuals were expressed as the logarithmic observed values minus logarithmic synthetic values. They vary around zero and have an average close to zero in the whole frequencies of 0.1-20 Hz. This shows that the residuals are independent on the hypocentre distance, indicating that the non-parametric inversion provides a good representation of the observed recordings considered in this study.

S O U RC E S P E C T R A
The bootstrap analysis proposed by Oth et al. (2008Oth et al. ( , 2011 was performed in this study to assess the stability of the inverted source spectra. 150 strong-motion recordings, accounting for approximately 16 per cent of the total recordings, were randomly removed from the data set, and the remaining ones were assembled as a new data set used in the inversion. We repeated this procedure 100 times to investigate the stability of the inverted source spectra. Fig. 4 shows the inverted source spectra resulting from 100 bootstrap inversions for four typical earthquakes representing four magnitude levels. The deviation from the source spectra obtained using the whole data set remains small, implying that the source spectra are stable. The cumulative attenuation within the reference distance is not included in the A(f,R) derived from the first-step inversion. The inverted source spectra from the second-step inversion absorb this cumulative attenuation when the trade-off between the site and source is solved using the known site response of the reference site. Therefore, the real source spectrum can be expressed as: where ψ(f) represents the cumulative attenuation within the reference distance and S inverted (f) is the inverted source spectrum. If the real source spectrum of an earthquake is known, the cumulative attenuation can be derived from eq. (4).
Assuming that the source spectrum follows the omega-square source model (Brune 1970), where R θ is the average radiation pattern over a suitable range of azimuths and take-off angles set to 0.55. V = 1/ √ 2 accounts for the portion of total S-wave energy in the horizontal components. ρ s and β s are the density and S-wave velocity in the vicinity of the source set to 2700 kg m −3 and 3.6 km s −1 , respectively. M 0 and f c are the seismic moment and corner frequency. We used the relationship proposed by Hanks & Kanamori (1979) to convert moment magnitude (M w ) to M 0 (unit: dynecm = 10 −7 N · m): If a small earthquake is regarded as the empirical Green's function (EGF) event of a large earthquake, the differences of the path attenuation in the strong-motion recordings at the same station from large and small earthquakes can be neglected. Fourier amplitude spectral ratio O L (f)/O S (f) can be approximately expressed as the theoretical source spectral ratio S L (f)/S S (f), where subscripts L and S represent the large and small earthquakes, respectively. According to eq. (7), the values of seismic moment and corner frequency for both large and small earthquakes could be achieved by minimizing the differences between the Fourier amplitude spectral ratio of the observed strong-motion recordings averaged over all stations triggered in both earthquakes and the theoretical source spectral ratio.
In this study, an M s 6.5 earthquake (No. 01) that occurred on 2008 August 5 at 17:49:16 (Beijing time) at the northeastern part of the Longmenshan fault was selected as a large event, and four other earthquakes (M s 4.4, 5.1, 5.3 and 5.7) were selected as its EGF events. The basic information for these earthquakes is listed in Table 1, and their epicentres and the recorded strong-motion stations are shown in Fig. 5.
The grid-searching method was adopted to determine the best-fit seismic moment and corner frequency in eq. (7). The best-fitting theoretical source spectral ratios between large and small earthquakes are in good agreement with the Fourier amplitude spectral ratios calculated using observed strong-motion recordings at frequencies of 0.1-20 Hz, as shown in Fig. 6. The obtained M w values range from 5.86 to 6.12 and the f c values range from 0.105 to 0.206 Hz for the M s 6.5 earthquake, as shown in Table 1. The values of M w are in good agreement with the one from the Global Centroid-Moment-Tensor (CMT) catalogue, that is, 6.0. According to eq. (5), the theoretical source spectra of the M s 6.5 earthquake were obtained, then the values of ψ(f) were calculated using eq. (4), as shown in Fig. 6. It shows that the ψ(f) is not strongly dependent on the selected EGF event, implying its stable estimation. In this study, we adopted the ψ(f) derived from the spectral ratio between the M s 6.5 and the M s 5.3 earthquakes (i.e. Nos. 01 and 04 in Table 1), which is approximately median of all four ψ(f). The source displacement spectra corrected using ψ(f) are shown in Fig. 7 for seven magnitude bins from 3.0 to 6.5 at an interval of 0.5 mag. Source spectra at high frequencies are close to the ω −2 decay.
Note that for a proper quantification of the stability of ψ(f), it would be useful to consider additional pairs of collocated large events/EGFs, in particular in the southwestern part of the fault. Unfortunately, such pairs are not available. Because the hypocentres of large and small earthquakes are not close enough to remove the difference of path attenuation from their sources to sites, or the strong-motion recordings obtained in both earthquakes are not enough to calculate the reliable spectral ratio between them.  Table 1, and strong-motion stations operating during these earthquakes.  Table 1, and the best-fit theoretical source spectral ratio (dashed line). (b) The cumulative attenuation within the reference distance of 30 km.

S O U RC E PA R A M E T E R S
The grid-searching method was adopted to obtain the best-fit seismic moment and corner frequency for each earthquake, making the theoretical source spectrum expressed by eq. (5) closest to the attenuation-corrected source spectrum. It can be represented as: M s − 1.0 ≤ M w ≤ M s + 1.0, the corresponding searching ranges of M 0 are derived from eq. (6). The values of stress drop ( σ ) for small-to-moderate earthquakes generally vary from 0.1 to 100.0 MPa (Kanamori 1994). Following Brune (1970), the corner frequency is expressed as f c = 4.9 × 10 6 β s ( σ /M 0 ) 1/3 . The searching ranges of f c are estimated according to the possible variation ranges of σ . Fig. 7 Figure 7. The attenuation-corrected source displacement spectra (left-hand panel), and the best-fitting theoretical source spectra to the inverted source spectra for four typical earthquakes representing four magnitude levels (right-hand panel).
shows some examples of the best-fitting theoretical source spectra. Then, the seismic moment and corner frequency were used to determine the stress drop and the source radius r according to the Brune (1970) source model: We also calculated the S-wave energy E s in the frequency range from 0.01 to 30 Hz according to the relationship proposed by Vassiliou & Kanamori (1982): where α s = 6.1 km s −1 represents the primary (P) wave velocity. The apparent stress σ a was calculated by the following relationship: where μ = 3.5 × 10 10 N m −2 represents the rigidity modulus. All of these source parameters for earthquakes considered in this study are shown in Table 2.

Seismic moment M 0 and corner frequency f c
The M w values determined in this study are in good agreement with those derived from the Global CMT catalogue, although they are slightly higher than measurements provided by Zheng et al. (2009), as shown in Fig. 8(a). We obtained the relationship between M w and M s measured by CENC by a least-squares regression analysis: There are linear deviations between M w and M s measured by CENC. M w is systematically lower than M s for M s = 3.5-6.5. This overestimation of M s is more severe in the case of larger earthquakes with the maximum close to 0.5. In fact, such phenomena have been commonly found in other studies, such as in the 2013 April 20 Lushan earthquake sequence (Lyu et al. 2013), large numbers of small-tomoderate earthquakes in mainland China (Zhao et al. 2011), some small earthquakes in the Tangshan area (Matsunami et al. 2003), and earthquakes with magnitude greater than 4.0 in the Sichuan-Yunnan region of China (Xu et al. 2010a). This deviation may result from the inaccurate calibration functions and the neglect of the base correction in the process of measuring magnitude by CENC (Zhao et al. 2011). Table 2. List of source parameters including moment magnitude (M w ), seismic moment (M 0 ), corner frequency (f c ), source radius (r), stress drop ( σ ), S-wave energy (E s ) and apparent stress (σ a ) determined in this study.     Fig. 8(b) shows the plots of seismic moment versus corner frequency for the earthquakes considered in this study. These are also compared with the constant stress drop relations corresponding to 0.1, 1 and 10 MPa. The M 0 and f c vary from 2.0 × 10 14 to 1.7 × 10 18 N · m and from 0.1 to 3.1 Hz, respectively. Corner frequencies in our study are significantly lower than those obtained in the 2009 L'Aquila earthquake sequence (Ameri et al. 2011) and earthquakes in central-eastern Iran (Hassani et al. 2011), which implies the lower stress drop. Some much smaller earthquakes in Kumaon Himalaya, India also provide a similar distribution of M 0 versus f c with a low stress drop (Sivaram et al. 2013). The seismic moment is approximately inversely proportional to the cube of the corner frequency in this study, that is, M 0 ∝ f c −3 , and the data regression yields: M 0 f c 3 is equal to 2.87 × 10 15 N · m · s −3 , which corresponds to a constant stress drop of 0.522 MPa according to the Brune (1970) model.

Stress drop
The stress drop values mainly vary from 0.1 to 1.0 MPa (Fig. 9), which are consistent with the results (≤1.0 MPa) given by Hua et al. (2009) for most of the Wenchuan aftershocks. They do not exhibit a significant dependence on the moment magnitude and hypocentre depth (Fig. 9), which indicates that the earthquakes considered in this study follow self-similarity with a constant stress drop. Studies from Allmann & Shearer (2009), Oth et al. (2010, Zhao et al. (2011), etc. all confirmed the self-similarity of global earthquakes. However, some other studies obtained conflicting results, in which the self-similarity is broken down in some specific earthquake sequences (e.g. Tusa et al. 2006;Drouet et al. 2010;Mandal & Dutta 2011;Pacor et al. 2016). The stress drops are nearly lognormal distributed with a logarithmic mean of 0.52 MPa, which is dramatically lower than the median value of 3.31 MPa for interplate earthquakes (Allmann & Shearer 2009). The Wenchuan aftershocks have small stress drop values in comparison to other large earthquake sequences, such as the 2010-2011 Canterbury, New Zealand earthquake sequence with stress release of 1-20 MPa (Oth & Kaiser. 2014), the 1983 M JMA 7.7 Japan sea earthquake sequence with stress release of 1-30 MPa (Iwate & Irikura 1988), etc. They are also much smaller than the value of the main shock, which is approximately 1-3 MPa for different finite fault-slip models (Bjerrum et al. 2010). However, the Wenchuan main shock has a stress drop value similar to other earthquakes of the same magnitude, for example, the 2001 M w 7.8 Kunlun, China earthquake and the 2002 M w 7.7 Denali, Alaska earthquake (Shaw 2013). Therefore, the Wenchuan aftershocks are characterized by obvious low values of stress drop. This characteristic was also observed in some other earthquakes, such as the 2010 JiaSian, Taiwan earthquake (Hwang 2012), and for several smaller earthquakes in the Garhwal Himalaya region (Sharma & Wason 1994). Shaw et al. (2015) proposed a physical model that shows reduced stress drops for nearby aftershocks compared to similar magnitude main shocks, because they rerupture part of the fault ruptured by the main shock which may have been partially healed. This model was supported by ground motion observations, showing smaller ground motions generated by nearby aftershocks (e.g. Abrahamson et al. 2014). Smaller values of aftershock stress drops have been also observed using corner frequency analysis of seismic sequences (e.g. Drouet et al. 2011).
In this study, an indicator ε proposed by Zuniga (1993) was used to investigate the stress drop mechanism of the Wenchuan earthquake sequence: ε < 1.0 implies a partial stress drop mechanism where the final stress is greater than the dynamic frictional stress (Brune 1970;Brune et al. 1986), whereas ε > 1.0 indicates that frictional overshoot has occurred with the final stress lower than the dynamic frictional stress (Savage & Wood 1971). The well-known Orowan's hypothesis is met when ε = 1.0 (Orowan 1960). In this study, ε equals 0.75-0.85, which indicates that the Wenchuan aftershocks can be interpreted by the partial stress drop mechanism. Sharma & Wason (1994) pointed out that such kind of aftershocks occur either when the fault locks (heals) itself soon after the rupture of the main shock passes, so the average dynamic frictional stress drops over the whole fault, or when the stress release is not uniform and not coherent over the whole fault plane, and behaves like a series of multiple events with parts of the fault remaining locked. The blank area of the seismic moment release in the ruptured area during the Wenchuan earthquake, as well as the absence of the larger aftershocks, indicates a possibility of fault lock at the unruptured areas on the fault plane (Chen et al. 2013). The low stress drop may be related to parts of the fault remaining locked on the fault plane. The apparent stress of M ≥ 3.0 earthquakes during 2000-2004 in the Sichuan province calculated by Cheng et al. (2006) is approximately proportional to 0.21 σ . This means that ε equals 1.4 (eq. 15), indicating frictional overshoot prevails over partial stress drop. The stress drop mechanism associated with earthquakes along the Longmenshan fault belt changed after the Wenchuan earthquake.
The stress drop spatial distribution was obtained by assembling and interpolating the values of all 132 aftershocks, compared with the slip distribution on the fault plane of the main shock, which was determined by Fielding et al. (2013), as shown in Fig. 10. Aftershocks were mainly concentrated on the southwest and northeast segments of the Beichuan fault, and less on the central part. Stress drop contours were generated in three segments from southwest to northeast, respectively. The higher slips emerged on the southwestern segment close to Wenchuan County. In the main shock, the Pengguan Massif began to rupture, and a large amount of stress was released on this segment ). As a result, smaller stress releases occurred for aftershocks here, with a logarithmic average of stress drop of 0.46 MPa. However, the logarithmic average of stress drop is higher, approximately 0.64 MPa for the northeastern segment near the Qingchuan County. This segment also consists of Precambrian quartzite or other stiff geological bodies. Slip on this segment is relatively smaller, and the released stress is lower. The logarithmic average of the stress drop on the central segment is close to 0.52 MPa, and the median slip value corresponds to the median stress drop. Therefore, we infer that the stress drop of the aftershocks may be related to the slip distribution on the fault plane of the main shock. Higher stress release for aftershocks occurred in areas with lower slip in the main shock.
For further verifying the above inference, we investigated the magnitude and hypocentre depth distribution between northeastern and southwestern segments, as shown in Fig. 9. The results show that both segments have a homogeneous distribution of magnitude ranging from 3.5 to 6.0, and a homogeneous distribution of depth ranging from 8 to 25 km. Furthermore, the stiffness of crustal structure shows few changes over the whole ruptured area of the Wenchuan main shock according to the CRUST1.0 model (Laske et al. 2013). Therefore, the magnitude, hypocentre depth and crust stiffness could be excluded from the cause of inhomogeneous distribution of stress drop between two segments.  This relationship means that the S-wave energy-to-moment ratio is approximately equals to 1.32 × 10 −5 , which is consistent with the result of 1.2 × 10 −5 for small earthquakes in Anchorage, Alaska derived by Dutta et al. (2003). As shown in Table 2, the apparent stress σ a varies from 0.077 to 1.606 MPa, which is directly proportional to 0.74 σ with a correlation coefficient of 0.998. The apparent stress is independent of the earthquake size, since σ is independent of M 0 , as mentioned above.

AT T E N UAT I O N C H A R A C T E R I S T I C S
The attenuation curve A(f, R) can be described in terms of anelastic attenuation and other factors ( ) related to seismic attenuation: where Q s stands for the S-wave quality factor dependent on the frequency. must be greater than A(f, R). Suppose only contains the geometrical spreading in this study. In general, the geometrical spreading can be a linear, hinged bilinear, or hinged trilinear model of R. In this study, geometrical spreading is a simple model expressed as (R 0 /R) n , where n is the geometrical spreading exponent. The greater the n value is, the stronger the geometrical spreading. According to the necessary condition of lnA(f,R) − ln(R 0 /R) n < 0, we seek out a maximum n to meet this condition for each frequency, indicating the strongest geometrical spreading. Then Q s can be evaluated from the slope of a linear least-squares fit of eq. (17) at each frequency. In this study, both geometrical spreading and anelastic attenuation were considered frequency dependent in order to deal with the trade-off between them. This strategy was also used in the study of Bindi et al. (2004). The geometrical spreading exponents at frequencies of 0.5-20 Hz for R = 30-150 km are shown in Fig. 12. The values of n vary from 0.35 to 0.75, increasing with increased frequency from 0.1 to 0.4 Hz at first, then overall decreasing until a critical frequencyaround 3.5 Hz, and finally increasing up to 20 Hz, which indicates frequency-dependent geometrical spreading in this region. Frequency-dependent geometrical spreading was also observed in North America by Babaie Mahani & Atkinson (2013) through studying response spectral amplitudes and PGAs of ground motions. Geometrical spreading in the Northeast, central United States (CUS), and the Pacific Northwest/southwestern British Columbia (PNW/BC) has a tendency to decrease at first and then increase with the increased frequency, which is very similar to what we observed in this study.
Based on the analyses of larger numbers of strong-motions recordings, previous studies have shown that n is not lower than 1.0 for local distances, while n is approximately equal to 0.5 for regional distances (Atkinson & Mereu 1992;Bora et al. 2015). The threshold for local and regional distance is related to the crustal thickness. Our study region is located at the southeast margin of the Tibetan Plateau where the crustal thickness is about 50 km. Therefore, we regard 75 km (i.e. 1.5 times of crustal thickness) as the boundary between the local and regional distances (Atkinson & Mereu 1992). A general geometrical spreading model (R 0 /R) 1.0 for R < 75 km, and (R 0 /75)(75/R) 0.5 for R ≥ 75 km is assumed. In this study, n is lower than 0.5 at frequencies ranging from 3 to 15 Hz, and 0.5-0.75 at frequencies lower than 3 Hz and greater than 15 Hz (Fig. 12). The average n value is 0.57 with a standard deviation of 0.11, obtained over frequencies ranging from 0.1 to 20 Hz. We compared the average geometrical spreading (R 0 /R) 0.57 with the general geometrical spreading mentioned above, as shown in Fig. 13. We also compared the geometrical spreading in Yunnan and southern Sichuan determined by Xu et al. (2010b), which reflects a weak attenuation of ground motion. The average geometrical spreading in this study is slightly stronger than the one given by Xu et al. (2010b), while much weaker than the general geometrical spreading. This result implies that regions near the ruptured fault of the Wenchuan earthquake show weak geometrical spreading. Boore et al. (2014) determined that the observed ground motions from China, mainly derived from the Wenchuan earthquake sequence, exhibit a weaker attenuation, which is ascribed to a 'high Q'. This may also be related to the weak geometrical spreading inferred from our study. As shown in eq. (17), the path attenuation mainly consists of geometricalspreading and anelastic attenuation (represented by Q), a potential trade-off is inherent between them.
The S-wave quality factor Q s versus frequency from 0.1 to 20 Hz is shown in Fig. 14. Q s (f) is regressed in the form of Q s0 f η , and the least-squares solution is given by 151.2f 1.06 . Other studies also provided the Q s values for the adjacent region (Fig. 14). Hua et al. (2009) obtained the Q s for the western mountains (274.6f 0.423 ) and eastern plains (206.7f 0.836 ) in northern Sichuan, separated by the Longmenshan fault belt. Zhao et al. (2011) also determined Q s = 191.8f 0.59 for western Sichuan. Compared with results from Hua et al. (2009) andZhao et al. (2011), Q s0 , representing the quality factor at 1.0 Hz, is lower in our study. However, the attenuation coefficient η is much greater than that at the mountains but close to that at the plains. Q s is closer to the results for the plains from Hua et al. (2009). The study region in this Figure 13. Comparison of the average geometrical spreading derived in this study with the general geometrical spreading, and results from Xu et al. (2010b). The averaged results of this study represent (R 0 /R) 0.57 over the hypocentre distance from 30 to 150 km. Two dashed lines represent the plus or minus one standard deviation of the average. The general results represent (R 0 /R) 1.0 for R < 75 km and (R 0 /75)(75/R) 0.5 for R ≥ 75 km.
paper is located on both sides of the Longmenshan fault belt, where the elevation suddenly drops from about 4500 m on the plateau to 500 m in the Sichuan Basin. The low Q s0 and high η may be related to the propagation path passing through the highly heterogeneous active fault belt.

S I T E R E S P O N S E
The calculated site response functions of the 43 strong-motion stations are shown in Fig. 15. Site responses for most stations are generally in good agreement with those determined by Ren et al. (2013). Compared with the site responses derived from the horizontal-to-vertical spectral ratio (HVSR) method (Fig. 15), predominant frequencies are approximately identical, while site amplifications from the non-parametric GIT are significantly higher, except for some stations (51SFB, 51SPA, 51QLY and L0021). That is because the HVSR method can approximately evaluate the predominant site frequency but underestimates the site amplification Hassani et al. 2011).
Since many analyses related to site effects in the Wenchuan earthquake sequence have been made in the study of Ren et al. (2013), our study only focused on the performance of a terrain effect array in the Wenchuan aftershocks. Stations L2009, L2002 and L2007 compose a terrain effect array, which were installed on the top (altitude 969 m), middle (altitude 960 m) and foot (altitude 927 m) of a hill (Wen et al. 2014). Fig. 16(a) shows the locations of the three stations on the hill, which share similar geological conditions. The site response functions of the three stations determined by the non-parametric GIT and HVSR method are shown in Fig. 16(b). Site responses from non-parametric GIT have significant discrepancies among the three stations, especially at frequencies of 2.0-8.0 Hz. Site amplification increases with the   increased elevation and is 1.5-2.0 times larger at L2009 than that at L2007. The site amplifications given by the HVSR method have no significant difference at the three stations, implying the HVSR method may not effectively reflect the local terrain effect. This result is in agreement with the conclusion from other studies (Parolai et al. 2004;Massa et al. 2013).

C O N C L U S I O N S
Nine hundred twenty-eight strong-motion recordings with hypocentre distances smaller than 150 km were used for separating the source spectra, path attenuation and site responses in the frequency domain using the two-step non-parametric GIT. These recordings were obtained at 43 permanent and temporary strong-motion stations during 132 earthquakes of M s 3.2-6.5, which occurred on or near the fault plane of the 2008 Wenchuan earthquake from 2008 May 12 to 2013 December 31.
We assumed that the path attenuation equals 1.0 at the reference distance of 30 km. As a result, the cumulative attenuation within this distance is transferred to the inverted source spectra when the trade-off between the source effect and site response is solved using a reference site. The cumulative attenuation was supposed as a ratio of the inverted source spectrum over the theoretical source spectrum for an M s 6.5 earthquake. Its theoretical source spectrum was determined using the Fourier amplitude spectral ratio method. Then the inverted source spectra of all 132 earthquakes were corrected by the cumulative attenuation to obtain the real source spectra, which show approximately close to ω −2 decay at high frequencies. Furthermore, a grid-searching method was used to determine the best-fit seismic moment and corner frequency. Moreover, the stress drop, source radius, S-wave energy and apparent stress were successively calculated. We investigated the scaling properties of these source parameters, and draw the following conclusions: (1) Moment magnitude M w has a linear deviation from the surface wave magnitude M s measured by CENC. M w is generally lower than M s , and is in agreement with previous studies. M 0 is approximately proportional to the f c −3 , and M 0 f c 3 = 2.87 × 10 15 N · m · s −3 . The average S-wave energy-to-moment ratio is close to 1.32 × 10 −5 . The apparent stress σ a is approximately equal to 0.74 σ , independent of the earthquake size.
(2) The value of stress drop σ for individual earthquakes varies mainly from 0.1 to 1.0 MPa, following an approximately lognormal distribution with an average of 0.52 MPa. The value is significantly smaller than the median stress drop of interplate earthquakes (Allmann & Shearer 2009), and some other large earthquake sequences. It is also much smaller than the stress drop of the Wenchuan main shock which is similar to some other large earthquakes with similar magnitude (∼8.0). This characteristic with low stress drop of Wenchuan aftershocks was investigated using the ε indicator. The results show that ε is less than 1.0, ranging from 0.75 to 0.85, indicating that the low stress drop may be interpreted by the partial stress drop mechanism. Explanations of the low stress drop in aftershocks may be related to the remaining locked parts on the fault plane of the main shock.
(3) The investigation shows that the stress drop σ has no significant dependence on the earthquake size and the hypocentre depth, indicating that the Wenchuan aftershocks follow self-similarity over the M w range of our data. The stress drop of aftershocks may be correlated to the slip distribution on the fault plane of the Wenchuan main shock. A relatively larger stress drop appeared at areas with relatively smaller slip.
The geometrical spreading is weak around the Wenchuan area within distances of R = 30-150 km, and is strongly dependent on the frequency. The S-wave quality factor Q s (f) is regressed by Q s (f) = 151.2f 1.06 . The quality factor shows strong dependence on frequency, which can be ascribed to the high heterogeneity of the crustal medium. Our study region is located on the southeast edge of the Tibet Plateau where the elevation suddenly drops from about 4500 m on the plateau to 500 m in the Sichuan Basin.
The inverted site responses of three stations from a terrain effect array show that the site amplification is strongest at the hilltop and smallest at the hillfoot, implying that the local topography considerably affects the ground motions. The site responses, calculated using the HVSR method, were not very different among the three stations. This suggests that the HVSR method may not be effectively used for analysing the local topography effect on ground motion.

A C K N O W L E D G E M E N T S
This work is supported by the Science Foundation of Institute of Engineering Mechanics, China Earthquake Administration under grant no. 2016A04, Nonprofit Industry Research Project of China Earthquake Administration under grant no. 201508005 and National Natural Science Foundation of China under grant no. 51308515. We are grateful to the editor Dr Ana Ferreira, Sylvia Hales and two anonymous reviewers for their valuable comments that helped improve our work.
All the strong-motion recordings in this paper were derived from the China Strong Motion Network Center at www.csmnc.net, (last accessed 2015 December). Earthquake parameters, including hypocentre location and measured magnitude M s , were obtained from China Earthquake Network Center at www.csndmc.ac.cn, (last accessed 2015 December). Moment magnitude M w in the Global Centroid-Moment-Tensor catalogue was obtained from http://www.globalcmt.org/CMTsearch.html, (last accessed 2015 December). ρ s , β s and α s were derived from the CRUST1.0 model at http://igppweb.ucsd.edu/∼gabi/crust1.html.