https://doi.org/10.20965/jdr.2021.p0658
Bhat, D. R. et al.

Paper:

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia
Deepak Raj Bhat∗1,† , Soichiro Osawa∗2 , Akihiko Wakai∗2 , Katsuo Sasahara∗3 ,Netra P. Bhandary∗4 , Fei Cai∗2 , Hirotaka Ochiai∗5 , and Norihiro Tanaka∗1∗1 Okuyama Boring Co., Ltd.
3-5-9 Higashi-nihonbashi, Chuo-ku, Tokyo 103-0004, Japan† Corresponding author, E-mail: deepakbhat@okuyama.co.jp∗2 Gunma University, Gunma, Japan∗3 Kochi University, Kochi, Japan∗4 Ehime University, Ehime, Japan∗5 Japan Forest Technology Association, Tokyo, Japan[Received December 7, 2020; accepted February 10, 2021]

In this study, novel finite element approaches are proposed for numerical analysis of stress-dependent landslide movement with groundwater fluctuation by rainfall. Two new constitutive parameters that are capable of directly controlling the relationship between theapparent factor of safety and sliding velocity are incorporated into a specific material formulation usedin finite element analysis for the first time. For thenumerical simulation of the measured time historyof the sliding displacement caused by the groundwater fluctuations, such required analytical parameters
can also approximately be determined by back analysis. The proposed models are applied to a landslide
field experiment on a natural slope caused by rainfall in real time in Futtsu City, Chiba Prefecture ofJapan to check its applicability. The predicted andmeasured time histories along the horizontal direction on the upper, middle, and lower slope are compared. In addition, the deformation pattern, shear
strain pattern, and possible failure mechanisms of thenatural slope of such a field experiment landslide are
discussed in detail based on the analysis results of thefinite element method (FEM)-based numerical simulation. Moreover, the creeping landslides and possible
landslide sites for further application of the proposedmodels are briefly discussed in the cases of Nepal andJapan as examples in Asia. It is believed that the proposed newly developed numerical models will help inunderstanding the secondary creep behavior of landslides triggered by extreme rainfall, and at the sametime, long-term management of such landslides will bemuch easier in monsoon Asia. Finally, it is expectedthat this study will be extended for simulation of thetertiary creep behavior of landslides induced by rainfall in the near future.

Keywords: monsoon, sliding velocity, factor of safety,

finite element method, constitutive model of soil

1. Introduction
1.1. Background
Soil creep is a time-dependent phenomenon in which
a deformation emerges under the application of constantstress. Theoretically, the different stages of creep are defined based on the ideal creep curve concept, in which
different stages of creep depend upon the strain rate. In
the first stage of creep, the strain rate decreases continuously with time; this is called the primary stage of creep.
After this, a constant strain rate emerges, which is called
the secondary stage of creep. Finally, the strain rate suddenly increases, which leads to a failure of the soil, whichis called the tertiary stage of creep. Terzaghi [1] was most
likely the first to consider the relationship between soil
creep and landslides. Ter-Stepanian [2] introduced thethreshold approach to explain soil creep in simple natural slopes by considering the zone of creep and its rate
as being dependent on the groundwater level. Eyring [3]introduced the theory of rate processes based on the ideathat the strain process consists of the mutual displacementof flow units by surmounting energy barriers and is quite
promising in prospect. Shimokawa [4] discusses the creep
deformation and creep strength of cohesive soils and theirrelationship to landslide based on experimental and theoretical studies using the rate process theory. Feda [5] hasalso interpreted the creeping behavior of soils using therate process theory.
Murayama and Shibata [6] and Christensen and Wu [7]have developed rheological models to represent the timedependent behavior of soil. Ter-Stepanian [8] proposed arheological model to understand the creeping behavior ofclayey soils during shear. In several publications about therheological creep behavior of soils to evaluate the range of
creep, Ter-Stepanian [8] explains that soil viscosity is pos-

658

Journal of Disaster Research Vol.16 No.4, 2021
© Fuji Technlogy Press Ltd.

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

itively related to a soil structure deformability coefficient.
Further investigation on Ter-Stepanian’s [2] creep model
has been done by Yen [9], and an apparent creep threshold line (yield strength line) is described, where the material is in a creep phase. In addition, Yen [9] concludes
that creep deformation is ensured when the shear stressdue to gravity exceeds the soil residual shear strength, andSuhaydu and Prior [10] and Iverson [11] also agree uponthis.
Sekiguchi and Ohta [12] and Zhou et al. [13] proposed
the time-dependent anisotropic model, which assumes
fixed anisotropic models. Desai et al. [14] proposed a
constitutive model for describing the creeping behaviorof natural slopes. In parallel, Pestana and Whittle [15],Wheeler et al. [16], and Dafalias et al. [17] proposed
time-independent anisotropic models accounting for both
initial and evolving anisotropy. A new anisotropic viscous model based on Wheeler et al.’s [16] model has
been proposed by Leoni et al. [18] to understand the ratedependent behavior of soft soils. Liingaard et al. [19] presented a comparative study on creep models. Calvelloet al. [20] presented an extensive study using a numerical model that focuses on active landslides controlled
by rainfall-induced pore pressure fluctuations with movements concentrated within a relatively narrow shear zone
above which the sliding mass moves essentially as a rigid
body. Leoni et al. [18] presented a new anisotropic modelfor the time-dependent behavior of soft soils using triaxial
compression test data. Huvaj and Maghsoudloo [21] useda simple Mohr-Coulomb constitutive model to model aslow-moving landslide in Italy using the computer program PLAXIS. Fernández-Merodo et al. [22] proposeda two-dimensional viscoplastic finite element model forslow-moving landslides and also applied this model tothe Portalet landslide of Spain as a case study. A newconcept of the residual-state creep test and a new testing setup and method for residual-state creep have beendeveloped by Bhat et al. [23, 24]. Based on that apparatus and method, the residual-state creep behavior of artificial clay and landslide clay has been studied [25, 26].
Bhat and Yatabe [27] proposed a new regression modelto understand the creeping behavior of clayey soils at theresidual-state of shear.
Creeping landslides are controlled by groundwaterfluctuations [2, 14, 28–31]. Therefore, groundwater fluctuations should be incorporated into the numerical simulation of such landslides. However, most of the previous numerical approaches [2, 5, 9, 18, 31–36] to soil creep
and associated problems have focused on laboratory creep
tests (i.e., consolidation/oedometer test and triaxial test),
which could not address the fluctuation of groundwaterlevel. Huvaj and Maghsoudloo [21] simulated the fluctuation of groundwater level in different phases to understand the creeping behavior of a slow-moving landslide,but the exact value of the deformation at any required
point (location) could not be captured perfectly. Recently,
a few researchers [22, 28, 29, 37] have proposed a twodimensional Elasto-viscoplastic constitutive model using
the finite element method based on the field instrumenJournal of Disaster Research Vol.16 No.4, 2021

tation and monitoring results, but they only considered asingle control constitutive parameter based on the trial anderror method, which could not control the displacement
rate of the landslide and also far from addressing the realistic field problem of the creeping behavior of landslides.
Therefore, if a new numerical approach, which can incorporate more than one control constitutive parameter for
directly controlling the displacement rate and total factor of safety of the landslide, can be developed, such approach may be useful for investigating the realistic fieldproblem of the creeping displacement behavior of landslides in the future.
The main objective of this study is to address the abovementioned problems to understand the creeping behavior of landslides using a newly developed finite elementbased numerical model and its implications for disasterprevention in monsoon Asia. The specific objectives include the following: i) to present the current creeping
landslide scenario in Nepal and Japan as examples of Asia
while searching for the possibility of application of the
proposed numerical models in the cases of Nepal andJapan in the near future; ii) to propose a simplified procedure for the determination of two new control constitutive parameters; iii) to estimate such new control constitutive parameters based on the relation between the apparent factor of safety and sliding velocity using the finite
element method (FEM); iv) to propose a new numericalapproach to evaluate the creeping behavior of landslides
owing to groundwater fluctuations due to rain water; v) tounderstand the failure mechanisms along the slip surfaceof creeping landslides triggered by extreme rainfall; vi) toapply the proposed models to understand the creeping behavior of a landslide field experiment on a natural slope
in Futtsu City of Chiba Prefecture, Japan, for the validation of the newly developed numerical models; and vii) topresent the simulation results of the deformation patternand shear strain pattern and discuss the possible failure
mechanisms of the slope based on variations in the local
factor of safety and mobilized shear strength along the slip
surface.

1.2. Creeping Landslides in Nepal and Japan:
Examples in Asia
Landslides are one of the major natural disasters in Asia
during the monsoon period. They have claimed thousands of human lives and extended economic losses ofbillions of dollars in individual property and infrastructuredamage each year, especially during the monsoon period.
The disaster records of Asian countries show that the frequency of fatal landslides was highest in China, followedby Indonesia, India, the Philippines, Japan, Pakistan,and Nepal between 1950–2009 [38]. Therefore, a heavy
amount of the national budget of Asian countries has beenspent on the prevention and mitigation of these landslides,
in which the major interest remains minimizing the creeping displacement of these landslides, thereby reducing thedamage risk to human settlements and structures over thelandslide mass or to the environment. However, the creep659

