Spatial Patterns and Scaling Distributions of Earthquake-Induced Landslides—A Case Study of Landslides in Watersheds Along Dujiangyan–Wenchuan Highway

Earthquake-induced landslide has various spatial characteristics that can be effectively described with the frequency–area curve. Nevertheless, the widely used power-law curve does not reflect well the spatial features of the distribution, and the power exponent does not show the association with the background factors. There is a lack of standards for building the relationship, and its implication on the spatial distribution of landslides has never been analyzed. In this study, we propose a new form of frequency distribution and explore the parameters in the typical watersheds along the highway from Dujiangyan to Wenchuan in the Wenchuan earthquake region. The obtained parameters are related to the landslide density and proportions of the large-scale landslides. Furthermore, a hot spot analysis of landslides in the watersheds is conducted to assess the relationship between the parameters and the spatial cluster patterns of landslides. The hot spots highlight the size and distance of landslide areas that cluster together, whereas the distribution parameters reflect the density and proportions of landslides. This research introduces a new method to analyze the distribution of landslides and their association with the spatial features, which can be applied to the landslide distribution in relation to other influential factors.


INTRODUCTION
Earthquake-induced landslides have different spatial distribution characteristics, which depend on the tectonic features of the earthquake and the geomorphological conditions of landslides. For example, the 2002 Alaska earthquake is a surface rupture earthquake with the induced landslides that were obviously distributed along the long axis in the NW-SE direction of the seismic Denali fault (Gorum et al., 2011). The 2013 Lushan earthquake occurred in a hidden fault and resulted in the relatively scattered landslides with distribution over the affected area (Xu et al., 2015).
Studies have observed that landslides follow the power-law frequency-area relationship (Malamud et al., 2004;Chong et al., 2013;Xu et al., 2014;Tanyaş et al., 2017), yet there is no standard criterion for determining the power-law exponent. Cumulative and noncumulative frequency-area curves are used for the research and creation of the landslide inventories that have rainfall or earthquake origin (Fujii, 1969;Hovius et al., 1997Hovius et al., , 2000. Landslide inventories in different regions seem to have different exponents, yet without comparability in relation to the landslide backgrounds. The landslide frequency-area curve is generally used to show the completeness of the landslide inventory. According to some researchers, the phenomenon that the tail of the frequencyarea curve obeys the power law cannot be attributed to the incompleteness of data. The overall variables of landforms, especially the different forms of the slope, control the two parts of the frequency-area curve. The frequency-area distribution is the combined result of the distance to slope and the actual terrain constraints of landslide (Pelletier et al., 1997;Hovius et al., 2000;Guthrie and Evans, 2004). Some scholars have quantified the spatial distribution characteristics of the earthquake-induced landslides (Harp and Jibson, 1996;Cardinali et al., 2000;Bucknam et al., 2001;Stark and Hovius, 2001;Guzzetti et al., 2002;Brardinoni and Church, 2004;Guthrie and Evans, 2004;Malamud et al., 2004;Korup, 2005;Van Den Eeckhaut et al., 2007;Corominas and Moya, 2008;Guan, 2018). In some cases, the special probability distributions (e.g., the three-parameter inverse-gamma distribution) fit the landslide area distribution well (Malamud et al., 2004), yet no physical implications are found in the parameters. Moreover, the landslide distribution varies with time (Fan et al., 2018), but there are still no specific associations with the utilized distribution parameters.
In summary, previous studies have used the power-law distribution for the entire landslide inventory in an earthquake area, often confusing landslides with different seismic and landform conditions or even including landslides from different historical events. Therefore, the power-law exponents have no relationship with any of the specific influential factors. Furthermore, it is frequently noticed that the landslides in the subregions of the earthquake-affected area do not follow the power-law distribution at large scales. This requires a more effective method to describe the size distribution and spatial patterns of landslides in an earthquake area.
This study investigates landslides in several watersheds in the Wenchuan earthquake area based on the latest complete inventory of seismic landslides (Xu et al., 2014). A new distribution form of landslide area is proposed, and then the distribution parameters in relation to the spatial characteristics of the landslides are discussed. This method is expected to be applicable for the landslide distribution characteristics under other influencing factors.

DISTRIBUTION OF LANDSLIDES INDUCED BY THE WENCHUAN EARTHQUAKE Data Source and the Frequency-Area Curve
The landslide inventory used in the study was published in 2014 (Xu et al., 2014). It is based on 86 high-resolution images with pre-and post-earthquake interpretations, including aerial photographs with 1-, 2-, 2. 4-, and 5-m resolutions, SPOT 5 with a 2.5-m resolution, CBERS02B with a 19.5-m resolution, IKONOS with a 1-m resolution, ASTER with a 15m resolution, IRS-P5 with a 2.5-m resolution, QuickBird with 0.6-and 2.4-m resolutions, and ALOS with a 2.5-m resolution. The inventory contains 197,481 landslides with a total area of 1,160 km 2 , and individual landslide areas range from 0.00003 to 6.972824 km 2 , covering an area of 110,000 km 2 (Figure 1). This inventory provides more detailed landslide data and is referenced according to Harp's earthquake landslide distribution cataloging and mapping rule. Accordingly, more small-scale landslides were identified, and the result is more complete and closer to the actual situation of the Wenchuan earthquake landslide (Harp and Jibson, 1996;Xu et al., 2014).
In the process of obtaining the frequency-area curve, binning is used to divide landslide areas, and both equivalent and inequivalent bin widths are used. The inequivalent bin width increased with the increase of landslide area, so the bin widths are approximately equal in logarithmic coordinates (Malamud et al., 2004). Each bin includes the proportion of the area occurrence, and both cumulative and noncumulative proportions are used. It is generally observed that the frequency-area curve has a straight line for medium and large landslides, which implies that the distribution satisfies the power law: where A is the landslide area, P(>A) is the cumulative proportion of landslides with areas greater than A, and c and α are parameters. Similarly, the noncumulative curve also applies to the distribution in many cases: where p(A) is the noncumulative proportion of landslides, and c and β are parameters. For the present case, we put forward the cumulative frequency-area curve of the entire landslides (Figure 2), which shows the power-law range from 0.1 to 1 km 2 , with c = 1.234 × 10 −5 , α = 2.46741, and R 2 = 0.9976. This differs remarkably in the power exponent compared with the result of Xu et al. (2014), who used the same data and obtained c = 13, α = 2.0745, and R 2 = 0.9931, ranging from 0.01 to 1 km 2 . It illustrates that the power-law curve strongly depends on the processing method and that the same data may result in different exponents. Moreover, even if the power-law curve holds, it covers only a very small portion of the landslide (about less than 10 or 1%). There is no detailed study on the power-law parameters, partly because in the process of obtaining the landslide frequency-area relationship, some researchers used an equivalent bin width, whereas others used inequivalent, some used cumulative frequency-area, whereas others used noncumulative (Fujii, 1969;Hovius et al., 1997Hovius et al., , 2000. There is no uniform standard for the bin width and proportion of area, and the valid range of fitting line is greatly influenced by subjective factors, which has greatly influenced the distribution of parameters. Therefore, it is difficult to analyze these parameters and to use them for the analysis of the spatial characteristics of landslides.

Study Area
The power-law curve is usually applied to the entire inventory of earthquake-induced landslides. This inevitably confuses landslides with different geological and geomorphological conditions and makes the power exponent meaningless in relation to the influencing factors. Therefore, it is necessary to consider landslides in certain control conditions. For example, we can consider landslides in watersheds under similar seismic and geological conditions, and it should reflect the landslide distribution under geomorphological conditions. Watershed is taken as a unit area because it is a natural site where landslides occur, and all influential factors can be conveniently determined in an individual watershed.
The suitable area for this study should have the following characteristics: (i) different landslide density distributions, (ii) different topographies, and (iii) greatly affected by earthquakes, namely with a similar background but local differences. Thus, for the study area, we selected the region along the highway from Dujiangyan to Wenchuan (DW) in the Wenchuan county in the Sichuan province. The study area is approximately 1,200 km 2 .
The Wenchuan earthquake created 1 × 10 9 m 3 of landslides along the highway, damaged 80% of the road, more than 10 km was completely covered by the collapsed material, more than 50 bridges were damaged, and traffic was disrupted, causing great difficulties for disaster relief (Zhuang et al., 2010). The epicenter is near the region and lies in its middle and lower part, with the Longmenshan (LMS) fault passing through the region from southwest to northeast. The fault broke from the epicenter to the northeast, so landslides were more severe in the northeast from the epicenter than in the southwest (Figure 1).
The distance to the fault in the study area is from 3 to 18 km, and the peak ground acceleration is from 0.56 to 1.32 g in the Wenchuan earthquake. Mountains and steep gorges characterize this region, the high terrain is in the northwest, and the low terrain is in the southeast, the elevation is from 1,014 to 4,450 m, and the average slope angle is from 30 to 35 • . Lithologies are mainly granitoid and intermediate with coal seam, pyroxene, and other sedimentary rocks. Belonging to the mountainous subtropical humid monsoon climate zone, it is one of the rainfall centers in the western Sichuan province, where heavy rain often occurs.
In this study, a digital elevation model with a resolution of 30 m was adopted, and a SWAT modular of ArcGIS was used to divide the study area into watersheds. This study included 11 watersheds along the DW highway and a total of 11,631 landslides with an area of 66.82 km 2 . The watersheds are the Longtan (LT) watershed, the Xingfu (XF) watershed, the Taiping (TP) watershed, the Yeniu (YN) watershed, the Chediguan (CDG) watershed, the Luoquanwan (LQW) watershed, the Dayin (DY) watershed, the Huanglian (HL) watershed, the Xiaogou (XG) watershed, the Zhuanjinglou (ZJL) watershed, and the Qicenglou (QL) watershed. Among them, the XG watershed has the largest total number of landslides (TNl) (2,441), whereas the smallest number of landslides (166) is in the ZJL watershed. The largest total landslide area (TAl) (11.11021 km 2 ) is in the XG watershed, whereas the smallest is in the ZJL watershed (0.421146 km 2 ), as shown in Figure 3 and Table 1.

Scaling Distribution
Considering the characteristics of landslide distribution in the study area, we used the equivalent bin width. The bin interval was set at 0.0001 km 2 , and by using the cumulative curves, we found that the curves do not follow the power-law form even on the tails (Figure 4). Rather, we found that the landslides of all scales are subject to a scaling distribution (Eq. 3) that combines the powerlaw with the exponential distribution (Yong et al., 2013(Yong et al., , 2017: where, C, ρ, and A c are the parameters listed in Table 2. We found that C and ρ are related by a logarithmic relation (Figure 5), which means C is not independent, and the distribution is determined by the parameters ρ and A c . It is expected to be derived from some underlying probability distributions. Accordingly, the distribution is reduced to the pair ρ and A c .
We found that ρ can well reflect the density of the landslide area (Dl), and A c coincides with Pl, the proportion of large landslides, with a lower limit of 0.01 km 2 as determined by A c (Figure 6).
As for ρ and A c , we found that they are linearly related in six watersheds and show a discrepancy in the other five watersheds. Two points (CDG watershed, LQW watershed) have deviated from the fitting line and are located on the left and above the line, whereas three of them (TP watershed, XG watershed, LT watershed) are on the right and below the line, meaning that  A c of the CDG watershed and the LQW watershed are higher than ρ, whereas A c of the TP watershed, the XG watershed, and the LT watershed are lower than ρ (Figure 7). Accordingly, we considered that the number of large-scale landslides is more than small-scale landslides in the CDG and LQW watersheds, but the small-scale landslides have a major part within the TP watershed and the XG watershed. The tail landslide distribution within watersheds did not satisfy the power law. Landslides of all scales are subject to the scaling distribution and not just the tail of the distribution. ρ and Ac can reflect the density of the landslide area and the proportion of the large landslides, respectively.

Spatial Analysis
After the analysis discussed earlier, we examine the spatial distribution characteristics of landslides within the watersheds. We have noticed that landslides inside the watersheds can be clarified in three categories, as shown in Figure 8. The first type includes the watersheds ZJL, QL, LT, and XF. Their ρ values are negative, their landslide densities and areas are small, and distributions are dispersed ( Figure 8A). The second type includes the watersheds HL, CDG, LQW, and XG. Their ρ values are from 0.0455 to 0.0642, and their landslide densities and the degree of concentration are higher than in the four watersheds mentioned before ( Figure 8B). The third type includes the watersheds TP, YN, and DY. Their ρ values (from 0.0817 to 0.1191) and densities are the highest among the three types of watersheds ( Figure 8C). From the analysis discussed earlier, it can be seen   Figure 9C).
In addition, we found that the watersheds with the lowest A c values are also those with the lowest ρ values, such as the watersheds ZJL, LT, QL, and XF (Figures 8A, 9A). However, the medium and large values of the two parameters are not consistent. For example, the CDG watershed and the LQW watershed have medium ρ values and large A c values. We then assessed their spatial distribution characteristics and found that the CDG watershed and the LQW watershed have more largescale landslides than other watersheds. These landslides are concentrated together, and despite their medium densities, their A c values are large, and ρ values are medium (Figures 8, 9). In contrast, more medium-and small-scale landslides are dispersed within the TP watershed despite large densities in the watershed; their A c values are medium, and ρ values are large (Figures 8, 9). This characteristic may well prove our conclusion that the proportion of the large-scale landslides is higher than the smallscale landslides within the CDG and LQW watersheds, but within the TP watershed is the opposite. We found that the CDG watershed has the highest A c value, and the landslide concentrate degree of this watershed is one of the greatest. Thus, we conjectured that A c could reflect the concentration degree of large landslides, which is expected to be confirmed by sufficient data.

Z-Scores for the Watersheds
To analyze the concentrated degree of landslides within the watersheds, we performed a hot spot analysis using the Getis-Ord Gi * tool of ArcGIS, which provides a value related to the clustering of landslides. Taking the landslide area as weight, this tool works by looking at each landslide in the context of a neighboring landslide to calculate which features with either high or low values are spatially clustered. In our research, high and low values correspond to large and small landslide areas, respectively. A landslide with a large area surrounded by landslides with a large area will be defined as a statistically significant hot spot. The local sum for one landslide and its neighbors is proportionally compared with the sum of all landslides in the region. When the local sum is very different from the expected local sum and when that difference is too large to result from a random chance, a statistically significant Z-score will be generated (Eq. 4). The larger the value of the statistically significant positive Z-scores, the more intensive is the clustering of high values (hot spot). For statistically significant negative Z-scores, the smaller the Z-score, the more intensive the clustering of low values (cold spot). There is no obvious spatial clustering if the Z-score is close to zero (not significant) [Copyright(C) 1995-2013 Esri]: where x j is the area of jth landslide, ω i,j is the weight of the distance between ith landslide and landslide jth, and n is the total number of landslides:X  To analyze the clustering degrees within the watersheds, we calculated the Z-score sum for hot spots (Zsh), the number of hotspots (Nh), the average Z-score for hot spots (Zah), the number of cold spots (Nc), the Z-score sum for cold spots (Zsc), and the average Z-score for cold spots (Zac), as shown in Table 3.

Distribution Parameters and Z-Scores
The sum of the Z-scores for hot spots (Zsh) and the sum of the Z-scores for cold spots (Zsc) have corresponding relations with ρ, as shown in Figures 10A,B. However, the Z-scores have no linear relationship with A c , as shown in Figures 10C,D. (Notes: The ZJL watershed does not have a cold spot, and there are only 10 watersheds for Zsc-ρ and Zsc-A c .) Although the CDG watershed has the highest A c , the sum of the Zsh is not the highest in this watershed, but it is the highest in the XG watershed, as shown in Tables 2, 3. We then observed the hot spots distribution of landslides with the two watersheds ( Figure 11) and found that the CDG watershed has many concentrated landslides, but their areas are not as large as the landslides inside the XG watershed. As the Z-score is affected by the landslide area, the larger area leads to a higher Z-score. That is why the Z-score for hot spots of the CDG watershed is lower than the XG watershed. In addition, despite the total number of landslides within the CDG watershed is smaller than in the XG watershed, the proportion of the largescale landslides with the CDG watershed is greater than within the XG watershed, so A c of the CDG watershed is greater. We then analyzed the watersheds with small-scale landslides to determine whether they conform to this rule. The XF watershed and the QL watershed have low landslide density, low A c and ρ values, low sums of Zsh, but the sum of the Zsc are high because of small-scale landslides within the two watersheds are in the majority (Table 3 and Figure 11). We consider that A c of the total landslide within the watershed could not reflect the local spatial clustering degree of landslides. Because the hot spots analysis is more focused on the size and distance of landslides that are clustered together, A c is more focused on the proportion of the large-area landslides within the watershed.  In addition, by analyzing the hot spots, we found that the large-scale landslides are located in the lower part of the watersheds and concentrated near the outlets of the watersheds or the sub-watersheds; the small-scale landslides are located in the middle or upper part of the watersheds and were always far from the outlets, as shown in Figure 11.

DISCUSSION AND CONCLUSION
Research generally considers that the tail of the landslide frequency-area relationship fits the power law, but the relationship between parameters and the spatial distributions has never been analyzed. There is no standard for building the frequency-area, and the range of fitting with power-law is greatly influenced by subjective factors. Therefore, it is difficult to relate the power exponent to the spatial distribution of landslides. We considered landslides in typical watersheds and found that the cumulative frequency satisfies the scaling distribution. The watershed with low density, small area of landslides, and low proportion of the large-scale landslides and clustering degree has small ρ, Ac, and Z-score values. If A c is greater than ρ, then the large-scale landslides are more frequent than the smallscale landslides, whereas when ρ is greater than A c , the landslide distribution characteristics are the opposite. Watershed can have a high Z-score and high density of landslides, yet not necessarily high A c because the Z-score is affected by the landslide area. On the contrary, a watershed with a large proportion of large-scale landslides and high A c yet the small landslides area will reduce the Z-score.
In our research, the cumulative frequency-area of landslides within several watersheds satisfies a scaling distribution that derives two parameters: the exponent ρ and the characteristic area A c . It was found that the two parameters reflect the landslide density and the proportion of the large-scale landslides, respectively, and that the relationships are determined by the data. We then performed a hot spot analysis to investigate the relationship between the parameters and the spatial clustering degree of landslides. From the analysis, we consider that the hot spot analysis emphasizes the size and distance of landslide area, which clusters together. A c emphasizes the proportion of the large-area landslide within the watershed, whereas ρ emphasizes the density of the landslides.
In addition, considering the watershed as a spatial cell rather than a grid is more in line with the natural characteristics of the landslide and can provide more details on the spatial distribution of landslides within a watershed. In this research, the watershed is considered as a factor of landslide division; it is expected that the same method can be easily used for landslide distribution under other influential factors.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.