Effects of finite source rupture on landslide triggering: The 2016MW 7.1 Kumamoto earthquake

The propagation of a seismic rupture on a fault introduces spatial variations in the seismic wavefield surrounding the fault. This directivity effect results in larger shaking amplitudes in the rupture propagation direction. Its seismic radiation pattern also causes amplitude variations between the strike-normal and strike-parallel components of horizontal ground motion. We investigated the landslide response to these effects during the 2016 Kumamoto earthquake (MW 7.1) in central Kyushu (Japan). Although the distribution of some 1,500 earthquake-triggered landslides as function of rupture distance is consistent 5 with the observed Arias intensity, the landslides were more concentrated to the northeast of the southwest-northeast striking rupture. We examined several landslide susceptibility factors: hillslope inclination, median amplification factor (MAF) of ground shaking, lithology, land cover, and topographic wetness. None of these factors sufficiently explains the landslide distribution or orientation (aspect), although the landslide headscarps have an elevated hillslope inclination and MAF. We propose a new physics-based ground motion model that accounts for the seismic rupture effects, and demonstrate that the low-frequency 10 seismic radiation pattern is consistent with the overall landslide distribution. Its spatial pattern is influenced by the rupture directivity effect, whereas landslide aspect is influenced by amplitude variations between the fault-normal and fault-parallel motion at frequencies < 2 Hz. This azimuth-dependence implies that comparable landslide concentrations can occur at different distances from the rupture. This quantitative link between the prevalent landslide aspect and the low-frequency seismic radiation pattern can improve coseismic landslide hazard assessment. 15


Introduction
Landslides are one of the most obvious and hazardous consequences of earthquakes. Acceleration of seismic waves alters the force balance in hillslopes and temporarily exceeds shear strength (Newmark, 1965;Dang et al., 2016). Greatly increased landslide rates have been reported on hillslopes close to earthquake rupture, mostly tied to ground acceleration (Gorum et al., 2011) and lithology (Chigira and Yagi, 2006). Substantial geomorphological and seismological data sets are required to assess 20 the response of landslides to ground motion and a growing number of studies has shed light on the underlying links (e.g. Lee, 2013;Allstadt et al., 2018;Roback et al., 2018;Fan et al., 2018). Several seismic measures such as vertical and horizontal peak   Kubo et al. (2016), the dashed box encompasses landslides related to the triggered event in Yufu (epicenter location after Uchide et al. (2016)). The inset map shows the station network within 150 km of the rupture. ground acceleration (PGA) (Miles and Keefer, 2009), root-mean square (RMS) acceleration or Arias intensity (I A ) (Arias, 1970;Keefer, 1984;Harp and Wilson, 1995;Jibson et al., 2000;Jibson, 2007;Torgoev and Havenith, 2016), seismic sourcemoment release, hypocentral depth, as well as rupture extent and propagation (Newmark, 1965) correlate with landslide density (Meunier et al., 2007). 5 Landslides concentrate in the area of strongest ground acceleration (Meunier et al., 2007), whereas total landslide area decreases from the earthquake rupture with the attenuation of peak ground acceleration (Dadson et al., 2004;Taylor et al., 1986).
This general pattern is modified by morphometrics (e.g. local hillslope inclination, curvature) and geological parameters (e.g. lithology, geological structure, land cover) (Gorum et al., 2011;Havenith et al., 2015) that influence landslide susceptibility (Pawluszek and Borkowski, 2017) on top of seismic amplification (Maufroy et al., 2015). For instance, Tang et al. (2018) found 10 that lithology, PGA, and distance from the rupture plane are important in assessing the distribution of landslides triggered by the 2008 Wenchuan earthquake (M W 7.9). Fan et al. (2018) found that hillslope aspect and slope were important determinants of the landslide distribution resulting from the 2017 Jiuzhaigou earthquake (M W 6.5).
On April 16, 2016 at 16:25 UTC central Kyushu was hit by a M W 7.1 earthquake (Fig. 1). The left-lateral dip-slip event ruptured along the Futagawa and Hinagu faults striking NW-SE with a hypocentral depth of 11 km (e.g. Kubo et al., 2016). 15 The rupture propagated northeastward and stopped at Mt. Aso. Fault source inversions show a northeast propagation of the rupture originating under Kumamoto city with highest slip near the surface at the western rim of the Aso caldera (e.g. Kubo et al., 2016;Asano and Iwata, 2016;Moore et al., 2017;Uchide et al., 2016;Yagi et al., 2016;Yoshida et al., 2017). The earthquake triggered approximately 1,500 landslides (National Research Institute for Earth Science and Disaster Prevention, 2016) that concentrated mainly inside the caldera and the flanks of Mt. Aso on the Pleistocene and Holocene lava flow deposits 20 (Paudel et al., 2008;Sidle and Chigira, 2004), although most of the terrain near the earthquake rupture is rugged (Fig. 1).
Thus, we hypothesize that rupture directivity causes an asymmetric distribution of landslides around the rupture plane, because of more severe ground motion along the propagating rupture (Somerville et al., 1997;Hovius and Meunier, 2012). Similarly asymmetric landslide distributions attributed to rupture directivity were repoprted for the 2002 Denali earthquake (M W 7.9) (Frankel, 2004;Gorum et al., 2014), and the 2015 Gorkha earthquake (M W 7.8) (Roback et al., 2018). In case of the 1999 25 Chi-Chi earthquake, Lee (2013) speculated that the prevalent landslide aspects were correlated to the fault movement direction (Ji et al., 2003;Meunier et al., 2008). These observations indicate that the rupture process introduces variations on the incoming energy on hillslopes.
Here we link those dominant near-surface seismic characteristics relevant to the pattern and orientation of coseismic landslides. We investigate the geological conditions (lithology, aspect, hillslope inclination, topographic amplification, soil wetness) 30 in central Kyushu as well as seismic waveform records from 240 seismic stations within 150 km of the rupture (Fig. 1). The two most prominent seismic effects-well founded in seismological theory (e.g. Aki and Richards, 2002) and documented in empirical relationships (e.g. Somerville et al., 1997)-are the rupture directivity and amplitude variations of fault-normal and fault-parallel motion. We examine whether the geomorphic characteristics around the Aso caldera made this area more susceptible to landslides than the surrounding topography near the earthquake rupture; or whether rupture effects control the asymmetric distribution of the landslides. We introduce a ground motion metric related to azimuth-dependent seismic energy (i.e. seismic velocity), because these effects attenuate with increasing frequency and are less captured by acceleration based metrics. We conclude by proposing a new ground-motion model that is consistent with the observed coeseismic landslide 5 pattern.

