Size distribution law of earthquake-triggered landslides in different seismic intensity zones

The Ms 8.0 Wenchuan earthquake in 2008 and Ms 7.0 Lushan earthquake in 2013 produced thousands of 10 landslides in the southern region of the Longmen Mountains in China. We conducted field investigations and analyzed remote sensing data to determine the distribution law of earthquake-triggered landslides. The results show a strong negative power-law relationship between the size and frequency of landslides in VII, VIII, and IX seismic intensity zones, a weak power law in the X seismic intensity zone, and a lognormal distribution in the XI seismic intensity zone. Landslide density increases with increasing seismic intensity. A sand pile cellular automata model was built under the conceptual framework of 15 self-organized criticality theory to simulate earthquake-induced landslides. Data from the simulations demonstrate that with increasing disturbance intensity, the dynamical mechanism of the sand pile model changes from a strong power law to a weak power law and then to a lognormal distribution. Results from shaking table experiments of a one-sided slope sand pile show that for peak ground acceleration (PGA) in the range of 0.075 g–0.125 g, the relation between the amount and frequency of sand follows a negative power law. For PGA between 0.15 g and 0.25 g, the relation obeys a lognormal 20 distribution. This verifies that the above-mentioned distribution of earthquake-induced landslides should be a universal law from a physical viewpoint and may apply to other areas. This new perspective may be used to guide development of an inventory of earthquake-triggered landslides and provide a scientific basis for their prediction.

