1-D inversion analysis of a shallow landslide triggered by the 2018 Eastern Iburi earthquake in Hokkaido, Japan

Destructive landslides were triggered by the 6.7 Mw Eastern Iburi earthquake that struck southern Hokkaido, Japan, on 6 September 2018. In this study, we carried out 1-D inversion analysis of one of the shallow landslides near the epicenter using a Bing debris-flow model. At this site, the slope failure comprised cover soil with an initial down-slope length of ~ 80 m and a thickness of ~ 7 m on a slope with < 20° dip. The landslide moved southeastward with a run-out distance of ~ 100 m. Inversion analysis of the post-failure deposit geometry was conducted with the Markov Chain Monte Carlo method (MCMC) to optimize the Bingham rheological parameters of the debris. The analysis reproduced several features of the deposit geometry with a yield stress of ~ 1500 Pa and dynamic viscosity of 800–3000 Pa s. The results suggest that the shallow landslide can be approximated by the flow of a viscoplastic fluid with high-mobility debris and a maximum frontal velocity of 6–9 m/s, with a flow duration of 2–4 min.


Introduction
The Eastern Iburi earthquake (M w = 6.7) occurred in southern Hokkaido, Japan, on 6 September 2018 (Fig. 1a). The resulting strong ground motion had a maximum intensity of 7 on the Japan Meteorology Agency (JMA) intensity scale and caused shallow landslides as well as a few large-scale deep-seated landslides near the epicenter (Yamagishi and Yamazaki 2018;Osanai et al. 2019). Most of the shallow landslides show long runouts of debris on relatively gentle slopes (< 30°; Kasai and Yamada 2019;Osanai et al. 2019) and are classified as earthflows or earth slides according to the scheme of Varnes (1978). The widespread hills in southern Hokkaido are covered by volcanic soils derived from material produced by several nearby volcanos. The soils are thought to have been weakened by heavy rainfall from a typhoon that passed through the area the day before the earthquake (Osanai et al. 2019).
Field observations of the shallow landslides revealed that the volcanic soil close to the slip zone is clay-rich and contains the clay mineral halloysite (Chigira et al. 2019;Kameda et al. 2019). Halloysite is a typical alteration product of volcanic glass (Wada 1977) and is commonly observed on the slip surfaces of other landslide sites in volcanic areas of Japan (Chigira and Yokoyama 2005;Chigira et al. 2012). The field survey also found evidence of liquefaction of volcanic soil (Kameda et al. 2019). Geotechnical experiments using a triaxial apparatus have demonstrated that the volcanic soils around the slip zone have a low shear resistance and are more susceptible to liquefaction than are other layers (Li and Wang 2020). The earthquake possibly caused liquefaction and fluidization of the water-saturated soils, resulting in multiple shallow landslides in this area.
In addition to the works mentioned above, other studies have analyzed the landslides by field surveys, laboratory experiments, and image analysis (e.g., Kawamura Okamoto Earth, Planets andSpace (2021) 73:116 et al. 2019;Wang et al. 2019;Zhang, et al. 2019;Ito et al. 2020), and they reported on the geological history and triggering mechanisms. To complement these studies, we conducted numerical modeling of the landslide and discuss its post-failure behavior. We selected one of the shallow landslides for our study, whose geometry was investigated in a previous field survey (Li and Wang 2020). We employ a Bing debris-flow model  for inversion analysis of the landslide geometry and discuss the rheological properties of the collapse debris to interpret its post-failure flow behavior.

Landslide
The studied landslide is located in the Asahi district of the Atsuma town in Hokkaido (42° 44′ 10.6″ N, 141° 52′ 39.8″ E; Fig. 1a). The volcanic soil slid southeastward on a slope that dips 19° at its highest point and decreases in dip down-slope. The run-out distance of the landslide is ~ 100 m (Fig. 1b, c;Li and Wang 2020). This landslide was selected for numerical analysis because of its simple geometry. Figure 1c shows a cross-section through the slope. Due to poor field exposure, we assume the depth of the slip plane on the slope based on the thickness of the post-failure debris, with the depth increasing down-slope. Despite making this assumption, our model successfully reproduces the whole geometry of the deposit including the slope.

