10th Asian Regional Conference of IAEG (2015)

Study for the Elucidation of the Sedimentation Structure Using the Numerical Analysis
Yusaku ISOBE(1), Tomoaki HACHIMURA(2), Hiroyuki OHNO(3), Kenji MATSUMOTO(2),Tetsuya MIYAHARA(2), and Shuichi NAGAOKA (2)(1) Geoscience Research Laboratory, Co., Ltd., Japan
E-mail: isobe@geolab.jp(2) Japan Environmental Sanitation Center, Japan, (3) Kankyo Chishitsu, Co., Ltd., Japan
Abstract
Sedimentation in the liquid (the water or the seas) had been thought based on the
Stokes' theorem, but the study of authors presented that it did not follow the Stokes'
theorem when a sedimentation particle has the large Reynolds number. Tsunami sedimentand slump deposit (landslide) are considered as the environment causing the large Reynolds
number, and understanding these sedimentation forms will leads to the valuable knowledge
for future disaster. Therefore, the authors tried to express the sedimentation phenomenon in
the liquid by the numerical analysis technique. It is confirmed that the numerical analysis isable to express the phenomenon that an enormous number of fine particles lose their energy
and sank in the liquid at uniform velocity.
Keywords: MPS method, Sedimentation, Numerical analysis
1.

Introduction

consideration of influence of the movement among
particles by analytical technique. In this approach, adifferential equation is discretized using a particle
interaction model shown in Fig.1 for a differential
operator.

A turbidite is the geologic deposit of a turbidity
current, which is a type of sediment gravity flow
responsible for distributing vast amounts of clastic
sediment into the deep ocean. The origin of theturbidity current is regarded as the sea landslide andthe tsunami with the earthquake, sudden vaporizationof the methane hydrate layer, submarine volcano
eruption at the bottom of sea. It is related tocomponents such as the methane hydrate tounderstand sequence of the turbidite formation. In
addition, it is very important in disaster prevention
and the environment to comprehend how a sea areachanges by the landslide. By the cause of such abackground, in this study, the numerical analysis has
been carried out to express the process of thelandslides under the sea water and collision of thesediment with the seafloor. Furthermore, the statethat seafloor is stirred up and the sediment lies
beneath the sea ground after collision is focused inthis research.
2.

reFig. 1 Concept of the particle interaction model in
MPS
A particle interacts with its neighboring particles
covered with a weight function w(r), where r is adistance between two particles. The weight functionin this study is defined as

Numerical Method
 re
 1
wr    r
 0

2.1 Governing equations
The MPS (Moving Particles Semi-implicit)
method (Koshizuka et al., 1996, 1998) is one of theparticle methods assuming a material as a particle andexpressing the behavior of the material in0  r  re 
re  r(1)

where re is a finite distance. Interactions are restricted
to this finite distance.

1

10th Asian Regional Conference of IAEG (2015)

When a particle i and its neighbors j are located atri and rj, particle number density is define asni   w r j  rij iIn addition, the different sizes of the weightfunction are used in this study. The size which is used
for the particle number density and the gradientmodel is re=2.1l0, where l0 is the distance between
two adjacent particles in the initial configuration. The
value of 2.1 was selected to avoid the concentrationof particles near the free surfaces. On the other hand,the size which is used for the Laplacian model is
re=4.0l0 (Koshizuka et al., 1996, 1998). The value of
4.0 was selected by the balance between computation
time and accuracy.
Governing equations are expressed byconservation laws of mass and momentum.
Kinematic viscous flows are considered in this study.

(2)

Because an incompressible non-viscous flows is
assumed here, the density of the fluid is constant.
Therefore, this constant value is assumed n0 because
the particles number density becomes constant.
The gradient operator is modelled using theweight function. A gradient vector is evaluatedbetween two neighboring particles i and j as

 j  i ij r j  rir  r 
j

2

(3)

i

