Effects of finite source rupture on landslide triggering: the 2016 M_w 7.1 Kumamoto earthquake

The propagation of a seismic rupture on a fault introduces spatial variations in the seismic wave field 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 (M_w 7.1) in central Kyushu (Japan). Although the distribution of some 1500 earthquake-triggered landslides as a function of rupture distance is consistent 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, the 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 head scarps have an elevated hillslope inclination and MAF. We propose a new physics-based ground-motion model (GMM) that accounts for the seismic rupture effects, and we demonstrate that the low-frequency 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.

1 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 the response of landslides to ground motion, and a growing number of studies have 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 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, and 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., 1986). This general pattern is modified by morphometrics (e.g. local hillslope inclination and curvature) and geological parameters (e.g. lithology, geological structure, and 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 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 16 April 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). 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 1500 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; Hovius and Meunier, 2012). Similarly asymmetric landslide distributions attributed to rupture directivity were reported 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 the 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 landslides. We investigate the geological conditions (lithology, aspect, hillslope inclination, topographic amplification, and soil wetness) 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 in 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 (GMM) that is consistent with the observed coseismic landslide pattern.

Figure 1 The area of Kyushu affected by coseismic landslides triggered by the 2016 M_w 7.1 Kumamoto earthquake. The coloured patch is the slip distribution of the rupture model by Kubo et al. (2016), and the dashed box encompasses landslides related to the triggered event in Yufu (epicentre location after Uchide et al. (2016)). The inset map shows the station network within 150 km of the rupture.

2 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).

Figure 2 Topographic and geological features of central Kyushu with landslides (black dots), landslide-affected area (outer black line), rupture area (inner black line), hypocentre (black diamond), and mountain peaks from Fig. 1 (triangles). (a) Hillslope inclination. (b) Median amplification factor (MAF). (c) Topographic wetness index (TWI). (d) Geology of central Kyushu. The most common geological units of the landslides are shown in (e). For the landslide-affected area the dominant geological units are in (f). (g) Land cover. Land cover in the landslide areas is shown in (h) and is shown for the entire landslide-affected area in (i).

2.1 Topographic data

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 and agriculture) and vegetation, while data on major geological units are from the Seamless Digital Geological Map of Japan (scale of 1:200 000) by the Geological Survey of Japan.

2.2 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., 2012, 2015). The 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 travelling at frequency f:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{MAF}\left(f\right)=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}{\displaystyle \frac{v_{\mathrm{S}}}{f}}{C}_{\mathrm{S}}\left({\displaystyle \frac{v_{\mathrm{S}}}{\mathrm{2}f}}\right)+\mathrm{1},\end{array}$$

where $C_{\mathrm{S}}\left(\frac{v_{\mathrm{S}}}{\mathrm{2}f}\right)$ is the topographic curvature convolved with a normalized smoothing kernel based on two 2-D boxcar 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 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 modelled for uniform conditions to the first order using the topographic wetness index (TWI) of Böhner and Selige (2006):

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{TWI}=\log {\displaystyle \frac{A_{\mathrm{c}}}{\tan \mathit{\beta }}},\end{array}$$

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

2.3 Ground-motion data

Ground-motion data are from the KiK-Net and K-Net of the National Research Institute for Earth Science and Disaster Prevention (NIED) of Japan. NIED operates both borehole and surface stations for KiK-Net, 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 Kyushu are available with complete azimuthal coverage within 150 km of 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 of 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 piecewise 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 (2019).

An earthquake was triggered approximately 80 km to the northeast in Yufu 32 s after the Kumamoto earthquake (Fig. 1; Uchide et al., 2016). 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 travel time differences between the P wave (v_P=5700 m s⁻¹) arrival of the Kumamoto earthquake and the S wave arrival (v_S=3300 m s⁻¹) of the triggered event. The respective travel paths to the stations are measured from the hypocentres. Since fewer instruments are located to the northeast and the triggered event close 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 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 the rupture-plane model, which 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) is available for 23 layers down to the mantle in $\approx \mathrm{0.1}{}^{\circ }$ resolution covering all of Japan; we only consider the layers of the upper 0.5 km to compute a velocity average for the MAF.

NIED provides data for the subsurface shear wave velocities (v_S30) as well as site amplification factors S_amp. Contrary to v_S by Koketsu et al. (2012), v_S30 is derived for the upper 30 m only and is more suitable for energy estimates, which require velocities at the surface (recording station). The site amplification factor S_amp describes how much seismic waves are amplified by, independent of their frequency.

2.4 Landslide data

Detailed landslide data are provided by NIED as polygons (Fig. 1), mapped from aerial imagery with sub-metre resolution at different times after the Kumamoto earthquake. The first data set contains landslides that were identified between 16 and 20 April, though the area close to the summit of Mt Aso was not covered. A second data set was collected on 29 April 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 landslides triggered by volcanic activity. We selectively combine the two data sets for this study, using only those landslides from the second database, which 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), which were 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 in landslides (≈50 km) between Yufu and Aso (Fig. 1) in otherwise steep topography.

Apart from the special release of landslide data for the 2016 Kumamoto earthquake, NIED hosts a landslide database for all Japan (National Research Institute for Earth Science and Disaster Prevention, 2014). This database covers unspecified 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 mapped 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 and because the landslide source area is generally not identical to the deposit area.

3 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² and is 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 area of 3914 km² (4406 km²) 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 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², i.e. A_topo=0.68.

4 Total landslide area

Total landslide area is linked to several earthquake parameters, mostly magnitude and hypocentre or average rupture-plane 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². 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 to be 3 km. However, Somerville et al. (1999) derived a relationship of asperity sizes based on the seismic moment that results in an average asperity length of 12.78 km for the 2016 Kumamoto earthquake. This length is consistent with the asperity sizes found by Yoshida et al. (2017) for their finite rupture model. The estimated landslide area with an asperity length of 3 km results in a predicted total landslide area of 12.90 km², while with the magnitude scaled asperity size of Somerville et al. (1999), the landslide area is 3.03 km². The landslide area estimates with constant asperity length and moment-dependent asperity length differ by a factor of 2 and 0.5 from the NIED data set, respectively.

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