Data
We combine data sets on the response of landslides to the earthquake, including topography, land cover, geology, seismic waveforms, velocity structure, near surface characteristics, and landslide location and planform (Fig .2).
2.1 Topographic data 10 Most topographic data used in this study are provided by the Japan Aerospace Exploration Agency (JAXA) and its Advanced Land Observing Satellite (ALOS) project with a horizontal resolution of 1" (≈ 30 m). This digital surface model (DSM) forms the basis for computing aspect, hillslope inclination, the median amplification factor (MAF, Maufroy et al., 2015), and the topographic wetness index (Böhner and Selige, 2006). The ALOS project also provides data on land cover including anthropogenic influence (sealing, agriculture) and vegetation. While data on major geological units are from the Seamless

15
Digital Geological Map of Japan (scale 1:200,000) by the Geological Survey of Japan.

Topographic amplification of ground motion
Topographic features, such as mountains and valleys, can amplify or attenuate seismic waves (Massa et al., 2014;Maufroy et al., 2012Maufroy et al., , 2015. Largest ground-motion variations occur on hillslopes and summits, whereas variations are intermediate on narrow ridges, and low on valley floors. Maufroy et al. (2015) introduced proxies for these topographic site effects, of which 20 we use the median amplification factor (MAF), based on the topographic curvature, and the S-wave velocity v S traveling at frequency f : 2f is the topographic curvature convolved with a normalized smoothing kernel based on two 2D box-car functions as a function of v S and f .

25
The curvature is estimated from the DSM (Zevenbergen and Thorne, 1987;Maufroy et al., 2015) and the seismic velocity v S is the average S-wave velocity of the uppermost 500 m from the model by Koketsu et al. (2012).
Another site effect that influences landslide potential is the local soil or groundwater content, which can be modeled for uniform conditions to first order using the topographic wetness index (TWI) of Böhner and Selige (2006):

30
where A c is the upslope catchment area and β is the hillslope inclination derived from the DSM with filled sinks (Planchon and Darboux, 2001).

Ground motion data
Ground-motion data are from the Kik-Net/K-Net of the National Research Institute for Earth Science and Disaster Prevention (NIED) of Japan. NIED operates for Kik-Net both borehole and surface stations, and we use the latter only. The Japan Me- 5 teorological Agency (JMA) also released seismic data from the municipal seismic network for the largest earthquakes of the Kumamoto sequence. In total, data from 240 stations in Kyushu are available with complete azimuthal coverage within 150 km from the earthquake rupture (Fig. 1).
The analysis of seismic waveforms is based on accelerometric data only. Both the NIED and JMA data are unprocessed and we follow the strong motion processing guidelines by Boore and Bommer (2005). We use both acceleration and velocity in 10 further processing, and integrate the accelerograms to obtain velocity records. We correct the data with the automated baseline correction routine by Wang et al. (2011). The JMA accelerometric data further require a piece-wise baseline correction prior to the displacement baseline correction due to abrupt (possibly instrument related) jumps (Boore and Bommer, 2005;Yamada et al., 2007). We use the automated correction for baseline jumps by von Specht, 2018 (in prep.).
An earthquake was triggered approximatley 80 km to the northeast in Yufu 32 s after the Kumamoto earthquake (Uchide 15 et al., 2016) ( Fig. 1). Due to the close succession of the two events, waveforms of the triggered event interfere with the coda of the Kumamoto earthquake. We taper the data to reduce signal contributions by the triggered event. The taper position is based on theoretical traveltime differences between the P wave (v P = 5700 m s -1 ) arrival of the Kumamoto earthquake and the S wave arrival (v S = 3300 m s -1 ) of the triggered event. The respective travel paths to the stations are measured from the hypocenters. Since fewer instruments are located to the northeast and the triggered event close to the sea, less than 10% of the 20 data are strongly contaminated by the triggered event.
NIED hosts the rupture plane model of Kubo et al. (2016), which describes the slip history on a curved rupture plane (based on the surface traces of the Futagawa and Hinagu faults) with a total length of 53.5 km and width of 24.0 km (Fig. 1). We use the extent and shape of the rupture plane to estimate the landslide affected area and to define the rupture plane distance r rup , the shortest distance from the rupture plane. We follow the approach of Somerville et al. (1999) to identify the asperity from 25 the rupture plane model, which is is the area with more than 1.5 times the average slip.
The underground structure in terms of seismic velocities (v P , v S ) and density (ρ) (Koketsu et al., 2012) are available for 23 layers down to the mantle in ≈ 0.1 degree resolution covering all of Japan; we only consider the layers of the upper 0.5 km to compute a velocity average for and the MAF.
NIED provides data for the subsurface shear wave velocities (v S30 ) as well as site amplifications factors S amp . Contrary to 30 v S by Koketsu et al. (2012), v S30 is derived for the upper 30 m only and more suitable for energy estimates, which requires velocities at the surface (recording station). The site amplification factor S amp describes by how much seismic waves are amplified independent of their frequency.

Landslide data
Detailed landslide data are provided by NIED as polygons ( Fig. 1), mapped from aerial imagery with sub-meter resolution at different times after the Kumamoto earthquake. The first data set contains landslides that were identified between April 16 to 20, though the area close to the summit of Mt. Aso was not covered. A second data set was collected on April 29, 2016 and covers those parts of Mt. Aso that remained unmapped. However, the second data set may contain rainfall induced landslides, since the rainy season in Kyushu starts in May (Matsumoto, 1989), and there was rainfall after the Kumamoto earthquake and 5 landslides triggered by volcanic activity. We selectively combine the two data sets for this study, using only those landslides from the second database, that are also partly present in the first data set. We exclude any rainfall triggered landslides with this approach, though possibly omitting seismically induced landslides exclusive to the second database. However, the area in question is comparatively small to the full extent of the study area, and the missing landslides are minor in terms of their area.
Several landslides cluster ≈ 80 km to the northeast of the mainshock in the municipalities Yufu and Beppu ( Fig. 1), that were 10 hit by a triggered earthquake (Uchide et al., 2016). We hypothesize that the distant northeastern landslides were induced by this triggered event. This also explains the considerable gap of landslides (≈ 50 km) between Yufu and Aso ( Fig. 1 with the deposit area.

Total area affected by landsliding
We define the landslide affected area, in which coseismic landsliding occurred, as the area spanned by the rupture plane distance covering 97.5 % of the total landslide area (Harp and Wilson, 1995;Marc et al., 2017). Thus the total landslide affected area is 3968.6 km 2 and within 22.9 km distance from the rupture plane.

25
An M W 7.1 event with a fault length of 53.5 km and an asperity length of 12.78 km (3 km) results in a landslide affected area of 3914 km 2 (4406 km 2 ) using parameters proposed by Marc et al. (2016). We derived the event depth of 11.1 km as the moment weighted average of the rupture model of Kubo et al. (2016). Both estimates are consistent with our area estimate. Marc et al. (2016) introduced a topographic constant, A topo , relating the total landslide area to the area that excludes basins and inundated areas. We estimate A topo from the ALOS land cover finding that 97 % of all landslides occurred in areas without 30 anthropogenic influence, i.e. land with urban and agricultural use, and water bodies. We exclude water bodies, urban areas-predominantly the metropolitan area of Kumamoto City, and rice paddies from the topographic analysis, obtaining an affected area of 3037 km 2 , i.e. A topo = 0.68.

Total landslide area
Total landslide area is linked to several earthquake parameters, mostly magnitude and hypocenter or average rupture-plane 5 depth (Keefer, 1984;Marc et al., 2016). We adopted the relation by Marc et al. (2016) to check for completeness of the total landslide area of 6.38 km 2 . The actual total landslide failure plane is likely smaller, as the NIED data set provides the combined area of depletion and accumulation. The modal hillslope inclination is estimated at 15 • . Instead of the earthquake magnitude scaling relation (Leonard, 2010) used by Marc et al. (2016), we use the rupture extent reported by Kubo et al. (2016). The area model requires the average length of the seismic asperities, which Marc et al. (2016) globally assumed as 3 km. However, Landslide concentration is defined as landslide area per area at a given distance band (Meunier et al., 2007). For the seismic processing, we consider the rupture plane distance r rup based on the rupture model, instead of the hypocentral distance (Meunier et al., 2007) or the Joyner-Boore distance (Harp and Wilson, 1995). 20

Coseismic landslide displacement
The sliding-block model of Newmark (1965) is widely used to estimate coseismic hillslope performance (e.g. Kramer, 1996;Jibson, 1993Jibson, , 2007. The model estimates the permanent displacement on a hillslope affected by ground motion. Newmark (1965) established a relation for hillslope displacement in terms of the maximum velocity at the hillslope for a single rectangular where A is the magnitude of the acceleration pulse and a y m s −2 the yield acceleration, which is the minimum pseudostatic acceleration required to produce instability. For downslope motion along a sliding plane, a y is related to the angle of internal friction, φ f and the hillslope inclination, δ, by with the average factor of safety F S. Chen et al. (2017) characterized unstable hillslopes-related to both rainfall and earthquakes-5 by a safety factor of F S < 1.5.
An upper bound for the displacement d s , is based on two ground motion parameters (Newmark, 1965;Kramer, 1996): where PGA m s −2 and PGV m s −1 are the peak ground acceleration and velocity, respectively. Thus, the coseismic hillslope performance can be characterized by velocity and acceleration. In the following sections, we derive a ground-motion 10 model based on the acceleration-related Arias intensity and the velocity-related radiated seismic energy.

Ground motion metrics
Though PGA is the most widely used ground-motion metric in geotechnical engineering, the Arias intensity I A (Arias, 1970) is widely used to characterize strong ground motion for landslides: 15 where g = 9.80665 m s −2 is standard gravity and T 1 and T 2 are the times where strong ground motion starts and cedes. The acceleration a(t) has units of m s −2 and the Arias intensity has units m s −1 . The Arias intensity captures both the duration and amplitude of strong motion. Empirical relationships between I A and d s in terms of earthquake magnitude and epicenter distance have been developed (e.g. Jibson, 1993;Bray and Travasarou, 2007;Jibson, 2007).
Since PGA and I A are related to each other (e.g. Romeo, 2000) and the hillslope displacement depends on both velocity and 20 acceleration (Eq. (3), (5)), it is reasonable to characterize velocity similar to Arias intensity. The velocity counterpart to I A is IV 2, the integrated squared velocity (Kanamori et al., 1993;Festa et al., 2008): The squared velocity is also used in radiated seismic energy estimates. The quantity j E is the radiated energy flux of an earthquake and estimated by (Choy and Cormier, 1986;Kanamori et al., 1993;Newman and Okal, 1998) 25 where ρ and c are the density and seismic wave velocity at the recording site and S amp is the site specific amplification factor.
The distance from the rupture is given by r rup and k is a term for path attenuation (Anderson and Richards, 1975), and effects of transmission and reflection (Kanamori et al., 1993). The attenuation constant k is also influenced by anisotropy and structure heterogeneity (Campillo and Plantet, 1991;Bora et al., 2015). The full definition of the energy flux includes two terms for compressional waves (c = v P ) and shear waves (c = v S ). The radiated energy of an earthquake, E S , results from the integral over the wavefront surface where A is the area of the surface through which the wave passes at the recording station and represents the geometrical 5 spreading.
The radiated seismic energy E S describes the energy leaving the rupture area and is related to the seismic moment (Hanks and Kanamori, 1979) where ∆σ is the stress drop, µ the shear modulus, and M 0 the seismic moment. We make use of this relation, when considering 10 the magnitude related term in the ground motion model. Since most seismic energy is released as shear waves, we apply the shear wave velocity at the recording site (v S ) to the entire waveform, i.e. we assume that all waves arrive with velocity v S at a site. This assumption has the advantage that it does not require a separation of the record into P-and S-waveforms, simplifying the computation. In appendix B we show from a theoretical perspective that using a uniform v S has only a small impact on the overall energy estimate. The site-specific correction term for the energy estimateÊ, based on Eq. (8) and (9) becomes While E S is the radiated seismic energy at the source,Ê is estimated from the velocity records at a station and only approximates E S . Therefore,Ê may differ from the true and unknown radiated energy E S (Kanamori et al., 1993). Several assumptions characterizeÊ -All energy is radiated as S-waves in an isotropic, homogeneous medium 20 -Geometrical spreading is corrected for an isotropic, homogeneous medium -Since IV 2 (Eq. (7)) depends on the radiation pattern,Ê depends on azimuth -Attenuation is homogeneous -Surface waves are not considered -Site amplification is frequency-independent 25 Below, we investigate the azimuthal variation of the energy estimates to characterize the radiation pattern.
The estimated wavefront areaÂ is related to the rupture extent and r rup , andÂ corresponds to a simplified version of the wave front area approximation by Schnabel and Bolton Seed (1973); Shoja-taheri and Anderson (1988): The extent of the rupture is assumed to be rectangular with length L and width W . Equation (12) describes a cuboid with rounded corners with only half of its surface considered, because no energy flux is assumed to be transmitted above ground. 5 While the geometrical spreading correction is expressed analytically as the wavefront areaÂ, we estimate the attenuation parameter k. Attenuation changes with distance as a power law at short distances (< 150 km) (Anderson and Richards, 1975) and longer distances are not considered. An empirical attenuation relationship is: where Y is i.e. the logarithm of the energy estimate without the attenuation term e −krrup from Eq. (11). The dummy variable C is only used for estimating k and not in the final correction for attenuation. A distance independent form of the Arias intensity, i.e.
corrected for geometrical spreading and attenuation, is defined by 15 where k is determined by Eq. (13) and setting Y = I AÂ . The corrected Arias intensity I A,A is the acceleration based counterpart toÊ.
Low-frequency effects, like directivity, are better captured with a velocity based metric (e.g. azimuth-dependent energy estimate), than an acceleration based metric (Arias intensity) alone.
In terms of the Fourier transform, the sensitivity of acceleration at higher frequencies becomes apparent, as the Fourier 20 transform of the time derivative of any function is and thus scales with frequency in the spectrum. The frequency sensitivity of IV 2 and I A is related to the squared spectrum given the the metrics. As an example, we show in Fig. 3 the different spectral sensitivities of IV 2 and I A for a theoretical seismic source spectrum (Brune, 1970). IV 2, and thusÊ has a higher sensitivity to lower frequencies than I A . The low-frequency part 25 of the spectrum can be accounted for by considering IV 2 in a ground-motion model.

Landslide related ground-motion models
The basic form of landslide related ground motion models for Arias intensity is based on earthquake magnitude M and distance from the earthquake rupture r (e.g. Harp and Wilson, 1995).   Brune (1970). The spectrum can be read as displacement (red), velocity (black) and acceleration (blue).
b) The squared Brune spectrum corresponds to the frequency sensitivity of velocity based IV 2 (blue) and the acceleration based IA (black).
This form is widely used (Keefer, 1984;Harp and Wilson, 1995). In engineering seismology, ground motion models usually have an additional distance term for anelastic attenuation This is a modified version of the model template by Kramer (1996). While Eq. 17 and Eq. 18 share some parameters (p 1 , c 1 and p 2 , c 2 ), the geometric spreading term includes not only distance dependence (p 3 , c 4 ) but also has a magnitude dependent 5 component (c 5 ). In addition, anelastic attenuation is included as well (related to c 3 ) in Eq. 18. The template of Kramer (1996) relates to the majority of GMPEs in engineering seismology. Models of this kind address strong motion in the context of landsliding (Travasarou et al., 2003;Bray and Rodriguez-Marek, 2004). The incorporation of anelastic attenuation is less common in landsliding GMMs and not mentioned in these studies but included in more recent studies (Meunier et al., 2007(Meunier et al., , 2013; Yuan et al., 2013). 10 We exchange the magnitude term from Eq. (18) with a site-dependent energy term, assuming that landsliding is more related to the energy of incoming seismic waves than to the moment at the source. We replace moment magnitude by the logarithm of energy (Eq. (11)), since energy is proportional to the seismic moment M 0 (Eq. (10)). Based on the site-dependent energy estimateÊ, we propose the model The five coefficients are inferred by non-linear least squares (e.g. Tarantola, 2005). We use the rupture plane distance (r rup ), i.e. the shortest distance between a site and the rupture plane.

Rupture directivity model
In the NGA-west2 guidelines (Spudich et al., 2013), the directivity effect is modeled by isochrone theory (Bernard and Madariaga, 1984;Spudich and Chiou, 2008) or the azimuth between epicenter and site (Somerville et al., 1997). We use 20 the latter approach and model directivity for estimated energy and corrected Arias intensity in a simplified way: ln whereÊ 0 and I A,A,0 are the offset (average), a E and a I the amplitude of variation with azimuth and θ E and θ I are the azimuths of the maximum. The definition of θ is similar to that of Somerville et al. (1997) as the angle measured between the epicenter and the recording site with the difference of being measured clockwise from north. The azimuths of the maximum, θ E and θ I , are free parameters because (1) the rupture is assumed to have occurred on two faults and has thus variable strike, and (2) the event is not pure strike-slip and has a normal faulting component. We therefore do not expect a match between the rupture strike and θ E and θ I . The geometrical spreading is already incorporated in the energy estimate as a distance term (Somerville et al., 1997;Spudich et al., 2013).

Model for fault-normal to fault-parallel ratio
The ratio of the response spectra of the horizontal sensor components is a function of oscillatory frequency f osc .
The north and east components (E, N ) of the sensor are rotated to be fault-normal (F N ) and fault-parallel (F P ).
The response spectra are calculated from accelerograms after Weber (2002) with a damping of ζ = 0.05.
The amplitudes of waves parallel to rupture propagation differ from waves normal to propagation on top of the directivity effect. This variation depends on the azimuth and is larger only at high periods. The fault-normal response amplitude is larger 10 than the fault-parallel response if directed (anti)parallel to the rupture. We model the ratio similar to Somerville et al. (1997) ln where parameters b i describe the relationship of the oscillatory frequency to the ratio, θ is the azimuth (Eq. (20)), and θ R is the azimuth of the maximal ratio. The ratio azimuth is subject to assumptions as is its counterpart θ E . The Heaviside function H(·) avoids negative values in the model, which would be equivalent to an undesired phase shift in the cosine term. 15 We introduced a functional form for oscillatory frequency dependence with four parameters in Eq. (25). We did not introduce a distance term and apply the model only to data with r rup ≤ 50 km.

Topographic analysis
Landslides occurred mostly in tephra layers (Fig. 2a,b) covered by forests (Fig. 2d,e) and predominantly along the NE rupture 20 segment. Nearly all landslides concentrated on hillslopes between 15 • and 45 • steep and a MAF ≈ 1 (Fig. 4a,c). Hillslope inclination and MAF were higher towards the landslide crown (Fig. 4b,d), indicating a progressive landslide failure starting from the crown, consistent with numerical simulations by Dang et al. (2016). TWI is linked to land cover and is highest in areas with rice paddies (Fig. 2i). Terrain with landslides has uniformly low TWI, thus we cannot evaluate the hydrological impact on the earthquake related landslides (e.g. Tang et al., 2018).

25
Most landslides originated at locations with amplified ground accelerations and steep hillslopes and ran out on flatter areas with less amplified ground acceleration. Landslides-interpreted as shear failure-start as mode II (in-plane shear) failure at the scarp and mode III (anti-plane shear) failure at the flanks (McClung, 1981;Fleming and Johnson, 1989;Martel, 2004). At later stages of the landslide rupture, mode I (widening) failure can also occur in the process (Martel, 2004). Simulations of elliptic landslides by Martel (2004) show that either the most compressive or the most tensile stresses are parallel to the major   Kernel density estimate Exponential distribution (λ = 1.72 × 10 -4 km -1 ) Figure 5. a) Landslide concentration with a) rupture distance rrup and b) asperity distance raspof the Kumamoto earthquake landslides.
The rate parameter of the the exponentially decaying landslide concentration is estimated by maximum likelihood. The distances to the four peaks shown in Fig. 1 are given. Densities change little with distance metric, as highlighted by the similiar kernel density estimates and the near-identical rate parameter estimatesλ. The landslide concentration for Mt. Aso depends more on the distance metric than for the other three locations. The more distant mountains have very similar concentrations despite differences in distances (in particular Mt. Otake).
However, when compared to Fig. 7 Mt. Shutendoji has a higher landslide concentration than Mt. Kinpo and Mt. Otake, despite being the farthest away.   Mt. Aso and its caldera and Mt. Shutendoji had a high density of landslides (Fig. 5), whereas Mt. Kinpo and Mt. Otake had none, despite being closer to the epicenter and at comparably close to the rupture (Fig. 5). All these locations have comparable rock type, land cover, comparable hillslope inclination, and MAF. Hence, lithology, land cover, and topographic characteristics are insufficient to explain the landslide distribution and concentration with respect to the hypocenter or the asperity.
The azimuthal density-with respect to the epicenter and asperity centroid-of the unspecified landslides follows to some 5 extent the distribution of hillslope inclinations > 19 • in the landslide affected area (Fig. 6). This similarity shows that the abundance of unspecified landslides mimics the steepness of topography in the region. Densities are higher towards Mt. Kinpo

