Prediction of the area affected by earthquake-induced landsliding based on seismological parameters

. We present an analytical, seismologically consistent expression for the surface area of the region within which most landslidestriggered by an earthquake are located (landslide distribution area). This expression is based on scaling laws relating seismic moment, source depth and focal mechanism with ground shaking and fault rupture length and assumes a globally constant threshold of acceleration for onset of systematic mass wasting. The seismological assumptions are identical to those recently used to propose a seismologically consistent expression for the total volume and area of landslides triggered by an 5 earthquake. To test the accuracy of the model we gathered geophysical information and estimates of the landslide distribution area for 83 earthquakes. To reduce uncertainties and inconsistencies in the estimation of the landslide distribution area, we propose an objective deﬁnition based on the shortest distance from the seimsic wave emission line containing 95% of the total landslide area. Without any empirical calibration the model explains 56% of the variance in our dataset, and predicts 35 to 49 out of 83 cases within a factor two, depending on how we account for uncertainties on the seismic source depth. For most 10 cases with comprehensive landslide inventories we show that our prediction compares well with the smallest region around the fault containing 95% of the total landslide area. Aspects ignored by the model that could explain the residuals include, local variations of the threshold of acceleration and processes modulating the surface ground shaking, such as the distribution of seismic energy release on the fault plane, the dynamic stress drop and rupture directivity. Nevertheless, its simplicity and

Abstract.We present an analytical, seismologically consistent expression for the surface area of the region within which most landslides triggered by an earthquake are located (landslide distribution area).This expression is based on scaling laws relating seismic moment, source depth, and focal mechanism with ground shaking and fault rupture length and assumes a globally constant threshold of acceleration for onset of systematic mass wasting.The seismological assumptions are identical to those recently used to propose a seismologically consistent expression for the total volume and area of landslides triggered by an earthquake.To test the accuracy of the model we gathered geophysical information and estimates of the landslide distribution area for 83 earthquakes.To reduce uncertainties and inconsistencies in the estimation of the landslide distribution area, we propose an objective definition based on the shortest distance from the seismic wave emission line containing 95 % of the total landslide area.Without any empirical calibration the model explains 56 % of the variance in our dataset, and predicts 35 to 49 out of 83 cases within a factor of 2, depending on how we account for uncertainties on the seismic source depth.For most cases with comprehensive landslide inventories we show that our prediction compares well with the smallest region around the fault containing 95 % of the total landslide area.Aspects ignored by the model that could explain the residuals include local variations of the threshold of acceleration and processes modulating the surface ground shaking, such as the distribution of seismic energy release on the fault plane, the dynamic stress drop, and rupture directivity.Nevertheless, its simplic-ity and first-order accuracy suggest that the model can yield plausible and useful estimates of the landslide distribution area in near-real time, with earthquake parameters issued by standard detection routines.

Introduction
Triggered landslides are a significant secondary hazard of earthquakes, and may be the dominant cause of damage to infrastructure and lifelines, especially roads (Bird and Bommer, 2004).The severity of this hazard and the associated risks is clear after most large earthquakes in steep landscapes, and was underlined by the devastation and fatalities caused by landsliding induced by recent large earthquakes in Sichuan (China) 2008 and central Nepal, 2015 (Yin et al., 2009;Kargel et al., 2015).The earthquake-induced landslide hazard is defined in the first instance by the number, size, and location of landslides.These variables are correlated with a combination of local factors, such as the peak ground acceleration (Meunier et al., 2007(Meunier et al., , 2008)), hillslope geometry (Parise and Jibson, 2000;Gorum et al., 2013), and the strength (Parise and Jibson, 2000;Gallen et al., 2015) and degree of hydrological saturation of near-surface materials, which are difficult to quantify due to their inhomogeneity across epicentral areas (Dreyfus et al., 2013).A simpler approach is to predict first-order variables such as the total volume and area of landsliding caused by an earthquake based on simple seismological considerations (e.g.O. Marc et al.: Prediction of the area affected by earthquake-induced landsliding Marc et al., 2016b).If the input parameters for such a model can be quantified in minutes to hours after an earthquake, then this could yield a quick insight into the scale of the affected area and the total amount of landsliding within.Here we focus on the landslide distribution area, A d , that is the surface area of the region within which landslides induced by a given earthquake are likely to be concentrated.This is an important risk parameter, as it defines the zone within which the landslide hazard is focused.It can be intersected with areas of vulnerability for an assessment of risk and defines the extent of the area experiencing seismically induced hillslope denudation, with geomorpohological (e.g.Marc et al., 2016a), geochemical (e.g.Jin et al., 2015), tectonic (e.g.Steer et al., 2014), and biological consequences (e.g.Garwood et al., 1979).Further, the area with concentrated landsliding should be more closely related to seismic forcing and less dependent on unusual hillslopes conditions.
Several global or regional compilations of earthquakeinduced landslide data have reported A d , and explored its scaling with earthquake magnitude (Keefer, 1984;Rodriguez et al., 1999;Hancox et al., 1997;Keefer, 2002).They found that log 10 (A d ) scales linearly with the moment magnitude M w of an earthquake: where M o is the earthquake seismic moment.Because the scatter around the central trend of this relation is substantial, a common approach is to base a prediction of A d on the envelope defining the maximum area affected by landsliding for a given earthquake magnitude.Further, some case studies have illustrated how specific seismic or geomorphic conditions can lead to landslide triggering over exceptionally long distances and therefore large areas (Keefer, 2002;Jibson and Harp, 2012).However, the definition of the relationship between A d and seismic moment and of the influence of other seismological and geomorphic parameters has a slender theoretical basis.Recently, expressions for the total area and volume of landslides triggered by an earthquake have been derived based on seismological scaling laws and simple topographic characterizations (Marc et al., 2016b).
Here we show that this same treatment may be used to predict the shape and size of the landslide distribution area.First, we present the basis for an expression of the landslide distribution area and the landslide maps and compilations of estimated landslide distribution areas against which model predictions can be compared.Then we assess the validity and accuracy of our theoretical approach and discuss its limitations, and finally we suggest directions for future improvements.