Forward modeling
We used a Bing debris-flow model  to conduct a numerical analysis of the landslide. This model assumes two layers of debris (i.e., an upper plug-layer and a lower shear-layer) that flow laminarly with rheological (b) (a) (c) Fig. 1 a Location of the epicenter of the Eastern Iburi earthquake (red cross) and the studied landslide in the Asahi district (42° 44′ 10.6″ N, 141° 52′ 39.8″ E; black arrow). b Aerial view of the landslide (taken by the Geospatial Information Authority of Japan 2018). c Cross-section of the landslide based on the field survey by Li and Wang (2020). Soils before and after the failure are marked by light blue and green, and yellow and green, respectively Page 3 of 9 Kameda and Okamoto Earth, Planets and Space (2021) 73:116 soil properties described by a viscoplastic Herschel-Bulkley model: where τ is the shear stress, τ y is the yield stress, and γ is the strain rate. A reference strain rate, denoted by γ r , is defined as follows: where µ is the dynamic viscosity. The volcanic soil of the studied landslide is assumed to behave as a Bingham fluid, which is a limiting case for the Herschel-Bulkley model with n = 1. The Bingham model is the most commonly used viscoplastic model that describes the behavior of debris or mud flows (Jiang and LeBlond 1993;Wan and Wang 1994;Julien 1995;Naruse 2016). The integral equations of the debris flow on a Lagrangian framework are described by the following three conservation equations for mass and momentum (Jiang and LeBlond 1993;Huang and Garcia 1998;Pratson et al. 2001;Imran et al. 2001;Naruse 2016): where U P is the plug-layer velocity, D P , D S , and D are the plug-layer, shear-layer, and total thickness of the debris, respectively, S is the slope angle, t is time, x is the position along the slope (downward is positive), g is gravitational acceleration, ρ is the density of air, and ρ m is the density of the debris material. According to Li and Wang (2020), the clay-rich volcanic soil near the base of the landslide (known as Ta-d pumice MS) has a specific gravity of 2.631 ( ρ m = 1177 kg/m 3 ) and a water content of 206% at a degree of saturation (S r ) of 92.8%. The upper volcanic soil (Ta-d pumice with volcanic ash) has a specific gravity of 2.553 ( ρ m = 1122 kg/m 3 ) and a water content of 121% at a degree of saturation (S r ) of 76.2%. If these soils were fully saturated at failure (i.e., S r = 100%), ρ m is estimated to be 1240 kg/m 3 for the upper layer and 1310 kg/m 3 for the lower layer. Therefore, we assume ρ m = 1300 kg/m 3 in our model. We developed a MATLAB code to numerically solve Eqs. (3)-(5). The equations were assembled in a staggered lattice using a finite difference method with a deformable moving grid system. Following Imran et al. (2001), a numerical viscosity of 0.001 was used to prevent the solution from becoming unstable. More detailed information on the solution procedure can be found in Imran et al. (2001). The calculation was stopped when the frontal velocity of the debris flow decreased to 1 cm/s.

Inversion analysis
Based on the geometry of the landslide deposit, inversion analysis was conducted to optimize the parameters τ y and γ r in Eq. (2). Since our analysis is 1-D, lateral flow is ignored. However, the volume of debris before and after failure is not equal along the constructed cross-section ( Fig. 1c), with the latter being larger than the former by 17%. This implies that soils from outside the survey line contributed to the final geometry of the deposit. For this reason, we use a volume of the post-failure deposit that is artificially reduced (shortened vertically), so that it is identical to the pre-failure volume. Due to these limitations, we aim at reproducing key features of the debris flow such as run-out distance and position of the hump of the deposit.
In the inversion analysis, we defined an objective function as a residual sum of squares between the observed (corrected) and modeled deposit thickness where n is the number of grid points and D oi is the observed thickness of the deposit at i grid point. Based on the result of the forward model, the thickness of the modeled deposit at each grid point ( D i ) was calculated using the one-dimensional data interpolation method. The analysis was done over the length of the deposit, from x = 76.5 to 184.68 m, which consisted of 58 grid points. The source area was excluded from the inversion analysis Page 4 of 9 Kameda and Okamoto Earth, Planets and Space (2021) 73:116 (x = 0 to 76.5 m), because we assumed the thickness of the deposit on this area as described above.
Since the F function is non-linearly dependent on the fitting parameters, a Markov Chain Monte Carlo (MCMC) method was applied to optimize the fitting parameters over a range of values. Details of this method are outlined by Metropolis et al. (1953). For an MCMC simulation, a candidate value for one of the unknown parameters is tested by adding a random number (n r = [− 1,1]) to the old candidate, namely log 10 (P i,can ) = log 10 (P i ) + n r , where P i and P i,can are the old and current candidates, respectively. Therefore, the candidate parameter vector for one local step can be described as ϕ can = P 1,can , P 2 , . . . , P N . F (ϕ can ) is then calculated from these parameters and F is evaluated The candidate vector ϕ can is accepted if the probability is min (1, exp (−�F β N )) , where β N is the "inverse temperature", which controls the acceptance or rejection of the candidate values. At low β N , the candidate value is updated when − F > 0. At high β N , this is not the case. The appropriate β N is problem-dependent, and in this work, β N was set to 0.5-1.0 by trial-and-error (β N was fixed during one series of calculations). After one local update, the next unknown parameter is tested in a similar way. One Monte Carlo Step (MCS) is defined as N trials (i.e., the number of unknown parameters examined) of the local update, and E sum is obtained by the arithmetic sum of F of each MCS. Schematic flowchart of the MCMC method for our inversion analysis is shown in Fig. 2.