5 Ground motion and seismically induced landsliding

5.1 Coseismic landslide displacement

The sliding-block model of Newmark (1965) is widely used to estimate coseismic hillslope performance (e.g. Kramer, 1996; Jibson, 1993, 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 pulse, v_max (m s⁻¹):

$$\begin{array}{}\text{(3)}& {\displaystyle }{d}_{\mathrm{s}}={\displaystyle \frac{v_{\max }^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{1}}{A{a}_{\mathrm{y}}}}-{\displaystyle \frac{\mathrm{1}}{A}}\right),\end{array}$$

where A is the magnitude of the acceleration pulse and a_y (m s⁻²) is 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

$$\begin{array}{}\text{(4)}& {\displaystyle }{a}_{\mathrm{y}}=g\left({\displaystyle \frac{\tan {\mathit{\varphi }}_{\mathrm{f}}}{\tan \mathit{\delta }}}\right)\sin \mathit{\delta }=g(\stackrel{\mathrm{‾}}{\mathrm{FS}}-\mathrm{1})\sin \mathit{\delta },\end{array}$$

with the average factor of safety $\stackrel{\mathrm{‾}}{\mathrm{FS}}$. Chen et al. (2017) characterized unstable hillslopes – related to both rainfall and earthquakes – by a safety factor of FS<1.5.

An upper bound for the displacement, d_s, is based on two ground-motion parameters (Newmark, 1965; Kramer, 1996):

$$\begin{array}{}\text{(5)}& {\displaystyle }{d}_{\max }={\displaystyle \frac{\mathrm{PGA}}{a_{\mathrm{y}}}}{\displaystyle \frac{{\mathrm{PGV}}^{\mathrm{2}}}{a_{\mathrm{y}}}},\end{array}$$

where PGA (m s⁻²) and PGV (m s⁻¹) 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 model based on the acceleration-related Arias intensity and the velocity-related radiated seismic energy.

5.2 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:

$$\begin{array}{}\text{(6)}& {\displaystyle }{I}_{\mathrm{A}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}g}}\underset{T_{\mathrm{1}}}{\overset{T_{\mathrm{2}}}{\int }}a(t{)}^{\mathrm{2}}\mathrm{d}t,\end{array}$$

where g=9.80665 m s⁻² is standard gravity and T₁ and T₂ are the times where strong ground motion starts and cedes (the acceleration a(t) has units of m s⁻² and the Arias intensity has units of m s⁻¹). 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 epicentre distance have been developed (e.g. Bray and Travasarou, 2007; Jibson, 1993, 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 (Eqs. 3 and 5), it is reasonable to characterize velocity similarly to Arias intensity. The velocity counterpart to I_A is IV2, the integrated squared velocity (Kanamori et al., 1993; Festa et al., 2008):

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathrm{IV}\mathrm{2}=\underset{T_{\mathrm{1}}}{\overset{T_{\mathrm{2}}}{\int }}v(t{)}^{\mathrm{2}}\mathrm{d}t.\end{array}$$

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), and Newman and Okal (1998):

$$\begin{array}{}\text{(8)}& {\displaystyle }{j}_{\mathrm{E}}={\displaystyle \frac{\mathit{\rho }c}{S_{\mathrm{amp}}^{\mathrm{2}}}}{e}^{-k{r}_{\mathrm{rup}}}\mathrm{IV}\mathrm{2},\end{array}$$

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:

$$\begin{array}{}\text{(9)}& {\displaystyle }{E}_{\mathrm{S}}=\iint j_{\mathrm{E}}\mathrm{d}A,\end{array}$$

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):

$$\begin{array}{}\text{(10)}& {\displaystyle }{E}_{\mathrm{S}}={\displaystyle \frac{\mathrm{\Delta }\mathit{\sigma }}{\mathrm{2}\mathit{\mu }}}{M}_{\mathrm{0}},\end{array}$$

where Δσ is the stress drop, μ the shear modulus, and M₀ the seismic moment. 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 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 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 $\widehat{E}$ based on Eqs. (8) and (9) becomes

$$\begin{array}{}\text{(11)}& {\displaystyle }\widehat{E}={\displaystyle \frac{\widehat{A}\mathit{\rho }{v}_{\mathrm{S}}}{S_{\mathrm{amp}}^{\mathrm{2}}}}{e}^{-k{r}_{\mathrm{rup}}}\mathrm{IV}\mathrm{2}.\end{array}$$

While E_S is the radiated seismic energy at the source, $\widehat{E}$ is estimated from the velocity records at a station and only approximates E_S. Therefore, $\widehat{E}$ may differ from the true and unknown radiated energy E_S (Kanamori et al., 1993). Several assumptions characterize $\widehat{E}$:

• All energy is radiated as S waves in an isotropic, homogeneous medium.

• Geometrical spreading is corrected for an isotropic, homogeneous medium.

• Since IV2 (Eq. 7) depends on the radiation pattern, $\widehat{E}$ depends on the azimuth.

• Attenuation is homogeneous.

• Surface waves are not considered.

• Site amplification is frequency-independent.

Below, we investigate the azimuthal variation in the energy estimates to characterize the radiation pattern.

The estimated wavefront area $\widehat{A}$ is related to the rupture extent and r_rup, and $\widehat{A}$ corresponds to a simplified version of the wavefront area approximation by Schnabel and Bolton Seed (1973) and Shoja-Taheri and Anderson (1988):

$$\begin{array}{}\text{(12)}& {\displaystyle }\widehat{A}=\mathrm{2}WL+\mathit{\pi }{r}_{\mathrm{rup}}(L+\mathrm{2}W)+\mathrm{2}\mathit{\pi }{r}_{\mathrm{rup}}^{\mathrm{2}}.\end{array}$$