A seismologically consistent expression for the landslide distribution area
Earthquakes trigger landsliding due to transient accelerations during ground shaking, which can shift the force balance in a slope and cause damage in the substrate, reducing its cohesion and resistance to failure (Newmark, 1965).This link between landsliding and ground shaking, specifically peak ground acceleration (Wilson and Keefer, 1985;Khazai and Sitar, 2004), is confirmed by a growing number of detailed observations (Meunier et al., 2007;Yuan et al., 2013).They constrain, for example, the statistical occurrence of landsliding within areas where the ground acceleration exceeds a threshold acceleration, a c ∼ 0.1-0.2(Meunier et al., 2007;Hovius and Meunier, 2012;Yuan et al., 2013).Note that, here, a c and all other accelerations are normalized by the gravitational acceleration and thus non-dimensional.Below this threshold, landsliding is rare or minor, and therefore a reasonable estimation of the size and shape of the landslide distribution area could be found by intersecting the region where the peak ground acceleration exceeds a c and the region with sufficiently steep topography for landsliding to occur.Marc et al. (2016b) have successfully modelled the total volume and area of the population of landslides due to a given earthquake, using seismological scaling laws to constrain the magnitude and extent of ground shaking.They assumed that a c is a threshold acceleration for damage above which an effective reduction of strength occurred within the hillslope materials, increasing the likelihood of failure of slopes steeper than the angle of internal friction and thus stable due to a significant cohesion.Therefore in this paper, a c is not the classical critical acceleration based on the safety factor and the slope gradient, but is a threshold acceleration for damage.With this definition we consider that a c is much more uniform across a landscape than cohesion which can vary greatly between soil and fractured and intact bedrock.This is consistent with modest variations of the minimum acceleration at which landsliding (Meunier et al., 2007;Hovius and Meunier, 2012;Yuan et al., 2013) and minor rockfall (Jibson and Harp, 2016) have been observed.The average value of a c across a landscape is an important constraint on A d .At the outset we assume a c = 0.15, which is conservative and allows a focus on areas with significant landsliding.Note that, in this approach, a c is independent of the slope gradient, but also that a slope experiencing an acceleration larger than a c will fail only if the resulting drop in cohesion renders it unstable.Thus the number and size of landslides on a hillslope will also depend on local strength, pore pressure, and topographic steepness, but we assume that where ground acceleration reaches a c some minor failures will initiate.Further, their model considers attenuation of seismic waves due to geometric spreading, with S-wave amplitude decreasing with distance between the earthquake source and the affected topography.The pattern and intensity of ground shaking are modelled as the superposition of Nat.Hazards Earth Syst.Sci., 17, 1159Sci., 17, -1175Sci., 17, , 2017 www.nat-hazards-earth-syst-sci.net/17/1159/2017/ point sources associated with asperities along the fault rupture.Therefore, the source of wave emission is considered to be distributed along the fault rupture length, while the wave acceleration at the source, that sets the distance over which waves can travel before becoming insufficient to trigger landslides, scales non-linearly with the magnitude of the earthquake.With the same assumptions we can predict the area affected by ground shaking exceeding a c , A s : where R HMAX is the horizontal distance from the surface projection of the wave source at which the ground shaking reaches a c .Assuming that wave attenuation is dominated by geometric spreading and that non-linear attenuation can be , with b the nearsource acceleration at a reference distance R ref = 1000 m from the seismic source and R 0 the mean depth of seismic wave emission.The assumptions are justified near the fault where most landslides occur (Marc et al., 2016b), but additional non-linear attenuation would result in a reduction of the predicted R HMAX and A s .The mean depth of emission is assumed to be the mean asperity depth because asperities emit most of the high-frequency waves (Ruiz et al., 2011;Avouac et al., 2015) and seem to explain best the observed patterns and amounts of landsliding (Meunier et al., 2013b;Marc et al., 2016b).Note that waves with frequencies of 0.5 to 10 Hz are often the most important for landslide triggering because they have wavelengths ranging from the landslide size to the hillslope size (Marc et al., 2016b).
Following Marc et al. (2016b), we use the scaling of fault rupture length, L, with seismic moment proposed by Leonard (2010): with µ the shear modulus, assumed to be 33 GPa, and C 1 and C 2 empirical constants.Although Leonard (2010) Leonard (2010) and reported above.
To construct a simplified ground shaking model that is consistent with first-order seismological observations and landslide triggering we use the scaling of the near-source acceleration, b, with earthquake magnitude proposed by Boore and Atkinson (2008) where M h = 6.75 is a hinge magnitude above which the acceleration carried by seismic waves saturates near b sat , and e 5 = 0.6788, e 6 = −0.1826,and e 7 = 0.054 are empirical constants for 1 Hz pseudo-spectral accelerations (Boore and Atkinson, 2008).We prefer to use the scaling of pseudospectral accelerations (rather than the one for peak ground accelerations for example) because it focuses on ground motion around a specific frequency, and therefore should better represent the frequency-dependent modulation of the ground shaking.The mechanical assumptions associated with pseudo-spectral accelerations (elastic oscillator with base fixed on the ground and 5 % damping) are probably not entirely adequate for hillslope failure, but not inconsistent with the resonance of topographic ridges themselves (Meunier et al., 2008).We refrain from using b sat values as derived by Boore and Atkinson (2008) because they should be included within a model accounting for site effects and nonlinear attenuation, attributes that are beyond the scope of this work.Therefore, to be consistent with the model of Marc et al. (2016b), we use b sat = 0.4, yielding normalized surface acceleration of about 0.4-0.8above asperities located at 5-10 km depth, consistent with observations during earthquakes with M w > M h .Fault type may influence both rupture length and ground shaking.For the shaking term, we follow the model choices of Marc et al. (2016b), attributing 30 % less shaking for normal fault earthquakes than for equivalent strike-slip or reverse slip events (Boore and Atkinson, 2008).Large earthquake ruptures on strike-slip faults often break the seismogenic layer over its entire depth, fixed to 17 km in our model (Leonard, 2010).Therefore, above a critical magnitude the increase of fault rupture area with seismic moment is only due to fault rupture length and the scaling exponent between moment and rupture length becomes 2/3.Although, a similar scaling break may exist in principle for dip-slip faults it is less clearly observed (Leonard, 2010) and therefore not prescribed in our model.With these considerations, both the fault rupture length and near-source wave amplitude can be estimated from the seismic moment and fault type.Thus Eq. ( 2) can be rewritten more explicitly and corrected for the presence of topography with a topographic index, C topo , to obtain the predicted area affected: Hence, the landslide distribution area can be predicted from the earthquake moment and mechanism, and the mean asperity depth for any earthquake.If a c and b sat are relatively constant as suggested by Marc et al. (2016b) then the model is fully constrained without any free parameter.The assumptions that a c and b sat are constant across all settings are discussed in Sect. 5. Note that C topo cannot be computed as the fraction of A s within which slopes are less than 10 • .This is because local flats will impede landsliding and change the landslide density locally, but not necessarily reduce the envelope containing all earthquake-induced landslides, which typically defines A d .Only large flat areas, extending beyond R HMAX and with a width similar to the fault length will matter.Typically the presence of a basin or inundated areas along the entire fault will halve the landslide distribution area (C topo = 0.5).
In the model, the critical seismic moment above which A s assumes a non-zero value is modulated by R 0 and ranges between 10 16 and 10 19 N m (Fig. 1).Above this critical moment, A s rises sharply, driven by the exponential increase of the source acceleration (b) with increasing moment (Eq.4).Upon reaching the hinge magnitude, M h = 6.75, b saturates and A s increases primarily due to increase of the fault rupture length (L) with moment.Therefore, for these large events A s scales as a power law of the moment, with an exponent of 2/5 for dip-slip events, and an exponent of 2/3 for strikeslip events (Fig. 1).For large events, b/a c is about 27 km and therefore R HMAX is almost independent of R 0 unless the latter reaches depths >∼ 20 km.Thus uncertainties on R 0 , which is the least well-constrained input variable of Eq. (5), will not substantially affect predictions of the landslide distribution area for most large, shallow earthquakes that rupture the upper crust (Fig. 1).However, for some intermediate-size earthquakes or earthquakes deeper than 25 km, the prediction may vary dramatically due to minor changes in the estimated source depth.