Impact of finite source on ground motion and landslides
The results of the seismic analysis are given for waveforms, the basis forÊ and I A , and response spectra, used for F N/F P .
To the northeast, signals with forward-directivity are shorter in duration with one or few strong pulses (Fig. 8,  in the fault-parallel direction (Fig. 8, top left).
We estimated energiesÊ from the three-component waveforms. For the Arias intensity, both horizontal components are used.
The geometrical spreading A is calculated according to Eq. (12) with a rupture length of L =53.5 km and width of W =24.0 5 km. Any remaining distance dependence has been corrected for by estimating and applying the attenuation parameter k (Eq. This study Uncertainty (±σ) Somerville et al. (1997) Uncertainty ( The azimuthal variation ofÊ and I A,A is modelled according to Eq. (20). We estimate parameters for two scenarios: directivity is assumed, resulting in azimuthal variations, where a E and a I are free parameters, directivity is not assumed, resulting in no azimuthal variations with a E = a I = 0.
The two models are compared with the Bayesian Information Criterion (BIC, Schwarz, 1978) for a least squares fit: Somerville et al. (1997) Uncertainty (±σ) This study Uncertainty (±σ) Azimuth Somerville et al. (1997) Uncertainty ( Figure 11. Kernel density estimate of F N/F P with azimuth obtained from response spectra for three different oscillatory frequency ranges: a) 0.1 -1 Hz, b) 1.0 -2.5 Hz, c) >2.5 Hz. For each plot, our F N/F P model and the model by Somerville et al. (1997)  where n is the number of estimated parameters (n =4 for first case, n =2 for second case), N is the number of data, andσ 2 is the variance of the model residuals. The model with the smaller BIC is preferable. The starting values of the parameters are the mean ofÊ and I A,A , no azimuthal variation (a d =0), and the azimuths of the maximum ofÊ θ and I A,A,θ are set to the strike of the fault (θ E = θ I = 225 • ). 10 The directivity model forÊ follows the trend of the data and the running average closer than the model without directivity (Fig. 9a). According to BIC, the model with directivity is preferable (BIC directivity = -110, BIC no directivity = -11). In case of the Arias intensity, the difference in BIC between the two models is less compared to the azimuth-dependent energy (Fig. 9b).
Here, the model without directivity is the preferred one (BIC directivity = 30, BIC no directivity = 22). In consequence, azimuthal variations in wave amplitudes and energy related to the directivity effect occur at lower frequencies.
The forward directivity waves contain a very strong low frequency pulse (Fig. 8). The pulse amplitude depends on the ratio of rupture and shear wave velocity and the length of the rupture (Spudich and Chiou, 2008). The forward directivity pulse is superimposed by high frequency signals in acceleration traces but becomes more prominent in velocity traces (Baker, 2007), due to its low frequency nature, i.e. below 1.6 Hz (Somerville et al., 1997).  (Fig. 11, Fig. 11). F N/F P is most variable with azimuth at low oscillatory frequencies (0.1 -1 Hz, Fig. 11a); variations are much smaller between 1 and 2.5 Hz (Fig. 11b); and nearly absent above 2.5 Hz (Fig. 11c). This decrease with frequency is captured by the F N/F P model (Eq. (25), Fig. 10). Since our model is an average over the covered distance, with an average rupture distance of 25.06 km, we compare it to the F N/F P model of Somerville et al. (1997) at 25 km (Fig. 10, Fig. 11).
Both models show a similar decay with frequency with our model predicting a slightly higher F N/F P . Therefore, the wave polarity ratio related to rupture directivity is pronounced at lower frequencies and dissipates with increasing frequency, similar 5 to the azimuthal variations observable in energy estimates (lower frequencies), but not in Arias intensity (higher frequencies).
The pattern of low-frequency ground motion is well reflected in that of the landslides. The azimuthal variation ofÊ coincides with that of landslide concentration (Fig. 9). Both azimuth-dependent energy and landslide concentration have a similar trend with the maximum parallel to rupture direction and the minimum strike anti-parallel. The orientation of maximum F N/F P is also reflected in the landslide aspect. The northwest and east directions show higher landslide density (Fig. 12a). The highest density of landslides has a northwestern aspect in agreement with maximum F N/F P , both perpendicular to the strike. The eastward increased density is mostly due to landslides very close to the rupture. A look at different distances reveals that the increased density of landslides facing east by southeast is at very short distances (r rup ≤ 2.5 km, Fig. 14), while the northwest facing landslides are further away (2.5 km < r rup ≤ 6 km). Only minor landslides are farther away with no specific pattern.
The distribution of aspect and hillslope inclination in the landslide affected area varies little with aspect (Fig. 12b).. The distinct northwest and east orientation of landslides is not an artefact of the orientation of the topography in the landslide affected area (12a,b). The unspecified landslides in the affected area have a near northward aspect and deviate by ≈ 30 • from the earthquake triggered landslides (Fig. 12c). This highlights that the earthquake affects landslide locations (Fig. 6, 7); and 5 will force failure on specific slopes facing in the direction of ground motion (Fig. 12, 14).

Ground motion model for Kumamoto
We derived two ground-motion models for Arias intensity from data with r rup ≤ 150 km (Tab. 1, Fig. 15). One model incorporates the azimuth-dependent seismic energy (Eq. (19)). The other is a conventional isotropic moment magnitude dependent model (Eq. (18)). The decay of Arias intensity with distance for both models fits well the running average and is proportional 10 to the decrease in landslide density with distance. Variation of estimated energy is well covered by the model and spans more than two orders of magnitude resulting in a variation of Arias intensity of nearly one order of magnitude.

•
A sp ec t Figure 14. Distribution of landslides with aspect and rupture distance. The rupture distance is measured from the model by Kubo et al. (2016). This model does not completely reach the surface, truncanting distances below 1 km. The distribution has been normalized by the distribution of aspect of the affected area.

Discussion
We provide a framework for characterizing coseismic landslides with an integrated approach of geomorphology and seismology emphasizing here the role of low-frequency seismic directivity and finite source. Given the observations of ground motion of the Kumamoto earthquake, two questions arise: (1) How specific is the observed ground motion, i.e. is the Kumamoto rupture 10 particularly distinct? (2) As a rupture very close to the surface, how much does seismic near-field motion contribute? The second question arises, because many landslides occurred very close to the rupture plane. However, it is not possible to separate the observed waveforms into near-, intermediate-, and far-field terms. To investigate both questions, we computed theoretical waveforms after Haskell (1964); Savage (1966); Aki and Richards (2002) from a circular rupture on an elliptic finite source with constant rupture velocity in a homogeneous, isotropic, and unbound medium (see Appendix).
Despite the simplified assumptions behind this waveform model, low-frequency ground motions capture the most prominent features of the observed waveforms. Simulated waveforms close to the rupture plane change in polarity orientation towards east-west, while a strong fault-normal polarity appears at larger distances. A decomposition into a near-field term and combined intermediate-and far-field term reveals that the former highly contributes to the ground-motion at close distances. The impact 5 of the near-field term may explain the dominance of east-facing landslides close to the rupture (Fig. 14).
The simulations also demonstrate the effect of directivity on estimates of radiated energy and Arias intensity. The azimuthal variations of simulatedÊ are similar to the observed variations. The Arias intensity of the simulations also has azimuthal variations with the same characteristics as the energy estimate. These variations in Arias intensity are absent in the observed data, indicating that Arias intensity is more influenced by local heterogeneities and scattering than the energy estimates as these 10 are ignored in the simulations.
The results show that the Arias intensity is not as susceptible to the directivity effect and variations in fault-normal to fault-parallel amplitudes as is the radiated energy: because of its higher sensitivity towards higher frequencies, these effects are masked by high-frequency effects such as wave scattering and a heterogeneous medium. We found that the radiation pattern related to the directivity effect is recoverable from energy estimates but not from Arias intensity. This low-frequency 15 dependence is also seen in the response spectra ratios for F N/F P where directivity related amplitude variations with azimuth have been identified only for frequencies <2 Hz, and in agreement with previous work (Spudich et al., 2004;Somerville et al., 1997). We introduced a modified model for Arias intensity using site-dependent seismic energy estimates instead of the source-dependent seismic magnitude to better capture the effects of low-frequency ground motion.
The conventional magnitude-based, isotropic model and the azimuth-dependent seismic energy model correlate with the 20 landslide concentration over distance (Fig. 15). As in Meunier et al. (2007) it is therefore feasible to use the ground-motion model to model the landslide concentration, P ls (I A ), by a linear relationship ln P ls (I A ) = a I + b I ln I A Azimuthal variations of landslide density correspond to azimuthal variations in seismic energy and can be described by a similar relationship For the Kumamoto earthquake data we estimate a I = 2.1, b I = 2.6 and a E = −31.5, b E = 2.3. The azimuth-dependent landslide concentration implies similar landslide concentrations at different distances from the rupture, thus partly explaining some of the deviation in Fig. 5 and Fig. 15.
Compared to the model of Harp and Wilson (1995) (Fig. 15) our model uses rupture-plane distance, as opposed to the Joyner- 30 Boore distance (r JB ). When using the hypocentral depth as pseudodepth, the model of Harp and Wilson (1995) overpredicts I A both at shorter and longer distances-irrespective of the pseudodepth at larger distances. This misestimate is most likely due to the lack of an additional distant dependent attenuation term in their model (Eq. (17)).
The use of MAF instead of curvature alone provides a proxy by how much a seismic wave is amplified (or attenuated) for a given wavelength and location. We show that both hillslope inclination and MAF tend to be lower towards the landslide toe ( Fig. 4). This effect is linked to the convention that landslide polygons cover both the zone of depletion and accumulation. (Sato et al., 2017) consider the tephra layers rich in halloysite to be the main sliding surfaces indicating shallow landslides (Song et al., 5 2017). When relating coseismic landsliding to the seismic rupture, only the failure plane of the landslide matters, because this is the hillslope portion that failed under seismic acceleration. Chen et al. (2017) noted, for example, that landslide susceptibility and safety factor calculation depends on whether the entire landslide or only parts-scarp area or area of dislocated massare considered. The reconstruction of the landslide failure planes is limited to statistical assessments of landslide inventories (Domej et al., 2017;Marc et al., 2019). However, failure may have likely originated close to the crown and then progressively 10 propagated downward the hillslope, because MAF > 1 indicates an amplification of ground motion towards the crown of the landslides.
Coseismic landslide locations have a uniformly low topographic wetness index, indicating that hydrology may have added little variability to the pattern of the earthquake triggered landslides; at least we could not trace any clear impact of soil moisture on the coseismic landslide pattern (Tang et al., 2018). 15

Conclusions
We investigated seismic waveforms and resulting landslide distribution of the 2016 Kumamoto earthquake, Japan. We demonstrate that ground motion at higher frequencies controls the isotropic (azimuth-independent) distance dependence of Arias intensity with landslide concentration. In addition, ground motion at lower frequencies influences landslide location and hillslope failure orientation, due to directivity and increased amplitudes normal to the fault, respectively. Topographic controls (hillslope 20 inclination and MAF) are limited predictors of coseismic landslide occurrence, because areas with similar topographic and geological properties at similar distances from the rupture had widely differing landslide activity Massey et al., 2018). Nonetheless landslides concentrated only to the northeast of the earthquake rupture, while unspecified landslides have been identified throughout the affected region.
We introduced a modified model for Arias intensity using site-dependent radiated seismic energy estimates instead of the source-dependent seismic magnitude to better model low-frequency ground-motion in addition to the ground-motion at higher frequencies covered by the Arias intensity.
Compared to previous models widely used in landslide related ground motion characterization our model is based on stateof-the-art ground-motion models used in engineering seismology, which have two different distance terms, one for geometrical 5 spreading and one for along-path attenuation. The latter is rare in landslide studies (e.g. Meunier et al., 2007;Massey et al., 2018). Our results emphasize that the attenuation term should be considered in ground-motion models, as the landslide concentration with distance mirrors such ground-motion models.
The effect of the earthquake rupture on the rupture process of the landslides results in landslide movements parallel to strongest ground-motion. Due to the surface proximity of the earthquake rupture plane, near-field ground motion influences Figure A1. Setup of the rupture model. Gray ellipse represents the rupture: light gray area is unruptured, medium gray area is slipping, and the dark gray area is after slip arrest.
the aspect of close landslides to be east-southeast. The intermediate-and far-field motion of the earthquake promoted more landslides on northwest exposed hillslopes, an effect that overrides those of steepness and orientation of hillslopes in the region.
We highlight that coseismic landslide hazard estimation requires an integrated approach of both detailed ground-motion and topographic characterization. While the latter is well established for landslide hazard, ground-motion characterization has been only incorporated by simple means, i.e. without any azimuth-dependent finite rupture effects. Our results for the Kumamoto earthquake demonstrate that seismic waveforms can be reproduced by established methods from seismology. We suggest that 5 these methods can improve landslide hazard assessment by including models for finite rupture effects.
Appendix A: Synthetic waveforms from displacement of a finite rupture We illustrate the generation of ground displacement as a discontinuity across a rupture fault (e.g. Haskell, 1964Haskell, , 1969Anderson and Richards, 1975;Aki and Richards, 2002). The displacement for any point x at time t is given by where c is the fourth order elasticity tensor from Hooke's law, G is the Green's function describing the response of the medium, D(ξ, t) is the displacement on the fault with area Σ and coordinates ξ, n is the fault normal vector. Summation over i,j,p,q is implied. While the surface integral is carried out numerically, the derivatives of the Green's function for an isotropic, homogeneous, and unbound medium can be solved analytically.
and δ ij is Kronecker's delta. The terms in Eq. (A2) are commonly separated in groups with respect to their distance r. In Eq.
(A2a) is the near-field (NF) term, as its amplitude decays with r −4 , it affects the immediate vicinty of a rupture only. Terms with a distance attenuation proportional to r −2 are called intermediate-field (IF) terms for P -waves (Eq. (A2c)) and S-waves 25 (Eq. (A2d)). The remaining two terms are the far-field (FF) terms for P -waves (Eq. (A2e)) and S-waves (Eq. (A2f)) with a decay proportional to r −1 . A major difference between the NF and IF terms, and the FF terms is that the former depend on the slip on the rupture and they are the cause for static and dynamic displacement; whereas the latter are functions of the time derivative of slip and result in dynamic displacement only.
The slip function in time is related to the Yoffe function Yoffe (1951); Tinti et al. (2005) with rise time T . We use the slip distribution of Savage (1966) to describe the amplitude distribution of the slip on the rupture, as well as the elliptical fault 5 shape and rupture propagation from Savage (1966). The slip amplitude is given by where D 0 is the maximum displacement at the center of the fault, L and W are the length and width of the fault, and p determines whether the rupture starts at the focus at the front of the rupture plane (strike-parallel, p = 1) or at the focus at the end (strike-anti-parallel, p = −1). The rupture originates in one of the two foci and propagates radially away from the source 10 with constant velocity ζ and terminates when it reaches the rupture boundary. The slip vectorŝ describes the orientation of the displacement D(ξ) on the fault plane. We follow the definition ofn andŝ in terms of fault strike φ s , dip δ, and rake λ from Aki and Richards (2002): The displacement vector D in Eq. (A2) is given by We consider an isotropic medium and the elasticity tensor c from Eq. (A1) is where λ and µ are the Lamé parameters We set λ = µ, resulting in the widely observed relation v P = v S √ 3.
With the assumptions outlined above it is possible to calculate the displacement of an earthquake at location x with 12 parameters (Fig. A1):