where  is a physical quantity. Here, the angle
brackets on the left side of equation indicates asymbol for expressing an interaction model between
particles. In MPS, this equation is expressed as below,
using the gradient vector at ri which is a weighted
average of these vectors,
 i  

 j  idr j  ri w r j  ri
0 
n j i r  r 2



j



(4)

i

u v

x y

u  u , v (5)

By this equation, the divergence model used in MPS
is defined as u i  

u j  uid
 r j  ri w r j  ri
0 
n j i r  r 2



j



i



2d

 





(7)

This model is conservative because the quantity lostby a particle i is gained by particles j. Using Equation(7), the Laplacian operator is discretized tosimultaneous linear equations with respect to i. A
parameter  is introduced so that the varianceincrease is equal to the analytical solution defined as2 r j  ri w r j  ri

  j i

 w r j  ri
j i



(10)

2.2 Stability of the numerical computation
Numerical computation is unstable by explicit
calculation of convective term or diffusion term in thesemi-implicit difference method. The convective term
is not explicitly calculated, and the time difference isconfined by CFL (Courant-Friedrich-Lewy) condition,
which is defined asi  j  i w r j  ri

n 0 j i1Du  P   2 u  g
Dt


(6)

Furthermore, the Laplacian model for diffusion in
MPS is given as
 2

(9)

where,  is the density, P is the pressure, u is thevelocity vector,  is kinematic viscosity, and g is
gravity vector.
The Equation (9) is the law of conservation ofmass, and the mass conservation equation isrepresented by density, while velocity divergence isusually used in the finite difference method. The
Equation (10) is called Navier-Stokes equations, andthe left side of Equation (10) is the Lagrangianderivative involving advection terms. In MPS, theadvection terms are directly incorporated into thecalculation by moving particles. Only gravity is
considered as the external force. Two-dimensional
problems are solved in this study.
Although these equations are expressed as single
phase flow, they are extended to multi-phase flow
combined with surface tension and rigid body modeland multi-phase flow is applied in this paper.

where d is the number of space dimensions. In
addition, the particles number density ni in thelocation of particle i should be originally adopted, but
n0 is assumed to simplify a calculation.
On the other hand, the divergence theorem isgenerally expressed in the case of two dimensions.
 u 

D
0Dt

CFL  C i 

u i  tl0

(11)where ui is an absolute value of the velocity and l0 is
the distance between two adjacent particles in theinitial configuration. In addition, this value isdifferent with respect to each particle in MPS. A limit
is established for the maximum value of Courant
number Ci as follows to keep the numerical stability.

(8)maxCi  

These models are built and an incompressible
non-viscous flows are calculated.

2

tu max
 C maxl0(12)10th Asian Regional Conference of IAEG (2015)

Because the time difference t and the distance
between two adjacent particles l0 are constant, theparticle which Courant number Ci is equal to the
maximum value Cmax has the characteristic that theabsolute value of the velocity of its particle becomesmaximum value umax. If the time difference t islarger, an efficient calculation is possible. Therefore,
it is desirable to use the time difference as large as
possible to keep the numerical stability. Then thetime difference is given as follows.
t 

C max l0
u max

that floating mud was stirred up when the sediments
arrived at the seafloor was also confirmed by in-situ
experiments. In this paper, even analysis proves that
this phenomenon is able to be expressed, and it is a
purpose to see whether sediments deposit under thefloating mud after it is stirred up. Considering such abackground, the parameter of a kinematic viscosity offloating mud is focused.

3.2 Analyses condition
Fig.2 shows the analysis model of the case that the
sediment begins to slide from the ground above thewater level assuming a coast landslide, which is case1.
Fig.3 describes the case of the sediment starting toslide under the water level to simulate a submarine
landslide, which is case2. The floating mud isassumed as seabed and the sediment is assumed as
sandy gravel in both cases. By these assumptions, thedensities of both cases are determined as shown in
Table.1 but the density of the floating mud wasquoted from a measured value (Miyahara et al., 2013).
Also as shown in Table.1, three cases of thekinematic viscosity of the floating mud is tried tograsp the tendency from small viscosity to large one.
The kinematic viscosity of the sediment is determinedaccording to the past study (Kawahara et al, 2012).
The distance between two adjacent particles and the
CFL condition are experimentally determined inconsideration of calculation time (Kawahara et al,2012).