A sand pile is a deceptively simple system that serves as a paradigm for SOC. Held and co-workers (1990) devised a precision apparatus that slowly and uniformly pours sand grain by grain onto a flat, circular surface. The sand pile stops 60 growing when the amount of sand added is balanced on average by the amount of sand that falls off the edge, at which point the system has reached a critical state. When a grain of sand is added to a pile in the critical state, it can initiate an avalanche of any size, including a catastrophic event, and the size and the frequency are related by a power law. SOC has explained the dynamics of many catastrophes such as earthquakes (Olami et al., 1992), forest fires (Drossel and Schwabl ., 1992;Malamud et al., 1998), mountain and rock slides, snow avalanches (Hergarten, 2002(Hergarten, ,2003, solar and stellar flares (Aschwanden, 65 2011a(Aschwanden, 65 , 2016, and stock market crashes (Sornette, 2003). This paper investigates SOC's application to landslides triggered by earthquakes.

Two great earthquakes struck the Longmen Mountains in 2008 and 2013 and triggered numerous landslides. The Longmen
Mountains, on the eastern edge of the Qinghai-Tibetan Plateau in Sichuan Province, China, have the greatest gradient change in the surrounding mountain system. At least 5-10 km of material in the Longmen Mountains has been eroded away since 70 the Miocene, rising at a rate of about 0.6 mm/y (Li Yong et al., 2006;Liu Shugen, 1993). The continuous effects of uplifting and denudation formed a deep and steep valley landscape with an average slope in the southern and Minjiang River areas above 27°. According to Davis's theory of erosion cycles (Davis, 1899), the Longmen Mountains are in the early maturity of geomorphologic evolution and their slopes have evolved to a critical stage. Common characteristic of slope disasters include dilapidation of natural terrain or cut slope, rock avalanches, debris flows, landslides, rock piles, and snowslides, all of which 75 are mainly composed of granular mixtures and associated with energy dissipation owing to material instability and slipping during slope accumulation. The sand pile model is ideal for reflecting the energy dissipation process within slope accumulation, and its dynamic characteristics should therefore be explainable under the conceptual framework of SOC Yan and Huang, 2016).
In this study, we aim to better understand the universal behavior of slope disasters from the point of view of SOC. We 80 studied the distribution of landslides induced by two earthquakes in the Longmen Mountains by field investigations and remote sensing interpretation. A sand pile cellular automata model was also built under the conceptual framework of SOC to simulate earthquake-induced landslides. Shaking table experiments of a one-sided slope sand pile under seismic excitations were conducted. The results of this study provide a physical explanation of the general distribution of earthquake-induced landslides.
2 Landslide magnitude-frequency distributions in different seismic intensity zones

Study area
The 2008 Wenchuan earthquake triggered the largest number of landslides observed in human history. On April 20, 2013, the Lushan earthquake also triggered a large number of landslides. Both events were caused by thrust faulting in the Longmen Mountains of China in areas under similar geographical, geological, and geomorphological conditions. 90

Landslide field investigation
The majority of our field work was accessed by car, landslides were recorded that were visible along the road, and then 95 followed by foot for detailed investigation. The landslide investigation of the Wenchuan earthquake began in the rescue time https://doi.org/10.5194/npg-2020-40 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License. from May 31, 2008 soon after the main shock until June 5, 2008. We collected preliminary survey data by recording the location of the landslide deposits that covered the road, estimating the landslide volume, and asking the highway department staff for the amount of cleaned landslide deposits. We chose to limit the survey to >10 m 3 landslides that blocked roads.
A relatively detailed reconnaissance field study was conducted from August 2 to 27, 2009. After two rainy seasons, the loose 100 material on the slope was nearly washed away, revealing a complete sliding bed or collapsed back wall. We measured the spatial position coordinates of the sliding bed and compared with a 1:2000 topographic map along the highway from Dujiangyan to Yingxiu, which is a section of Chinese national highway 213. A three-dimensional digital figure of the landslide was generated to obtain the landslide volume and depth. According to the slope engineering design scheme of the highway, the landslide depth can be estimated more accurately by measuring the length of exposed anchors after the 105 earthquake.
After the Lushan earthquake, we conducted landslide field surveys along the road from April 26 to August 26, 2013. We accumulated survey routes exceeding 500 km in length, including the China national highway 318, provincial highways 210 and 211, country town roads from Leyin to Rongjin, Longmen to Taiping, Renjia to Shuangshi, Shuangshi to Linguan, and Taiping to Dachuan. We focus on landslide volume, which is mainly obtained from field measurements and the amount of 110 cleaned landslide deposits by the highway department during the rescue time.

Landslide inventory mapping
Landslides triggered by earthquake are usually large in number and widespread in distribution, and detailed coseismic landslide data cannot be based on field investigations alone. At present, the visual interpretation of high-resolution remote sensing images is the primary method to obtain large-area earthquake-induced landslide data. 115 We collected pre-and post-earthquake high-resolution remote sensing data to facilitate the landslide analysis. The Wenchuan earthquake data included TM satellite images (pre-earthquake), ALOS satellite images (10-m spatial resolution, taken on The color and shape characteristics of landslides (e.g., morphology, hue, shadow, texture) of the optical remote sensing images is clearly distinguishable from the surrounding area, especially shortly after the earthquake. Following the criteria proposed by Harp et al. (2011) andXu (2015), we manually mapped the earthquake-triggered landslides using the GIS platform.
Coseismic landslides are easily detected on high-resolutions images (~1-10 m or higher). Many coseismic landslides may overlap in areas of high landslide density, which makes them difficult to uniquely distinguish and subject to the interpreter's discretion. Densely vegetated areas also pose a problem because vegetation is often destroyed during landsliding and 130 "tadpole" traces appear on the image. Special care must be taken to distinguish the landslide from the vegetation damage range and the identified landslide inventory map must be validated to improve the interpretation accuracy. To address this problem, we first removed landslides in areas with slopes less than 20° because landslide events require certain terrain conditions. Because some farmland, bare land, quarries, sand and gravel yards, old landslides, and other traces of human activities are difficult to distinguish from coseismic landslides on the images, we performed detailed field investigations and 135 combined pre-earthquake satellite images to exclude inaccurate assessments. Finally, we converted the projected landslide area identified in the remote sensing image with the slope gradient to obtain the final landslide area.

Analysis and results
Although the distribution of earthquake-induced landslides is controlled by a variety of factors (e.g., slope aspect, faulting, topography, lithology) that have been discussed in several studies, we believe that the most important control factor is 140 seismic intensity. Therefore, in this study, we analyze landslide distribution according to seismic intensity zone. The range of seismic intensity zones is based on official seismic intensity distribution maps published by the China Earthquake Administration.
Frequency-area (or volume) distribution of a landslide event often exhibits power-law scaling over a limited scale range (Guzzetti et al. 2002;Stark 2001). The relationship can be represented by Eq. (1): 145 where a , b is constants, b is the exponent power that describes the statistical features of the landslide distribution. is landslide size, characterized by area and volume, even depth. ( ) is the number of landslides beyond a given .

Landslide volume (depth)-frequency distributions in different seismic intensity zones
We performed detailed field investigations of 10 5 landslides triggered by the Wenchuan earthquake and 261 landslides 150 triggered by the Lushan earthquake. The data were fitted by Eq. (1) where is landslide volume and ( ) is the number of landslides beyond a given volume , ℎ is landslide depth and (ℎ) is the number of landslides beyond a given depth ℎ. The power-law formula is fitted by a least-squares regression where R 2 represents the goodness of fit (i.e., larger R 2 reflect a better fit). The results are listed in Table 1. https://doi.org/10.5194/npg-2020-40 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License. In the VII, VIII, and IX seismic intensity zones of the Lushan earthquake, the landslides volume-frequency distribution follows a strong power-law relationship. In the IX seismic intensity zone of the Wenchuan earthquake, the volume-frequency distribution and depth-frequency distribution all follow a power-law relationship. The X seismic intensity zone of the Wenchuan earthquake shows similar yet weak power-law characteristics, even though the number of landslides is small (29 sites). A substantial number of large landslides are observed in the XI seismic intensity zone of the Wenchuan earthquake 160 and are connected in a single mass, which is not easy to distinguish into individual sites. Almost two-thirds of the mountains show evidence of mountain peeling. Of the 15 landslides surveyed, the minimum volume is 20 m 3 and the maximum is 16800 m 3 . The minimum depth is 0.3 m and the maximum is 5 m. Although the number of samples is insufficient to draw statistical conclusions, large landslide events appear to dominate. A power-law relationship is not observed between the volume (or depth) and frequency in the XI seismic intensity zone. 165

Landslide frequency-area distributions in different seismic intensity zones
We identified 20,254 coseismic landslides of the Wenchuan earthquake and 1608 coseismic landslides of Lushan earthquake by manual image interpretation, as shown in Figs. 2 and 3, respectively. The statistics of the landslide area were fitted by Eq.
(1), where is landslide area and ( ) is the number of landslides beyond a given area . To eliminate the influence of map resolution and undersampling of smaller landslides, we statistically analyze landslides over 1000 m 2 . The power-law formula 170 is fitted by a least-squares regression. We also test the lognormal relationship of the data by the chi-square test method. The results are listed in Table 2. https://doi.org/10.5194/npg-2020-40 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License.  In the VII, VIII, and IX seismic intensity zones of the Lushan earthquake, the landslides area-frequency distribution follows a power-law relationship (R 2 > 0.9), although none pass the lognormal distribution test. The landslides area-frequency distribution follows a power-law relationship (R 2 > 0.9) in the IX seismic intensity zone of the Wenchuan earthquake, a 180 weak power law in the X seismic intensity zone (R 2 =0.873), while a lognormal distribution in the XI seismic intensity zone.

Landslide density in different seismic intensity zones
We define landslide density as the number of landslides per square kilometer. The landslide numbers and calculated densities in the different seismic intensity zones of the two earthquakes are listed in Table 3. Landslide density clearly increases with seismic intensity, in agreement with previous studies (Qi et al., 2010).To test if the statistical results obtained from landslides generated by the Wenchuan and Lushan earthquakes may also apply to landslides in other areas, we combine cellular automaton simulations with laboratory experiments under the conceptual framework of SOC to interpret the distribution mechanism of earthquake-triggered landslides from a physical viewpoint.

Cellular automata simulations 190
Computer models are an integral component of SOC research, mostly from numerical simulations of a sand pile cellular automaton. The Bak, Tang, and Wiesenfeld (1987) model (BTW) is the earliest and most classic model of sand pile cellular automata and many other models have been developed based on the BTW model for different physical systems that display SOC (e.g., earthquakes, forest fires, magnetic vortex motion, interface growth, biological evolution). In this study, we introduce a sand pile cellular automata model to simulate earthquake-induced landslides. 195

Model procedure
According to the physical features of earthquake-induced landslides, the model has following aspects: (1) The BTW model is used to study the avalanche size distribution over time of a sand pile under disturbances. Because earthquake-induced landslides occur simultaneously, the avalanche size distribution of many sand piles must be studied when they are disturbed by a single event. 200 (2) In the BTW model, a system may be partially disturbed and remain strictly energetically conservative. During earthquake-induced landslides, the slope is disturbed as a whole and its instability must overcome its self-stabilization ability, which requires input energy. The disturbed disseminate is therefore energetically non-conservative.
(3) The BTW model is used to study the dynamic behavior of a system under the same perturbation conditions. For earthquake-induced landslides, the disturbance force on the slope is sensitive to the seismic intensity zone and may exceed 205 the perturbation level. We must therefore modify the disturbance intensity to simulate this physical process.
According to the above characteristics, we constructed a sand pile cellular automata model to simulate earthquake-induced landslides. The model is defined on a two-dimensional × lattice. The sites were numbered with a pair of sub-indexes ( , )(1 ≤ , ≤L), and each site has four nearest neighbors located in the upper, lower, left, and right directions. The state of each site is characterized by a non-negative integer variable , , which is the state value that reflects the stability of site 210 ( , ) (equivalent to the site energy) and each site has a threshold ℎ . We introduce a as the disturbed transmission parameter for the four neighbors, which are not larger than 0.25 owing to non-conservative factors. The model is described by the following algorithm.
Step 1. N sand piles of equal size are built simultaneously but with different initial states. For each sand pile, all sites are initialized to a random value between 0 and ℎ , and Fmax is the maximum value of all sites in the sand pile. 215 Step 2. The N sand piles are simultaneously disturbed. The state value of all sites increases uniformly by the disturbance intensity ′ . , → , + ′ Step 3. The N sand piles are investigated individually. If all sites in a given sand pile remain less than ℎ , nothing happens.
Conversely, if , ≥ ℎ , the site ( , ) becomes unstable and relaxes according to the rule: 220 The relaxation may cause some of the neighbors to become unstable. If so, step 3 is repeated all sites are less than ℎ ( , < ℎ ). 225 By changing ′ , we can determine the relationship between avalanche magnitude and occurrence frequency.
It should be noted that when a disturbance is applied, some sand piles will react while others will not depending on ′ . When ′ = ℎ − , at least one sand pile will react. When ′ is small, sometimes only a few sites will be triggered. But as ′ increases, a batch of sites may be triggered, each of which may trigger chain reactions that may ultimately cross paths in space. Parallel processing is therefore adopted in the algorithm and all disturbed sites react simultaneously in a parallel-230 updating manner. We measured avalanche size in terms of the number of sites participating in the relaxation. This property is called cluster size and is as a measure of the area affected by the avalanche.

Results
The model parameters are = 0.2, ℎ = 1, and = 50. One million sand piles (N = 10 6 ) were generated. Each sand pile was continuously reacted 10 5 times with a disturbance intensity of ′ = 1 − in succession to ensure that each sand pile 235 evolves to a critical state before the formal experiment. Eight groups of simulation experiments were then carried out by increasing ′ from 0.00001 to 0.01, with each group reacted only once. Let the number of avalanches be and the frequency of the avalanche size equal to be ( ). The avalanche density is equal to the number of sand piles with avalanche events divided by the total number of sand piles. The statistical results are shown in Table 4.  Table 4 shows that the disturbance intensity 1 − can be divided into two intervals and the dynamic characteristics of the sand pile model exhibit different properties. When ′ < (1 − ), the avalanche scale and occurrence frequency basically obey the same power-law distribution but the avalanche density increases monotonously with increasing ′ (Fig. 4).
When ′ > (1 − ), the dynamics of the sand pile model exhibit a strong to weak power-law relationship and then to a lognormal distribution with increasing ′ (Fig. 5). 245 physical systems that show SOC. To better understand the physical phenomena of earthquake-induced landslides, we 255 performed shaking table tests to study the dynamic behavior of a sand pile under different earthquake forces.

Experimental procedure
A landslide triggered by an earthquake is a natural phenomenon that occurs over a tremendously large size range (~10 2 -10 8 m 3 ). The purpose of the experiment is to study the dynamic behavior of a sand pile and the model sand piles need not simulate a certain prototype. Previous sand pile experiments have shown that the gradation of model material, physical and 260 mechanical parameters, and model size may influence the collapse size, but there is no influence on the relationship between collapse size and its occurrence frequency. We therefore did not consider similarity relations in the tests.
Large-scale shaking table experiments were conducted in the Key Laboratory of High-speed Railway Engineering at Southwest Jiaotong University in China (Fig. 6). The shaking table is a single-direction table with a size of 2 × 4 m, capacity of 25,000 kg, and loading frequency range of 0.4-15 Hz. In the absence of an applied load, the maximum acceleration is 1.2 265 g and the displacement ranges from −100 to 100 mm.
The one-side slope sand pile was built in a steel model box with a 3.75-m length, 1.75-m width, and 2.75-m height placed on the shaking table. The sand pile material was a dried natural sand gravel collected from earthquake-triggered landslides in the Longmen Mountains. Particles larger than 50 mm were removed and the gradation was measured (Fig. 7). When the sand pile reaches its natural angle (e.g., soil angle of internal friction), it is in a critical stable state. The sand pile has a length of 270 2.58 m, width of 1.5 m, height of 1.95 m, and total weight of 6800 kg (Fig. 8).
Slope responses under the excitation of field seismic waves were recorded at the Wolong seismic station of the 2008 Wenchuan earthquake (WL wave, Fig. 9). The input WL wave was proportionally scaled to its peak value.
To study the variation of sand pile dynamic characteristics with increasing disturbance intensity, we designed five sets of tests with input peak accelerations of 0.075 g to 0.250 g. After inputting the excitation, the weight of the sand gravel collapse 275 was measured as a test value. After each test, sand gravel was added from the top of the slope to ensure that the slope remained in a critical stable state. Each set of tests was run no less than 60 times.

Data analysis of result
Let the collapse weight be and the collapse weight frequency equal to be ( ). The analysis results are shown in Table 5.
The collapse density is equal to the number of tests with collapse events divided by the total number of tests. Note: When the PGA is 0.075 g, 0.1 g, and 0.125 g, some tests occurred without collapse events and the number of tests was increased. The other groups were repeated 60 times.
When the peak acceleration input was between 0.075 g and 0.125 g, some tests occurred without collapses. When collapses did occur, small collapses were significantly more common than large collapses and the results obey a power-law distribution. When the PGA was between 0.15 g and 0.25 g, no tests occurred with zero collapses, the power-law 295 relationship was weakened, and the results followed a lognormal distribution. The collapse density increases with increasing peak acceleration.
Studying the dynamic response of sand pile with seismic waves as disturbance sources is a unique experimental method to study SOC characteristics. Previous results have shown that the disturbance mode does not affect the sand pile dynamics characteristics, but does influence the power-law relationship parameters (Yao et al.,1998Yang et al., 2007). The 300 shaking table sand pile model tests show that changes in disturbance intensity lead to a shift in system dynamics. The physical process of the shaking table sand pile test is close to the prototype problem of earthquake-induced landslides, even though the number of experiments remains limited, and provides good support for the universality of earthquake-triggered landslides in different intensity zones.

Conclusions and discussion 305
(1) We analyzed data from landslides triggered by the 2013 Lushan and 2008 Wenchuan earthquakes. The results show a negative power-law relationship between landslide size and frequency in the VII, VIII, and IX seismic intensity zones. The relationship becomes a weak power law in the X seismic intensity zone and changes into a lognormal distribution in the XI seismic intensity zone. Landslide density increases gradually with increasing seismic intensity. Cellular automaton simulations reveal that with increasing disturbance intensity, the dynamical mechanism of the sand pile model changes from 310 a strong power law to a weak power law and then to a lognormal distribution, and the avalanche density increases. The results of the shaking table sand pile model tests verify these findings. The overall landslide distribution law is therefore constrained, even though these landslides are complex and very random, and the evolution mechanism of the distribution law in different intensity zones is clarified. The distribution probability model of earthquake-triggered landslides and evolution model with increasing seismic intensity presented here exceed the statistical relation level of a typical sample, which is 315 therefore a universal law.
(2) SOC was founded in 1987 as a branch of non-equilibrium thermodynamics and study of its phenomenology and precise definition continues. SOC has highlighted that thresholds, metastability, and large-scale fluctuations play a decisive role in the spatiotemporal behavior of a large class of multi-body systems. However, the influence of disturbance intensity on the system dynamics behavior of SOC has not received much attention. The disturbance intensity range of a catastrophic event 320 may span several orders of magnitude (e.g. the energy difference between a magnitude 9 earthquake is 32,768 times more powerful than a magnitude 6 earthquake). In the case of the 2008 Wenchuan and 2013 Lushan earthquakes, the dynamic characteristics of the SOC system undergo a strong power law, weak power law, and lognormal distribution evolution process with increasing disturbance intensity. This constraint on the evolution pattern of a SOC system behavior goes beyond the traditional field of SOC and makes a strong impetus to further develop basic SOC theory. 325 (3) Compared with chaos theory, SOC does not greatly emphasize the initial conditions and system details, which facilitates the experimental and analytical procedures. The sand pile example tends to explain nearly everything from mountain formation to stock market volatility. But if many of the unique details of mountain systems can be understood by simple cellular automata numerical simulations, it may not be realistic for most geographers and requires further validation. For example, slope structures (e.g., joints, fracture surfaces) in nature are non-uniformly distributed and natural granular 330 materials often exhibit a wide gradation. The nonuniformity of components is one of the important characteristics of a mountain system, however, the influence of non-uniform components on the dynamics is not considered here. We briefly discuss the effect of nonuniformity on dynamics in a previous study (Guo et al., 2017). However, the heterogeneity of cell geometry, arrangement randomness, and interaction anisotropy can be considered to investigate the unique details of how nonuniformity property affects system dynamics. Further research will therefore aim to determine the deep-seated law of 335 earthquake-induced landslides.
(4) In the IX seismic intensity zone, the cumulative number-area distribution of landslides triggered by the earthquakes exhibits a negative power-law relationship but with different power exponents: 1.543 for the Wenchuan earthquake and 0.955 for the Lushan earthquake. Whether the observed variability in the power-law exponents is incidental or applies only under certain conditions should be investigated in future studies. 340