Arias intensity attenuation relationship in Sichuan–Yunnan region, China

Arias intensity is an essential ground motion measure correlating with the potential for earthquake-induced landslides. The Sichuan–Yunnan region, which is primarily mountainous, is a high incidence region of earthquake-induced landslides in China. However, there is no available attenuation relationship for this intensity measure due to the backward construction of the stations. In this study, we developed a region-specific Arias intensity attenuation relationship using the China Strong-Motion Networks Center (CSMNC) database which was established in 2008. We recommend this relationship be applied in the Sichuan–Yunnan region for moment magnitudes ranging between 4.2 and 7.0, distances ranging between 0 and 300 km and with Vs30 (the average shear-wave velocity in the upper 30 m of a soil profile) ranging between 128 and 760 m/s. The current study finds that this relationship’s intra-event, inter-event, and total standard deviations are greater than for other regions. This is likely caused by the complicated seismotectonic activities, nonlinear site effects, error from inferring Vs30, basin effects, etc. However, this relationship has the best performance in fitting and predicting the data from the Sichuan–Yunnan region.

The Sichuan-Yunnan region located in the southwest of China is one of China's most seismically active regions. During the Cenozoic, the Sichuan basin was hindering the Indian plate from squeezing the Eurasian plate, resulting in many fault zones such as the Longmenshan Thrusts zone (Fig. 1). The east boundary of the Sichuan-Yunnan faultedblock (Tapponnier et al. 1982;Zhang et al. 2003) formed a substantial left-lateral strikeslip fault system with a total length of more than 1100 km named the "Kangding fault" system (Tapponnier and Molnar 1977). Allen et al. (1991) and Wen et al. (2000) investigated the history of earthquakes in the Sichuan-Yunnan region that had occurred over a period of hundreds of years and the active strike-slip fault zones in the intra-plate area of China and found that intensive seismotectonic activities often appear around the "Kangding fault" system.
The Sichuan-Yunnan region, primarily mountainous, is a high incidence region of earthquake-induced landslides leading to massive casualties and economic losses (e.g. 197,481 landslides in the 2008 Wenchuan M w 7.9 earthquake ; 22,528 landslides in the 2013 Lushan M w 6.8 earthquake (Xu et al. 2015); 10,559 landslides in the 2014 Ludian M w 6.2 earthquake (Wu and Xu 2018); 4800 landslides in the 2017 Jiuzhaigou M w 6.5 earthquake ). Due to the backward construction of the station, there was a lack of strong-motion recordings in this region before 2000. The China Strong-Motion Networks Center (CSMNC) system completed in March 2008 has collected a significant amount of high-quality digital strong-motion recordings in this region, including numerous near-field recordings of large earthquakes such as the Wenchuan M w 7.9 earthquake. The collection of these useful recordings is benefiting from a dense strong motion station distribution. Those recordings have greatly enriched the strong-motion database in China and enabled the establishment of the attenuation relationship in the Sichuan-Yunnan region. Li et al. (2020) developed the ground motion prediction model for horizontal PGA, 5% damped response spectrum in the Sichuan-Yunnan region. Xu et al. (2012) fitted the empirical relationship of Newmark displacement with PGA using the Wenchuan M w 7.9 earthquake recordings. Jin et al. (2018) regressed a new empirical relationship of Newmark displacement with Arias intensity, with the smallest error and the maximum correlation coefficient compared to other relationships with Arias intensity, and recommended for southwest China. For predicting Newmark displacement based on Arias intensity, there would be a lack of region-specific attenuation relationship in the Sichuan-Yunnan region.
In this study, we developed an Arias intensity attenuation relationship in the Sichuan-Yunnan region, considering source effects (magnitude and focal mechanism), path effects (source to site distance), and site effects (V s30 , the average shear-wave velocity in the upper 30 m of a soil profile). Through model comparison, uncertainty analysis, model testing, and prediction comparison of Newmark displacement we explored the accuracies of estimation and prediction of this model. This study is of great significance to the seismic risk analysis of earthquake-induced landslides in southwest China, which is the intention of our work.