Data and methods
Our model predicts the landslide distribution area A d , and its first-order shape.However, for many earthquakes we have poor constraints on the shape and definition of landslide distribution area.Therefore, we test the model in two different ways.First, we use a large database of earthquakes containing geophysical information and a, sometimes crude, estimate of A d , to assess whether the scaling in the model matches the data and appears to yield a correct first-order prediction of A d .Secondly, we focus on a limited number of cases for which we have a detailed landslide inventory that allow us to define a simple and objective parameter to characterize A d and to compare it quantitatively to the model prediction.  1) are shown for each earthquake, followed by S or N for strike-slip and normal faults, while all other cases are reverse fault earthquakes.R 0 is the mean depth of wave emission.The prediction of the seismologically consistent model is shown for reverse, normal, and strike-slip faults at different depth with solid black, solid grey, and dashed lines, respectively.Inset: distribution of the uncertainty factor (upper estimate/best estimate or best estimate/lower estimate).

Landslide maps and compilations of landslide distribution area
To assess if our theoretical framework captures the observed scaling between A d and Mo, we test it against a large database of 83 crustal earthquakes, for which magnitude and location can be reasonably well constrained.These 83 cases have been harvested from published compilations (Keefer, 1984;Hancox et al., 1997;Rodriguez et al., 1999;Bommer and Rodriguez, 2002;Martino et al., 2014) and from recent landslide maps (Table 1).They include the 10 cases with comprehensive landslide inventories described separately below, 36 inventories for which we could access one or several maps with isolines of landslides density or point inventories to check the values reported in published compilations (Bonilla, 1960;Keefer et al., 1980;Harp et al., 1984;Harp and Keefer, 1990;Jibson et al., 1994;Tibaldi et al., 1995;Hancox et al., 1997;Keefer and Manson, 1998;Hancox et al., 2003Hancox et al., , 2004;;Jibson and Harp, 2006;Mahdavifar et al., 2006;Sato et al., 2007;Kamp et al., 2008;Mosquera-Machado et al., 2009;Alfaro et al., 2012;Collins et al., 2012;Jibson and Harp, 2012;Gorum et al., 2014;Martino et al., 2014;Xu et al., 2014aXu et al., , b, 2015;;Martha et al., 2016;Zhou et al., 2016), and a further 37 cases for which we could not access any raw data to evaluate the reported values (Table 1).For 10 earthquakes, detailed landslide inventories with comprehensive maps of the landslide as polygons are available, allowing an objective characterization of A d (as discussed below): the 1976 Guatemala, 1991Limon, 1993Finisterre, 1994Northridge, 1999Chi-Chi, 2004Niigata, 2007Aysen, 2008Iwate, 2008 Wenchuan and 2010 Haiti earthquakes, ranging from M w 6.2 to 7.9, and with hypocentre depths between 3 and 25 km (Table 1), (Harp et al., 1981;Harp and Jibson, 1996;Liao and Lee, 2000;Yagi et al., 2007;Meunier et al., 2008;Yagi et al., 2009;Gorum et al., 2013Gorum et al., , 2014;;Xu et al., 2014c;Marc et al., 2016b).Most of these inventories were produced from extensive imagery and include all landslides that could be detected at the full image resolution.Two inventories deviate from this.The 1976 Guatemala inventory is based on high-resolution air photos, but only covers a limited area containing the most intense landsliding.Published values of A d for this case are larger than our estimate from the landslide inventory, considered as a lower bound.However, the inventory for Guatemala is considered to contain more than 90 % of the landslides triggered by the 1976 earthquake (R. W. Jibson, personal communication, 2013).Landslides triggered by the 1991 Limon earthquake were mapped across a wide-swathe Landsat-5 image and the limit of the disturbed areas could be constrained, but the low image resolution (30 m) did not allow the delineation of all individual landslides in the most intensely af-fected area.Therefore, the distribution area A d is probably not significantly underestimated, but the cumulative surface area of landslides within it is, and any calculations based on that measure may be biased (see Marc et al., 2016b).We note that for some earthquakes such as the Wenchuan and Haiti events, several inventories are available.For the Wenchuan earthquake the inventory by Xu et al. (2014c) seems the most comprehensive and robust, compared to earlier mapping.For Haiti we analyse and compare the inventories of Gorum et al. (2013) and Harp et al. (2016), which have different interpretations in some areas, likely due to differences in the preand post-earthquake imagery used.Most of the 10 landslide inventories are not complete because they are limited by the resolution of the available imagery, and do not include very small landslides and rockfalls that can be detected in the field (e.g.Jibson and Harp, 2012).As a result these inventories can yield A d estimates smaller than those from previous compilations, but our estimates may be more representative of the area affected by dense landsliding, which is of primary interest in terms of both hazard and erosion.For most earthquakes we have information about moment magnitude, hypocentral depth, and focal mechanism from the international seismological centre catalogue (Storchak et al., 2013), but no estimate of the mean asperity depth most relevant to describe the mean wave emission depth (Table 1).However, we have shown that for large earthquakes, A dp is not very sensitive to depth, while for small earthquakes we expect the hypocentral depth and mean asperity depth to be close to each other.
C topo was crudely estimated based on the topography and the fault position, and is about 1 for most cases (68 out of 83), about 0.5 for 10 events (coastal/basin geometry) and less than 0.5 for five cases (Table 1).For a small fraction of cases in our database we could find a stress drop estimate (Allmann and Shearer, 2009;Baltay et al., 2011) (Table 1), which must correlate with the wave emission at the source and thus the ground shaking intensity (Hanks and McGuire, 1981;Baltay and Hanks, 2014).All stress drop estimates were converted to an equivalent dynamic stress drop as defined by Brune (1970).Subduction earthquakes and distant offshore earthquakes were ignored because the area affected by strong shaking is mostly submerged and hillslopes are only present at large distance where the shaking intensity may not be well approximated by our simple model (see Marc et al., 2016b).