Results
The results of MCMC inversion are shown on a yield stress-viscosity ( τ y − µ ) plot (Fig. 3). The parameters in the simulation are treated in logarithmic form (log 10 (τ y ) and log 10 (γ r )), but to clearly illustrate the results, E sum values are plotted against τ y and µ . We used several initial value conditions indicated by the red crosses in Fig. 3. However, after several MCSs, the resulting viscosity and yield stress values converged to two distinct domains: a high-viscosity domain and a low-viscosity domain. E sum values for these domains are similar, ranging from ~ 17 to 20, which corresponds to an F value of ~ 9 to 10 and a trial number N = 2. Figure 4 shows how τ y and µ are updated with the progress of the MCMC simulation after a converged state is achieved. In a converged state, the accepted yield stress values vary steadily between 1450 and 1800 Pa, while the viscosity values fluctuate between two distinct ranges of values: ~ 1000 Pa s and 1500-3000 Pa s (Fig. 4). Of the > 8000 MCSs, the ~ 6200 MCSs yield viscosity values in Based on these results, the probability distributions of log 10 (τ y ) and log 10 (γ r ) values in the two domains are individually estimated (Fig. 5). The mean and standard deviation of log 10 (τ y ) and log 10 (γ r ) in the high-viscosity domain are 3.192 (1σ = 0.0169) and − 0.151 (1σ = 0.0772), respectively, corresponding to a yield stress of 1554 Pa and viscosity of 2200 Pa s (Fig. 5a, b). The mean and standard deviation of log 10 (τ y ) in the low-viscosity domain are comparable to the values in the high-viscosity domain at 3.201 (1σ = 0.0166), which gives a similar yield stress of 1588 Pa (Fig. 5c). However, log 10 (γ r ) in the low-viscosity domain has a larger mean value of 0.225 and lower standard deviation (1σ = 0.0347) than that in the high-viscosity domain and, consequently, a lower viscosity of 945 Pa s (Fig. 5d). Figure 6 shows the results of the forward modeling based on the two optimized parameter sets obtained from the MCMC simulation. The deposit thickness and frontal velocity are plotted against the distance from the top of the source area for both the high-viscosity and low-viscosity parameter sets. In the case of the high-viscosity parameter set (Fig. 6a), the debris flow accelerates to 6.5 m/s before quickly decelerating to 1 m/s. It then continues to flow with a velocity of < 1 m/s until it stops moving, 252 s after initial failure, at a distance of 175 m (F = 9.185). In the lowviscosity case (Fig. 6b), the debris flow accelerates to a frontal velocity of 9 m/s and keeps this high-velocity state for a certain period. The debris flow then quickly decelerates and stops flowing after 115 s (F = 9.033). Although the low-viscosity debris flow is completed in a shorter timeframe than the high-viscosity debris flow, the final deposit profiles are similar in both cases.

