MATHEMATICAL MODELLING OF SUSPENDED SEDIMENT MOVEMENT

USING FINITE DIFFERENT SCHEME

Edy Sulistyono, Ir.,MT, Profesional Utama SDA
KTA No : 074865

The Objective.

The objective of this project is to simulate movement and distribution of dredging

disposal whic is damped in the coastal area and to simulate movement of sediment
transport in the rivermouth area. Simulation result will be used to dicide the location of
disposal dumping point with minimum impact to coastal line environment and to estimate
time period of the entrance channel re excavation (maintenance).

Tujuan Kegiatan

Maksud kegiatan ini adalah untuk mensimulasikan pergerakan dan distribusi buangan
hasik pengerukan yang akan berdampak pada pantai dan mensimulasikan pergerakan
sedimen transport di area muara sungai. Sedang tujuannya adalah menentukan titik
lokasi terbaik untuk membuang hasil dreging yang memberikan dampak lingkungan
minimum pada pada lingkungan pantai dan muara, serta untuk memperkirakan berapa
kala ulang yang diperlukan untuk melakukan penggalian ulang dalam rangka
pemeliharaan saluran masuk.

Project Location.

Location of this study is in the Flores sea, offshore area of South Lasalimu district, Buton
regency, south east Sulawesi, see figure 1.1. , 1.2. and 1.3

Lokasi Kegiatan.

Study ini berlokasi di perairan flores, di oofshore kecamatan Lasalimu Selatan,

kabupaten Buton, provinsi Sulawesi selatan, lihat gambar 1.1, 1.2, dan 1.3
 PROJECT LOCATION

Figure 1.1. Buton Regency

Figure 1.2. South Lasalimu District

 Figure 1.3. Entrance Channel.

Theoritical Background.

Hydrodynamic Model
The Hydrodynamic model is described by the following equations. These equations, the
conservation of mass and momentum, are derived from vertical integration of Navier-
Stokes formulae.

Latar Belakang Teori.

Model Hydrodinamic
Model dynamic di formulasikan melalui peramaan berikut. Persamaan ini merupakan
gabungan antara persamaan konservasi massa dan konservasi momentum, yang
diturunkan dari persamaan Navier-Stokes.

∂η ∂ P ∂ Q
+ + =0
∂t ∂ x ∂ y

( ) ( )
∂ P ∂ P2 ∂ PQ
+ +
∂t ∂ x D ∂ y D
∂η gPn 2 √ P2 +Q2
+gD +
∂x D7/3

(
− ∂ ED ∂
∂x
P
∂x D ( ))
− ∂ ED ∂
∂y
P
∂y D ( ( ))
−Ω Q−fVV x=0

( ) ( )
∂Q ∂ PQ ∂ Q 2
+ +
∂t ∂ x D ∂ y D
+gD
∂ η gQn2 √ P2 +Q 2
∂y
+
D7/3
 (
− ∂ ED ∂
∂x
Q
∂x D ( ))
− ∂ ED ∂
∂y
Q
∂y D (
+Ω P−fVV y =0 ( ))
The following symbols are used in the equations:

P,Q : flux densities in x- and y-directions

D : water depth（=η+ h ）
η : surface elevation
n : Manning roughness
E : eddy viscosity
f : wind friction factor
V ,V x ,V y : wind speed and components in x- and y- directions
Ω( x , y ) : Coriolis parameter
g : gravity acceleration
ρw : density of water

These equations are described by the Finite Difference Method and solved by using the
Alternating Direction Implicit (ADI) technique.

Silt Transport Model

The silt transportation is calculated by solving the so-called advection-dispersion
equation:

Model Transportasi Sedimen

Transportasi sedimen di analisis berdasarkan persamaan yang dikenal sebagai persamaan dispersi
- advecsi Navier Stokes berikut :

∂C ∂C ∂C
+u +v =
∂t ∂x ∂y
1 ∂
h ∂x (
hD x
∂C 1 ∂
+ )
∂x h ∂ y
hD y
∂C
∂y ( 1
+Q L C L −S
h )
Symbol List

C : depth averaged mass concentration（kg/m3）

u,v : depth averaged flow velocities (m/s)
Dx,Dy : dispersion coefficients (m3/s)
h : water depth (m)
S : accretion/erosion term（kg/m3/s）
QL : source discharge per unit horizontal area (㎥/s/㎡)
CL : concentration of source discharge (㎏/㎥)
 The advection-dispersion equation is solved using an explicit, third-order Finite
difference scheme, known as the ULTIMATE scheme, Leonard (1991).

Persamaan dispersi – adveksi diselesaikan dengan methode explisit, skeme