An objective definition of the landslide distribution area
It is important to secure consistency between estimates of A d from different sources and to constrain the degree of uncertainty associated with these estimates.Commonly, the landslide distribution area is estimated by locating all landslides caused by an earthquake as accurately and comprehensively as possible, and drawing a single, smooth envelope containing all landslides, regardless of possible variations of landslide density within it (Keefer, 1984;Hancox et al., 1997;Rodriguez et al., 1999).However, many published inventories are limited by the spatial extent and quality and resolution of the available imagery, and may not include the small landslides in the far field, which, if taken into account, would give rise to a much greater A d .For example, in the case of the 2008 Wenchuan earthquake, accounting for sparse landsliding in the far field would give A d ∼ 200 000 km 2 instead of A d ∼ 44 000 km 2 (Xu et al., 2014c).It is also likely that for large earthquakes with widespread landsliding, there is a tendency to focus on the most intensely affected areas, while for smaller earthquakes, the extent over which small rockfalls occur may be investigated in greater detail through field investigations (Jibson andHarp, 2012, 2016).For small to medium earthquakes triggering a relatively limited number of landslides erroneous inclusion of landslides triggered by other processes just before or after the main shock may also cause a significantly upward bias of A d estimates.In most cases we lack the information required to assess the accuracy and consistency of A d estimates between different events, but we note that, using the common method described just above to determine A d from published maps or detailed inventories, we could reproduce within ∼ 20 % A d estimates reported in global compilations and citing the same source study (Keefer, 1984;Rodriguez et al., 1999).For 27 earthquakes, we found different estimates of A d , from different publications, methods, and/or source imagery.These include different A d reported when considering only the area affected by intense landsliding, or including more distant, sparse landsliding (e.g.Hancox et al., 1997;Xu et al., 2014c).This is not an adequate quantification of the uncertainties in the dataset as a whole, but serves to illustrate how estimates of A d may vary or be biased (Fig. 1

inset).
We propose a robust, alternative approach to define A d , based on the fault rupture and the landslide inventory.In this approach, A d is defined as the surface area of the region within a distance R 95 from the seismic wave emission line projected at the surface.This region is set to contain 95 % of all landslides triggered by an earthquake, as measured by their surface area.In this definition, R 95 is a 1-D measure of the spread of landsliding away from the seismic source.The source is modelled as a wave emission line (or series of lines), defined by the location and length of the earthquake rupture, as determined by geophysical inversion of the slip distribution (or moment distribution for the Guatemala earthquake) on the seismogenic fault plane (Figs. 3, 4), (Kikuchi and Kanamori, 1991;Wald et al., 1996;Zeng and Chen, 2001;Hikima and Koketsu, 2005;Hayes et al., 2010;Suzuki et al., 2010;Fielding et al., 2013).The rupture length of such inversions agrees (within 30 % or less) with the a priori predicted length, based on the scaling proposed by Leonard (2010), for the Chi-Chi, Haiti, and Northridge earthquakes, but is smaller (by 40-50 %) for the Niigata and Iwate earthquakes, and larger for the Wenchuan and Guatemala earthquakes (Table 1).This results in uncertainties when attempting to measure or predict R 95 without knowledge of the rupture length and position, for example for old earthquake or just after an earthquake has occurred.For the Finisterre case we only have the position of the epicentre, the fault strike and the aftershock locations.They define a long fault rupture with the main shock, M w 6.9, being closer of the northwestern fault tip and a large secondary shock, M w 6.5 farther east (Stevens et al., 1998).Accordingly, we defined two separate emission lines, both with an epicentre located one-third of the rupture length from the respective tips (Figs. 3, 4).For the Limon and Aysen earthquakes we placed the emission line centred on the maximum of moment emission and epicentre, respectively (Goes et al., 1993;Legrand et al., 2011), and oriented according to the focal mechanism.In the latter case, choosing the alternative focal mechanism (90 • rotation) would not change significantly our results.
Our treatment differs from previous definitions of the maximum distance for landsliding (e.g.Keefer, 1984;Hancox et al., 1997;Rodriguez et al., 1999), because we consider a seismic line source that may be offset from the surface rupture of the seismogenic fault, and because it is not based on individual detected landslides but on the distribution of landsliding.Advantages of using the R 95 criteria as compared to previous approaches are its objectivity and reproducibility, and its robustness to the accidental addition or omission of minor landsliding in the far field.Thus, R 95 is strongly related to the seismic forcing and still representative of the area where hazard and erosion are likely to be most signifi- cant.A drawback is that this approach requires polygon inventories.Using 95 % of the landslide number as a threshold for point-based inventories would be an adequate solution but this definition would still be quite sensitive to effects affecting the number of identified landslides such as the imagery resolution and/or amalgamation of adjacent smaller landslides (Marc and Hovius, 2015).Another drawback is that R 95 assumes equal rates of decrease of landslide density with distance from the emission line in all directions, which may not always be the case.