The extent of the rupture is assumed to be rectangular with length L and width W. Equation (12) describes a cuboid with rounded corners and with only half of its surface considered because no energy flux is assumed to be transmitted above ground.

While the geometrical spreading correction is expressed analytically as the wavefront area $\widehat{A}$, 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

$$\begin{array}{}\text{(13)}& {\displaystyle }\ln Y=C+k{r}_{\mathrm{rup}},\end{array}$$

where Y is

$$\begin{array}{}\text{(14)}& {\displaystyle }Y={\displaystyle \frac{\widehat{A}\mathit{\rho }{v}_{\mathrm{S}}}{S_{\mathrm{amp}}^{\mathrm{2}}}}\underset{T_{\mathrm{1}}}{\overset{T_{\mathrm{2}}}{\int }}\mathrm{IV}\mathrm{2},\end{array}$$

i.e. the logarithm of the energy estimate without the attenuation term $e^{-k{r}_{\mathrm{rup}}}$ 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

$$\begin{array}{}\text{(15)}& {\displaystyle }{I}_{\mathrm{A},\mathrm{A}}={\displaystyle \frac{\widehat{A}}{S_{\mathrm{amp}}^{\mathrm{2}}}}{e}^{-k{r}_{\mathrm{rup}}}{I}_{\mathrm{A}},\end{array}$$

where k is determined by Eq. (13) and setting $Y=I_{\mathrm{A}}\widehat{A}$. The corrected Arias intensity I_A,A is the acceleration-based counterpart to $\widehat{E}$.

Figure 3 (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 IV2 (blue) and the acceleration-based I_A (black).

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 transform of the time derivative of any function is

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathcal{F}\left(\dot{f}\right(t\left)\right)=i\mathit{\omega }\mathcal{F}\left(f\right(t\left)\right),\end{array}$$

and thus scales with frequency in the spectrum. The frequency sensitivity of IV2 and I_A is related to the squared spectrum given the metrics. For example, in Fig. 3 we show the different spectral sensitivities of IV2 and I_A for a theoretical seismic source spectrum (Brune, 1970). IV2, and thus $\widehat{E}$, has a higher sensitivity to lower frequencies than I_A. The low-frequency part of the spectrum can be accounted for by considering IV2 in a ground-motion model.

5.3 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):

$$\begin{array}{}\text{(17)}& {\displaystyle }\ln I_{\mathrm{A}}=p_{\mathrm{1}}+p_{\mathrm{2}}M+p_{\mathrm{3}}\ln r.\end{array}$$

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:

$$\begin{array}{}\text{(18)}& {\displaystyle }\ln I_{\mathrm{A}}=c_{\mathrm{1}}+c_{\mathrm{2}}M+c_{\mathrm{3}}r+(c_{\mathrm{4}}+c_{\mathrm{5}}M)\ln .\end{array}$$

This is a modified version of the model template by Kramer (1996). While Eqs. (17) and (18) share some parameters (p₁, c₁ and p₂, c₂), the geometric spreading term includes not only distance dependence (p₃, c₄) but also has a magnitude-dependent component (c₅). In addition, anelastic attenuation is included as well (related to c₃) in Eq. (18). The template of Kramer (1996) relates to the majority of GMMs 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 is included in more recent studies (Meunier et al., 2007, 2013; Yuan et al., 2013).

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₀ (Eq. 10). Based on the site-dependent energy estimate $\widehat{E}$, we propose the model

$$\begin{array}{}\text{(19)}& {\displaystyle }\ln I_{\mathrm{A}}=c_{\mathrm{1}}+c_{\mathrm{2}}\ln \widehat{E}+c_{\mathrm{3}}r+(c_{\mathrm{4}}+c_{\mathrm{5}}\ln \widehat{E})\ln r.\end{array}$$

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.4 Rupture directivity model

In the NGA-west2 guidelines (Spudich et al., 2013), the directivity effect is modelled by isochrone theory (Bernard and Madariaga, 1984; Spudich and Chiou, 2008) or the azimuth between epicentre 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:

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }\ln {\widehat{E}}_{\mathit{\theta }}=\ln {\widehat{E}}_{\mathrm{0}}+a_E\cos (\mathit{\theta }-{\mathit{\theta }}_E),\text{(21)}& {\displaystyle }& {\displaystyle }\ln I_{\mathrm{A},\mathrm{A},\mathit{\theta }}=\ln I_{\mathrm{A},\mathrm{A},\mathrm{0}}+a_I\cos (\mathit{\theta }-{\mathit{\theta }}_I),\end{array}$$

where ${\widehat{E}}_{\mathrm{0}}$ and $I_{\mathrm{A},\mathrm{A},\mathrm{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) for the angle measured between the epicentre and the recording site, with the difference of being measured clockwise from the 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).

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

Figure 4 Distribution of hillslope inclination and MAF. The left column shows (a) hillslope inclinations and (c) MAF within the landslide-affected area (green) and within the landslide areas (black). The right column presents (b) hillslope inclinations and (d) MAF in different segments of the landslide areas which is expressed as relative height. Segments towards to the toe (relative height 0.0–0.5) are in green, and segments towards the crown are in red (relative height 0.5–1). The solid line is the mean, and the dashed lines enclose the 95 % uncertainty range. The concept of relative height is illustrated for the Aso Ohashi landslide in (e). MAF<1 indicates attenuation and MAF>1 amplification of seismic waves due to topography. The cyan line in (d) highlights MAF=1, i.e. no amplification or attenuation.

Figure 5 (a) Landslide concentration with (a) rupture distance r_rup and (b) asperity distance r_asp of the Kumamoto earthquake landslides. The rate parameter of 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 similar kernel density estimates and the near-identical rate parameter estimates $\widehat{\mathit{\lambda }}$. 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.

Figure 6 Kernel density estimates of azimuth and distance of (a) landslide concentration of the coseismic landslides, (b) concentration of landslide-susceptible terrain with hillslope inclinations >19^∘, and (c) landslide concentration of unspecified landslides. Azimuth and distance are with respect to the asperity centroid. The marginal densities with respect to azimuth are shown in blue as outer ring. The densities are normalized to their maxima.