25
fault size and orientation: length L, width W , strike φ, dip δ material: 1st and 2nd Lamé parameters λ and µ, density ρ (alternatively: compressional and shear wave velocities v P and v S and density ρ) rupture and slip: rupture velocity ζ, slip D 0 , rise time T , rake λ, rupture orientation with respect to strike, p The fault size and displacement of earthquakes are correlated with each other and are scaled to the magnitude. The number of parameters reduces to ten (nine if the Lamé constants are equal), when scaling relations (e.g. Leonard, 2010;Strasser et al., 2010) are used in combination with the seismic moment M 0 . The moment can be decomposed in with shear modulus (2nd Lamé constant) µ, the rupture area-here an ellipse-A = π 4 LW , and average displacement,D which follows from Eq. (A4) asD = 2 3 D 0 . α = 2 km s −1 α = 3 km s −1 α = 4 km s −1 α = 5 km s −1 α = 6 km s −1 α = 7 km s −1 Figure B1. Ratio between the approximate and exact energy estimates for different P-wave velocities in the medium. The exact estimate assumes that P-and S-waves arrive at different velocities at the recording site, while the approximate estimate assumes that all waves arrive with shear wave velocity at the site. This approximation introduces only a minor underestimation, since most radiated energy is released as S-waves. The distance variation arises from the different distance and velocity dependencies of the intermediate-field terms and the far-field terms.
The results are not strictly comparable to observed data due to the models simplicity. The computed amplitudes will be smaller than observed values, because no free surface is assumed. Assuming a free surface would nearly double the amplitudes from wave reflection, as well as the amplifications from wave transmissions (from high to low velocity zones). Only 10 direct waves are computed, and effects of reflections of different layers are not covered due to the isotropy and homogeneity.
Corresponding waveforms-in particular surface waves-are not exhibited. However, the purpose of this model is to show (1) the general behavior of waveforms in the vicinity of a rupture, which is dominated by direct waves, and (2) how amplitudes distribute relatively in space.