Comparison of observed and predicted landslide distribution areas
Most of our landslide distribution area data fall within the range of theoretical predictions from Eq. ( 5), based only on fixed global parameters and variable earthquake settings (Fig. 1).In Fig. 2, the predicted distribution areas, A dp , calculated accounting for the hypocentral depth and adjusting for the abundance of hillslopes where landslides may occur (> 10 • ) (Marc et al., 2016b), are plotted against values estimated from observations, A d , for the earthquakes in our database.For 42 % of 83 earthquake cases, Eq. ( 5) yields A dp within a factor of 2 of A d when considering the hypocentral data as the exact emission depth (R 2 = 0.56).This increases to 59 % of the events when setting the emission depth within 25 % of the hypocentral depth (Fig. 2, inset).Landslide distribution areas vary in size between 10 and 10 5 km 2 , but half of the earthquakes have A d between 10 3 and 10 4 km 2 .In this range, the predictions are mostly within a factor of 2 from the estimated area.When constrained, uncertainties on A d estimates are within a factor of 2 and 4 for about two-fifths and four-fifths of the well-constrained cases, respectively, but occasionally the ratio between A d estimates from different sources may reach up to a factor of 10 (Fig. 1 inset).For a number of small-to moderate-magnitude earthquakes with a small A d , the model predicts no landsliding.This may be due to uncertainties on the hypocentral depth estimation that we assume to be the depth of wave emission, R 0 , and assigning a < 25 % uncertainty on R 0 allows for the prediction to match A d within a factor of 2, for five out of seven cases in this category (Fig. 2).Most earthquakes with A d > 10 4 km 2 have M w > 7.5 and are shallow enough (R 0 < 20 km) to be relatively insensitive to depth.Nevertheless, they are often underpredicted by our model by a factor 1.5 to 3 (Fig. 2).Notwithstanding these observations, the global distribution of errors does not correlate significantly with the seismic moment nor with the hypocentral depth or the focal mechanism of the earthquakes (Fig. S1 in the Supplement).
For comparison, an empirical fit of A d against the seismic moment for the earthquakes in our database has an accuracy of and overall scatter similar to that of our model predictions (Fig. 2, inset).Nevertheless, the fact that Eq. ( 5), based on physical considerations and computed without any free parameter, has this accuracy for a global catalogue suggests that our approach captures essential aspects of earthquakeinduced landsliding.Further validation of the model is inherently limited by the uncertainties associated with the estimation of A d and the inconsistency of the reported values for individual cases.However, additional insights into the validity and limitations of the model may be gained by comparing its prediction to objective landslide distribution areas obtained from well-constrained earthquakes.This is done in the next section for a limited number of comprehensive inventories where the landslide distribution is constrained in detail.

Model comparison against an objective measure of the landslide distribution area
To quantify the error of the model we evaluate the proportion of the total area affected by landsliding located within the A dp perimeter predicted by the model.We also consider the difference between the radius of the area affected by landsliding in the model, R HMAX , and R 95 defined earlier as the distance from the seismic wave emission line within which 95 % of the area affected by observed landsliding is located (Figs. 3, 4).For the 10 earthquakes for which we have comprehensive landslide inventories, the model distribution area always contains between 88 and 100 % of the cumulative surface area of all mapped landslides (Table 2).These numbers indicate that the model always captures the region of most intense landsliding, but that it sometimes overpredicts the affected area, with its limit well beyond the outermost mapped landslide in most directions, as for the Niigata and Iwate cases  (Hikima and Koketsu, 2005;Suzuki et al., 2010).Additionally, the distribution area proposed by Harp et al. (1981) for the Guatemala earthquake is shown in purple.For the Limon earthquake, the total area of landslides is likely underestimated in the most affected area and therefore R 95 is likely reduced and more similar to the dashed red contour.
(Fig. 4).The difference in 1-D radius gives a more accurate view of the merits and limits of the model (Table 2).For eight cases, the Northridge, Limon, Haiti, Aysen, Finisterre, Guatemala, Wenchuan, and Chi-Chi earthquakes, R 95 is well predicted, with an absolute error < 6 km, that is within ∼ 20 % of R 95 in all cases (Fig. 5).However, we note that the Limon inventory is incomplete in the most affected area, sug- Table 2. Summary of the emission length L, the maximum horizontal radius for landsliding R HMAX , the landsliding included in the prediction (% of A T ), and the distance from the emission line containing 95 % of the total landslide area, R 95 .In parentheses we give the values of L derived from rupture length scaling with moment (Eq. 3) that were used when we could not access a rupture slip model.For the Haiti earthquake we show the results from the inventory of Gorum et al. (2013) and of Harp et al. (2016) (in parentheses).We use the two largest shocks to model the Finisterre event.* For the Limon earthquake some landslides in the most affected area could not be delineated and the total area is underestimated.Therefore we likely underestimate the amount of landsliding within the model prediction and overestimate R 95 , that may actually be smaller than R HMAX .(Gorum, et al. 2013) Figure 5. Relative and absolute errors in the prediction of the distance from the wave emission line containing 95 % of the total landslide area plotted against the seismic moment.Horizontal black lines delimit cases where the relative error is within < 20 %.The absolute error (red circles) indicates the difference in kilometres between the prediction and the boundary of 95 % of the landsliding.