Figure 7 Spatial distribution of landslides. (a) Coseismic landslides. The total landslide area at a location is shown as a colour-coded smooth function in the background; (b) same as in (a) but for unspecified landslides within the landslide-affected area of the Kumamoto earthquake.

The north and east components (E, N) of the sensor are rotated to be fault-normal (FN) and fault-parallel (FP) with fault strike ϕ.

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }\mathrm{FN}=E\cos \mathit{\varphi }-N\sin \mathit{\varphi },\text{(23)}& {\displaystyle }& {\displaystyle }\mathrm{FP}=E\sin \mathit{\varphi }+N\cos \mathit{\varphi },\text{(24)}& {\displaystyle }& {\displaystyle }\mathrm{FN}/\mathrm{FP}={\displaystyle \frac{{\mathrm{SA}}_{\mathrm{FN}}\left(f_{\mathrm{osc}}\right)}{{\mathrm{SA}}_{\mathrm{FP}}\left(f_{\mathrm{osc}}\right)}}.\end{array}$$

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 than the fault-parallel response if directed parallel or antiparallel to the rupture. We model the ratio similar to Somerville et al. (1997):

$$\begin{array}{}\text{(25)}& {\displaystyle }\ln (\mathrm{FN}/\mathrm{FP})=(b_{\mathrm{1}}+b_{\mathrm{2}}{f}_{\mathrm{osc}}^{b_{\mathrm{3}}}\cos (\mathrm{2}(\mathit{\theta }-{\mathit{\theta }}_R)\left)H\right(b_{\mathrm{1}}+b_{\mathrm{2}}{f}_{\mathrm{osc}}^{b_{\mathrm{3}}}),\end{array}$$

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

We introduced a functional form for oscillatory frequency dependence with four parameters in Eq. (20). We did not introduce a distance term and apply the model only to data with r_rup≤ 50 km.

6 Results

6.1 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 with a steepness between 15 and 45^∘ and an MAF ≈1 (Fig. 4a, c). Hillslope inclination and the 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). The TWI is linked to land cover and is highest in areas with rice paddies (Fig. 2i). Terrain with landslides has a 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 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 axis of the landslide, coinciding with the average landslide aspect. Yamada et al. (2013, 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 long side of the landslide source area.

Figure 8 Characteristic waveforms observed in the vicinity of the rupture exhibiting rupture directivity effects. Station 93002 is in the forward-directivity region with a large amplitude pulse on the fault-normal component. Station 93039 is also in the forward-directivity region but with an offset to the rupture. In this region the fault-parallel component has higher amplitudes. The station KMM012 is in the backward-directivity region, and waveforms have longer duration without large amplitudes. The waveforms are low-pass filtered at 1.2 Hz.

Figure 9 Energy estimates ($\widehat{E}$) over azimuth; (b) same as in (a) but for Arias intensity with correction for geometrical spreading (I_A,A).

Figure 10 Kernel density estimate of the amplitude ratio of response spectra of fault-normal and fault-parallel components (FN∕FP) with respect to oscillatory frequency. Beyond 2–3 Hz FN∕FP variations cease, as highlighted by our model and the model by Somerville et al. (1997).

Figure 11 Kernel density estimate of FN∕FP with azimuth obtained from response spectra for three different oscillatory frequency ranges: (a) 0.1–1 Hz, (b) 1.0–2.5 Hz, and (c) >2.5 Hz. For each plot, our FN∕FP model and the model by Somerville et al. (1997) are shown for (a) 0.55 Hz, (b) 1.75 Hz, and (c) 4 Hz. As in Fig. 10, amplitudes decrease with increasing oscillatory frequency.

Figure 12 (a) Aspect and hillslope inclination distribution within areas of the earthquake-triggered landslides. This distribution is normalized by the distribution of the aspect of all hillslopes in the landslide-affected area. The black line denotes the strike of the Kumamoto earthquake (225^∘). (b) Distribution of aspect and hillslope inclination in the landslide-affected area; (c) same as in (a) but for unspecified landslides.

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 epicentre and being comparably close to the rupture (Fig. 5). All these locations have comparable rock type, land cover, hillslope inclination, and MAFs. Hence, lithology, land cover, and topographic characteristics are insufficient in explaining the landslide distribution and concentration with respect to the hypocentre or the asperity.

The azimuthal density – with respect to the asperity centroid – of the unspecified landslides follows to some extent the distribution of hillslope inclinations $>\mathrm{19}{}^{\circ }$ in the landslide-affected area (Fig. 6b, c). This similarity shows that the abundance of unspecified landslides mimics the steepness of topography in the region. Densities are higher towards Mt Kinpo (W), Mt Otake (WSW), Mt Shutendoji (N), Mt Aso (E), and the Kyushu mountains (SE). The coseismic landslide distribution differs completely from the distributions of unspecified landslides and their surrounding topography (Fig. 6), respectively, as nearly all landslides happened to the northeast of the epicentre, close to the rupture plane (Fig. 7). Chen et al. (2017) identified only 29 landslide reactivations during the Kumamoto earthquake. The contrast between the distributions of unspecified landslides and earthquake-related landslides indicates a contribution by the rupture process.

6.2 Impact of finite source on ground motion and landslides

The results of the seismic analysis are given for waveforms, the basis for $\widehat{E}$ and I_A, and response spectra, used for FN∕FP. To the northeast, signals with forward directivity are shorter in duration, with one or a few strong pulses (Fig. 8, top right). Waveforms 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 an intermediate duration. Waveforms in either a 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 $\widehat{E}$ 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 km. Any remaining distance dependence has been corrected for by estimating and applying the attenuation parameter k (Eq. 13)

After the determination of k, $\widehat{E}$ and I_A,A are considered distance-independent and can be investigated for azimuthal variations. With a reference point for the azimuth at the epicentre, $\widehat{E}$ 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 von Mises kernel (κ_vM=50) of $\widehat{E}$ and I_A,A shows increased $\widehat{E}$ between 45 and 135^∘, i.e. approximately parallel to the strike. Minimal values of $\widehat{E}$ occur in the opposite direction (200–300^∘). The running average of I_A,A has several fluctuations, but these are not as wide and large as those of $\widehat{E}$. The azimuthal variation in $\widehat{E}$ indicates the rupture directivity, and the absence of large variations in I_A,A indicates that the directivity effect is only evident at lower frequencies (compare with Fig. 3).

The azimuthal variation in $\widehat{E}$ 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=\mathrm{0}$.

The two models are compared with the Bayesian information criterion (BIC; Schwarz, 1978) for a least-squares fit:

$$\begin{array}{}\text{(26)}& {\displaystyle }\mathrm{BIC}=n\ln N+N\ln {\widehat{\mathit{\sigma }}}^{\mathrm{2}},\end{array}$$

where n is the number of estimated parameters (n=4 for the first case and n=2 for the second case), N is the number of data, and ${\widehat{\mathit{\sigma }}}^{\mathrm{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 $\widehat{E}$ and I_A,A, no azimuthal variation (a_d=0), and the azimuths of the maximum of ${\widehat{E}}_{\mathit{\theta }}$ and $I_{\mathrm{A},\mathrm{A},\mathit{\theta }}$ are set to the strike of the fault (${\mathit{\theta }}_E={\mathit{\theta }}_I=\mathrm{225}{}^{\circ }$).

The directivity model for $\widehat{E}$ 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${}_{\mathrm{directivity}}=-\mathrm{110}$, BIC${}_{\mathrm{no}\mathrm{directivity}}=-\mathrm{11}$). In the 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).

The low-frequency azimuthal variations are also reflected in the spectral response of the waveforms. Spectral accelerations of stations with r_rup≤50 km were computed from 0.1 to 5 Hz at intervals of 0.01 Hz for the fault-normal and fault-parallel component. The distribution of FN∕FP shows decreasing azimuthal variability with increasing oscillatory frequency (Fig. 11). FN∕FP 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 FN∕FP model (Eq. 20; 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 FN∕FP model of Somerville et al. (1997) at 25 km (Figs. 10, 11). Both models show a similar decay with frequency, with our model predicting a slightly higher FN∕FP. 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).

Figure 13 Orientation of horizontal peak ground acceleration for the simulated waveforms. The arrow length scales with magnitude of acceleration. The simulated rupture plane is oriented as the rupture plane of the Kumamoto earthquake (strike: 225^∘, dip: 70^∘) and of elliptic shape (grey). The upper side is denoted by the green line, and the lower half is denoted by black. The rupture process originated at the hypocentre (red dot) with circular propagation outwards (white arrow).

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, truncating distances below 1 km. The distribution has been normalized by the distribution of aspect of the affected area.

Figure 15 Ground-motion model for I_A. The solid lines are the model with energy estimates for three different energy levels as in Fig. 9a. The inset figure shows, for comparison, the ground-motion model of Harp and Wilson (1995) (green) and landslide concentration density (red).

The pattern of low-frequency ground motion is well reflected in that of the landslides. The azimuthal variation in $\widehat{E}$ coincides with that of landslide concentration (Fig. 9). Both azimuth-dependent energy and landslide concentration have a similar trend, with the maximum being parallel to rupture direction and the minimum strike being antiparallel. The orientation of maximum FN∕FP 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 FN∕FP, 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 ($\mathrm{2.5}\mathrm{km}<r_{\mathrm{rup}}\le \mathrm{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 (Fig. 12a, b). The unspecified landslides in the affected area have a near-northward aspect and deviate by $\approx \mathrm{30}{}^{\circ }$ from the earthquake-triggered landslides (Fig. 12c). This highlights that the earthquake affects landslide locations (Fig. 6) and will force failure on specific slopes facing in the direction of ground motion (Figs. 12, 14).

6.3 Ground-motion model for Kumamoto

Table 1 Parameters for ground-motion models.

We derived two ground-motion models for Arias intensity from data with r_rup≤150 km (Table 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 the running average well and is proportional to the decrease in landslide density with distance. Variation in estimated energy is well covered by the model and spans more than 2 orders of magnitude, resulting in a variation in Arias intensity of nearly 1 order of magnitude.

The magnitude-based model is nearly equivalent to the energy-based model with $\widehat{E}=\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{15}}$ J. This value is close to the average energy estimate found from energy estimates of the directivity model from Eq. (20) ($\widehat{E}=\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{15}}$ J). The closeness of the two values implies that the magnitude-based model can be seen as an average over the azimuth of the energy-based model.

7 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 a 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), and 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 motion captures 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 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 in simulated $\widehat{E}$ 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 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 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 dependence is also seen in the response spectra ratios for FN∕FP where directivity-related amplitude variations with azimuth have been identified only for frequencies <2 Hz, 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 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:

$$\begin{array}{}\text{(27)}& {\displaystyle }\ln P_{\mathrm{ls}}\left(I_{\mathrm{A}}\right)=a_I+b_I\ln I_{\mathrm{A}}.\end{array}$$

Azimuthal variations in landslide density correspond to azimuthal variations in seismic energy and can be described by a similar relationship:

$$\begin{array}{}\text{(28)}& {\displaystyle }\ln P_{\mathrm{ls}}\left(E\right)=a_E+b_E\cos (\mathit{\theta }-{\mathit{\theta }}_E).\end{array}$$

For the Kumamoto earthquake data, we estimate a_I=2.1, b_I=2.6, and $a_E=-\mathrm{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 Figs. 5 and 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 pseudo-depth, the model of Harp and Wilson (1995) overpredicts I_A both at shorter and longer distances – irrespective of the pseudo-depth 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 the 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 the 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., 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 depend on whether the entire landslide or only parts – scarp area or area of dislocated mass – are 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 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).

8 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 inclination and the 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 (Havenith et al., 2016; 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 state-of-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 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 the strongest 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 northwestern 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 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, 1964, 1969; Anderson and Richards, 1975; Aki and Richards, 2002). The displacement for any point x at time t is given by

$$\begin{array}{}\text{(A1)}& {\displaystyle }{u}_i(x,t)=\underset{\mathrm{\Sigma }}{\iint }{c}_{jkpq}{\displaystyle \frac{\partial G_{ip}\left(D_j\right(\mathit{\xi },t\left)\right)}{\partial x_q}}{n}_k\mathrm{d}\mathrm{\Sigma },\end{array}$$

where c is the fourth-order elasticity tensor from Hooke's law, G is Green's function describing the response of the medium, D(ξ,t) is the displacement on the fault with area Σ and coordinates ξ, and n is the fault-normal vector. Summation over i, j, p, and q is implied.

While the surface integral is carried out numerically, the derivatives of Green's function for an isotropic, homogeneous, and unbound medium can be solved analytically:

$$\begin{array}{ll}\text{(A2)}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial x_q}}{G}_{ip}\left(D_j\left(\mathit{\xi },t\right)\right)={\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{15}{\mathit{\gamma }}_i{\mathit{\gamma }}_p{\mathit{\gamma }}_q-\mathrm{3}\left({\mathit{\delta }}_{ip}{\mathit{\gamma }}_q+{\mathit{\delta }}_{iq}{\mathit{\gamma }}_p+{\mathit{\delta }}_{pq}{\mathit{\gamma }}_i\right)}{\mathrm{4}\mathit{\pi }\mathit{\rho }{r}^{\mathrm{4}}}}\end{array}$$