persamaan beda hingga, yang dikenal sebagai skema ultimate Leonard (1991)

Deposition

In this model, a stochastic model for flow and sediment interaction is applied. This
approach was first developed by Krone (1962).
Krone suggested that the deposition rate is expressed by the equation:

Pengendapan
Di model ini, untuk interaksi antara aliran dan sedimen diterapkan prinsip model acak
atau stocastik. Pendekatan ini untuk pertamakalinya dikembangkan oleh Krone (1962).
Krone menyarankan bahwa laju proses pengendapan di tuangkan dalam persamaan
sebagai berikut :

S D=w s c b pd ）

ws : settling velocity（m/s）
cb : near bed concentration（kg/m3）

( )
τb
Pd ¿ 1− , τ ≤τ cd
: probability of deposition τ cd b
τb : the bed shear stress（N/m2）
τ cd : critical bed shear stress for deposition（N/m2）

Erosion

Erosion process can be described in two ways depending upon whether the bed is
dense and consolidated or soft and partly consolidated, Mehta et al. (1989).

For dense, consolidated bed:

Erosi

Proses erosi dapat dijelaslan melalui du acara, tergantung pada apakah dasar saluran
padat dan terkonsolidasi atau lunak dan setengah terkonsolidasi .

Untuk kondisi terkonsolidasi :

S E =E (τ b /τ ce -1)n, τ b >τ ce

where
 E : erodibility of bed（kg/m2/s）
τ ce : critical bed shear stress for erosion（N/m2）
n : power of erosion

For soft, partly consolidated bed:

Untuk kondisi dasar lunak atau setengah terkonsolidasi :

S E =Eα exp [α (τ b -τ ce -1)1/2], τ b >τ ce

where

α : coefficient (m/N1/2)

Deposition in Entrance Channel Model.

In order to estimate maintenance cost, it is necessary to estimate return period of

redredging work. Considering design concept of the entrance channel, excavation will be
done when the deposition in the channel more than tolerance depth.

Based on tolerance depth and rate of sediment deposition, return period of excavation
work for maintenance purpose can be estimated.

Formula

Untuk mengestimate baiayan pemeliharaan, diperlukan perkiraan kalua ulang pekerjaan

penggalian sedimen. Mempertimbangkan konsep disain saluran masuk atau entrance, ,
penggalian akan dilakukan setelah sedimentasi di saluran sudah melebihi kedalaman
toleransi.

Berdasarkan kedalaman toleransi dan laju sedimentasi, kala ulang pekerjaan penggalian
untuk keperluan pemeliharaan dapat di estimas menggunakan formula sebagai berikut.

Formula

k s =2. 50∗D50

C=18 *log
( )
12∗h
ks

1
C =18 *log
( )
12∗h
D 90
 ( )
1. 5
C
μ=
C1

( )
ρs −ρw
Δρ=
ρw

U
q b =5∗D 50∗ ∗μ∗√ g . exp
C μ∗U 2 (
−0 . 27∗Δρ∗C2∗D50
)
With :

D50 = sediment grain size 50% screen pass

D90 = sediment grain size 90% scren pass
h = water depth
U = current velocity
ρs = sediment specific gravity
ρw = specific gravity of sea water
qb = sediment rate/m width

Defining the Hydrodynamic Model

Identify the hydraulic problems

When preparing to do a study of, for example. The currents in the Bay of X-town or
Extreme water levels along a certain stretch of coast, we have to assess the following
before we start to setup the model:

What are the important hydraulic phenomena ? (are tide, wind, wave breaking,
stratified flow etc. important ?).

What is the area of influence of these phenomena ? (for example how big an area has
to be included in the model so that the surge from a hurricane can be computed
correctly at the point of interest and what area can be neglected ?)

What is the time scale of these phenomena ? (e.g. are we looking for eddies in an
entrance to a harbor or are we going to compute a larger number of tidal constituents in
an area ?).

Perencanaan Model Hydraulic

1. Identifikasi persoalan hydraulic yang terjadi

2. Ketika menyiapkan pekerjaan studi, disiapkan antara lain sebagai contoh. Arus di
teluk di desa X, atau elevasi maksimum sepanjang ruas pantai, kita harus
memperkirakan sebelumnya untuk setup model sebagai berikut :
 a. Apakah phenomena hidraulik yang dominan terjadi (apakah pasang surut,
angin, gelombang pecah, aliran berlapis, apakah penting di pertimbangkan?)
b. Apakah pengaruh dari phenomena diatas ? (sebagai contoh, seberapa besar
area yang terpengaruh dan harus dimodelkan, sehingga pengaruh surge dari
badai dapat dihitung dengan benar pada titik tertentu, dan dimana saja area
yang tidak perlu dipertimbangkan).
c. Berapa skala yang diperlukan untuk menggambarkan phenomena yang terjadi?
(sebagai contoh apakah kita perlu mempertimbangkan terjadinya eddies current
pada pintu masuk Pelabuhan? Atau apakah kita perlu menghitung banyak
parameter pasang surut diarea model ?)

