Effects of finite source rupture on landslide triggering: The 2016 <i>M<sub>W</sub></i> 7.1 Kumamoto earthquake

Abstract. The propagation of a seismic rupture on a fault introduces spatial variations in the seismic wavefield surrounding the fault during an earthquake. 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 Kyūshū (Japan). Although the distribution of some 1,500 earthquake-triggered landslides as function of rupture distance is consistent with the observed Arias intensity, the landslides are 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 can sufficiently explain the landslide distribution or orientation (aspect), although the landslide headscarps coincide with 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 seismic radiation pattern consistent with the overall landslide distribution. The spatial landslide distribution is primarily influenced by the rupture directivity effect, whereas landslide aspect is influenced by amplitude variations between the fault-normal and fault-parallel motion at frequencies 



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 exceed cohesion and friction (Newmark, 1965;Dang et al., 2016). 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). Several seismic measures such as vertical and horizontal peak 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 source-moment release, hypocentral depth, as well as rupture extent and propagation (Newmark, 1965) correlate with landslide density (Meunier et al., 2007).
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., 5 1986). Those metrics are combined with 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), since they could alter landslide susceptibility (Pawluszek and Borkowski, 2017), in addition to seismic amplification (Maufroy et al., 2015). For instance, Tang et al. (2018) found that lithology, PGA, and distance from the rupture plane are important in assessing coseismic landslide distribution that were triggered by the 2008 Wenchuan earthquake (M W 7.9). 10 On April 16, 2016 at 16:25 UTC central Kyūshū 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).
The rupture propagated northeastward and stopped at Mt. Aso. Fault source inversion studies 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., 15 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 (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). Similarly asymmetric landslide 20 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 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 land- 25 slides. To this end we investigate the geological conditions (lithology, aspect, hillslope inclination, topographic amplification, soil wetness) in central Kyūshū 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 30 area more susceptible to landslides than the surrounding topography near the earthquake rupture; or instead rupture effects control the asymmetric distribution of the landslides. We introduce a ground motion metric related to azimuth-dependent seismic energy (i.e. related to seismic velocity), because these effects attenuate with increasing frequency and are less captured by acceleration based metrics. From this we propose a new ground-motion model that is consistent with the observed coeseismic landslide pattern. 35 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.

Topographic data
Most topographic data used in this study are provided by the Japan Aerospace Exploration Agency (JAXA) and its Advanced

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 we use the median amplification factor (MAF), based on the topographic curvature, and the S-wave velocity v S traveling at 15 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 . The curvature is estimated from the DSM (Zevenbergen and Thorne, 1987;Maufroy et al., 2015) and the seismic velocity 20 v S is the average S-wave velocity of the uppermost 500 m from the model of Koketsu et al. (2012). The frequency f of the seismic wave is the fundamental frequency of the hillslope section on which landsliding occurred (Massa et al., 2014).
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): 25 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 Meteorological 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 Kyūshū are available with complete azimuthal coverage within 150 5 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 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 10 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 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 15 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 lose to the sea, less than 10 % of the 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 20 on the surface traces of the Futagawa and Hinagu faults) with a total length of 53.5 km and 24.0 km width ( 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.
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 25 compute a velocity average for the computation of MAF.
NIED provides data for the subsurface shear wave velocities (v S30 ) as well as site amplifications factors S amp . Contrary to v S from Koketsu et al. (2012), v S30 is derived for the upper 30 m only and more suitable for energy estimates, which require 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 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 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 considered to be 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), As mentioned above, this area was hit by a triggered earthquake (Uchide et al., 2016). We believe that the distant northeastern landslides of any origin. We extract a subset from this landslide database to compare it with the landslides triggered by the Kumamoto earthquake. Contrary to the special Kumamoto release, only the landslide deposits are given as polygons, whereas the scarps are mapped as lines. We manually define polygons representing the total landslide area bound by the scarp line and covering the deposit area to make both data sets comparable. This step is essential, because the landslide source area is generally not identical 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.
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 25 area of 3914 km 2 (4406 km 2 ) using parameters proposed 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. 97 % of all landslides occurred in areas without anthropogenic influence, i.e. land with urban and agricultural use, and water bodies. We exclude water bodies, urban areas-predominantly 30 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 depth. For instance, Keefer (1984) developed empirical relations, and Marc et al. (2017) introduced a model with physical assumptions. We adopted the relation by Marc et al. (2017) to check for completeness of the total landslide area of 6.38 km 2 . The actual landslide failure plane is likely smaller, as the NIED data set provides the combined area of depletion and accumulation. Landslide concentration is defined as landslide area per area at a given distance band (Meunier et al., 2007). For the seis-15 mic 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).