$$\begin{array}{}\text{(A2a)}& {\displaystyle }& {\displaystyle }\cdot \underset{\frac{r}{\mathit{\alpha }}}{\overset{\frac{r}{\mathit{\beta }}}{\int }}{D}_j\left(\mathit{\xi }t-\mathit{\tau }\right)\mathit{\tau }\mathrm{d}\mathit{\tau }\text{(A2b)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{6}{\mathit{\gamma }}_i{\mathit{\gamma }}_p{\mathit{\gamma }}_q-\left({\mathit{\delta }}_{ip}{\mathit{\gamma }}_q+{\mathit{\delta }}_{iq}{\mathit{\gamma }}_p+{\mathit{\delta }}_{pq}{\mathit{\gamma }}_i\right)}{\mathrm{4}\mathit{\pi }\mathit{\rho }{\mathit{\alpha }}^{\mathrm{2}}{r}^{\mathrm{2}}}}{D}_j\left(\mathit{\xi },t-{\displaystyle \frac{r}{\mathit{\alpha }}}\right)\text{(A2c)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{6}{\mathit{\gamma }}_i{\mathit{\gamma }}_p{\mathit{\gamma }}_q-\left(\mathrm{2}{\mathit{\delta }}_{ip}{\mathit{\gamma }}_q+{\mathit{\delta }}_{iq}{\mathit{\gamma }}_p+{\mathit{\delta }}_{pq}{\mathit{\gamma }}_i\right)}{\mathrm{4}\mathit{\pi }\mathit{\rho }{\mathit{\beta }}^{\mathrm{2}}{r}^{\mathrm{2}}}}{D}_j\left(\mathit{\xi },t-{\displaystyle \frac{r}{\mathit{\beta }}}\right)\end{array}$$