Data acquisition and processing
The strong-motion data used in this study were collected from CSMNC (see Data Availability). The selection of strong-motion recordings mainly adhered to the following criteria: (1) Selecting the recordings collected from free-field stations, and non-free-field station recordings were removed. (2) Excluding the recordings not in the Sichuan-Yunnan region.
(3) Recordings of obviously poor-quality were excluded following the visual inspection of the strong-motion recordings (for example, recordings with multiple wave packets, severe drop tail, spikes, noise-dominated recordings, and unreliable recordings). (4) Recordings with incomplete two horizontal components were excluded. (5) We did not use recordings with the source to site distance greater than 400 km (rupture distance for M w > 6, and hypocentral distance for M w ≤ 6).
The V s30 (the average shear-wave velocity in the upper 30 m of a soil profile) is a popular parameter for expression of site effects in attenuation relationships (e.g. Anderson et al. 1996;Castro et al. 1997;Park and Elrick 1998;Anderson et al. 1996;Lee et al. 2012). The V s30 values of all stations in this dataset were not directly measured from the shear wave velocity profile. Instead, the V s30 values of approximately half of the stations were derived from the literature (Yu and Li 2015;Zhang et al. 2020) and the NGA-West2 database flat file (https:// peer. berke ley. edu/, last accessed on January 5, 2021). For other sites where the shear wave velocity is available (Li et al. 2013;Zhao et al. 2019) according to the GB 50,011-2010 seismic design code (Chinese Standard 2016), the V s30 values were calculated by Eq. 1 modified from Wang et al. (2010). For sites without borehole information, the V s30 values were inferred from geological units, considering such factors as geotechnical category, terrain-based categories, etc. where d is the depth to the rock, in m and V se is the equivalent shear wave velocity (the average shear-wave velocity of the soil layers in the top 20 m or the soil above rock for sites with d less than 20 m), in m/s. Figure 3 shows the distribution of magnitudes, distances, and V s30 values of these recordings. The dataset includes recordings having moment magnitudes ranging between 4.2 and 7.9, distances ranging between 0.481 and 396.29 km, and V s30 values ranging between 128 and 760 m/s (according to the National Earthquake Hazards Reduction Program (NEHRP) site classification standard Building Seismic Safety Council (BSSC) (1994)). Most of the recording stations belong to class C and class D sites, and 6 stations belong to class E sites.
(1) The Arias intensity formula selected for this study is the average of the Arias intensity from two perpendicular horizontal components, as expressed in Eq. 2: where I xx and I yy are the Arias intensities from two perpendicular horizontal components.

Predictive attenuation relation for Arias intensity
The Arias intensity attenuation relationship form used in this study is modified from Lee et al. (2012) as follows: where I a is the average of the Arias intensity from two horizontal components, in m/s, M is the moment magnitude, R is the source to site distance (the closest distance to the rupture plane (rupture distance) for large earthquakes (M w > 6) and hypocentral distance for others (M w ≤ 6), in km), V s30 is the average shear wave velocity of 30 m on the soil profile, in m/s, F N and F R are dummy variables indicating the fault types ((1, 0) for normal and normaloblique faults, (0, 1) for reverse and reverse-oblique faults, (0,0) for other (strike-slip and unknown faults)), η are the inter-event residuals conforming to the normal distribution of N(0, ϕ 2 ); ε are the intra-event residuals conforming to the normal distribution of N(0, τ 2 ); c 1 , c 2 , …, c 8 are regression coefficients. Equation 4 determine the total standard deviation σ.

Results and evaluations
After using the regression analysis mentioned above, the regression coefficients are listed in Table 2. The inter-event standard deviation ϕ is 0.852 (in natural log format), the intraevent standard deviation τ is 1.270 (in natural log format), and the total standard deviation σ consisting of ϕ and τ is 1.529 (in natural log format).
In Eq. 3, the negative value for F N means that the Arias intensity of the normal (including normal-oblique) event is 0.456 (in natural log format) smaller than that of the strikeslip event. The positive value for F R means that the Arias intensity of the reverse (including reverse-oblique) event is 0.901 (in natural log format) larger than that of the strike-slip event. The Arias intensity attenuation curves for different magnitudes and fault styles (Fig. 4) show the same trends with Eq. 3, indicating effects of focal mechanism, which is consistent with other previous research (e.g. Travasarou et al. 2003;Stafford et al. 2009;   Fig. 6f). It is observed that the attenuation curve obtained in this study fits the distribution of the Arias intensity of earthquakes in the Sichuan-Yunnan region well. It is worth noting that the far-field data (distance > 300 km) seems to be overestimated for the Wenchuan M w 7.9 earthquake (Fig. 6a). We can clearly see from Fig. 3 the limited data for events of M > 7, especially in the farfield. As a result, it is hard to constrain our model in M > 7 and long-distance range without more data from big earthquakes in this region. Thus, our model should be used carefully for events of M > 7.
Several models were selected for shallow earthquakes with the same consideration factors (e.g. fault types, site effect, etc.) as the comparative models for this study, including Travasarou et al. 2003 (Scherbaum et al. 2004). Figure 7 compares this study's median attenuation curve (strike-slip event, V s30 = 280 m/s or class D sites) with the comparative models. For the large event (M w = 7.5), the estimated value at near sites (R < 10 km) for this study is slightly larger than that of FPS12 and much larger than those of other comparative models. The estimated value in the sites of R > 10 km of this study is smaller than those of Hwa04 and SA17 (the gap increasing with distance) and minimally different from those of other comparative models (Fig. 7a). For the middle event (M w = 6), the estimated value of this study model at the near sites (R < 10 km) is larger than those of other models (except Sta09, FPS12, SA17 and CB19 in sites of about 1 km < R < 10 km). The estimated value in the sites of R > 10 km of this study is smaller than those of Tra03 and Hwa04 (the gap increasing with distance) and also smaller than those of Sta09, SA17 and CB19 (the gap decreasing with distance). The estimated value of this study is minimally different than those of FPS12 and Lee12 (Fig. 7b). For the small event (M w = 4.5), there is an intersection at the distance of around 70 km in the curves of this study, Tra03, Hwa04, Sta09 and SA17, and the estimated value of this study is always larger than those of FPS12, Lee12 and CB19 for the entire distance range (Fig. 7c). The above comparison indicates that there are noticeable differences between this study model and the comparative models (the comparison of normal and reverse events for class D sites as Figs. 15,16 in Appendix,and class C sites comparsion as Figs. 17,18,19 in Appendix). The differences in the near sites may primarily be due to, the datasets being derived from differing regions. The differences in the middle-and far-fields may be partly due to the difference in characteristics of the datasets, and partly, the difference in slopes of the attenuation curves related to Q values. Figure 8 shows the distributions of intra-event residuals with V s30 , distance, and earthquake magnitude. The V s30 -binned means of intra-event residuals show a slightly decreasing trend with V s30 (Fig. 8a), probably representing the nonlinear site effects. There is no systematic trend in intra-event residuals with distance (Fig. 8b) and earthquake magnitude (Fig. 8c). The intra-event residuals obey a log-normal distribution (Fig. 8d).

Uncertainty analysis
The distribution of inter-event residuals with earthquake magnitude is shown in Fig. 9. There is no obvious trend of inter-event residuals for reverse (including reverse-oblique), normal (including normal-oblique), or strike-slip earthquakes. The positive values of interevent residuals for most unknown earthquakes indicate the possibility of reverse events. The dispersion of inter-event residuals for small and moderate events is noticeably stronger than for large events. As a result, the inter-event standard deviation of this study is 0.852, which is larger than those of the other comparative models noted earlier. One possible explanation may be the different source effects derived from the complicated seismotectonic activities in the Sichuan-Yunnan region (Fig. 1).
The total standard deviation σ in this study (1.529) is larger than that of Tra03 (1.328), larger than that of Hwa04 (1.290, 1.230, 1.250, and 0.820 in class B, C, D, and E sites, respectively), larger than that of Sta09 (0.973 (form 1), 0.952 (form 2) and 1.019 (form 3) in soil sites; 1.153 (form 1), 1.138 (form 2) and 1.170 (form 3) in rock sites), larger than that of FPS12 (1.171), larger than that of Lee12 (0.994), larger than that of SA17 (1.290-1.308), and larger than that of CB19 (0.919-1.32). In addition to the inter-event standard deviation analyzed above, the intra-event standard deviation of this study is 1.270, contributing significantly to the total standard deviation.
The intra-event residuals ε could be further decomposed into inter-site (site-to-site) residuals ε s and remain residuals ε r . The intra-event standard deviation τ could also be decomposed into inter-site standard deviation τ s and the remainder τ r (Chen and Tsai 2002). The ε s reflects the site effects in station s can be estimated as: where n es is the number of earthquakes recorded by station s, ε i is the intra-event residual of the i th earthquake in station s.
The inter-site standard deviation τ s can be expressed as follows: where n s is the number of stations. This dataset consists of 1605 recordings obtained at 491 stations, and most stations have only a few recordings. The inter-site residuals of those stations have no statistical meaning. Nineteen stations with more than 10 recordings each were selected, with the distribution of inter-site residuals calculated by Eq. 5 with V s30 is shown in Fig. 10. The decreasing trend of the inter-site residuals with estimated V s30 indicates substantial nonlinear site effects. Due to the error from inferring V s30 , the inter-site residuals of some sites are biased very heavily (e.g. 2.040 in 53YPX station, − 0.635 in 53DFD station, etc.). As a result, there is an insignificantly decreasing trend of inter-site residuals in inferred V s30 sites. The inter-site standard deviation, which was calculated according to Eq. 6, peaked at 0.603, indicating an outstanding contribution to intra-event standard deviation.
The latitude and longitude values of the stations provided by the CSMNC retain one decimal place. Therefore, we could not accurately locate them by coordinates. For stations located in the Sichuan basin, we can confirm that they are in the basin because of the massive area of the Sichuan basin (Fig. 1). Almost stations in mountainous areas are constructed in townships. And the township in China in mountainous areas are generally located on river terraces in intermountain basins because of convincing transportation. Thus, the basin effects are also contributing to intra-event standard deviation. However, due to the lack of overburden thickness informaiton such as Z 1.5 or Z 2.5 , we cannot characterize the basin effect in our model.
After the Wenchuan M w 7.9 earthquake, the damage investigation showed more earthquake-induced landslides in the hanging wall than the foot wall (Xu and Li 2010), likely indicating the hanging wall effect on Arias intensity for large earthquakes. Following the criteri of the hanging wall effect found by Abrahamson and Silva 2008 (M w > 6, dip < 70, |R x |< 50 km (R x , the horizontal distance to the top edge of the rupture measured perpendicular to strike, positive value for hanging wall and negative value for foot wall), and R jb < 30 km (R jb , Horizontal distance to the surface projection of the rupture)), hanging wall effect recordings were selected. Only some near-field recordings of the Wenchuan (M w 7.9) and Lushan (M w 6.8) earthquakes can meet the criteria. Figure 11 shows the distribution of intra-event residuals of selected recordings with R x . There is no apparent difference between the hanging and foot wall, seemly indicates  Table 1 are potentially mainshock-aftershock sequences (e.g. 2 # to 3 # , 5 # to 9 # , 24 # to 26 # ). Studies (e.g. Douglas and Halldórsson 2010;Lee et al. 2020) have observed some differences in the characteristics of spectral accelerations between mainshock and aftershock ground motions. After examining the residuals, we find intra-event residuals of aftershocks are not sensitive to distance (Fig. 12a), magnitude (Fig. 12b), and V s30 (Fig. 12c). The inter-event residuals of aftershocks are sensitive to fault styles (Fig. 12d), and show positive values for reverse events and negative values for strike-slip events. Therefore, our model would have overestimated to aftershocks of strike-slip events and underestimated to aftershocks of reverse events. Bulletin of Earthquake Engineering Fig. 8 Distribution of intra-event residuals with different parameters: a V s30 , b distance, and c earthquake magnitude. d Histogram of intra-residuals and its fit to a normal distribution curve. Box and error bar denotes the mean ± one standard deviation of binned intra-event residuals Fig. 9 Distribution of interevent residuals with earthquake magnitude Fig. 10 The inter-site residuals against V s30 . The big circles and curve represent inter-site residuals of estimated V s30 sites (having shear wave velocity information in Data Acquisition and Processing) and fit curve of them; the small circles represent inter-site residuals of inferred V s30 sites (no shear wave velocity information in Data Acquisition and Processing) Fig. 11 Distribution of intra-event residuals with R x and binned means and standrad deviations of them for a Wenchuan M w 7.9 earthquake, and b Lushan M w 6.8 earthquake Fig. 12 Distribution of intra-event residuals of aftershocks with a distance, b Vs30, and c magnitude. d Distribution of inter-event residuals of aftershocks with magnitude of this study model was tested using the strong-motion data of the Yangbi earthquake (25 recordings).
The goodness-of-prediction of the median attenuation curves of this study to the Yangbi earthquake data is shown in Fig. 13a. Compare to the current study (Fig. 13b), other models have the tendency of overprediction (Fig. 13c-e) or underprediction (Fig. 13f) of the data from the Yangbi earthquake.
For comparison of the prediction results quantitatively, we introduced the root mean square error (RMSE) as follows: where Δ s is the predicted residuals for site s as defined in Eq. 8, and n is the number of recordings.
where ln I as is the natural log of observed Arias intensity for site s, ln I as is the natural log of predicted Arias intensity for site s.
∑ n 1 Δ 2 s n (8) Δ s = ln I as − ln I as Fig. 13 a Prediction of the median and one sigma Arias intensity attenuation curve of this study and observed data in Yangbi earthquake. Distribution of the predicted residuals (Eq. 8) of b this study, c Tra03, Hwa04, d Sta09, e SA17, CB19, f FPS12, and Lee12 with distance A smaller RMSE denotes a higher predictive accuracy. The RMSE values of different models based on the Yangbi earthquake data are shown in Table 3. The RMSE of this study model is 1.203, the lowest of all models, indicating that this study model has the highest accuracy in predicting the data of the Yangbi earthquake. 7 Prediction of Newmark displacement of the Lushan M w 6.8 earthquake A dataset of 8 stations with PGA > 0.2 g from the Lushan M w 6.8 earthquake is used to test the prediction of Newmark displacement. Basic information concerning these stations and actual Newmark displacement for different critical acceleration A c (0.02 g, 0.04 g, 0.06 g, 0.08 g, 0.10 g and 0.12 g) are listed in Table 4. First, according to the magnitude, distances, and V s30 in Table 4, we estimated the Arias intensity of each station using our model. We then calculated the Newmark displacement with different A c using Jin's empirical relationship with Arias intensity (Jin et al. 2018). As a comparison, we used the PGA attenuation relationship proposed by Li et al. (2020) and Xu's empirical relationship (Xu et al. 2012) of Newmark displacement with PGA to estimate another group of Newmark displacement. The predicted residuals of the two methods are shown in Fig. 14. The discreteness of residuals based on Arias intensity is slightly fewer baizes than those based on PGA, especially in the range of Newmark displacement > 10 cm. As a result, the predicted RMSE using Arias intensity (1.054) is also less than that using PGA (1.118). Thus, the Arias intensity attenuation relationship in the Sichuan-Yunnan region would be a new choice for estimating Newmark displacement and further risk analysis work of earthquake-induced landslides in southwest China.

Conclusions
This study has developed a region-specific Arias intensity attenuation relationship for the shallow earthquakes in the Sichuan-Yunnan region using 1605 recordings from 78 earthquakes recorded by 491 stations from 2008 to 2020. The relationship formula is modified from Lee et al. (2012), considering source (magnitude and focal mechanism), path, and site effects. As with previous studies, the strong effects of focal mechanism and site category on Arias intensity were found in this study. However, hanging wall effects were found to be seemly insignificant, and the aftershock effect was sensitive to fault style. The total standard deviation of this relationship is 1.529, which is larger than found in previous studies and likely caused by the complicated seismotectonic activities, the nonlinear site effects, the error from inferring V s30 , the basin effects, and other possible factors. Development of a robust relationship would be proposed by using a constantly updated and more reliable database in the future. Bulletin of Earthquake Engineering Table 4 The selected Lushan M w 6.8 earthquake data used for testing of Newmark displacement Prediction The actual Newmark displacement was calculated by double integrating those parts of the strong-motion record that lie above the critical acceleration

Bulletin of Earthquake Engineering
Our comparisons with similar models for other regions indicate significant differences in Arias intensity (Figs. 7,15,16,17,18,19). The accuracy of the output of this relationship relative to the data in the Sichuan-Yunnan region and the highest predictive accuracy of this relationship to the data of Yangbi earthquake in Yunnan province demonstrate the best performance in modeling and prediction for this region. However, it is hard to constrain our model in M > 7 and long-distance range due to the data limited. Therefore, we consider this relationship to be valid in the Sichuan-Yunnan region for moment magnitudes ranging between 4.2 and 7.0, distances ranging between 0 and 300 km, and V s30 ranging between 128 and 760 m/s. This relationship may only provide a reference for the shallow earthquakes in the Sichuan-Yunnan region and the other regions with similar seismotectonic settings. Thus, it should be used carefully in other settings.