EQ
gesting that R 95 may actually be overestimated and R HMAX may exceed it (Fig. 4).Further, for the Haiti earthquake, the model underpredicts R 95 by about 35 % (9 km) if we consider Harp's 2016 inventory which extends far into little-affected areas (Fig. 3).For the two remaining cases the error is much larger, with an overprediction of about a factor of 2, and an absolute error of approximately 10 km, for the Niigata and Iwate earthquakes (Figs. 4,5).Additionally, we note that the agreement between R 95 and our model may sometimes hide important along-strike variations or trends, as for the Guatemala and Northridge earthquakes (Figs. 3 and 4).The possible reasons for such trends and other limitations of the model that could explain why some cases are substantially under-or overpredicted are explored in the discussion, below.
We stress that the 5 % of the total landsliding outside of R 95 may entail a significant hazard which can extend much further, especially for large earthquakes like the Wenchuan, Chi-Chi, or Guatemala cases which triggered huge amounts of landslides (total area of 1160, 128, and 61 km 2 , respectively).In the Wenchuan case, the distance from the emission line required to encompass 97.5 % of the total landslide area is 48 km instead of 34 km for R 95 .These landslides are more sparsely distributed and most of the time smaller than the ones close to the fault (Valagussa et al., 2017), but they remain difficult to predict in many cases.

Discussion
We have shown that a first-order a priori seismic shake map coupled with a universal shaking threshold for landsliding can reproduce reasonably well the landslide distribution areas in a compilation of 83 cases, and that it matches the surface area encompassing 95 % of the total landslide area for most of the cases for which we have comprehensive landslide inventories.In this section, we identify and try to quantify the different sources of uncertainties and potential ways to improve the model.

Asperity length and wave emission along the fault
For the Niigata, Iwate, and Nagano earthquakes, the difference between the observed and predicted area of landsliding would be greatly reduced if we treated the earthquake as a single point source centred on the largest slip patch (Fig. 4).This illustrates the problem with the implicit assumption of Eq. ( 5) that there is an equal emission of waves at the relevant frequencies with source amplitude b, along the entire rupture.The amount of wave emission can vary along the rupture and may concentrate in a zone much smaller than the rupture length.In the model of Marc et al. (2016b) the number of sources contributing to wave emission and landsliding along a ruptured fault was given as L/l asp , with l asp the characteristic length scale of an asperity, set arbitrarily to 3 km.Thus, a 20 km long rupture is represented by six sources, distributed along most of the rupture length.The efficacy of their model did not depend much on the value of l asp because their seismological model was calibrated empir-ically, and l asp may well be larger.If we consider asperities as circular patches with diameter ∼ l asp , located strictly within the fault rupture and behaving as point sources with emission from their centre, then the relevant length along which waves are emitted is L − l asp and Eq.(5) becomes with H the Heaviside function, to represent that there should always be at least a point source in the middle of the fault even if it requires an asperity smaller than l asp .This has no significance for large earthquakes and long faults where L l asp , but may significantly reduce the predicted value of A s for smaller earthquakes.This modification of the model, although plausible, would not improve much the global residuals because Eq. ( 5) yields similar numbers of under-predicted and overpredicted small earthquakes, for which the model prediction is systematically reduced when using Eq. ( 6).Moreover, cases such as the Niigata or Iwate earthquakes, are still overpredicted when modelled with a single point source.This suggests that for these cases, with well-constrained source depth, a better prediction of R HMAX is needed, and therefore of either the source term b sat , or the threshold of acceleration for damage a c .

Threshold acceleration for landsliding
Equation (2) can be rewritten as a second-degree equation and solved for R HMAX and a c .Thus, assuming that Eqs. ( 4) and (3) hold and that b sat = 0.4 is constant, we can use A d , C topo , Mo and R 0 to invert for R HMAX and a c (Fig. 6).Inverted values of a c cluster around 0.15, which is the global value we have used here for the direct prediction, and about 50 % of the inverted values are between 0.1 and 0.2, consistent with values reported for case studies (Meunier et al., 2007;Hovius and Meunier, 2012;Yuan et al., 2013).The rest of the inverted values are mostly within the range 0.05-0.3and it is difficult to assess whether the inverted values of a c are representative of specifically weak or strong areas or whether other processes, not described by the model, affect the inversion.Defined as the threshold acceleration at which significant cohesion loss occurs, a c should be independent of hillslope steepness and depend only on material properties.Consistent with this view we find no correlation between a c and the landscape steepness as described by the modal slope of the landscape (see Marc et al., 2016b).To define and obtain quantitative estimates of substrate strength or of the ground pore pressure at the landscape scale is an outstanding challenge and lacking relevant constraints, we cannot assess further their influence on the variability of a c and A d .Nevertheless, it is interesting to focus on those earthquakes that have the lowest inverted values of a c < 0.05.A striking example is the 1988 Saguenay (Canada) earthquake, which caused landsliding over an exceptionally large area (A d = 45 000 km 2 ), despite its moderate magnitude (M w = 5.8) and large depth (R 0 = 28 km).Equation ( 5) predicts no landsliding for this event, and the inverted value of a c is as low as 0.01, more than an order of magnitude below the global threshold of 0.15.The Virginia 2011 (USA) earthquake is a very similar example.One explanation is that the furthest landslides defining A d were indeed occurring at very low shaking levels (see Jibson and Harp, 2012) but seismological processes may also be the cause.The Saguenay earthquake was very peculiar from a seismological point of view, with a larger than expected stress drop (Boore and Atkinson, 1992), and therefore probably larger strong motions (e.g.Baltay et al., 2011;Baltay and Hanks, 2014).Also, its occurrence close to the Moho discontinuity may have led to reflected waves reaching the surface with a stronger amplitude than the direct S-waves (Somerville et al., 1990), causing strong shaking up to ∼ 100 km from the epicentre.These effects have likely contributed to an exceptional pattern of strong ground motion and a significantly extended landslide distribution area.Evidently, such effects and mechanisms are not captured by our simple model.Therefore, inverting for a c and finding anomalous values, for example a c < 0.05, which is more than 3 times smaller than the global average, may suggest that an earthquake had an anomalous or complex seismic process.Hence, low a c values suggest that the Bihar, and possibly the Whittier Narrow, Alaska, and Wenchuan earthquakes had some mechanistic specificities, either in wave emission or wave propagation, that have affected the distribution of landsliding.For the Wenchuan case we have already highlighted the complex rupture and the importance of variable initial stress state and rupture velocity (Wen et al., 2012a, b).