Coseismic landslide displacement
The sliding-block model of Newmark (1965) is widely used to estimate coseismic hillslope performance (e.g. Kramer, 1996;20 Jibson, 199320 Jibson, , 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 25 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 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.

Ground motion metrics
Though PGA is the most widely used ground-motion metric in geotechnical engineering, the Arias intensity I A (Arias, 1970) 5 is widely used to characterize strong ground motion for landslides.

The Arias intensity is defined by
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) is given in units of m s −2 and the Arias intensity in m s −1 . The Arias intensity captures both the duration 10 and amplitude of strong motion. Empirical relationships between I A and d s 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 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) where ρ and c are the density and seismic wave velocity at the recording site and S amp is the site specific amplification factor. 20 The distance from the rupture is given by r rup and k is a term for along-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 ), respectively. 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 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 and µ the shear modulus. We make use of this relation, when considering the magnitude related term in the ground motion model. Since most seismic energy is released as shear waves, we apply the shear wave velocity to the entire waveform and the site-specific correction term for the energy estimate, and the energy estimate,Ê, based on Eq. (8) and (9) becomeŝ 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 -Geometrical spreading is corrected for an isotropic, homogeneous medium 15 -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 Below, we investigate the azimuthal variation of the energy estimates to characterize the radiation pattern. 20 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 of 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. 25 While the geometrical spreading correction is expressed analytically, in form of 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: 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 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. 10 In terms of the Fourier transform, the sensitivity of acceleration at higher frequencies becomes apparent, as the Fourier 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 squared nature of the metrics. As an example, we show in Fig. 3 the different spectral sensitivities of IV 2 and I A for 15 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 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 r (e.g. Harp and Wilson, 1995). 20 ln This form is widely used (Keefer, 1984;Harp and Wilson, 1995). In engineering seismology, a typical ground motion model has an additional distance term for anelastic attenuation This is a modified version of the model template by Kramer (1996). 25 We exchange the magnitude term from Eq. (18) with a site-dependent energy term, as 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  . a) Far-field spectrum after 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).  (11)), since energy is proportional to the seismic moment M 0 (Eq. (9)). 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. 5

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 the latter approach and model directivity for estimated energy and corrected Arias intensity in a simplified way: E 0 and I A,A,0 are the offset (average), a E and a I the amplitude of variation with azimuth and θ E and θ E 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 for two reasons: (1) the rupture is assumed to have occurred on two faults and has thus variable strike, (2) 15 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 .

20
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.

25
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 like the directivity effect. The fault-normal response amplitude is larger than the fault-parallel response if directed (anti)parallel to the rupture. We model the ratio similar to Somerville et al. (1997)  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 like 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.
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 segment. Nearly all landslides concentrated on hillslopes between 15 • and 45 • inclination and MAF ≈ 1 (Fig. 4a,c). Hillslope inclination and MAF were higher towards the landslide crown (Fig. 4b,d), indicating a landslide failure process starting from 10 the crown and according to simulations by Dang et al. (2016). They showed a progressive failure for two major landslides of Kumamoto earthquake by numerical simulation. 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). Most landslides originated at localities with amplified ground accelerations and steep hillslopes and propagated progressively to flatter areas with less amplified ground accelerations. 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 co-occur the process (Martel, 2004).
Simulations of elliptic landslides by Martel (2004) show that the most compressive and the most tensile stresses are parallel to 5 the major axis of the landslide, coinciding with the average landslide aspect. Yamada et al. (2013) and Yamada et al. (2018) show for several japanese landslides that peak forces were aligned parallel to the long side of the landslides; Allstadt (2013) shows from waveform inversion for the Mt. Meager landslide that force and acceleration were parallel to the longer side of the landslide source area.
Mt. Aso and its caldera and Mt. Shutendōji had a high density of landslides (Fig. 5), whereas Mt. Kinpō and Mt.Ōtake lack 10 landslides, though these locations are closer to the epicenter and at comparable distances from the rupture (Fig. 5). However, all these locations have the same rock type and 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, because both differ widely despite near-identical conditions.    The azimuthal density-with respect to the epicenter-of the unspecified landslides follows to some extent the distribution of hillslope inclinations ≥ 19 • in the landslide affected area (Fig. 6). This similarity shows that the abundance of unspecified  (Fig. 7). The contrast between the distributions of unspecified landslides and earthquake related landslides indicates a contribution by the rupture process.

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, top right). Wave-10 forms with backward-directivity to the southwest of the rupture are longer with no dominant pulse (Fig. 8, bottom left). Waveforms parallel to the rupture have intermediate duration. Waveforms in either forward-or backward-direction have stronger amplitudes in the fault-normal direction, whereas waveforms outside the directivity-affected regions have stronger amplitudes 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. 15 The geometrical spreading A is calculated according to Eq. (12)   km. Any remaining distance dependence has been corrected for by estimating and applying the attenuation parameter k (Eq. (13)) After the determination of k,Ê and I A,A are considered distance independent and can be investigated for azimuthal variations. With a reference point for the azimuth at the epicenter,Ê shows oscillating variations in amplitude with azimuth ( Fig.   9a), while I A,A exhibits a similar amplitude variations over the entire azimuthal range (Fig.9b). The running average based on a but not as wide and large as that ofÊ. The azimuthal variation ofÊ indicates the rupture directivity and the absence of large variations in I A,A shows that the directivity effect is only evident at lower frequencies (compare with Fig. 3).
The azimuthal variation ofÊ and I A,A is modelled according to Eq. (20). We estimate parameters for two scenarios:

10
directivity is assumed, resulting in azimuthal variations and 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: 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 preferred. 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 5 of the fault (θ E = θ I = 225 • ).
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 10 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).
The low-frequency azimuthal variations are also reflected in the spectral response of the waveforms. Spectral accelerations of 5 stations with r rup ≤ 50 km were computed from 0.1 Hz -5 Hz with an interval of 0.01 Hz for the fault-normal and fault-parallel component. The distribution of F N/F P shows decreasing azimuthal variability with increasing oscillatory frequency (Fig.   11a). F N/F P shows highest variation 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 10 distance of 25.06 km, we compare it to the F N/F P model of Somerville et al. (1997) at 25 km (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 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 landslides. The azimuthal variation ofÊ coincides 15 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. Two directions show higher landslide density, one to the northwest and one to the east (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 with aspect east by southeast occurs at very short distances (r rup ≤ 2.5 km, Fig. 14), while the northwest oriented landslides are further away (2.5 km < r rup ≤ 6 km). Only minor landslides have larger distances with no specific pattern.
The distribution of aspect and hillslope inclination in the landslide affected area shows little variation over the entire aspect 5 range (Fig. 12b) with slightly more hillslopes facing westward. Neither the highly preferred orientation of landslides to the northwest nor to the east is reflected by the topography of the landslide affected area (12a,b). The unspecified landslides within the affected area have a near northward aspect and deviate by ≈ 30 • from the earthquake triggered landslides (Fig. 12c). This highlights that the earthquake not only affects the landslide locations (Fig. 6, 7), but also their orientation (Fig. 12). 10 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 model is a conventional isotropic moment magnitude

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 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 5 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 show a change in polarity orientation towards east-west, while the 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 short distances. 10 The impact of the near-field term may explain the increase in eastward aspect (Fig. 14) of landslides close to the rupture. indicating that Arias intensity is more influenced by local heterogeneities and scattering than the energy estimates as these 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 faultparallel amplitudes as is the radiated energy: Because of its higher sensitivity towards higher frequencies, these effects are 5 masked by high-frequency effects, e.g. 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 dependence has also been demonstrated by the response spectra ratios for F N/F P where directivity related amplitude variations with azimuth have been identified only for frequencies < 2 Hz which is 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 10 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 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 (r rup ), by a linear relationship ln P ls (r rup ) = a r + b r ln r rup (27) 15 Azimuthal variations of landslide density, P ls (θ), correspond to azimuthal variations in seismic energy and can be described by a similar relationship ln P ls (θ) = a θ + b θ cos(θ − θ 0 ) For the Kumamoto earthquake data we estimate a r = −5.4, b r = 2.6 and a θ = −46.8, b θ = 3.0. The azimuth-dependent landslide concentration implies that similar landslide concentrations can occur at different distances from the rupture, thus partly 20 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-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 is most likely due to the lack of an additional distant dependent attenuation term in their model (Eq. (17)). 25 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 showed 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.
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. For instance, Chen et al. (2017) found that landslide susceptibility 30 and safety factor calculation depends on whether the entire landslide or only parts-scarp area or area of dislocated massare considered. It is intractable to reconstruct each individual landslide failure plane from the mapped data. However, failure may have likely originated close to the crown and then progressively propagated downward the hillslope, because MAF > 1 indicates an amplification of ground motion towards the crown of the landslides. Additionally, (Sato et al., 2017) consider the tephra layers rich in halloysite to be the main sliding surfaces indicating shallow landslides (Song et al., 2017).
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, i.e. we could not trace any clear impact of soil moisture.

8 Conclusions
We investigated the Kumamoto earthquake and its triggered landslides. We demonstrated that the pattern of coseismic landslides is consistent with that of low-frequency ground motion. At low frequencies it exhibits two aspects: rupture directivity and increased amplitudes normal to the fault. While the former effect influences landslide locations, the latter effect relates to the landslide orientation (aspect). Topographic controls (hillslope inclination and MAF) are limited predictors of coseismic 10 landslide occurrence, because areas with similar topographic and geological properties at similar distances from the rupture showed 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. 15 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 spreading and one for along-path attenuation. The latter is not commonly incorporated 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 shows similarities to such ground-motion models. 20 The effect of low-frequency ground motion on the rupture process of the landslides results in landslide movements parallel to highest ground-motion. Due to the surface proximity of the earthquake rupture plane, near-field ground motion influences 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 demonstrated that earthquake triggered landslide hazard estimation requires an integrated approach of both detailed 25 ground-motion and topographic characterization. While the latter is well established for landslide hazard, ground-motion char-  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.
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, neous, 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.

10
(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 (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 15 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 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 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 We consider an isotropic medium and the elasticity tensor c from Eq. (A1) is c jkpq = δ jk δ pq λ + (δ jp δ kq + δ jq δ kp )µ (A8) where λ and µ are the Lamé parameters λ = ρ(v 2 P + 2µ) µ = ρv 2 S (A9) 10 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): 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

15
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 .