$$\begin{array}{}\text{(A2d)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{\gamma }}_i{\mathit{\gamma }}_p{\mathit{\gamma }}_q}{\mathrm{4}\mathit{\pi }\mathit{\rho }{\mathit{\alpha }}^{\mathrm{3}}r}}{\dot{D}}_j\left(\mathit{\xi },t-{\displaystyle \frac{r}{\mathit{\alpha }}}\right)\text{(A2e)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{{\mathit{\gamma }}_i{\mathit{\gamma }}_p{\mathit{\gamma }}_q-{\mathit{\delta }}_{ip}{\mathit{\gamma }}_q}{\mathrm{4}\mathit{\pi }\mathit{\rho }{\mathit{\beta }}^{\mathrm{3}}r}}{\dot{D}}_j\left(\mathit{\xi },t-{\displaystyle \frac{r}{\mathit{\beta }}}\right),\end{array}$$

where

$$\begin{array}{}\text{(A3)}& {\displaystyle }r=\left|\mathbf{x}-\mathit{\xi }\right|\mathrm{and}{\mathit{\gamma }}_i={\displaystyle \frac{x_i-{\mathit{\xi }}_i}{r}},\end{array}$$

and δ_ij is Kronecker's delta. The terms in Eq. (A2) are commonly separated into groups with respect to their distance r. In Eq. (A2) is the near-field (NF) term; as its amplitude decays with r⁻⁴, it affects the immediate vicinity of a rupture only. Terms with a distance attenuation proportional to r⁻² are called intermediate-field (IF) terms for P waves (Eq. A2a) and S waves (Eq. A2a). The remaining two terms are the far-field (FF) terms for P waves (Eq. A2b) and S waves (Eq. A2b) with a decay proportional to r⁻¹. 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.

Figure A1 Set-up of the rupture model. Grey ellipse represents the rupture: light grey area is unruptured, medium grey area is slipping, and the dark grey area is after slip arrest.

The slip function in time is related to the Yoffe function of Yoffe (1951) and 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

$$\begin{array}{}\text{(A4)}& {\displaystyle }D\left(\mathit{\xi }\right)=D_{\mathrm{0}}\sqrt{\mathrm{1}-{\left({\displaystyle \frac{{\mathit{\xi }}_{\mathrm{1}}-p\mathit{ϵ}\frac{L}{\mathrm{2}}}{\frac{L}{\mathrm{2}}}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{{\mathit{\xi }}_{\mathrm{2}}}{\frac{W}{\mathrm{2}}}}\right)}^{\mathrm{2}}},\end{array}$$