Source term and earthquake stress drop
The near-source wave amplitude, b sat , is the only explicit parameter representing the earthquake source in our model.It was kept constant in our analyses.However, as was the case in the Saguenay earthquake, various seismological processes may increase or decrease the amplitude of waves emitted at the source of an earthquake.For example, it has been established that earthquakes with larger dynamic stress drops must have stronger surface ground motions, especially at high frequency (Hanks and McGuire, 1981;Baltay et al., 2011;Baltay and Hanks, 2014).For the 22 earthquakes with a constrained dynamic stress drop (Table 1) we do not find a clear correlation between the residuals of our model and the magnitude of the stress drop (Fig. S2).Except for three substantially under-predicted earthquakes (Saguenay, Arthur's Pass 1995, Erzincan) for which the large A d may be due to relatively large stress drops, the 19 other residuals do not seem to be controlled by the stress drops.Rupture speed and rupture speed variability may also influence the source emission and its characteristics (Wen et al., 2012a;Causse and Song, 2015).It has been argued that rupture speed could correlate negatively with stress drop, blurring the relation between earthquake stress drop and strong ground motion (Causse and Song, 2015).Accurate measurements of the rupture speed and its variations are now made for large earthquakes (e.g.Wen et al., 2012a), but they are lacking for most historic events or smaller events, impeding further exploration of their effects on landsliding.

Directivity and asymmetry in landslide distribution
We have shown that our model can predict the general distance from the emission line to a contour containing 95 % of all the landslides as measured by their cumulative surface area, with a reasonable reliability.However, this 1-D parameter does not describe potential 2-D complexities in the shape of the landslide distribution area.For example, across-strike asymmetry often exists for thrust faults where the hanging wall may experience larger shaking (Oglesby et al., 2000) and therefore has more intense landsliding over a larger distance (Gorum et al., 2011).This hanging wall effect is difficult to isolate because for dip-slip faults the topography in the hanging wall is usually closer to the earthquake source and also more prone to landsliding due to larger relief and steepness.Along-strike asymmetry of the landslide distribution may arise from asymmetry of the ground shaking pattern due to seismic directivity.Directivity is the result of interference between waves emitted at different times and locations along the rupture.It can enhance strongly the ground shaking in the direction of rupture propagation, but reduces shaking at the fault tip opposite the rupture direction (Wallace and Lay, 1995).These effects were clearly articulated in the 1976 Guatemala earthquake, with both high seismic intensities and dense landsliding limited to a narrow band along the fault near the rupture initiation, and spreading over a wider area at the other extremity of the fault (Harp et al., 1981) (Fig. 4).
To further explore these effects, we attempt to quantify asymmetry of the landslide distribution and in the rupture mechanism with very simple parameters available for most cases.We define a landslide asymmetry indicator L as = 100• (A op −A ep )/A T , with A T the total landslide area and A ep and A op the total landslide area beyond a line perpendicular to the fault strike and located at the fault tip closest to and furthest from the epicentre, respectively.This indicator effectively compares both fault extremities and quantifies the amount of the asymmetry as a percentage of the total landsliding.This ratio is corrected for the relative abundance of flat or gentle topography below the 10 • threshold in the relevant areas.It can be compared to the distance between the earthquake epicentre (projected on the fault trace) and the middle point of the fault rupture, normalized by the half length of the fault rupture, Epi as .This latter measure tends to zero for a centred epicentre and to one for an epicentre at the tip of the fault, suggesting low and high directivity, respectively.These indices cannot be computed for the Finisterre, Aysen, and Limon earthquakes for which we do not have constraints on the location of the actual portion of the fault that ruptured and its position relative to the epicentre.For the other cases, we find a degree of correlation between the asymmetry of the earthquake epicentre and the landslide distribution, at least where Epi as > 0.4, but results for the 10 cases with comprehensive landslide inventories are not straightforward (Fig. 7).
For Epi as < 0.4 we observed negligible landslide asymmetry.For larger Epi as values we observe a strong asymmetry oriented with the directivity for the Northridge and Guatemala earthquakes, and a moderate asymmetry for the Iwate case.For the Wenchuan case, we find little (10 %) asymmetry in the landslide distribution, even though the earthquake hypocentre was at one extremity of the fault rupture.In this earthquake, the rupture propagated over multiple fault segments with different initial stresses, with important rupture speed changes at the transition between segments (Wen et al., 2012a).This is likely to have complicated the ground shaking pattern (Wen et al., 2012b;Meunier et al., 2013a), blurring any directivity effect.Similar effects may be at play for the Gorkha and Denali earthquakes where published landslide maps (Gorum et al., 2014;Martha et al., 2016) indicate little (10 %) to no asymmetry (< 1 %), respectively (Fig. 7).These cases also had complex ruptures, with the Denali earthquake rupturing into three segments, the last one in super-shear (Frankel, 2004), and the Gorkha earth- against along-strike landsliding asymmetry, defined by the difference between the total landslide area beyond the tip of the fault opposed to the epicentre and the total landslide area beyond the other tip normalized by the total landslide area.quake propagating along a complex fault geometry with laterally varying mechanism (Kumar et al., 2016).Finally, we note that the amplitude of ground shaking is reduced by directivity but the shaking duration is longer, possibly balancing any effects on landslide triggering.In summary, the landslide asymmetry proxy proposed here varies between earthquakes and is not simply related to the epicentre position relative to the fault trace.Analysis of a larger number of wellconstrained cases is necessary to constrain and quantify the effect of seismic directivity on landslide patterns.

Conclusions
We have presented an analytical expression for the distribution area of earthquake-induced landslides.It shares its derivation with the model of Marc et al. (2016b) that has been shown to predict total landslide volume and area with good accuracy in most cases.The expression is based on scaling relationships between essential seismic parameters such as the seismic moment, the source depth, and the rupture length, and aims at predicting where the shaking level is likely to exceed a universal threshold for rock damage and onset of landsliding.Compared to a large compilation of estimates of landslide distribution areas from observational constraints, the analytical expression is shown to explain 56 % of the variance without any adjusted parameters, and to predict 35 or 49 of 83 cases within a factor of 2, when taking the hypocentral depth as the emission depth, or when allowing the emission depth to be within 25 % of the hypocentral depth, respectively.Analysis of outliers is not easy with such compilations as many cases are poorly constrained and because definition and measurement of the landslide distribution area is not uniform for all cases.For detailed inventories we see that the prediction of our model agrees relatively well with the region in which 95 % of the total landslide area is concentrated.However, some earthquakes are significantly overpredicted and others have important along-strike asymmetry not captured by our model.These discrepancies may arise from variability in wave emission along a complex rupture (e.g.Wen et al., 2012a), as well as wave interference leading to directivity.Thus, modelling of the ground shaking pattern must be improved to aid better understanding and forecasting of landslide hazard associated with earthquakes.Nevertheless, the overall agreement of our prediction is striking given that it does not use any geomorphological calibration.It suggests that 0.15 g is an appropriate threshold for ground acceleration to predict the area of concentrated landsliding, while the outer limits of landsliding may be controlled by a smaller acceleration at about 0.05 g (Jibson and Harp, 2016).However, at lower ground acceleration thresholds, site effect, non-linear attenuation, and other secondary controls become increasingly important, disproportionately complicating any prediction of landsliding.Still, understanding how the threshold of acceleration varies from a landscape to another is probably key to increasing the accuracy of our prediction as much as better understanding control on the ground shaking.The prediction of landslide distribution area, together with the recent expressions for total volume and area of landslides (Marc et al., 2016b) and constraints on the spatial pattern of landslide density and size (Meunier et al., 2007;Valagussa et al., 2017), defines a consistent framework for the evaluation of seismological parameters controlling ground shaking for quantitative landslide hazard.Finally, we note that the simplicity and limited number of parameters of the expression presented here makes it well suited for hazard assessment in earthquake scenarios as well as in the immediate aftermath of an earthquake.

Figure 1 .
Figure1.Total area affected by landsliding against seismic moment and moment magnitude for 83 earthquakes.Vertical error bars represent different estimate of A d for the 27 cases where they could be obtained.Name codes (defined in Table1) are shown for each earthquake, followed by S or N for strike-slip and normal faults, while all other cases are reverse fault earthquakes.R 0 is the mean depth of wave emission.The prediction of the seismologically consistent model is shown for reverse, normal, and strike-slip faults at different depth with solid black, solid grey, and dashed lines, respectively.Inset: distribution of the uncertainty factor (upper estimate/best estimate or best estimate/lower estimate).

Figure 2 .
Figure 2. Predicted landslide distribution area plotted against estimated landslide distribution areas for 83 earthquakes.For visibility, cases where the predicted area is 0 are set to 1 km 2 and only the name codes of earthquakes outside of a factor of 2 from the 1 : 1 line are shown.Vertical error bars represent the range of predicted values when R 0 is varied between 75 and 125 % of the best estimate of the hypocentral depth.Inset: histograms of model residuals (A d /A dp ).Histogram for the best empirical fit of A d against Mo (Emp, grey line) for the model prediction (Pre, dashed line) and for the prediction accounting for 25 % of uncertainty on R 0 (Pre R0 , solid black line) are shown, with the number of earthquakes correctly predicted within a factor of 2 and 4.

Figure 3 .
Figure 3. Distribution of landsliding (yellow polygons) and predicted landslide distribution area (green circles) for the Northridge (a), Haiti (b), Chi-Chi (c), Aysen (d), and Finisterre (e) earthquakes.The background is a topographic slope gradient map derived from 30 m Aster GDEM, provided by NASA, allowing us to locate flatlands where no landslides are expected even if the shaking threshold is exceeded.Emission line source and the area where the ground shaking is expected to be larger than a c are represented with green lines and circles, respectively.Red perimeters show the area encompassing 95 % of the total area of landsliding defined by a uniform distance away from the wave emission line.For reference rupture slip distribution maps are shown for the Northridge and Haiti earthquakes(Wald et al., 1996;Hayes et al., 2010).For Haiti, in addition to the landslide inventory from Gorum et al. (2013), we also show the one ofHarp et al. (2016) (cyan polygons) and its associated R 95 with red dashed line.Note that the Aysen strike-slip fault is located only based on the focal mechanism and epicentre, and that the other solution (north-oriented) would give a similar radius for R 95 .

Figure 4 .
Figure 4. Same as Fig. 3 for the Wenchuan (a), Niigata (b), Iwate (c), Limon (d), and Guatemala (e) earthquakes.For reference, rupture slip distribution maps are shown for the Niigata and Iwate earthquakes(Hikima and Koketsu, 2005;Suzuki et al., 2010).Additionally, the distribution area proposed byHarp et al. (1981) for the Guatemala earthquake is shown in purple.For the Limon earthquake, the total area of landslides is likely underestimated in the most affected area and therefore R 95 is likely reduced and more similar to the dashed red contour.

Figure 6 .
Figure6.Inverted value of the threshold of acceleration for damage, a c , against hypocentral depth.Most values cluster between 0.1 and 0.2, but some events (with their name tags displayed) have exceptionally low inverted values.

Figure 7 .
Figure7.Hypocentre asymmetry (0 = centre, 1 = tip of the fault) against along-strike landsliding asymmetry, defined by the difference between the total landslide area beyond the tip of the fault opposed to the epicentre and the total landslide area beyond the other tip normalized by the total landslide area.