Bhat, D. R. et al.

Fig. 1. Location of selected creeping landslides for further application of the proposed model in Nepal.

behaviors of such landslides and the associated geological and geotechnical issues are not fully understood, so it
is difficult to predict the creep displacement behavior oflandslides, which may be useful for the design of countermeasure works for these landslides in the future. On
the other hand, although the creep displacement behaviorof landslides is very important, such creep displacement
behavior and its associated problems of landslides havenot been considered during the stability analysis of slopes
and landslide countermeasure works until now. In thispaper, creeping landslides in Nepal and Japan are briefly
discussed as examples in Asian countries.
Nepal is a highly mountainous country that is coveredby mountainous and hilly regions in more than 80% of
the area. Nepal has been formed by the collision of the
Indian and Eurasian plates, which occupy approximately
800 km along the central part of the Himalayan arc.
A huge variation of topography from less than 60 m in
the Southern Terai plain to 8848 m at Mount Everest in
the norther Higher Himalayas within an average lateral
distance of less than 200 km characterizes the steep topography of Nepal [39, 40]. This kind of topography isprone to landslides in Nepal [40]. Large-scale landslides
are very common in the Lesser and Higher zones of the
Nepal Himalayas due to the steep mountainous slopes anddynamic geological conditions [40]. Based on a databasestudy of landslides from 1978–2005, Petley et al. [41]state that here is a high level of variability in the occurrence of landslides from year to year in Nepal.
In Nepal, monsoons are a major source of rainfallin summer, and approximately 80% of the total annual precipitation takes place during the monsoon period(i.e., June to September) [39]. The annual average rainfall
is about the 1600 mm and approximately 2000–3000 mm
in the central part of Nepal. During the rainy season, daily
660

distribution of the precipitation is also uneven. For example, sometimes 10% of the total annual precipitation
occurs per day and 50% of the total annual precipitation per 10 days occurs during the monsoon period [39].
Therefore, such type of uneven distribution of precipitation may play an important role in triggering creeping
landslides in Nepal. According to Petley et al. [41], there
is a strong increasing trend in the number of fatalities occurring in Nepal at any point in the monsoon cycle.
In this context, the first and third authors visited thedifferent landslide areas along the major highway andBarkpark area of Gorkha in August 2017 for detailed
observation and sampling of the creeping landslide sites
of Nepal. During this visit, they observed a few old
activated creeping landslides and also identified a few
newly activated creeping landslides [42]. Fig. 1 shows
the location of the few representative creeping landslides, such as the Narayanghat-Mugling creeping landslide (Block I) and the Narayanghat-Mugling creeping
landslide (Block II) along the Narayanghat-Mugling roadsection, the Jalkini creeping landslide and Okharpuwa
creeping landslide along Trishuli highway, as well as theBalefi creeping landslide, Jure creeping landslide, Hindi
creeping landslide, and Tusare-Batish creeping landslide
along the Araniko highway, which seems applicable forthe further study using this newly developed FEM-based
numerical model in the near future. The soil samples fromthese landslide sites were collected, and further laboratory investigations have also been performed. Such field
and lab investigation results will be useful in estimating
the landslide simulation parameters for the application of
the proposed numerical model for numerical analysis of
those landslides in the near future. Figs. 2–4 show representative creeping landslides in Nepal. The scenario of
creeping landslides in Nepal has been discussed based on
Journal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

Fig. 2. View of Narayanghat-Mugling creeping landslide(Block II), Chitwan, Nepal.

Fig. 3. View of Jure creeping landslide, Sindupalchok, Nepal.

field observation and laboratory investigation of eight representative creeping landslides of Nepal [42].
During the field observation in Nepal, it was observedthat a large number of human settlements and rapid infrastructure developments are situated either on old or
newly activated creeping landslide masses or on or near
creeping landslide-prone areas, especially in the mountainous regions of Nepal. Landslides are a major problem
in Nepal, causing high levels of economic loss and substantial numbers of fatalities each year [41]. There is asimilar situation in many other Asian counties. Therefore, the damage to lives, properties, infrastructure, andthe environment are gradually increasing day by day due
to rainfall-induced landslides and the associated creeping
phenomenon in Asian countries during monsoon season.
Similarly, Japan is also a mountainous country,
where 75% of the total land area consists of mountainous terrain, and the remaining 25% is flat and low lying
with plateaus. Therefore, Japan has suffered from numerous landslide disasters since ancient times. On the other
hand, many people have lived in the mountainous region
since ancient times. Therefore, damage caused by landslide disasters, especially creeping landslides, has beenplaying a major role in influencing the quality of life [43].
Wakai et al. [44] and Bhat et al. [45] reported that many
creeping landslide sites accommodate human settlements,
agricultural fields, roads and highways, bridges and tunnels, and nature conservation sites in Japan. When the
Journal of Disaster Research Vol.16 No.4, 2021

Fig. 4. View of Okharpuwa creeping landslide, Nuwakot, Nepal.

displacement rate of such landslides is suddenly increased
and accelerated, it leads to significant mass failure, which
damages human lives, property, nature, and the environment. If some approaches to predicting the creep failureof landslides are possible, all damage from such landslidedisasters can be prevented [25, 46]. Therefore, a study ofthe creeping behavior of landslides and the associated geological and geotechnical engineering issues toward creep
failure prediction seems very important in Japan.
The geology of Japan is fragile due to the plate-tectonic
activities of the four major plates (i.e., North American
plate, Pacific plate, Eurasian plate, and Philippine plate)along the Japanese archipelago in the north-south direction, which is the one of the major causes of the activation and reactivation of creeping landslides in Japan [43,45]. Creeping landslides are one of the major natural disasters in mountainous regions of Japan and aremainly triggered by tectonic activities, occasional volcanic activities, heavy rainfall events, or a combinationof themes [44, 47, 48]. Continuous high-intensity rainfall
has been found to play a key role in provoking landslide
movement [49]. Sasaki et al. [50] revealed that soil creeps
usually occur during and after rain on the north side of the
Tsukuba Mountain in Japan. Ogita et al. [51] discuss theshapes and mechanisms of large-scale creeping landslides
in Japan. Bhat et al. [45] investigated the preexisting shear
surfaces of reactivated landslides in Nepal and Japan andhighlight the importance of the study of creep displacement behaviors of landslides in Nepal and Japan. Wakaiet al. [44] discuss the finite element-based numerical simulation of creeping landslides in Japan and its implication
for long-term monitoring and management of such landslides in the future.
The Japan Landslide Society [43] has introduced
27 representative landslides and slope failures all overJapan. The Japan Landslide Society [43] has also highlighted that several landslide mitigation measures are conducted in these landslide areas in order to stop or reduce
the landslide movement so that the resulting damage canbe minimized. However, the movement of such land661

Bhat, D. R. et al.

Fig. 7. View of Yui creeping landslide, Shizuoka, Japan [53].
Fig. 5.
View of Jinnosuke-dani creeping landslide,
Haku-san Mountain area, Japan [52].

seems almost impossible to complete landslide displacement predictions based on field monitoring results and, atthe same time, long-term management of such landslide
sites in the near future. Therefore, if a numerical toolor model that will be capable of better understanding
the creep displacement behavior of landslides inducedby rainfall and also predicting landslide displacement isneeded. Then, the significant damage to lives, properties,
and the environment from such landslide problems canbe mitigated. Therefore, we apply our newly developed
numerical models to a few landslide sites as case studies
of Nepal and Japan for landslide displacement predictionin the near future (Figs. 2–7). Later, we will apply suchmodels in cases of other Asian countries.
Fig. 6. View of Tozawa creeping landslide, Akita, Japan(Source: Yuri Forestry Office, Akita Prefectural Government, Japan).

2. Proposed Finite Element Analysis Method
2.1. Factor of Safety and Sliding Velocity

slides is not totally controlled or stopped even now. It
is also well understood that most of such landslides arereactivated due to the fluctuation of groundwater level,
mainly in rainy seasons. Therefore, it is believed that the
proposed models in this study would be useful for better
understanding stress-dependent landslide movement withgroundwater fluctuations for such representative landslides and slope failure, and at the same time, mitigation
measures for long-term management of such landslides
will be made easier in the near future. Figs. 5–7 show
a few representatives creeping landslides in Japan, where
the Okuyama Boring Co., Ltd. (the first and last authors)
is currently working. Recently, the authors (first, second,and last) have also been searching for the possibility of application of the proposed numerical models in these representative landslide sites as case studies in Japan.
For most Asian countries, regular monitoring(i.e., measuring the movement rate with groundwater fluctuations with respect to time) results of creeping
landslides is still lacking. On the other hand, most creeping landslides are in mountainous regions and remoteareas, where regular field monitoring of such landslides
is also difficult and time consuming [44, 45]. Hence, it
662