where D₀ is the maximum displacement at the centre of the fault, L and W are the length and width of the fault with eccentricity $\mathit{ϵ}=\sqrt{\mathrm{1}-{\left(\frac{W}{L}\right)}^{\mathrm{2}}}$, 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-antiparallel, $p=-\mathrm{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 $\widehat{\mathit{s}}$ describes the orientation of the displacement D(ξ) on the fault plane. We follow the definition of $\widehat{\mathbf{n}}$ and $\widehat{\mathbf{s}}$ in terms of fault strike ϕ_s, dip δ, and rake λ from Aki and Richards (2002):

$$\begin{array}{}\text{(A5)}& {\displaystyle }\widehat{\mathbf{n}}=\left(\begin{array}{c}-\sin \mathit{\delta }\sin {\mathit{\varphi }}_s\\ \sin \mathit{\delta }\cos {\mathit{\varphi }}_s\\ -\cos \mathit{\delta }\end{array}\right),\end{array}$$

$$\begin{array}{}\text{(A6)}& {\displaystyle }\widehat{\mathbf{s}}=\left(\begin{array}{c}\cos \mathit{\lambda }\cos {\mathit{\varphi }}_s+\cos \mathit{\delta }\sin \mathit{\lambda }\sin {\mathit{\varphi }}_s\\ \cos \mathit{\lambda }\sin {\mathit{\varphi }}_s-\cos \mathit{\delta }\sin \mathit{\lambda }\cos {\mathit{\varphi }}_s\\ -\sin \mathit{\lambda }\sin \mathit{\delta }\end{array}\right).\end{array}$$

The displacement vector D in Eq. (A2) is given by

$$\begin{array}{}\text{(A7)}& {\displaystyle }\mathit{D}=D\left(\mathit{\xi }\right)\widehat{\mathbf{s}}.\end{array}$$

We consider an isotropic medium, and the elasticity tensor c from Eq. (A1) is

$$\begin{array}{}\text{(A8)}& {\displaystyle }{c}_{jkpq}={\mathit{\delta }}_{jk}{\mathit{\delta }}_{pq}{\mathit{\lambda }}_{\mathrm{M}}+({\mathit{\delta }}_{jp}{\mathit{\delta }}_{kq}+{\mathit{\delta }}_{jq}{\mathit{\delta }}_{kp}){\mathit{\mu }}_{\mathrm{M}},\end{array}$$

where λ_M and μ_M are the Lamé constants of the isotropic medium:

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{M}}=\mathit{\rho }(v_{\mathrm{P}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{M}}),{\mathit{\mu }}_{\mathrm{M}}=\mathit{\rho }{v}_{\mathrm{S}}^{\mathrm{2}}.\end{array}$$

We set λ_M=μ_M, resulting in the widely observed relation $v_{\mathrm{P}}=v_{\mathrm{S}}\sqrt{\mathrm{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, including length L, width W, strike ϕ_s, and dip δ;

• material, including first and second Lamé constants λ_M and μ_M and density ρ (alternatively: compressional and shear wave velocities v_P and v_S and density ρ);

• rupture and slip, including rupture velocity ζ, slip D₀, rise time T, rake λ and rupture orientation with respect to strike ϕ_s and rupture orientation parameter 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 10 (9 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₀. The moment can be expressed as

$$\begin{array}{}\text{(A10)}& {\displaystyle }{M}_{\mathrm{0}}={\mathit{\mu }}_{\mathrm{M}}A\stackrel{\mathrm{‾}}{D},\end{array}$$

with shear modulus (second Lamé constant) μ_M, the rupture area – here an ellipse – $A=\frac{\mathit{\pi }}{\mathrm{4}}LW$, and average displacement $\stackrel{\mathrm{‾}}{D}$, which follows from Eq. (A3) as $\stackrel{\mathrm{‾}}{D}=\frac{\mathrm{2}}{\mathrm{3}}{D}_{\mathrm{0}}$.

The results are not strictly comparable to observed data due to the model 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 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 behaviour 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

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

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\mathit{\alpha }}{\mathit{\beta }}}=\sqrt{\mathrm{3}}.\end{array}$$

From this and Eqs. (A2b) and (A2b), it follows that the amplitude of compressional waves is $\approx \frac{\mathrm{1}}{{\sqrt{\mathrm{3}}}^{\mathrm{3}}}$ of the shear wave amplitude. If we say that the P wave train has a similar duration as the S wave train, then the energy contribution of the P waves with respect to the S waves becomes $(\frac{\mathrm{1}}{{\sqrt{\mathrm{3}}}^{\mathrm{3}}}{)}^{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{27}}$. The total energy of a signal is (Rudnicki and Freund, 1981)

$$\begin{array}{}\text{(B2)}& {\displaystyle }{E}_{\mathrm{total}}=E_{\mathrm{P}}+E_{\mathrm{S}},\end{array}$$

and can be estimated by

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\widehat{E}}_{\mathrm{total}}={\mathit{\alpha }}_{\mathrm{S}}\mathrm{aIV}{\mathrm{2}}_{\mathit{\alpha }}+{\mathit{\beta }}_S\mathrm{aIV}{\mathrm{2}}_{\mathit{\beta }},\end{array}$$

with the integrated squared velocity (IV2) 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 to

$$\begin{array}{}\text{(B4)}& {\displaystyle }& {\displaystyle }{\widehat{E}}_{\mathrm{total}}=\mathrm{aIV}{\mathrm{2}}_{\mathit{\beta }}\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{P}}}{\mathrm{27}}}+{\mathit{\beta }}_{\mathrm{S}}\right)\text{(B5)}& {\displaystyle }& {\displaystyle }=\mathrm{aIV}{\mathrm{2}}_{\mathit{\beta }}\left({\displaystyle \frac{{\mathit{\beta }}_{\mathrm{S}}\sqrt{\mathrm{3}}}{\mathrm{27}}}+{\mathit{\beta }}_{\mathrm{S}}\right)\text{(B6)}& {\displaystyle }& {\displaystyle }\approx \mathrm{aIV}{\mathrm{2}}_{\mathit{\beta }}\left({\displaystyle \frac{{\mathit{\beta }}_{\mathrm{S}}}{\mathrm{27}}}+{\mathit{\beta }}_{\mathrm{S}}\right).\end{array}$$