(13)

The Equation (13) is calculated every time step in thisstudy and then the size of the time difference is
Cmax=0.2
isautomaticallycalculated.
But
experimentally used for an upper limitation of
Courant number in anticipation of security because
various factors are related to numerical stability.

3.

Outline of Analyses

3.1 Purpose of analyses
In our past study, the results that the falling sediments
were accelerated in the air and were sank at constantvelocity were obtained by in-situ experiments(Matsumoto et al, 2011; Matsumoto et al, 2012), andthe analysis results using CIP and MPS methods wereagreed with the results of in-situ experiments (Isobeet al, 2012; Kawahara et al, 2012). The phenomenon

Fig. 2 Configuration of particles at 0.0 sec of case1

Fig. 3 Configuration of particles at 0.0 sec of case23

10th Asian Regional Conference of IAEG (2015)

Table.1 Analyses parameters of both cases
parameters

unit
3

fluid

sediment

floating mud

density



kg/m

1000

2100

1300

kinematic viscosity



m /s

2

1.0E-06

1.0E-03three cases; 1.0E-02, 1.0E-03, 1.0E-4

distance between two
adjacent particles

l0m

gravitygm/s

CFL condition

0.1
2

9.80665
0.2-(a) kinematic viscosity of floating mud1E-02 (m2/s)(b) kinematic viscosity of floating mud1E-04 (m2/s)

Fig. 4 Configuration of particles from 0.0 seconds to 60.0 seconds of case14

10th Asian Regional Conference of IAEG (2015)

3.3 Results
Fig.4 shows the results of two cases within case1
to clarify the difference of kinematic viscosity offloating mud (), which are two order difference
between 1E-02 (m2/s) and 1E-04 (m2/s). The
sediments raise up fluid and dive into fluid at 1.0 sec
in both cases. After the sediments are arrived at theseabed at 5.0 sec, the floating mud is stirred up to
nearby surface of fluid by collision of the sedimentsfrom 7.0 to 10.0 sec in the case of 1E-04 (m2/s).
On the other hand, the state of floating mud being
stirred up is not seen clearly in the case of 1E-02(m2/s).
Fig.5 shows the results of two cases within case2
as well as Fig.4. Although the collision with the fluid
does not occur because of starting to slide in the fluid,

the phenomenon such as case1 was similarly obtained
after the sediments are arrived at the seabed. The
sediments lie under the floating mud in all study
cases. It is a matter of course that this is caused bydensity.
To focus on the configuration at 60 sec, if thekinematic viscosity is large (Fig.4,5(a)), the floating
mud is not stirred up and the sediments subside into
the floating mud. In contrast, if the kinematic
viscosity is small (Fig.4,5(b)), the floating mud isdecreased as comparison with the initial configuration
because massive floating mud is stirred up. It isadequately considered that the state of the floating
mud corresponding to these kinematic viscosity exists.
Finally the sediments are sandwiched by the floating(a) kinematic viscosity of floating mud(b) kinematic viscosity of floating mud1E-02 (m2/s)
1E-04 (m2/s)
Fig. 5 Configurations of particles from 0.0 seconds to 60.0 seconds of case2

5

10th Asian Regional Conference of IAEG (2015)

in 2011 off the Pacific coast of Tohoku Earthquake
and its relations with the tsunami are discussed(Kawamura, 2015). We think that it is possible fornumerical analysis technique to be made use of forthe elucidation of these phenomena.