Discussion and conclusions
Our 1-D inversion analysis yields rheological parameters that closely reproduce the observed deposit profile (Fig. 6). The parameters are separated into two domains that both lie outside the field of sand-rich soils as Yield stress (Pa) Viscosity (Pa·s) Fig. 4 Behavior of yield stress and viscosity over 8500 MCSs. Left: fitting parameter log 10 (τ y ) recalculated as τ y . Right: fitting parameter log 10 (γ r ) recalculated as µ Page 6 of 9 Kameda and Okamoto Earth, Planets and Space (2021) 73:116 classified by Jeong et al. (2010) (Fig. 3). The boundaries between the fields of sandy and silty soils, and between silty and clay-rich soils occur at lower viscosities than those obtained for the studied landslide (Locat 1997;Jeong et al. 2010). Although this suggests that the rheology of the volcanic soil is more comparable to that of a gravel-rich soil, care should be taken in making this comparison. The basal volcanic soil of a similar landslide (Ta-d) triggered by the same earthquake contains the clay mineral halloysite (Kameda et al. 2019), and suspensions of halloysite generally exhibit shear-ratedependent resistance as assumed in this study (Theng and Wells 1995), but the rheological properties of halloysite-bearing soils are poorly defined. Moreover, the parameters obtained in this study likely represent the bulk rheological properties of the collapsed soil rather than of the sorted fine-grained materials that are generally used in laboratory experiments (Locat 1997;Jeong et al. 2010). On the other hand, several studies have reported the rheological properties of soils involved in flow-like landslides based on numerical modeling (Malet et al. 2004;Carrière et al. 2018). Using a Bing model, Carrière et al. (2018)  Probability distributions of log 10 (τ y ) and log 10 (γ r ) obtained by the MCMC simulation. Mean values (dashed line) and standard deviations (grey shading) are shown for the high-viscosity (a, b) and low-viscosity (c, d) domains Kameda and Okamoto Earth, Planets and Space (2021) 73:116 not directly comparable to our study, since they used the Herschel-Bulkley model with n = ~ 0.3 and the debris is composed mainly of quartz and mica (Carrière et al. 2018), which are absent in the volcanic soil of the Atsuma region (Kameda et al. 2019). Nevertheless, the obtained yield stress is in the same order of magnitude as our results. Further experiments are necessary to test the validity of our results and assess which viscosity value gives the best fit to the flow of the volcanic soil.
Although the observed deposit profile is closely reproduced by our numerical model, there are some minor differences. For instance, the observed debris profile shows ~ 1 m of uplift at a distance of 130 m (Figs. 1 and 6). In our numerical model, this is observed as a transient phenomenon that gradually disappears with the progression of the debris flow. Figure 7 shows time snapshots of the debris profile to show temporal changes in the morphology of the debris flow in the high-viscosity case. During the first 2 min (Fig. 7a, b), uplift of the debris occurs at a distance of 130 m, which resembles the uplift feature in the observed debris profile (Fig. 6). During the period after 2 min, the uplift at 130 m ceases to occur (Fig. 7c, d). At 130 m from the top of the source area, the slide slope is steeper (14° dip) over a distance of ~ 5 m (Fig. 1c). This area corresponds to the sloping ground between the road and the field, which might explain the transient uplift (i.e., a stack of debris behind the slope). The conservation of this feature in the debris-flow deposit may be a result of the local variation in rheological soil properties and/or water saturation state. Such variation could arise from a vertical layering of different volcaniclastic strata, which is not incorporated in our model.
Another discrepancy between our modeled profile and the observed profile is the V-shaped depression directly behind the uplift at 130-145 m (Fig. 6). One explanation for this depression is a lateral crack in the debris deposit. To discuss such features, it would be necessary to use numerical modeling techniques that are capable of reproducing flow processes in discontinuous media (e.g., the discrete element method). However, another interpretation is a local erosion of a small drainage channel located south of the road. In fact, it seems there is a small impoundment on the right and it is possible to see a little stream on the left (Fig. 1b). If the latter interpretation is correct, then our model seems to reproduce the observed deposit profile quite well, including the microtopographic features.
In recent years, 3D simulations of landslides considering various rheological models have been conducted (Kelfoun and Druitt 2005;Pirulli and Mangeney 2008), but in most of these studies, the model parameters are calibrated by trial-and-error. In addition, seismic signals associated with landslides rather than deposit profiles are often used in inversions analysis of landslides Moretti et al. 2020). Although there are still some uncertainties in our analysis, the model successfully fits the observed deposit profile with reasonable yield stress and viscosity values, and therefore, the MCMC-based inversion analysis presented here is useful as a new methodology for reproducing landslide processes and inferring material properties. The results of our analysis suggest that the shallow landslide can be approximated by the flow of a viscoplastic fluid over a period of several minutes. The morphological features, geological background, and rheological properties of the landslide are similar to those of other shallow landslides triggered by the 2018 earthquake (Osanai et al. 2019;Li and Wang 2020). The many slope failures in the region might therefore have collapsed in a similar manner to that outlined in the present study.