Appendix B: Radiated seismic energy estimation
The exact calculation of radiated seismic energy is challenging. One simplifying assumption is that all waves arrive at the site with shear wave speed, which is a very good approximation for the far-field term. The reasoning can be justified from a theoretical perspective: For most earth media the ratio between the P-wave velocity α and S-wave velocity β is From this and Eq. A2e and A2f follows that the amplitude of compressional waves is ≈ 1 √ 3 the S-waves becomes ( 1 √ 3 3 ) 2 = 1 27 . The total energy of a signal is (Rudnicki and Freund, 1981) E total = E P + E S (B2) and can be estimated by 10Ê total = α S aIV 2 α + β S aIV 2 β , with the integrated squared velocity (IV 2) for P-and S-waves from Eq. 7, the P-and S-wave velocities α P and α S at the recording site, and a constant a covering the remaining factors which are identical for both terms (compare with Eq. 8). If we express the energy contribution of P-waves in terms of S-waves, we can summarize the above relation tô The last expression is differs only by 2.6 % from the exact term. While slightly underestimating the energy, this aproximate definition of using α S instead of β S does not require the identification of P-and S-waves. This is useful, since at short distances the S-wave train is usually inseparable from the P-wave train. 20 At shorter distances, the intermediate-field term needs also to be taken into consideration. The amplitude of the intermediate term decays with r 2 (Eq. A2c, A2d), while the far-field amplitude decays with r (Eq. A2e, A2f). That is, the amplitude scales by distance and velocities and thus the IV 2 are IV 2 α = α −4 r −2 (r −1 + α −1 ) 2 , IV 2 β = β −4 r −2 (r −1 + β −1 ) 2 .