Select Model area and Simulation Periods

Draw up our model area on a sea chart showing the area of interest and the area of
influence (as described above). This is normally an iterative process as, on the one
hand we should keep the model area as small as possible while, on the other hand we
have to include the total area of influence etc. Guidelines for Selecting the model area
is given under Bathymetry. We also have to consider the Boundary Conditions when
deciding on wer model area and the grid spacing (see Selecting the grid spacing under
Bathymetry).

The next step is to determine the time step to be used in the model. A description of
how to select this is found under Time step in the Reference Manual. We also have to
decide how long a simulation period is required. For example, to obtain a steady state
current pattern around a harbor we may only need some hours whereas we may need
to simulate 14 days to be able to distinguish certain tidal constituents. In the simula tion
period we should also allow time for the model to "warm up". For a small model this
might take only a few hours while, for a large tidal model it may take up to a couple of
days. We can determine the required "warm up" period during the first calibration
simulations. If we have a very long "warm up" period we may wish to use the "hot start"
facility.

Collecting Data

This task may take a long time if, for example, we have to initiate a monitoring
program. Alternatively it may be carried out very quickly if we are able to use existing
data which are immediately available. In all cases the following data should be
collected:

Bathymetric data such as sea charts from local surveys or, Boundary data, which
might be measurements (existing or planned specifically for wer model), observations,
tidal constants etc.

• Wind data
• Source and sink data (i.e. flow magnitudes and directions from rivers, sewage
outfalls, cooling water intakes etc.)
• Information on the bed resistance

Calibration and validation data: these might be measured water levels or current speed
and current direction at a number of points in the area of interest. Tidal constituents
 can also be used if a tidal analysis program is available. Float trackings or aerial
photos of eddies are other possibilities.

Setting up the Model

What does it mean? "Setting up the model" is actually another way of saying
"transforming real world events and data into a format which can be understood by the
numerical model MIKE 21". Thus, generally speaking, all the data collected have to be
resolved on the spatial grid selected and resolved in time with the time step selected.

Bathymetry

As mentioned under Bathymetry the use of the sea charts for the Hydrodynamic
Model is the most important task in the modelling process. We should therefore, to as
great an extent as possible, try to adhere to the rules listed under Bathymetry
(Selecting the model area and Selecting the grid spacing). A description of how the
depth values are entered into MIKE 21 format is given in the section Specifying the
bathymetry (also under Bathymetry).

Depending on which combination of boundary input we are going to use (see General
description under Boundary Conditions) we have to prepare one of the following types
of data for each boundary:

a. A few parameters describing, firstly, a constant or sinusoidal variation of the

water level or the flux and, secondly, the flux calculation procedure.

b. A type 0 data file (for example if we have measurements of the water levels at a
point on were boundary). If we have the data as an ASCII file we can read them
into the MIKE 21 data file format as described in the Pre- and Postprocessing User
Guide

c. A type 1 data file containing either water level or flux variations along the open
boundary during time (e.g. boundary conditions generated by the tidal analysis
and prediction program). A description of how these data can be entered into
NIIKE 21 can be found in the Pre- and Postprocessing User Guide.

d. Transfer data file (a type 1 data file) containing water levels and fluxes from an
encompassing model. A description of how these data are extracted is given in
the Pre- and Postprocessing.

Wind data

Wind data can be given in different ways. Please see Specifying the wind conditions

Source and sink data

Discharges etc. can also be given in various ways. Please see Specifying source and
sink under Sources and Sinks.
Bed resistance
 Bed resistance is normally not specified on the basis of measurements but used as a
calibration parameter. However, areas with high resistance can be specified initially; see
Specifying the bed resistance under Bed Resistance for a description of how to do this.

Calibration and verification data

If we have time series of measurements we can enter them into the MIKE 21 data file
format as described in the Pre- and Post processing User Guide (Reading ASCII files).
We can then plot them for comparison with simulation results.

Calibrating and Verifying the Model

Purpose

Having completed all the tasks listed above we are ready to do the first hydrodynamic
simulation and to start on the calibration of the model.

The purpose of the calibration is to tune the model in order to reproduce satisfactorily
known/measured conditions for a particular period known as the calibration period. The
calibrated/tuned model is then verified by running one or more simulations for which
measurements are available without changing any tuning parameters. This should
ensure that simulations can be made for any period similar to the calibration and verifica-
tion periods with satisfactory results. However, we should never use simulation results,
whether verified or not, without checking that they are reasonable.