Based on an extensive literature survey, Vulliet and
Hutter [54] proposed an empirical relation between sliding velocity and stress state. Following the Vulliet and
Hutter [54], Sugawara [55] proposed a relation between
the overall factor safety (Fs ) and the sliding velocity (v) ofa slope on a power function, as shown in Eq. (1).
 v  −1nFs =, . . . . . . . . . . . . . (1)A
where, A (velocity dimension) and n (dimensionless) areconstants.
For the application of such a concept to the twodimensional finite element method, the sliding velocity (v) assumes a translational motion and is replaced by
the maximum displacement or strain rate (γ̇max ) of eachelement. Additionally, for the elasto-plastic finite element
method, the overall factor safety (Fs ) of a slope has been
replaced by the local factor of safety (Fs,local ) to reach theplastic state of each element. After redefining each constant to balance dimensions on both sides of the equation,
the following Eq. (2) is obtained.

γ̇max =α̇, . . . . . . . . . . . . . (2)nFs,localJournal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

where, α̇ (dimension of time−1 ) and n (dimensionless) areconstants.
c and ϕ  are the cohesion force and internal friction angle of soil, respectively. x − y represents the direct component along the x and y direction of the Cartesian coordinate system.
γmax = (εx − εy )2 + γxy 2 , . . . . . . . (3)σx + σy
sin ϕ  + c cos ϕ 
2Fs,local = .
σx − σy 2+ τxy 2
2

. . . .

(4)

The details will be described later in Section 2.3. The
magnitude of the maximum displacement rate under anarbitrary stress condition is defined based on Eq. (2). Similarly, the maximum displacement (γmax ) at each node andlocal factor of safety (Fs,local ) at each element are definedby following Eqs. (3) and (4), respectively.
On the other hand, Sugawara [55] proposed the empirical relation (a linear function) by direct plotting of bothfactors of safety and landslide displacement rate in a semilog graph. Following this concept of Sugawara [55], we
also propose Eq. (5), as shown below.

γ̇max = β̇ × 10m(1−Fs,local ) ,. . . . . . . . (5)where, β̇ (dimension of time−1 ) and m (dimensionless)are constants.
A series of tests have been performed to check the applicability of both proposed models based on Eqs. (2)and (5) in real field problems by Bhat et al. [47, 56]. They
have concluded that there are slight changes in the analysis results obtained from both proposed methods. Therefore, the other proposed model based on Eq. (5) is notconsidered in this study.

2.2. Outline of Nonlinear Analysis
In this study, the authors have simplified the elastoplastic finite element analysis code GA3D [57] into a twodimensional plane strain condition and a calculation algorithm that can reproduce the time-dependent behavior.
First, in order to generate an initial stress field under
the stable conditions with low groundwater level, the selfweight of the entire system is used as the acting external
force vector ( f 0 ) by balance calculating the overall stiffness matrix ([K]) of the system and the initial displacement vector (uu0 ), as shown below in Eq. (6).
[K]uu0 = f 0 . . . . . . . . . . . . . . . (6)Let u 0 be the initial value of the displacement vector (uu)at each node. u 0 has no physical meaning because it isan intrinsic displacement that occurs when the slope isformed. To obtain the stress increment vector (Δσ  ) for
each element, Eq. (7) is used based on the element displacement increment vector (Δuuelem ), where Δuuelem is obtained by only extracting the nodal components of each
Journal of Disaster Research Vol.16 No.4, 2021

element from the overall node displacement vector (uu0 ).
Δσ  = [D][B]Δuuelem .

. . . . . . . . . . (7)[B] and [D] are the B and D matrices of each element.
In the elastic-state, [D] is obtained from Young’s modulus
of elasticity (E) and Poisson’s ratio (ν ) based on Hooke’s
law (i.e., elastic modulus matrix [De ]), but in the plasticstate, the elastoplastic coefficient matrix [Dep ] is used, as
shown in Eq. (8).
[Dep ] =





∂g
∂f T e
[D ]
∂σ
∂σ
,
[De ] −

T 
 

∂ f ∂h∂g∂ f T e ∂g−[D ] 
+∂h ∂ε p∂σ∂σ∂σ
. . . . . . . . . . . . . . . . . . . (8)[De ]where f is the yield function, and g is the plastic potential. In this analysis method, strain softening is not considered; hence, the hardening parameter is assumed to be
zero (i.e., h = 0). Both the yield function ( f ) and the plastic potential (g) are defined based on the Mohr-Coulomb
type equations, as shown in Eqs. (9) and (10).
f = σ1 − σ3 − 2c cos ϕ  − (σ1 + σ3 ) sin ϕ  , . (9)g = σ1 − σ3 − (σ1 + σ3 ) sin ψ , . . . . . . (10)where σ1 and σ3 are the maximum and minimum principal stresses, and ψ is the dilatancy angle. Under nondrained (i.e., equal volume) conditions, ψ = 0.
Each stress component will be updated by cumulatively
adding the stress increment vector (Δσ  ) from Eq. (7) tothe stress vector (σ  ) (initially zero) of each element. At
this time, a correction is needed to fix and pull back theconstant J2 up to the yield surface to prevent the driftingof stress points from outside of the yield plane [58].
The equivalent nodal force vector ( f ) is back calculated(the sum of the total stress component vector (σ ), which
is obtained by adding the pore water pressure to the normal stress component of the effective stress (σ  ) of eachelement) from the stress distribution obtained as above.
f =



V

[B]T σ dV.

. . . . . . . . . . . (11)

If all elements are in the elastic state, then, f = f 0 .
Therefore, an iterative calculation due to material nonlinearity is not necessary. However, the plastic-state isnormally occurring in some elements.
In general, the external loading force ( f 0 ), which actsin the first stage of forming the initial stress field, is extremely large. Therefore, it will not converge if the generated stress increments are collectively subjected to theelasto-plastic processing using Eqs. (7) and (8). Hence,the shear strength reduction method (SSRM) [59] is used
for the convergence of it in the given initial condition.
In this method, a large shear strength is assumed in theelastic-state at first, and then, it is gradually decreased andreaches the plastic region. Finally, the imbalance force isbalanced by the number of iterative calculations.

663

Bhat, D. R. et al.

The process of nonlinear iterative calculation within
one calculation step of the shear strength reduction
method is described below. When any element reaches
the plastic-state, the residual force vector, Δrr (= f 0 − f )is obtained, and it is distributed to the entire system by themodified Newton-Raphson method.
[K]Δuu = Δrr . . . . . . . . . . . . . . . (12)As in Eq. (12) above, the initial stiffness matrix ([K]) is
used to obtain the residual displacement increment vector (Δuu) for the external force (Δrr ). It is update byadding Δuu instead of u . In addition, Δσ  is obtained
from Eq. (7) using the element displacement increment
vector (Δuuelem ), which is obtained by only extracting thenodal components of each element from Δuu. Moreover, a
new residual force (Δrr ) is obtained based on Eq. (11), andthe next Δuu is obtained from Eq. (12) in the next.
While repeating such calculations and updating σ
and u , the norm of the residual force (Δrr ) becomes sufficiently small, and the iterative calculations will converge.
After the convergence, the procedure proceeds to the next
step in the shear strength reduction method (slightly reducing the shear strength), and the convergence calculation may be repeated.
After the initial stress distribution is obtained, it will
move to the time history analysis considering the groundwater level fluctuations for each time step. It is notedthat the groundwater level has greatly varied at different locations of the slope. Here, the hydrostatic pressure is only considered, and the fluctuation of the pressure head due to permeability is ignored. Each effective
stress component of each element is obtained based on
the pore water pressure, which is obtained from the corresponding groundwater depth; hence, the plasticity analysis is performed using Eq. (9). Similarly, in the case ofelastic-state, [D] = [De ], and [D] = [Dep ] for the case of
the plastic-state, and the stress increment vector (Δσ  ) is
calculated from Eq. (7). After that, the balance of the entire system from the modified Newton-Raphson method
using Eqs. (11) and (12) is obtained the same as the initialstress analysis. After the calculation is converged, it will
move to the calculation for the next time step (i.e., input
a new groundwater level). By repeating this process, thesequential time history analysis is performed by inputting
the groundwater level fluctuations up to the predeterminedtime.