mud after accumulation of the sediments, and thenthe alternation of beds as shown in Photo.1 isappeared.
Although it is assumed that the density of
sediments is larger than that of floating mud in thisstudy, the study that the densities of sediments andfloating mud are equal was carried out in advance. If
the kinematic viscosity of the floating mud is small, it
is not mixed up as shown in Fig.6. But when thevalue is large, like 1E-02, it is stirred up. This isindicated in the past study (Miyahara et al., 2013).
Therefore, stirring up the floating mud and thesediments sinking into the floating mud are caused bylarge difference of the density between them (floating
mud < sediment) and small viscosity of the floating
mud.

References
Koshizuka, S. and Oka, Y. (1996): Moving-particle
semi-implicit method for fragmentation ofincompressible fluid, Nuclear Science andEngineering, 123, pp. 421-434.
Koshizuka, S., Nobe, A., and Oka, Y. (1998):
Numerical analysis of breaking waves using themoving particle semi-implicit method, International
Journal for Numerical Methods in Fluids, 26, pp.
751-769.
Matsumoto, K., Owatari, T., Ohno, H., Miyahara, T.,
Hachimura, T., Buma, N., and Oshima, M. (2011):
Measurement and survey in water by using sonar
exploration, Journal of the Japan Society ofEngineering Geology, Vol. 52, No.4, pp. 128-136.
Matsumoto, K., Nagaoka, S., Ohno, H., Isobe, Y.,
Kawahara, H., Hachimura, T., Buma, M., Oshima,M., and Tosaka, H. (2012): Free-setting
phenomenon of particles at large Reynolds numbers– free-setting experiment, Journal of the Japan
Society of Engineering Geology, Vol. 53, No.2, pp.
70-79.
Isobe, Y., Kawahara, H., Matsumoto, K., Ohno, H.,
Tosaka, H., Nagaoka, S., Hachimura, T., Buma, M.,
and Oshima, M. (2012): Analysis of free-setting
phenomenon of reclaimed solids at offshore
disposal site by CIP method, Journal of the Japan
Society of Engineering Geology, Vol. 53, No.5, pp.
218-226.
Kawahara, H., Isobe, Y., Matsumoto, K., Ohno, H.,
Tosaka, H., Nagaoka, S., Hachimura, T., Buma, M.,
and Oshima, M. (2012): Analysis of free-setting
phenomenon of reclaimed solids at offshore
disposal site by MPS method, Journal of the Japan
Society of Engineering Geology, Vol. 53, No.5, pp.
227-234.
Miyahara, T., Nagaoka, S., Ohno, H., Hachimura, T.,
Kawahara, H., Isobe, Y., Hida, Y., Tomiyta, Y., andTosaka, H. (2013): For impact phenomenon to waterbottom by settling of particles in water on offshore
disposal site, Journal of the Japan Society ofEngineering Geology, Vol. 54, No.3, pp. 104-113.
Lee, H. J. (2009): Timing occurrence of large
submarine landslides on the Atlantic Ocean margin,
Marine Geology, Vol. 264, pp. 53-64.
Kawamura, K. (2015): Submarine landslides andmethane hydrate in cases of Norway and Japan,
Geotechnical Engineering Magazine, Vol. 63, No.2,
Ser.No.685, pp. 14-17.

Fig. 6 Configuration of particles at 60.0 seconds when
the density of the sediment is equal to the floating mud

Photo. 1 Sediments

4.

Conclusions

The phenomenon that the floating mud is stirred
up by collision of the sediments with the seabed andthe sediments lie under the floating mud after that
was able to be caught well by MPS. However, thekinematic viscosity of not only floating mud but also
sediments cannot be identified because those are
depended on the materials properties. The various
case studies such as the depth of sea, slope angle andthe amount of landslide are executed as future work.
In the occurrence of the submarine landslides, thecausation with the methane hydrate is doubted (Lee,2009). It is indicated that the landslide of the seabedwhich is width of 2-3km and length of 10km occurred6