Calibration periods

The periods which we select for calibration of the model should include all the different
situations we wish to investigate in the production runs. If, for example, we have a
complicated tidal pattern in were model, we should calibrate the model for neap tide as
well as for spring tide.

However, as we must have some measurements/observations against which to

calibrate and, as these measurements are often only available for short periods, we
may only have a few periods from which to choose.

Calibration parameters

When we run wer calibration run for the first time and compare the simulation results to
wer measurements we will, in many cases, see differences between the two. The
purpose of the calibration is then to tune the model so that these differ ences become
negligible. We can change the :

a. Bed resistance
b. Eddy viscosity
c. Bathymetry
d. Boundary conditions

Recommendations on how the specifications can be changed are given in the

Reference Manual
 Please note that the currents computed by Mike 21 are depth averaged. Thus when we
compare them with measurements we need measurements at several levels (depths)
or we must assume a theoretical current profile (logarithmic, for example).

Verifying the model

When we have finished the calibration we can run one more period for which we have
measurements (or other data) without changing the calibration parameters. If we then
get a satisfactory agreement between the simulation results and the measurements we
can consider were calibration to be successful:

Running the Production Simulations

As we now have a calibrated model we can get on with the "real" work: doing wer actual
investigation. This will, in some cases, only include a few runs if, for example, the current
patterns for two different lowest of a harbor are to be computed. In other cases hundreds
of runs will have to be made. An example of such a study is an investigation of extreme
water levels. This might involve the simulation of the 200 most severe storm events
during a 30 year period. Based on these simulations the extreme statistics can be made.

If we carry out many production simulations the short cuts for running the simulations, as
described in chapter 2 in the Preand Postprocessing User Guide, will be especially
useful. We might also use several catalog directories.

Presenting the Results

Throughout a modelling study we are working with large amounts of data and the best
way of checking them is, therefore, to look at them graphically. In only a few cases, such
as when we check were bathymetry along a boundary, should we look at the individual
number. Much emphasis has, therefore, been placed on the capabilities for graphical
presentation in Mike 21, and it is an area which will be expanded and focused on even
further in future versions.

Essentially, one plot gives more information than scores of tables and, if we can present
it in colors, were message will be even more easily understood.
Model result

Option no 1

Position : Coordinat UTM 51 M 525170 9421017 ( 2 mil from the coast line),
bottom elevation -301.57 m , at Grid ( 227 , 209 ), discharge 5 m3/s , V = 7.4
m/s .

Source Position Option 1

Suspended load distribution

Option no 2

a. Source

Position of the source 4 mil ( 6.44 km ) from the coast line, coordinate UTM 51 M
( 528420, 9421096), bottom elevation -458.25 m , and grid ( 292 , 209 ),
discharge 5 m3/s , V = 7.4 m/s .

Coastal Line

Suspended load distribution

Conclusions.

a. Discharging disposal material of channel excavation in the position approximate

2 mile away from coastal line and influence of wind force whic is flowing from
 east to west , show that suspended concentration in the coastal line very low, just
0.0000000061 ppm much lower than 20 ppm ( mension in Kepmen LH no
51/2004.

b. Wind direction during rainy season flow from west to east do not give significant
influence to the movement and distribution of the suspended load, otherwise
during dry season when wind flow from east to west will potentialy drive
suspended load flow to the west and coastal line

c. During simulation, show that influence of wind force do not significant in the
offshore area. Slope of sea bottom very steep, therefore the depth increase
significanly. Sea bottom slope give advantage significanly because dumping
point can be selected nearer to the coast line.

d. From the simulation result of the mud model can be cncluded tha concentration
of the suspendel load during rainy season more than 20 ppm (100 ppm – 1000
ppm), it is mean greather than limitation as mension in the Kepmen LH no
51/2004. In this simulation wind flow do not influence significanly.

e. From the simulation result of the suspended load model in the entrance channel,
can be shown that deposition rate of the suspenden load do not influence
significantly to the redustion of channel depth ( less deposition ). Reduction of
channel depth, is influenced significantly by deposition process of bed load
during flood ( rainy season, 76m3/m/year ).

Recommendations

a. Dumping activity of the disposal material can be carried out during rainy season
and dry season, but it is recommended carried out during rainy season, in the
position at least 2 miles from the coast line to the offshore side.

b. If rate of deposition in the entrance channel 76 m3/m/year, and permitted bed

load deposition less than 0.5 m ( deposited bed load shoud be excavated 73,124
m3/m/year), then maintenance excavation shoud be done at least every year.