The last expression is differs only by 2.6 % from the exact term. While slightly underestimating the energy, this approximate 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.

At shorter distances, the intermediate-field term needs also to be taken into consideration. The amplitude of the intermediate term decays with r² (Eqs. A2a, A2a), while the far-field amplitude decays with r (Eqs. A2b, A2b). That is, the amplitude scales by distance and velocities and thus the IV2 are

$$\begin{array}{}\text{(B7)}& {\displaystyle }& {\displaystyle }\mathrm{IV}{\mathrm{2}}_{\mathit{\alpha }}={\mathit{\alpha }}^{-\mathrm{4}}{r}^{-\mathrm{2}}{\left(r^{-\mathrm{1}}+{\mathit{\alpha }}^{-\mathrm{1}}\right)}^{\mathrm{2}},\text{(B8)}& {\displaystyle }& {\displaystyle }\mathrm{aIV}{\mathrm{2}}_{\mathit{\beta }}={\mathit{\beta }}^{-\mathrm{4}}{r}^{-\mathrm{2}}{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}.\end{array}$$

Again by replacing all P wave terms by S wave terms, the total energy becomes

$$\begin{array}{ll}\text{(B9)}& {\displaystyle }& {\displaystyle }{\widehat{E}}_{\mathrm{total}}={\mathit{\alpha }}_S\mathrm{aIV}{\mathrm{2}}_{\mathit{\alpha }}+{\mathit{\beta }}_S\mathrm{aIV}{\mathrm{2}}_{\mathit{\beta }}{\displaystyle }& {\displaystyle }=a{r}^{-\mathrm{2}}\left({\mathit{\alpha }}_S{\mathit{\alpha }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\mathit{\alpha }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right.\\ \text{(B10)}& {\displaystyle }& {\displaystyle }\left.+{\mathit{\beta }}_S{\mathit{\beta }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right),{\displaystyle }& {\displaystyle }=a{r}^{-\mathrm{2}}\left({\sqrt{\mathrm{3}}}^{-\mathrm{3}}{\mathit{\beta }}_S{\mathit{\beta }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\sqrt{\mathrm{3}}}^{-\mathrm{1}}{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right.\\ \text{(B11)}& {\displaystyle }& {\displaystyle }\left.+{\mathit{\beta }}_S{\mathit{\beta }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right).\end{array}$$

With the assumption that α_S=β_S, Eq. (B4) becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{E}}_{\mathrm{total}}^{\mathrm{appr}}\approx a{r}^{-\mathrm{2}}\left({\mathit{\beta }}_S{\mathit{\alpha }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\mathit{\alpha }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right.\\ \text{(B12)}& {\displaystyle }& {\displaystyle }\left.+{\mathit{\beta }}_S{\mathit{\beta }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right),{\displaystyle }& {\displaystyle }=a{r}^{-\mathrm{2}}\left({\sqrt{\mathrm{3}}}^{-\mathrm{4}}{\mathit{\beta }}_S{\mathit{\beta }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\sqrt{\mathrm{3}}}^{-\mathrm{1}}{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right.\\ \text{(B13)}& {\displaystyle }& {\displaystyle }\left.+{\mathit{\beta }}_S{\mathit{\beta }}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}\right).\end{array}$$

The ratio between the approximation and the exact solution is

$$\begin{array}{}\text{(B14)}& {\displaystyle }{\displaystyle \frac{{\widehat{E}}_{\mathrm{total}}^{\mathrm{appr}}}{{\widehat{E}}_{\mathrm{total}}}}={\displaystyle \frac{{\sqrt{\mathrm{3}}}^{-\mathrm{4}}{\left(r^{-\mathrm{1}}+{\sqrt{\mathrm{3}}}^{-\mathrm{1}}{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}+{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}}{{\sqrt{\mathrm{3}}}^{-\mathrm{3}}(r^{-\mathrm{1}}+{\sqrt{\mathrm{3}}}^{-\mathrm{1}}{\mathit{\beta }}^{-\mathrm{1}}{)}^{\mathrm{2}}+{\left(r^{-\mathrm{1}}+{\mathit{\beta }}^{-\mathrm{1}}\right)}^{\mathrm{2}}}}.\end{array}$$

The two limits with respect to distance are

$$\begin{array}{}\text{(B15)}& {\displaystyle }& {\displaystyle }\underset{r\to \mathrm{0}}{lim}{\displaystyle \frac{{\widehat{E}}_{\mathrm{total}}^{\mathrm{appr}}}{{\widehat{E}}_{\mathrm{total}}}}={\displaystyle \frac{{\sqrt{\mathrm{3}}}^{-\mathrm{4}}+\mathrm{1}}{{\sqrt{\mathrm{3}}}^{-\mathrm{3}}+\mathrm{1}}}\text{(B16)}& {\displaystyle }& {\displaystyle }\approx \mathrm{0.932},\text{(B17)}& {\displaystyle }& {\displaystyle }\underset{r\to \mathrm{\infty }}{lim}{\displaystyle \frac{{\widehat{E}}_{\mathrm{total}}^{\mathrm{appr}}}{{\widehat{E}}_{\mathrm{total}}}}={\displaystyle \frac{{\sqrt{\mathrm{3}}}^{-\mathrm{5}}+\mathrm{1}}{{\sqrt{\mathrm{3}}}^{-\mathrm{6}}+\mathrm{1}}}\text{(B18)}& {\displaystyle }& {\displaystyle }\approx \mathrm{0.974}.\end{array}$$

The second limit is identical to the far-field case derived above. The two limits show that even in the range of the intermediate-field term, the energy estimate deviates little when assuming that all waves arrive with β_S at the recording site. A comparison of the approximate energy estimate and the exact estimate as a function of distance and velocity is shown in Fig. B1.