2.3. Constitutive Law Considering the New Viscosity Formulation
In the nonlinear calculation algorithm described in Section 2.2, the constitutive law of soil is necessary for considering Eq. (2). The interested material is a group ofthin slip/sliding elements that are formed along the sliding surface of the slope, and all these elements are simply
sheared uniformly (i.e., the local safety factor (Fs,local ) of
all elements is equivalent or equal. Let us assume that the
local safety factor (Fs,local ) of such elements is equal to664

the overall safety factor (Fs ) (i.e., Fs,local = Fs ).

γ̇max =α̇. . . . . . . . . . . . . . . (13)Fsn

However, the acting stress and strain rate (γ̇max ) along
the slip surface are not uniform in the actual case and often exhibit a progressive failure pattern. In other words,
the magnitudes of shear stress and shear strain assigned toeach element are different, but the sum of the shear resistance of all elements along the sliding surface is always inbalance with the driving force of the sliding mass. Since
the magnitude of the strain rate generated in each element
is obtained as a result of such a balance calculation of thewhole system, the strain rate is at least controlled by the
local safety factor (Fs,local ) of all elements. Hence, thestrain rate (γ̇max ) should be less than or equal to the rightside of Eq. (2). Therefore,γ̇max α̇. .
nFs,local. . . . . . . . . . . . (14)On the other hand, let us consider a case where the local factor of safety (Fs,local ) of all elements is equal forsatisfying Eq. (14); then, the strain rate (γ̇max ) of all elements is equal, and it is obvious that it should not contradict Eq. (13). Based on the consideration above, Eq. (14)should be considered in terms of the constitutive law forsatisfying the internal force balance conditions of the entire system and matching the strain rate (γ̇max ) at each element based on the magnitude estimated from the local
safety factor (Fs,local ).
In addition, when a strain rate (γ̇max ) on the right side ofEq. (14) is exceeded in a certain time step, it is necessary
to make a new action by incorporating an internal force(viscous resistance force), in which such forces retain theshear deformation of the elements, in order to keep it below the regulation value.
In a previous study, Sakuramoto et al. [60] proposed amodel based on the nonlinear viscoelastic visco-plastic finite element method using an iterative method to model ofthe polymer materials for artificial joints, in which elasticviscoelastic-viscoelastic elements are connected in series.
The concept of the proposed models in this study is similar to that of Sakuramoto et al. [60] in the sense of viscoelasticity components as a visco-plasticity. However,
the elements are not connected in the series of such threeelements. Moreover, we have also prepared the iterative
calculation algorithm for it.
Figure 8 shows a schematic diagram of the mechanical structure of the constitutive law, including the two
types of viscous resistance. The strain increment (Δε )based on the stress increment (Δσ ) given to the element
for each time step is divided into the elastic strain component (Δε e ) and other components Δε vp , where Δε vp occurs
as the visco-plastic strain component (Δε vp ) at first. However, the occurring magnitude of strain (strain rate) for
each time step is regulated by Eq. (14); hence, there is nochance of further excessive shear strain for the elements inwhich γ̇max will exceed on the right side of Eq. (14). Thus,

Journal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

Fig. 8. Conceptual figure of the proposed constitutive modelin this study.

the elastic response is restored, and the increase in strain
is suppressed. This is the viscoelasticity, which is shownin the upper part of Fig. 8. Including this, the equivalent nodal force ( f ) is calculated using Eq. (11), and theresidual force (Δrr ) is distributed throughout the entire system according to the modified Newton-Raphson method.
Here, the stiffness of the viscoelastic spring is the same as
the normal elastic spring, which contributes to the elastic
stress increment before the yield limit.
Since the restoring force of this viscoelastic spring is
temporarily generated as a resistance to sliding force forthe time step, the restoring force is released for the nexttime step when the calculation of γ̇max is initialized. Thismeans that when the occurring viscoelastic force is larger
in the immediately previous time step, the load factor of
the sliding force will also be greater in the next time step.
Thus, even if the variation of Fs,local greatly depends on
the location on the slip surface, the shear stress is gradually redistributed from the concentration point to the adjacent element because the strain rate may exceed the upper
limit of γ̇max . Hence, the sliding velocity tends to be uniform.
Here, it should be noted that the cumulative amount ofshear strain for each interested time step in Eq. (14) is
the sum of the elastic component and visco-plastic component, and it does not include the strain caused by viscoelasticity after that. To reduce the occurrence of the exceeded strain on the right side of Eq. (14), the viscoelastic
spring strain should be sufficiently smaller than the viscoelastic strain. Therefore, Young’s modulus (E) of thematerial should not be too small. These points are thesame as the discussions on how to give the rigidity when
joint elements are used for modeling of the discontinuous
surface characteristics [61].

3. Application of the Proposed Model
3.1. Landslide Field Experiment on a Natural Slope
The proposed model is applied to a landslide field
experiment on a natural slope triggered by a rainfall
Journal of Disaster Research Vol.16 No.4, 2021

Fig. 9. Location of the study area and simplified topographical map of field test site in Futtsu City, Chiba Prefecture, Japan, showing the location of piezometers, water level
gauges, extensometers, and tensiometers.

simulator in real time in Futtsu City, Chiba Prefecture,
Japan (Fig. 9). Fig. 9 shows the location of the study
area and the simplified topographic map of the field test
site, showing the location of the piezometers, water level
gauges, extensometers, and tensiometers along line 1 of
the X-section (A–A ) and line 2. The monitoring data
along line 2 are not considered in this study. The size ofthe field experimental landslide measured approximately
20 m × 10 m, which spread by the precipitation intensityof 140–300 mm/h over an area of approximately 100 m2
on the slope (Fig. 10). Figs. 11–12 show the photographsof the experimental sites before and after the collapse.
The Standard Penetration Test (SPT) (Japanese Geotechnical Society, 2017) was performed at five locationsof Nos. 6–8 to understand the tentative soil type of the testsite and to locate the possible slip surface of the rainfall
area (Fig. 9). At the location point of No. 6, a slip surfaceof 1.13 m from the ground level was observed. Similarly,
a sliding surface depth of 0.84 m from the ground surface
was recorded at No. 7, but there was not any sign of failure
at No. 8. Based on the SPT test results and several videos
captured during the experiment, the possible slip surfaceof the landslide was confirmed at approximately 1 m fromthe ground surface, as shown in Fig. 10.
The landslide field experiment using a rainfall simulator in real time was conducted on December 12, 2018.
665

Bhat, D. R. et al.

Fig. 10. Cross-sectional profile of the experimental slope
along the line 1, A–A .

Fig. 11. Photograph of the field experimental site before
collapse in Futtsu City, Chiba Prefecture, Japan, showing thelocation of piezometers, water level gauges, and extensometers.

However, the pre-experiment was also conducted one day
before by spreading a precipitation intensity of 140 mm/h
for approximately one hour to adjust the distribution of
precipitation on December 11, 2018. The real field experiment of the landslide was conducted at 10:53 a.m. on
December 12, 2018, by spreading a maximum precipitation intensity of 140 mm/h approximately 10 m × 10 m
(Fig. 10). The precipitation was increased to 188 mm/h
at 12:00 p.m. Similarly, the precipitation was increased
to 200 mm/h at 1:00 p.m. and 240 mm/h at 2:00 p.m. Finally, precipitation of 300 mm/h was fixed at 3:18 p.m.,where the slope starts to apparently move after rainfallof 4 hours and 25 minutes (Fig. 13). The major objective of this study is to observe the creeping behaviors
of a landslide induced by rainfall using finite elementbased proposed methods, so the highlighted block in
Fig. 10 (i.e., 17 m × 12.63 m) was only considered for
the FEM-based numerical simulation and analysis in thisstudy. Six water level gauges and extensometers were installed at various locations of the slope along line 1 (close
to A–A ) and line 2 for the measurement of the groundwater fluctuations and displacement rate with respect to time
4
(Fig. 9). However, the recorded data along line 2 (i.e., ,
5 and )6 are not considered in this study.
,
666

Fig. 12. A representative photograph of the field experimental site after collapse.

Fig. 13. Precipitation and groundwater fluctuation on theupper, middle, and lower water level gauges.

Figure 13 shows the precipitation and groundwaterfluctuations on the upper slope (i.e., upper water level
3
gauge or ),
middle slope (i.e., middle water level gauge2or ),and lower slope (i.e., lower water level gauge
1
The groundwater level in the upper water levelor ).
gauge began to rise at 10:57 a.m. after 3 minutes of precipitation of 140 mm/h, followed by the middle water level
gauge and the lower water level gauge, respectively. The
groundwater level in the upper water level gauge gradually increased by 0.27 m from 10:57 a.m. to 11:05 a.m.
and remained almost constant until 2:09 p.m. Then, it
suddenly decreased until 2:59 p.m. After spreading the
precipitation of 300 mm/h, the groundwater level againrose (Fig. 13). In the case of the lower water level gauge,
the groundwater level began to rise at 11:16 a.m. and gradually increased by 0.31 m until 11:26 a.m. Then, it also
followed a similar pattern of fluctuations of groundwaterin the upper water level gauge. However, the groundwaterlevel began to rise at 11:15 a.m. in the middle water level
gauge, and groundwater of 0.24 m gradually increased after about 10 minutes.
After this, the change in groundwater level was almost constant in all cases. Although the precipitation increased from 240 mm/h to 300 mm/h, the groundwaterlevel suddenly decreased after 2:09 p.m. and 2:16 p.m. in
the case of the upper water level gauge and lower water
Journal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

Fig. 14. Variations of the pore water pressures at the different depths on the upper, middle, and lower tensiometers.

Fig. 15. Variations of the monitored cumulative displacement on the upper, middle, and lower extensometers.

level gauge, respectively. The reason behind the decreasing groundwater level was not clearly understood. However, the soil structure at different locations on the slope
may change with changes in the slope. This change insoil structure may be due to the formation of soil blocks
and changes in the distribution of large pores in the soil.
During the movement of the landslide mass, there may
be heterogeneous distribution of the groundwater in thesoil surface, which may create a new discharge systemwith newly formed soil blocks. Moreover, the flow rateshould also change with the change in slope. However,the change in flow rate was not measured during the experiments in this study.
Six tensiometers were buried at depths of 20 cm, 50 cm,
80 cm, and 120 cm along line 1 and line 2 for automatic measurement of the pore water pressure during
the field experiment. Fig. 14 shows the relation between the pore water pressure and time at the different depths on the lower, middle, and upper slope alongline 1. After about 11:00 a.m., most of the tensiometers showed an increase in water pressure, and the fluctuations of the pore water pressure reading was almost constant until about 2:13 p.m., although the recorded valueof the pore water pressure at the different depths were
varying (Fig. 14). After that, most of the tensiometers
recorded decreasing trends. However, the tensiometers at
the depth of 120 cm in the middle followed the trend ofincreasing. After 2:00 p.m., 2:23 p.m., and 2:18 p.m.,the increase/decrease in pore water pressure at a depthof 80 cm on the lower, middle, and upper slope were notable (Fig. 14). The reason behind the change in pore water pressure at the different depths should be related to thechange in soil structures at the different depths.
Similarly, the corresponding extensometers recorded
cumulative displacement with time, as presented inFig. 15. From the analysis of the recorded displacement of the extensometer reading, it was understood thatthe slope began to move at 11:25 a.m. after gradually
increasing the groundwater level (Fig. 15). The lower
part (lower extensometer) of the landslide slope began tomove first, followed by the middle extensometer and upper extensometer. The lower was slightly moving withan average constant displacement rate of 0.16 mm/minfrom 11:25 a.m. to 2:10 p.m. After that, the average constant displacement rate was slightly increased up
to 2:47 p.m., and then, the displacement rate remained
almost 0 after that. Similarly, an average constant displacement rate of 0.37 mm/min was recorded in the middle up to 2:40 p.m., but change in displacement remained
almost 0 after that.
For recording and observing the overall phenomenon of
the movement of the landslide mass, multiple cameras andvideo cameras were installed at different locations along
the directions of line 1 and line 2 and fixed the time according to the Internet immediately before the beginningof the experiments to synchronize the time code on thevideos. Based on the detailed observation and analysis of
the recorded videos with time codes, it was found that thelower extensometer and middle extensometer were also
moving along the same direction with the movement of
the landslide mass. Therefore, those lower and middle
extensometers along line 1 may not have measured theexact displacement of the landslide mass. In other words,
the data recorded by the lower extensometer and middle
extensometer may have some errors. Therefore, the lower
extensometer and middle extensometer could not record
the data of tertiary creep (Fig. 15). Hence, the change indisplacement rate after 3:03 p.m. remains almost constantin the case of the lower extensometer and middle extensometer reading (Fig. 15).
In the case of the upper extensometer reading, a constant displacement rate of 0.1 mm/min was recorded in
the beginning at 2:09 p.m., where the displacement rateof the lower extensometer was also increased. However,
the average displacement rate gradually increased, and adisplacement rate of 2.1 mm/min was reached at 3:03 p.m.
After that, the displacement rate suddenly increased, leading to failure of the slope. Finally, the slope was totally collapsed after 3:18 p.m., with a displacement rateof 56.8 mm/min (i.e., last recorded data) (Fig. 15). Here,

Journal of Disaster Research Vol.16 No.4, 2021

667

Bhat, D. R. et al.

Fig. 16. Discretized finite element meshes with geologicalstructure and groundwater level in the slope based on field
investigations.

it seems that soil materials along the slip surface show
the tertiary creep behaviors after 3:03 p.m. In the other
word, the slip soil materials behave as the secondary stageof creep before 3:03 p.m. Therefore, the point (i.e., attime 3:03 p.m. along x-axis and equivalent cumulative displacement) is a critical point, from where the tertiary stageof creep began or, alternatively, the end point of the secondary stage of creep. As we have already discussed, theproposed models in this study are only capable of simulating the secondary creep behavior of soil materials along
the sliding surface. Therefore, we have only considered
up to the end point of the secondary stage of creep inthis study. Therefore, the groundwater fluctuation datafrom 9:00 a.m. to 3:03 p.m. were only considered for numerical analysis in this study.

3.2. Numerical Modelling
The proposed numerical model is applied to an analysisof the creeping behavior of a landslide field experiment
triggered by rainfall in real time of the natural slope inFuttsu City, Chiba Prefecture of Japan. Fig. 16 shows thediscretized finite element meshes with geological structure and groundwater level in the slope based on field
investigations. The three major representative materials(layers) are assumed from the SPT test results of suchlandslide for simplification of the model. The weakest
material (i.e., indicated by red in Fig. 16) was named“Sliding Surface.” The thickness of this layer is assumed
to be 1.0 cm. The soil materials are above the sliding
surface, which were named “Sliding Mass” or “Landslide
Mass.” Similarly, it is assumed that the strongest material is below the sliding surface, which was referred as
the “Bedrock.” Here, the landslide mass and bed rockwere modelled as solid elements and elastic bodies. However, the sliding surface was modelled with joint elements
and assumed to be an elasto-viscoplastic body. The mesh
668

adopted in the calculations consists of a rectangular element with eight nodes. The total number of 4800 nodes,
1531 elements, and 175 boundary conditions exist in thefinite element model of the landslide body (Fig. 16). Thebase of the model (i.e., bottom boundary) was assumed
to be fully impervious and fixed, and the lateral sides(right boundary and left boundary) were constrained byrollers. The hydraulic head of the maximum and minimum groundwater was imposed at the lateral boundaries
based on the field monitoring results of groundwater fluctuation. In this study, the results of groundwater fluctuations from 9:00 a.m. to 3:03 p.m. were only considered(i.e., up to the secondary creep). The lower, middle, andupper were the nearest locations of points, where the displacement of the landslide body was measured during thefield experiment (Fig. 16).
A parametric study was done to obtain the two new unknown control constitutive parameters (∝,˙ n). Initially, thetotal factor of safety (Fs ) of the landslide block is calculated using FEM. In FEM method, the shear strength
reduction method (SSRM) is used to estimate the totalfactor of safety (Fs ) of the field experimental landslidein Futtsu City. This method was originally proposed by
Zienkiewicz et al. [59] for elasto-plastic finite element
slope stability analysis. The details of the SSRM by FEM
have been discussed by Zienkiewicz et al. [59], Ugai andLeshchinsky [62], Cai et al. [63], Cheng et al. [64], andMurakami et al. [65]. Using the abovementioned method,
the total factor of safety owing to groundwater fluctuations caused by precipitation has been estimated for every
minute. From the detailed analysis of the slope stability
data using the FEM method, it was found that Fs were varied in a range of 0.98–1.020. After the calculation of Fs ,the relations between γ̇max and Fs was established basedon Eq. (14). After that, the general equations were obtained based on the well fitted curve between γ̇max and Fs .
Then, the unknown two new control constitutive parameters (∝,˙ n) were estimated by solving these general equations. On the other hand, the detailed field and laboratory
studied were not performed on this landslide site. Therefore, effective strength parameters (c , ϕ  ) of the sliding
surface were estimated based on the back analysis. The
other requires simulation input parameters for each soilmaterial (layers) are estimated based on different related
reports and guidelines of textbooks. The estimated valueof the unknown two new control constitutive parameters
and the summary of the material parameters for landslide
simulation are tabulated in Table 1.

3.3. Numerical Simulation Results and Discussion
Figure 17 shows the comparison of predicted time histories of displacement obtained from the numerical analysis and field monitoring data from the experiment. Here,the predicted time histories of displacement were measured at the same point in the upper, middle, and lower,
where the landslide displacement was also measured during the field experiment. Moreover, the horizontal component of displacement was considered for predicting thetime histories of displacement in the numerical model.
Journal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

Table 1. Material parameters for each material.
Materials (→)
Parameters (↓)
Young’s modulus, E [kN/m2 ]Poisson’s ratio, ν
Cohesion force, c [kN/m2 ]Internal friction angle, φ [◦ ]Dilatancy angle, ψ [◦ ]Unit weight, γ [kN/m3 ]α̇ [min−1 ]nLandslide Mass

Sliding Surface

Bed Rock

10000
0.45
7
20
0
18
0.00031
25

10000
0.45
7
20
0
18
0.0003125

100000
0.3
100
40
0
22
–
–

Fig. 17. Comparison of time histories of calculated sliding
displacement of the sliding mass and observed extension on
the upper, middle, and lower slope along line 1.
Fig. 18.
Calculated residual deformation (exaggerated
20 times) after the final calculation step.

Figure 18 shows the results of the deformation patternat the end (i.e., time = 3:03 p.m.) as a representative result. The red dotted line shows the result of a maximum
deformation pattern of each node at the end of the numerical simulation with comparison to without the deformation (i.e., initial condition). Similarly, Fig. 19 shows theresults of the shear strain pattern at the end step of thecalculation (i.e., time = 3:03 p.m.). From the comparative study of the field monitoring and numerical analysis
results, it is found that the horizontal displacement at thesame point of upper, middle, and lower was not perfectly
matched to each other (Fig. 17).
In this study, the detail field and lab test to estimate thevarious soil parameters of each soil layer for the simulation were not performed. Therefore, the landslide simulation parameters of sliding mass, sliding surface, andbedrock were assumed based on the back calculation as
well as trial methods. Moreover, material parameters of
sliding mass and slip surface for numerical simulation
were also assumed to be the same. The failure mechanisms along the slip surface due to the fluctuations ofgroundwater are not clearly known. From the videos thatwere taken by the NHK from the multiple cameras installed in different locations, we observed that the slip surface contained much water immediately after the occurJournal of Disaster Research Vol.16 No.4, 2021

rence of the landslides. It also pulled out and brooked theroots along the slip surface after the landslide (Fig. 12).
Therefore, the failure mechanisms along the slip surfaceof the natural slope should be known to obtain better
simulation results. However, we briefly discuss the possible failure mechanisms along the slip surface of sucha slope based on the simulation analysis results in theparagraphs below. On the other hand, the groundwaterlevel decreased in the upper and middle after 2:09 p.m.
and 2:16 p.m., respectively, whereas the movement rate
after that time was almost constant (i.e., secondary stageof creep) until 3:03 p.m. The interesting point is the major
causes of suddenly decreasing groundwater level despite
the constant movement (i.e., secondary creep behaviors)of the landslide mass. Therefore, the measured and predicted time histories were not good agreements with eachother (Fig. 17). Thus, the numerical simulation up to thetertiary stage of creep has been demanded to obtain better simulation results of such type of landslides. Hence,
we believe that such proposed models will be extended or
modified to simulate the tertiary creep behaviors of slope
or landslide triggered by rainfall in the near future. It is669

Bhat, D. R. et al.

Fig. 20. Change of calculated shear strength mobilizedalong slip surface in time.
Fig. 19. Calculated maximum shear strain distributions after
the final calculation step.

expected that such extended study will be useful for predicting the creep displacement behaviors of slope failure
or landslides triggered by rainfall, and at the same time,
countermeasure works for long-term prevention and mitigation from such type of landslide in monsoon Asia will
be easier in the near future.
Skempton [66, 67] has suggested that the drained residual strength is mobilized along pre-existing shear surfacesof a creeping landslide. Picarelli et al. [31] have suggested
that continuous translative movements of slopes may be
the result of cyclic pore pressure changes along the slip
surface, causing shear strength to be mobilized along theslide/slip surface. Vulliet and Hutter [54] have also agreed
with Picarelli et al. [31]. The main cause of landslide displacement is due to the increase of the mobilized shear
stress, in which weight of surcharge above the potential
failure may be increased because of rising of groundwater levels. This means that when the groundwater level israised, the pore water pressure is also increased along theslip surface. Hence, the effective stress will be decreased,
and the mobilized shear stress will be greater than the effective stress. Then, the soil mass begins to move slowly
at first, finally leading to failure. In other words, the normal stress is decreasing with the increase in groundwaterlevel, and the driving force becomes equal to or slightly
greater than the resisting force, which leads to failure or
moving of the landslide body. Picarelli et al. [31] havereported that increasing pore water pressure decreases theshear strength of the soil, which may lead to slope failure.
To check the details of the failure mechanisms, theshear resistance mobilized and variation of local factorof safety along the sliding surface have been also studied. Fig. 20 shows the change of calculated shear strength
mobilized along the slip surface in time. Similarly, thecalculated time histories of the local factor of safety ateach element along the slip surface are shown in Fig. 21.
The distance of the interested points was selected in an
equal horizontal distance (i.e., 2-m interval) from the toe
670

Fig. 21. Calculated time histories of the local factor ofsafety at each element along slip surface.

of slope (i.e., the coordinate of the toe was assumed zero).
The shear resistance mobilized along the slip surface
was slightly increased from the beginning (at 10:53 a.m.)
to 11:15 a.m. for all the observed distances, but after that,
it was suddenly increasing (Fig. 20). On the other hand,the change in the local factor of safety for all observed
distances along the slip surface was found almost constant up to 11:20 a.m. After that, the local factor of safety
along the slip surface was suddenly decreased, which maylead to movement of the landslide mass (Fig. 21). Basedon the comparative study of the numerical analysis results of the local factor of safety and shear strength mobilized along the slip surface in time, it was understoodthat the mobilized shear strength along the slip surface
near to toe (i.e., 2 m) was high, whereas the local factorof safety at the same location was also low, and vice versa(Figs. 20–21). From an analysis of the field monitoring
results, the movement of the slope was beginning near thetoe (i.e., lower) at first at 11:25 a.m. Therefore, such numerical analysis results were also good agreements withthe field monitoring results of displacement in the case of
initiating the movement of the slope. Thus, it can be concluded that the shear strength mobilized along the slip surface has a major role in triggering the landslide. Hence,
Journal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia

the failure mechanism of a creeping landslide is strongly
dependent on shear stress mobilized along the slip surface due to the change in groundwater level. Previous researchers (e.g., [2, 30, 31, 34]) have also agreed on this.

4. Concluding Remarks
This paper has primarily addressed the stressdependent landslide movements with groundwater fluctuations with the help of newly developed FEM-based numerical models and numerical analysis results and has
provided an insight for disaster prevention in monsoon
Asia. A two-dimensional elasto-viscoplastic constitutive
model has been used to simulate the secondary creeping
behavior of the sliding surface soil materials of rainfallinduced landslides. The major findings of this study include the following:(1) As an example, the current scenario of creeping landslides in Nepal and Japan have been briefly discussed, and possible landslide sites for further application of the proposed numerical models in the near
future have also been presented.
(2) The novel numerical approaches have been proposed
to understand the secondary creep behavior of thelandslide owing to groundwater fluctuations by rainfall.
(3) A simplified procedure for the determination of newcontrol constitutive parameters has been developed.
Two control constitutive parameters have been incorporated for the first time to assess the realistic fieldproblem of landslides triggered by rainfall.
(4) Such constitutive parameters can easily be estimatedbased on the relation between the total factor ofsafety (Fs ) and the sliding velocity (γ̇max ).
(5) In previous studies, understanding the creeping behavior of landslides seemed almost impossible without detailed field monitoring data and detailed field
and laboratory investigations. However, the proposed models are useful for better understanding thecreep displacement behaviors of landslides by directly changing the control constitutive parameters
and other simulation parameters for each soil material. Moreover, this study may be useful for predicting the displacement rate of landslides caused byextreme rainfall.
(6) The proposed models have been applied to a landslide field experiment on a natural slope in the Futtsu
City of Chiba Prefecture, Japan, to check the applicability of the novel models in real field problems ofslope failure and landslide.
(7) The measured (experimental) and predicted (model)
time histories along the horizontal direction havebeen compared, and slightly different results were
found with each. The possible causes of such differences have also been discussed.
(8) Based on the numerical analysis results of the field
experiment landslide, the details of the deformation
Journal of Disaster Research Vol.16 No.4, 2021

pattern, shear strain pattern, and possible causes ofslope failure mechanisms have been discussed.
The proposed numerical methods have been applied to
study the creep displacement behavior of few landslides
case studies until now. However, it will be continuously
applied to many other rainfall-induced landslides, withvariation of groundwater level to check for further applicability of such proposed models. Moreover, such models are only prepared for the numerical simulation of secondary stage of creep. Therefore, further modification or
extension of such models is necessary for the numericalanalysis of tertiary creep behavior of landslides in the near
future.
Acknowledgements
The authors would like to acknowledge Mr. Yosuke Koyama(NHK: Japan Broadcasting Corporation) for the overall planningand preparation of the field experiments. The authors would liketo thank the Japan Forest Technics Co. for conducting the preliminary survey of the experimental site. The authors wish to
thank the OSASI Technos Inc. and Daiki Rika Kogyo Co., Ltd.
for providing the necessary measuring equipment for the experiment and their cooperation during the field measurements. The
authors are also grateful to Futtsu City, a Federation of Five Construction Firms in Futtsu City, the Chiba Prefectural Police, andthe Department of Civil Engineering, Chiba Prefecture, for theirsupport and cooperation during the field experiments. The authors
would also like to thank the editor and two reviewers for their critical and helpful reviews to improve the quality of this paper. This
work has also been supported by JSPS-KAKENHI Grant Number 18H01674.

References:[1] K. Terzaghi, “Mechanism of Landslide,” S. Paige (Ed.), “Application of Geology to Engineering Practice: Berkey Volume,” pp. 83123, Geological Society of America, 1950.
[2] G. Ter-Stepanian, “On the Long Term Stability of Slopes,” Norwegian Geotechnical Institute (NGI) Publication No.52, pp. 1-14,
1963.
[3] H. Eyring, “Viscosity, Plasticity, and Diffusion as Example of Absolute Reaction Rates,” J. Chem. Phys., Vol.4, pp. 283-291, 1936.
[4] E. Shimokawa, “Creep Deformation of Cohesive Soils and its Relationship to Landslide,” Mem. Fac. Agr. Kagoshima. Univ., Vol.16,
pp. 129-156, 1980.
[5] J. Feda, “Interpretation of Creep of Soil by Rate Process Theory,”Géotechnique, Vol.39, No.4, pp. 667-677, 1989.
[6] S. Murayama and T. Shibata, “Rheological Properties of Clay,”
Proc. of the 5th Int. Conf. on Soil Mechanics and Foundation Engineering, Part 1, pp. 269-273, 1961.
[7] R. W. Christensen and T. H. Wu, “Analysis of Clay Deformationas Rate Process,” J. Soil Mech. Found. Div., Vol.90, No.6, pp. 125157, 1964.
[8] G. Ter-Stepanian, “Creep of a Clay During Shear and its Rheological Model,” Géotechnique, Vol.25, No.2, pp. 299-320, 1975.
[9] B. C. Yen, “Stability of Slopes Undergoing Creep Deformation,” J.
Soil. Mech. Found. Div., Vol.95, No.4, pp. 1075-1096, 1969.
[10] J. N. Suhaydu and D. B. Prior, “Explanation of Submarine Landslide Morphology by Stability Analysis and Rheological Models,”
Proc. of 10th Offshore Technol. Conf., pp. 1075-1082, 1978.
[11] R. M. Iverson, “A Constitutive Equation for Mass-Movement Behavior,” J. Geol., Vol.9, No.3, pp. 143-160, 1985.
[12] H. Sekiguchi and H. Ohta, “Induced Anisotropy and Time Dependency in Clays,” Proc. of the 9th Int. Conf. on Soil Mechanics andFoundation Engineering (ICSMFE), pp. 229-238, 1977.
[13] C. Zhou, J.-H. Yin, J.-G. Zhu, and C.-M. Cheng, “Elastic
Anisotropic Viscoplastic Modeling of the Strain-Rate-Dependent
Stress-Strain Behaviour of K0 -Consolidated Natural Marine Claysin Triaxial Shear Test,” Int. J. Geomech., Vol.5, No.3, pp. 218-232,
2006.

671

Bhat, D. R. et al.

[14] C. S. Desai, N. C. Samtani, and L. Vulliet, “Constitutive Modelingand Analysis of Creeping Slopes,” J. Geotech. Eng., Vol.121, No.1,
pp. 43-56, 1995.
[15] J. M. Pestana and A. J. Whittle, “Formulation of a Unified Constitutive Model for Clays and Sands,” Int. J. Numer. Anal. Meth.
Geomech., Vol.23, No.12, pp. 1215-1243, 1999.
[16] S. J. Wheeler, A. Näätänen, M. Karstunen, and M. Lojander, “AnAnisotropic Elastoplastic Model for Soft Clays,” Can. Geotech. J.,
Vol.40, No.2, pp. 403-418, 2003.
[17] Y. F. Dafalias, M. T. Manzari, and A. G. Papadimitriou,“SANICLAY: Simple Anisotropic Clay Plasticity Model,” Int. J.
Numer. Anal. Meth. Geomech., Vol.30, No.12, pp. 1231-1257,
2006.
[18] M. Leoni, M. Karstunen, and P. A. Vermeer, “Anisotropic Creep
Model for Soft Soils,” Géotechnique, Vol.58, No.3, pp. 215-226,
2008.
[19] M. Liingaard, A. Augustesen, and P. V. Lade, “Characterization ofModels for Time-Dependent Behavior of Soils,” Int. J. Geomech.,
Vol.4, No.3, pp. 157-177, 2004.
[20] M. Calvello, L. Cascini, and G. Sorbino, “A Numerical Procedure for Predicting Rainfall-Induced Movements of Active Landslides Along Pre-Existing Slip Surfaces,” Int. J. Numer. Anal. Meth.
Geomech., Vol.32, No.4, pp. 327-351, 2008.
[21] N. Huvaj and A. Maghsoudloo, “Finite Element Modeling of
Displacement Behavior of a Slow-Moving Landslide,” Proc. of
Geo-Congress 2013, pp. 670-679, 2013.
[22] J. A. Fernández-Merodo, J. C. Garcı́a-Davalillo, G. Herrera,
P. Mira, and M. Pastor, “2D Viscoplastic Finite Element Modellingof Slow Landslides: The Portalet Case Study (Spain),” Landslides,
Vol.11, No.1, pp. 29-42, 2014.
[23] D. R. Bhat, N. P. Bhandary, R. Yatabe, and R. C. Tiwari, “ResidualState Creep Test in Modified Torsional Ring Shear Machine: Methods and Implications,” Int. J. of GEOMATE, Vol.1, No.1, pp. 39-43,
2011.
[24] D. R. Bhat, N. P. Bhandary, R. Yatabe, and R. C. Tiwari, “A New
Concept of Residual-State Creep Test to Understand the Creeping
Behavior of Clayey Soils,” Proc. of GeoCongress 2012, pp. 683692, 2012.
[25] D. R. Bhat, N. P. Bhandary, and R. Yatabe, “Residual-State Creep
Behavior of Typical Clayey Soils,” Nat. Hazards, Vol.69, No.3,
pp. 2161-2178, 2013.
[26] D. R. Bhat, R. Yatabe, and N. P. Bhandary, “Creeping Displacement
Behavior of Clayey Soil in a New Creep Test Apparatus,” Proc. of
Geo-Shanghai 2014, pp. 275-285, 2014.
[27] D. R. Bhat and R. Yatabe, “A Regression Model for Residual State
Creep Failure,” Proc. of 18th Int. Conf. on Soil Mechanics and Geotechnical Engineering, pp. 707-711, 2016.
[28] E. Conte, A. Donato, and A. Troncone, “A Finite Element Approachfor the Analysis of Active Slow-Moving Landslides,” Landslides,
Vol.11, No.4, pp. 723-731, 2014.
[29] Y. Ishii, K. Ota, S. Kuraoka, and R. Tsunaki, “Evaluation of Slope
Stability by Finite Element Method Using Observed Displacementof Landslide,” Landslides, Vol.9, No.3, pp. 335-348, 2012.
[30] E. Eberhardt, L. Bonzanigo, and S. Loew, “Long-Term Investigation of a Deep-Seated Creeping Landslide in Crystalline Rock.
Part II. Mitigation Measures and Numerical Modelling of Deep
Drainage at Campo Vallemaggia,” Can. Geotech. J., Vol.44, No.10,
pp. 1181-1199, 2007.
[31] L. Picarelli, G. Urciuoli, and C. Russo, “Effect of Groundwater
Regime on the Behaviour of Clayey Slopes,” Can. Geotech. J.,
Vol.41, No.3, pp. 467-484, 2004.
[32] H. G. Brandes and D. D. Nakayama, “Creep, Strength and Other
Characteristics of Hawaiian Volcanic Soils,” Géotechnique, Vol.60,
No.4, pp. 235-245, 2010.
[33] J. D. Nelson and E. G. Thompson, “A Theory of Creep Failure
in Overconsolidated Clay,” J. Geotech. Eng. Div., Vol.103, No.11,
pp. 1281-1294, 1977.
[34] F. D. Patton, “Groundwater Pressure and Stability Analyses ofLandslides,” Proc. of the 4th Int. Symp. on Landslides, Vol.3,
pp. 43-60, 1984.
[35] L. K. Walker, “Undrained Creep in a Sensitive Clay,” Géotechnique,Vol.19, No.4, pp. 515-529, 1969.
[36] Z.-Y. Yin, C. S. Chang, M. Karstunen, and P.-Y. Hicher, “AnAnisotropic Elastic–Viscoplastic Model for Soft Clays,” Int. J.
Solids and Struct., Vol.47, No.5, pp. 665-677, 2010.
[37] W. Z. Savage and A. F. Chleborad, “A Model for Creeping Flowin Landslides,” Bull. Assoc. Eng. Geol., Vol.19, No.4, pp. 333-338,
1982.
[38] K. Forbes and J. Broadhead, “Forests and Landslides: The Role ofTrees and Forests in the Prevention of Landslides and Rehabilitationof Landslide-Affected Areas in Asia,” p. 3, Rap Publication, 2011.

672[39] R. K. Dahal, “Rainfall-Induced Landslides in Nepal,” Int. J. Eros.
Control Eng., Vol.5, No.1, pp. 1-8, 2012.
[40] S. Hasegawa, R. K. Dahal, M. Yamanaka, N. P. Bhandary,R. Yatabe, and H. Inagaki, “Causes of Large-Scale Landslides in theLesser Himalaya of Central Nepal,” Environ. Geol., Vol.57, No.6,
pp. 1423-1434, 2009.
[41] D. N. Petley, G. J. Hearn, A. Hart, N. J. Rosser, S. A. Dunning,
K. Oven, and W. A. Mitchell, “Trends in Landslide Occurrence in
Nepal,” Nat. Hazards, Vol.43, pp. 23-44, 2007.
[42] D. R. Bhat and A. Wakai, “Investigation of Creeping Landslides
Along the Major Highway of Nepal Triggered by the 2015 Gorkha,
Nepal Earthquake,” Proc. of the Research Committee on Strong
Linear Phenomena of Surface Ground and its Effects at the Great
Earthquake Symp., Ground Society Association of GeoKanto,
pp. 77-80, 2018.
[43] Japan Landslide Society and National Conference of Landslide
Control, “Landslides in Japan,” 6th Revision, 2002.
[44] A. Wakai, D. R. Bhat, K. Kotani, and S. Osawa, “Numerical Simulation of a Creeping Landslide Case in Japan,” B. Tiwari, K. Sassa,
P. T. Bobrowsky, K. Takara (Eds.), “Understanding and Reducing
Landslide Disaster Risk: Volume 4 – Testing, Modeling and Risk
Assessment,” pp. 273-280, Springer, 2021.
[45] D. R. Bhat, R. Yatabe, and N. P. Bhandary “Study of Preexisting
Shear Surfaces of Reactivated Landslides from a Strength Recovery
Perspective,” J. Asian Earth Sci., Vol.77, pp. 243-253, 2013.
[46] D. R. Bhat, A. Wakai, and K. Kotani, “A Finite Element Approachto Understand the Creeping Behaviour of Large-Scale Landslides,”Proc. of the 19th Int. Summer Symp., pp. 9-10, 2017.
[47] D. R. Bhat and A. Wakai, “Numerical Simulation of a Creeping
Landslide Induced by a Snow Melt Water,” Technical J., Vol.1,
No.1, pp. 71-78, 2019.
[48] H. Saito, O. Korup, T. Uchida, S. Hayashi, and T. Oguchi, “Rainfall Conditions, Typhoon Frequency, and Contemporary Landslide
Erosion in Japan,” Geology, Vol.42, No.11, pp. 999-1002, 2014.
[49] Y. Hong, H. Hiura, K. Shino, K. Sassa, A. Suemine, H. Fukuoka,
and G. Wang, “The Influence of Intense Rainfall on the Activity of Large-Scale Crystalline Schist Landslides in Shikoku Island,
Japan,” Landslides, Vol.2, No.2, pp. 97-105, 2005.
[50] Y. Sasaki, A. Fujii, and K. Asai, “Soil Creep Process and its Role in
Debris Slide Generation – Field Measurements on the North Side of
Tsukuba Mountain in Japan,” Eng. Geol., Vol.56, Nos.1-2, pp. 163183, 2000.
[51] S. Ogita, W. Sagara, D. Higaki, and Research Committee on Elucidating Mechanisms of Large-Scale Landslides, “Shapes and Mechanisms of Large-Scale Landslides in Japan: Forecasting Analysis from an Inventory (WCoE 2014–2017),” K. Sassa, M. Mikoš,and Y. Yin (Eds.), “Advancing Culture of Living with Landslides:
Vol.1 – ISDR-ICL Sendai Partnerships 2015–2025,” pp. 315-324,
Springer, 2017.
[52] H. Yagi and T. Inokuchi, “Aerial Watching of Landslides in Japan-26- Jin-Nosuke-Dani Landslide, Mt. Hakusan, Central Japan,” J.
Jpn. Landslide Soc., Vol.49, No.6, pp. 348-349, 2012.
[53] H. Yagi and T. Inokuchi, “Aerial Watching of Landslides in Japan-50- Epilogue of the Series “Aerial Watching of Landslides inJapan,”” J. Jpn. Landslide Soc., Vol.54, No.6, pp. 289-293, 2017.
[54] L. Vulliet and K. Hutter, “Viscous-Type Sliding Laws for Landslides,” Can. Geotech. J., Vol.25, No.3, pp. 467-477, 1988.
[55] N. Sugawara, “Variation of Surface Displacement Rate with Factor
of Safety for Creep-Type Landslides (As a Consequence from Reviewing Literatures),” Oyo Technical Report, No.23, pp. 1-18, 2003
(in Japanese with English Abstract).
[56] D. R. Bhat, A. Wakai, and K. Kotani, “A Comparative Study of Two
Newly Developed Numerical Models to Understand the Creeping
Behaviour of Landslides,” Proc. of the 20th Int. Summer Symp.,
pp. 101-102, 2018.
[57] The Japanese Geotechnical Society, “FEM Series 3 for Ground Engineers: Using the Elasto-Plastic Finite Element Method,” 2003 (inJapanese).
[58] T. Tanaka, “Deformation and Stability Analysis by Finite Element
Method,” The Japanese Society of Soil Mechanics and Foundation
Engineering, “Principle of Soil Mechanics,” 1st Revision, pp. 109154, 1992 (in Japanese).
[59] O. C. Zienkiewicz, C. Humpheson, and R. W. Lewis, “Associatedand Non-Associated Visco-Plasticity and Plasticity in Soil Mechanics,” Géotechnique, Vol.25, No.4, pp. 671-689, 1975.
[60] I. Sakuramoto, T. Tsuchida, K. Kuramoto, and S. Kawano, “A Studyon Nonlinear Viscoelastic Viscoplastic Finite Element Method Using Iterative Method,” Trans. Jpn. Soc. Mech. Eng. A, Vol.66,
No.651, pp. 1984-1989, 2000 (in Japanese with English Abstract).
[61] The Japanese Geotechnical Society, “Japanese Geotechnical Society Q & A,” Vol.2, pp. 276-278, 2009 (in Japanese).

Journal of Disaster Research Vol.16 No.4, 2021

Rigorous Analysis of Stress-Dependent Landslide Movementswith Groundwater Fluctuations Applicable
to Disaster Prevention in Monsoon Asia[62] K. Ugai and D. Leshchinsky, “Three-Dimensional Limit Equilibrium and Finite Element Analysis: A Comparison of Results,” Soils
Found., Vol.35, No.4, pp. 1-7, 1995.
[63] F. Cai, K. Ugai, A. Wakai, and Q. Li, “Effects of Horizontal Drains
on Slope Stability Under Rainfall by Three-Dimensional Finite Element Analysis,” Comput. Geotech., Vol.23, No.4, pp. 255-275,
1998.
[64] Y. M. Cheng, T. Lansivaara, and W. B. Wei, “Two-Dimensional
Slope Stability Analysis by Limit Equilibrium and Strength Reduction Methods,” Comput. Geotech., Vol.34, No.3, pp. 137-150, 2007.
[65] A. Murakami, A. Wakai, and K. Fujisawa, “Numerical Methods,”
Soils Found., Vol.50, No.6, pp. 877-892, 2010.
[66] A. W. Skempton, “Long-Term Stability of Clay Slopes,”Géotechnique, Vol.14, No.2, pp. 77-102, 1964.
[67] A. W. Skempton, “Residual Strength of Clays in Landslides, Folded
Strata and the Laboratory,” Géotechnique, Vol.35, No.1, pp. 3-18,
1985.

Name:
Katsuo Sasahara

Affiliation:
Professor, Science and Technology Unit, Natural Sciences Cluster,
Research and Education Faculty, Kochi University

Address:
2-5-1 Akebonocho, Kochi, Kochi 780-8520, Japan

Name:
Netra P. Bhandary

Name:
Deepak Raj Bhat

Affiliation:
Professor, Faculty of Collaborative Regional Innovation, Ehime University

Address:
Affiliation:

3 Bunkyo-cho, Matsuyama, Ehime 790-8577, Japan

Chief (Research) Engineer, Okuyama Boring
Co., Ltd.

Address:
5F Goto Building, 3-5-9 Higashi-Nihonbashi, Chuo-ku, Tokyo 103-0004,
Japan

Selected Publications:

Name:
Fei Cai• “Residual-state creep behavior of typical clayey soils,” Nat. Hazards,
Vol.69, No.3, pp. 2161-2178, 2013.
• “Study of preexisting shear surfaces of reactivated landslides from astrength recovery perspective,” J. Asian Earth Sci., Vol.77, pp. 243-253,
2013.

Affiliation:

Name:

Name:

Soichiro Osawa

Hirotaka Ochiai

Affiliation:

Affiliation:

Graduate Student, Geotechnical Engineering Laboratory, Department ofEnvironmental Engineering Science, Gunma University

Technical Adviser, Japan Forest Technology Association (JAFTA)

Address:

7 Rokubancho, Chiyoda-ku, Tokyo 102-0085, Japan

Associate Professor, Gunma University

Address:
1-5-1 Tenjincho, Kiryu, Gunma 376-8515, Japan

Address:

1-5-1 Tenjincho, Kiryu, Gunma 376-8515, Japan

Name:
Name:

Norihiro Tanaka

Akihiko Wakai

Affiliation:
Affiliation:

Tokyo Office Chief, Chief (Research) Engineer, Okuyama Boring Co., Ltd.

Professor, Gunma University

Address:

Address:

5F Goto Building, 3-5-9 Higashi-Nihonbashi, Chuo-ku, Tokyo 103-0004,
Japan

1-5-1 Tenjincho, Kiryu, Gunma 376-8515, Japan

Journal of Disaster Research Vol.16 No.4, 2021

Powered by TCPDF (www.tcpdf.org)

673