SEDIMENT TRANSPORT AND DEPOSITION IN THE OK TEDI-FLY

RIVER SYSTEM, PAPUA NEW GUINEA:

THE MODELING OF 1998 - 1999
Report 7
Environment Department, Ok Tedi Mining Ltd., Tabubil, Papua New Guinea

Yantao Cui and Gary Parker

St. Anthony Falls Laboratory
University of Minnesota

June, 1999

1. INTRODUCTION

The Ok Tedi-Fly river system is located in the country of Papua New Guinea.
The Ok Tedi copper mine is presently disposing both tailings and rock waste into the
river system near the headwaters of the Ok Tedi. This sediment flows from the Ok Tedi
to the Fly River, eventually reaching the Gulf of Papua.

This document provides a brief report concerning the results of the modeling
effort of 1998 – 1999, which were obtained as part of an overall risk assessment of mine
activities. The results themselves are included in a set of two Zip discs which are
described later in this report.

Mine
Reaches 1 Reach 3a
Kiunga
Reach 2
Konkonda
Reach 3 D'Albertis Junction

flow Reach 4

Manda Reach 5 Strickland

River

Obo Everill
Junction

Figure 1: Definition of reaches.

 It is useful to begin by defining the relevant reaches of the river itself, as shown in
Figure 1. They are defined as follows for the purpose of this report:
• Reach 1: the Ok Mani, Ok Tedi, Ok Mabiong etc. upstream of the Ok Tedi and the
Ok Mani near Tabubil;
• Reach 2: the coarse-bed Ok Tedi from the confluence with the Ok Mani to
Konkonda, i.e. the point of transition from gravel-bed to sand-bed stream;
• Reach 3: the sand-bed Ok Tedi from Konkonda to D’Albertis Junction, i.e. the
confluence with the Fly River;
• Reach 3a: the upper Fly River from Kiunga to D’Albertis Junction;
• Reach 4: the middle Fly River from D’Albertis Junction to Manda;
• Reach 5: the middle Fly River from Manda to Obo.
The downstream end of Reach 5 is near, but not at Everill Junction, where the Fly River
joins the Strickland River.

The second author of this report has been served as the sediment transport
consultant for Ok Tedi Mining Ltd. since December of 1986. During this period a number
of sediment transport models have been developed for the company. Only those of
direct relevance to the present report are discussed here. The model Okgrav5 was used
to model sediment transport from the mine to a point near Konkonda, i.e. Reaches 1 and
2. The model OkFly was used to model sediment transport from the upstream points of
Konkonda and Kiunga to Obo, i.e. Reaches 3, 3a, 4 and 5. Okgrav5 is a program
running from DOS written in the Pascal language. OkFly runs from DOS in Fortran, but
has been provided with a Visual Basic interface.

2. GOALS OF THE MODELING

The primary focus of the modeling was the sand-bed Ok Tedi-Fly river system,
i.e. Reaches 3, 3a, 4 and 5. The modeling was designed to answer the following
questions over the 70 year period 1985 - 2055. The mine commenced sediment
disposal in 1985, and will likely close no later than 2010.
a. What will be the pattern of bed aggradation, and subsequent degradation over the
lifetime of the mine and beyond?
b. What will be the loadings of suspended sediment?
c. How will mine-derived sediment mix with natural sediment in the bed of the river?
d. What rates of overbank deposition will be observed?
e. What will be the provenance and size distributions of the bed and overbank sediment
deposits?
f. How will flood frequency respond to the sediment loading?
The answers to these questions are to be used in turn to develop a risk assessment
including such factors as dieback of riparian forest, potential for acid drainage, effect of
sediment on aquatic fauna, fate of copper etc.

3. ORGANIZATION OF THE RESULTS

The modeling results are reported in two pairs of Zip discs. The first pair is dated
1/29/99; these results are referred to as the January results. The second pair is dated
5/10/99; these results are referred to as the May results. As will be seen below, the
January results are more pessimistic and the May results are more optimistic. The May
results are recommended as the likelier of the two to be quantitatively accurate.

Results are presented for seven options, as outlined below. These options are
based on possible schemes for sediment control.
No Mine scheme: in this case it is assumed that the mine was never built. The
predictions should show negligible changes over the 70 years of modeling.
Null scheme: the mine is assumed to continue until year 2010 with no sediment
remediation.
Scheme A: dredge above Konkonda at 10 Mt/a for 1998, 1999, close mine on May 31,
2000.
Scheme B: dredge for three years, then tailings disposal to piped storage area for mine
lifetime.
Scheme C: dredge for two years at 10 Mt/a and stop dredging for mine lifetime.
Scheme DL: dredge low (15 Mt/a) for mine lifetime.
Scheme DH: dredge high (19 Mt/a) for mine lifetime.

The two pairs of discs are organized according to the same structure. This
organization is explained in the Word file “[Link]” on the root directory of the first
of each pair of discs.

The data needed to approach the questions of Section 2 were generated using
OkFly. In order to generate the input for that model, however, it was first necessary to
run Okgrav5. With this in mind, Okgrav5 is discussed first.

4. DESCRIPTION OF MODEL Okgrav5

Model geometry The model Okgrav5 covers Reaches 1 and 2. The model has
the following subreaches:
a. a debris flow runout zone from the Northern Dumps to a point on the Ok Gilor;
b. a subsequent fluvial reach extending to the Ok Tedi where it joins the Ok Mani;
c. a debris flow runout zone from the Southern Dumps to the confluence of the Ok Mani
and the Ok Kumiup;
d. a fluvial zone extending down the Ok Mani from the Ok Kumiup to the confluence
with the Ok Tedi; and
e. a fluvial zone extending down the Ok Tedi from the confluence with the Ok Mani to
the gravel-sand transition (near Konkonda).
Subreaches a, b, c and d are in Reach 1 and subreach e constitutes Reach 2.

Flow calculation Reaches 1 and 2 are assumed to be sufficiently steep to allow

for the assumption of quasinormal flow. As a result the full backwater equations are not
used to predict the flow. Rather than use daily discharges as input, the flows are
computed for a number of points in a flow duration curve. The flow duration curve is
specified for each node along the system, but is assumed to be invariant in time. The
model thus cannot distinguish between “wet” and “dry” years.

Sediment feed The model focuses on gravel transport. It assumes three

sources of gravel: natural gravel, gravel derived from rock waste from the mine and
gravel derived from slide material from either the walls of Harvey Creek (southern
dumps) or the Vancouver Slide (northern dumps). The model also tracks mine-derived
sand and silt. Sand and silt derive from rock waste, slide material and tailings. An
additional source of sand and silt is abrasion, as discussed below.

Rock waste is fed in at the northern and southern dumps based on observed or
projected annual dumping rates. Observed and predicted annual disposal rates are
updated based on information from the mine. The rock waste is initially partitioned into
gravel an coarser, sand and silt fractions based on expected breakdown from the dump
site to the base of the debris flow runout zones. Much of the gravel, and some of the
sand and silt are captured in the debris flow runout zones as they aggrade with constant
slope. The rest of the gravel is delivered to fluvial subreaches b and d.

Sediment is also produced indirectly in the form of slide material as a

consequence of the dumping of rock waste. This material takes the form of the
Vancouver slide mass of the northern dumps and erosion of the walls of upper Harvey
Creek in response to the disposal of waste sediment. Observed and predicted
schedules for this delivery are occasionally updated based on new information from
consultants working for the mine. The initial breakdown into gravel, sand and silt
fractions is assumed to be identical to that for rock waste.

Fluvial gravel transport Gravel transport is computed using the Parker (1980)
surface-based transport relation. It treats sediment mixtures ranging from 2 mm to 512
mm on a fractional basis. Three layers are assumed in the model: a bedload layer, a
bed surface (active) layer and a substrate layer. The average transport rate and
transport grain size distribution for a year is based on calculations from the flow duration
curve.

Abrasion of gravel The gravel-sized rock waste is relatively easily abraded.

Abrasion is included in the Exner relation for sediment continuity. The products of
abrasion are computed to be silt and sand, mostly the former. The coefficient of
abrasion is assumed to be relatively large near the mine and relatively smaller farther
downstream, to reflect the gradual dominance of abrasion-resistant gravel as the friable
material is ground out. Although the model allows for different abrasion rates for several
lithologies, only two lithologies have been used; harder natural gravel and softer mine-
derived gravel.

Sediment continuity: The debris flow runout zones are assumed to aggrade at a
constant slope of 8 percent, based on field measurements. The fractions of gravel, sand
and silt captured in this zone are also based on field measurements. The deposits are
assumed to be matrix-supported. The area of deposition is assumed to expand as the
bed aggrades in accordance with the geometric expansion of the respective river
valleys. The downstream ends of the debris flow runout zones are constrained by
confluence with a river with sufficient discharge to move the sediment fluvially.

Within the fluvial zones b and d (Ok Gilor-Ok Mabiong-Ok Tedi and Ok Mani) the
bed is allowed to aggrade or degrade between any pair of nodes according to whether or
not the gravel supply exceeds the sediment delivery. The width of deposition is divided
into a channel width, which is used to compute sediment transport, and a wider valley
width, over which aggraded sediment is allowed to deposit uniformly. This configuration
approximates the reworking of the valley flat by a braided or wandering gravel-bed
stream. The channel widths were determined from aerial photography and adjusted in
the “zeroing” procedure discussed below. Valley geometry was assumed to be V-
shaped, so that valley width increases with bed aggradation. This was particularly
important for the Ok Mani, which is now an order of magnitude wider than the pre-mine
state. Valley shape was obtained from topographic maps.

Within zones b and d the deposit is assumed to be clast-supported, but sand and
some silt was allowed to fill the pore space as the bed aggrades. This material is
released as the bed later degrades into its own deposits.

Exactly the same algorithm is used for zone e (Reach 2, i.e. the Ok Tedi from the
confluence with the Ok Mani to Konkonda), with the exception that a) valley width is not
allowed to change according to bed height and b) sand and silt are not stored in the
interstices of the gravel as it aggrades. Because the river later degrades into its own
deposit upon mine closure, it is necessary to store in memory both the thickness of the
deposit and the grain size distribution. The stored grain size distributions consist of
averages of the deposit over a specified thickness.

Fluvial transport of sand and silt Sand and silt deriving directly from rock waste,
slide material and tailings are allowed to deposit in or be eroded from zones a, b, c and
d. They are neither deposited in nor eroded from zone e. Sand and silt is be produced
by abrasion in zones b, d and e. The sand and silt transport at the confluence of the Ok
Tedi and the Ok Mani (upstream end of Reach 2) is delivered directly to Konkonda
(downstream end of Reach 2) with no intervening storage. The sand and silt produced
by abrasion along Reach 2 is also delivered directly to Konkonda with no storage. This
delivery is specified in terms of a mean annual rate. The assumption of no storage of
sand and silt in Reach 2 is based on estimates which suggest that the error in doing so
is small.

Source partitioning Okgrav5 partitions the sand and silt arriving at Konkonda into
six types: sand from rock waste, sand from slide material, sand from tailings and silt from
rock waste, silt from slide material and silt from tailings. This partitioning is for the
purpose of tracking copper, which is present in different amounts according to the
provenance of the material.

Model zeroing Models of bed evolution are prone to wild errors if used without
prudent adjustment. The model was adjusted in the following way. Assuming that no
mine was ever built, gravel infeed and channel width were adjusted so as to predict a
river system that changes a negligible amount over the planned lifetime of the mine.
The natural gravel infeed rates were determined in this way. They are very small
compared to rates of loading from the mine.

Model output The output for Okgrav5, along with sediment budgets determined
from the output, are given in the Excel workbook [Link] in directory General of the first
of the two discs. In running Okgrav5 to determine the fate of the gravel, only two
schemes needed to be considered; Scheme A, according to which the mine closes in
year 2000, and all the other schemes.

The workbook includes worksheets specifying production schedules, gravel

budgets for all relevant reaches, and schedules for the delivery of sand and silt to
Konkonda. The computed delivery schedules for sand and silt to Konkonda are given
for every scheme in worksheet SchemeResults. The sand and silt delivery rates are
broken down according to provenance. The assumed grain size distributions of the sand
and silt delivered to Konkonda (and Kiunga for the case of natural sediment) are given in
worksheet GrainSizeDist. Estimated natural loadings of sand and silt at Konkonda and
Kiunga are given in worksheet NatLoad. How these were determined is explained in
more detail below.

Worksheets SchemeResults, NatLoad and GrainSizeDist provide information

concerning sediment inputs used in the model OkFly to model Reaches 3, 3a, 4 and 5.
The flow of the models is as follows. First Okgrav5 is run to determine the fate of the
gravel in Reaches 1 and 2, and the sand and silt delivery rates from the mine at
Konkonda. The output of Okgrav5 is then manipulated with the spreadsheet Ok99.
Then OkFly is run to determine the fate of the sand and silt in Reaches 3, 3a, 4 and 5.

5. DESCRIPTION OF MODEL OkFly

Model network Although the model was designed to be implemented for the
case of Reaches 3 (lower Ok Tedi) and 3a (upper Fly River) flowing into Reaches 4 and
5 (middle Fly River), it was written for a river network of more complexity. A schematic
diagram of an arbitrary set of river segments that the model could handle is illustrated in
Figure 2.

9
5
8 10
6
4 11
7
2
3

Figure 2. Illustration of the river network scheme of Okfly.

Note that the numbering system actually used by OkFly does not correspond to the
numbering scheme adopted here for data presentation. Reaches 4 and 5 of the
modeling exercise correspond to Reach 1 of OkFly, Reach 3 of the modeling exercise
corresponds to Reach 2 of OkFly and Reach 3a of the modeling exercise corresponds to
Reach 3 of OkFly.

Nodes and channel geometry. The locations and characteristics of all 47 nodes
used in the study are specified in worksheet Geometry of workbook okfly99 in directory
General of the first of the two Zip discs. Information there includes node number,
upchannel distance from Obo, channel width, floodplain width, channel area between the
given node and the node just upstream and floodplain area between the given node and
the node just downstream. The model also requires information concerning initial
channel and floodplain elevation.

The channel and floodplain are assumed to have the simple geometries
illustrated in Figure 3. The floodplain width used in the numerical model is not the width
of the entire floodplain, but some “effective” width somewhat wider than the width of the
meander belt. This “effective” width was used to account for the fact that much of the
water in the floodplain is locally derived, and does not come from the channel. This fact
is shown clearly in aerial photographs of the river in flood, showing dirty water in a band
near the channel and ponded clear water away from the channel. The model also
requires information concerning initial channel and floodplain elevation.

floodplain

channel

Figure 3. Simplified channel and floodplain geometry.

Discharges and flow calculation As opposed to Okgrav5, OkFly uses daily

discharges based on flow records (including discharge and water surface elevation) at
Konkonda, Kiunga, Kuambit, Manda and Obo. The daily discharges are based on two
cycles of “dry,” “medium,” “wet” and “El Nino” years at each site. The flow records for
each such year are based on a real record. The specification of flows from year to year
is given in worksheet Hydrology of workbook okfly99 in directory General of the first of
the two discs. It should be noted that “dry” years, for example, are not necessarily dry
everywhere. An interpolation scheme was used to determine the daily discharges at
nodes for which flow records were not available.

The flow calculation is based on steady, gradually varied flow computed

upstream from Obo, where water surface elevations were specified. In the event that
the Froude number of the flow is in excess of 0.75 (an event that rarely if ever occurred)
the flow is treated as quasinormal. The technique does not assume constant water
discharge between nodes, as this is allowed to vary according to an interpolation
scheme, reflecting the interaction between channel and floodplain hydrology. The flow is
partitioned between channel and floodplain in a standard way.

Hydraulic resistance OkFly uses the Brownlie (1981) model of hydraulic

resistance to compute in-channel resistance. Expressed in terms of friction slope Sf, it
takes the form

−1.361
⎛ H ⎞
S f = 0.02054R 1.286
Fg
2.572
⎜⎜ ⎟⎟ σg
0.4130
for lower flow regime (1a)
⎝ D 50 ⎠
−1.304
⎛ H ⎞
S f = 0.01252R 1.086
Fg
2.172
⎜⎜ ⎟⎟ σg
0.2785
for upper flow regime (1b)
⎝ D 50 ⎠
where H denotes in-channel depth, D50 denotes the median grain size of the sediment in
the surface (active) layer of the bed, R denotes the submerged specific gravity of the
sediment and σg denotes the geometric standard deviation of the sediment in the surface
layer. In addition, Fg is a grain Froude number given by

Qw
Fg = (2)
BH RgD 50

where Qw denotes in-channel water discharge, B denotes channel width and g denotes
the acceleration of gravity.

Flow regime is defined as follows. If Sf > 0.006 then the flow is assumed to be in
upper regime. For Sf < 0.006 the following parameters are computed;

−1 / 3
Fg ' = 1.74S f (3a)

u * ' D 50
δ= (3b)
11.6ν

where ν denotes the kinematic viscosity of water and

u * ' = gHS f (4)

where Sf is computed from (1b).

The lower limit of upper regime flow is given by the criterion

⎧ ⎛ D 50 ⎞ ⎛D ⎞ D
⎪− 0.02469 + 0.1517lg⎜ ⎟ + 0.8381lg 2 ⎜ 50 ⎟, for 50 < 2
⎛ Fg ⎞ ⎪ ⎝ δ ⎠ ⎝ δ ⎠ δ
⎜ ⎟ =
⎜F '⎟ ⎨
lg
⎝ g ⎠ ⎪ D
⎪ lg(1.25 ), for 50 ≥ 2
⎩ δ

(5a)

where lg is shorthand for the logarithm to the base 10. If the value of lg(Fg / Fg ' ) is
grreater than that given in (5a), the flow is in upper regime. The upper limit of lower
regime flow is given by the criterion

⎧ ⎛ D 50 ⎞ ⎛D ⎞ D
⎪− 0.2026 + 0.07026 lg⎜ ⎟ + 0.9330lg 2 ⎜ 50 ⎟, for 50 < 2
⎛ Fg ⎞ ⎪ ⎝ δ ⎠ ⎝ δ ⎠ δ
lg⎜ ⎟ = ⎨
⎜F '⎟
⎝ g ⎠ ⎪ D
⎪ lg(0.8 ), for 50 ≥ 2
⎩ δ
 (5b)

If the value of lg(Fg / Fg ' ) is less than that given by (5b) the flow is in lower regime.
If the value of lg(Fg / Fg ' ) is less than that given by (5a) but greater than that
given by (5b), the flow is in transition, and the flow regime is decided depending on
whether or not the flow is rising (in which case it is assumed to be in lower regime) or
falling ( in which case it is assumed to be in upper regime).

Floodplain resistance was computed using a Manning formulation, such that

5/3 1/ 2
Hf S ff Bf
Q wf = (6)
n

where Hf denotes floodplain depth, Sff denotes floodplain friction slope, Bf denotes
floodplain width and n denotes the Manning coefficient of resistance. The value of n was
inferred by means of calibration. The relation between floodplain friction slope Sff and in-
channel friction slope Sf is as follows;

S ff ∆x f = S f ∆x (7)

where ∆x denotes along-channel distance between nodes, and ∆xf denotes the
corresponding along-valley distance.

There is reason to believe that floodplain discharge does not exceed 10% of total
discharge. In OkFly a maximum of 30% of flow is allowed on the floodplain.

Computation of sediment transport In order to compute the sediment transport

rate, the sediment in motion is first divided into wash load and bed material load. In the
calculations of January, 1999, bed material load was considered to be material in excess
of 20 mµ and wash load was considered to be material finer than 20 mµ. In the
calculations of May, 1999 the cutoff size was raised to 43 mµ.

Wash load was assumed to continue downstream with no storage in the bed of
the river. As such, no sediment transport equation was used to treat wash load.
Deposition of wash load on the floodplain was, however, allowed, as described below.
For this reason wash load was non-conservative, such that the transport rate of wash
load decreased down the channe due to loss to the floodplain during floods.

The basis for the calculation of bed material load was the Brownlie (1981)
formulation. This relation computes total bed material load Qs as a function of a single
grain size D50, here taken to be the median size of the sediment in the active layer of the
bed. The relation takes the following form:

−0.3301
0.6601 ⎛ H ⎞
Q w (Fg − Fgo )
cF
Q s = 7.115 × 10 −3 1.978
Sf ⎜⎜ ⎟⎟ (8)
R +1 ⎝ D 50 ⎠

where
 ⎧c F = 1 for laboratory data
⎨ (9)
⎩ c F = 1.268 for field data
In equation (6), the parameter Fgo is given as

−0.1405 −0.1606
Fgo = 4.596 τ *o σg
0.5293
Sf
(10)

where

τ *o = 0.22 Y + 0.06e −17.73 Y (11)

and

Y= ( RR g )
−0 . 6
(12)

where grain Reynolds number Rg is defined as

3
gD 50
Rg = (13)
ν

An important aspect of the present analysis is an accounting for differential

transport of grain sizes in a sand mixture. Under pre-mine conditions, for example, the
bed of the Fly River showed a consistent tendency for downstream fining, from about 0.3
mm near D’Albertis Junction to about 0.1 mm near Obo. The mine-derived sand
likewise shows a wide range in grain sizes. A review of the literature revealed no
transport relations for sand that specifically address the issue of differential transport of
sand in mixtures. It was thus necessary to devise a new technique for the case at hand.

A first attempt to treat the problem was based on an adjustment of the parameter
Fgo in (6) so that it varied as a function of grain size in a mixture. This corresponds to
the introduction of a “hiding function” into the condition defining the threshold of motion.
This procedure was used in the January calculations. It became apparent that most of
the bed material transport takes place under conditions that are so far above the
threshold of motion that the adjustment produced very little differential transport.

With this in mind, the following procedure was implemented in the May
calculations. For simplicity all of the bed material load was assumed to be in
suspension. While this statement is unlikely to be completely accurate, under the high
flow conditions that transport the bulk of bed material load it can be expected that
suspended load dominates bedload by a factor of 3:1 or more. The material of the
surface (active) layer was divided into a number of grain size ranges, each with a
characteristic grain size Dj. A near-bed volume concentration of suspended sediment cbj
was assumed for this range, and the vertical distribution of suspended load in this range
was computed from the Rousean distribution as
 v sj
⎛ (1 − ζ ) / ζ ⎞ κu *
c j = c bj ⎜⎜ ⎟⎟
⎝ (1 − ζ b ) / ζ b ⎠ (14)

where cj denotes the volume concentration of the jth fraction a distance z above the bed,
ζ = z/H, ζb = zb/H where the near-bed elevation zb is here chosen as 5% of the depth, vsj
= the fall velocity associated with size Dj, κ is the Karman constant and u∗ is the shear
velocity of the flow.

Assuming the vertical distribution of streamwise flow velocity to follow a 1/6

power law, it can be shown after some work that the volume transport rate of the jth
grain size range Qsj is given by

v sj
⎛ (1 − ζ ) / ζ
1 ⎞ κu * 1 / 6
Q sj = α r c bjUHB∫ ⎜⎜ ⎟⎟ ζ dζ (15)
ς b (1 − ζ ) / ζ
⎝ b b ⎠

where U is the average flow velocity and αr is a coefficient associated with in-channel
resistance. As long as bedforms keep the surface layer well mixed, it is reasonable to
assume a simple linear relationship between the fraction Fj of the jth grain size in the
active layer and cbj, such that

c bj = α sFj (16)

where αs is another constant. Substituting (16) into (15), it is found that

v sj
⎛ (1 − ζ ) / ζ
1 ⎞ κu * 1 / 6
Q sj = α r α sUHBFj ∫ ⎜⎜ ⎟⎟ ζ dζ (17)
ς b (1 − ζ ) / ζ
⎝ b b ⎠

The constant αrαs can be evaluated by summing (17) and equating the result to that
predicted by (8). Further reducing, the fraction pj of bed material load in the jth size
range is found to be given by the relation

v sj
1 ⎛ (1 − ζ ) / ζ ⎞ κu *
Fj ∫ ⎜⎜ ⎟ ζ 1 / 6 dζ
Q sj ς b (1 − ζ ) / ζ ⎟
pj = = ⎝ b b ⎠
(18)
v sj
Qs
1 ⎛ (1 − ζ ) / ζ ⎞ κu *
∑j Fj ∫ςb ⎜⎜ (1 − ζ ) / ζ ⎟⎟ ζ 1/ 6 dζ
⎝ b b ⎠

Thus in the present formulation the total transport rate of bed material load is
given by (8) and the size distribution of the bed material load is given by (18). Sorting is
associated with the variation of fall velocity vs with grain size D; finer grains have a lower
fall velocity, a more uniform concentration distribution in the vertical and thus a higher
transport rate relative to content in the surface (active) layer.
 Sediment continuity The Exner formulation for sediment continuity for sediment
mixtures can be decomposed into the following two relations.

∂η 1 ∂Q S
B + = −q l (19a)
∂t 1 − λ p ∂x

∂ (L aFj ) ∂ (η − L a ) 1 ∂Q sp j
B + fI jB + = −q l p l j (19b)
∂t ∂t 1 − λ p ∂x

In the above relations x denotes in-channel distance downstream, t denotes time, η

denotes bed height, La denotes the thickness of the active layer, λp denotes the porosity
of the surface (active) layer of the bed, ql denotes the volume loss of sediment per unit
distance downstream per unit time to the floodplain and p l j denotes the fraction in the
jth size range of sediment lost to the floodplain. The parameter fIj denotes the grain size
fractions that transfer across the interface as the bed aggrades or degrades. When the
bed degrades, fractions fIj are set equal to the corresponding fractions in the substrate
stored immediately below the surface layer. When the bed aggrades, the fractions fIj are
set equal to a weighted average of the surface and load fractions. The following relation
based on the work of Toro-Escobar et al (1996) has been used for this purpose;

fI j = 0.7p j + 0.3Fj (20)

Equations (18a) and (18b) are solved as follows for the evolution of the bed. A
solution to the shallow water equations allows for an evaluation of the parameters
necessary to compute both the transport rate Qs and grain size fractions of the load pj
everywhere. Equation (19a) is then solved to determine the change in bed elevation at
each node. Then (19b) is solved to determine the change in the composition of the
surface (active) layer. The model used to compute the deposition of sediment on the
floodplain is given below.

Floodplain deposition It was incumbent upon the authors to devise a simple but
useful model for floodplain deposition, in so far as nothing else seems available in the
literature. According to the convective model of Parker et al. (1996), the rate of loss to
the floodplain ql is specified as

q l ∂η f Q
= = α f sf2 (21)
Bf ∂t Bf

where Qsf denotes the portion of the total in-channel sediment transport rate that rides
above the level of the floodplain, ηf denotes floodplain elevation and αf is a
dimensionless order-one coefficient that must be obtained by calibration. Here Qsf is
taken to include both wash load and bed material load that reside in the water column
above the level of the floodplain during floods, as illustrated in Figure 4. The above
calculation is applied at each time: when the river is not in flood Qsf vanishes and there is
no deposition on the floodplain.
 portion of water column that
contributes sediment overbank

channel
floodplain

Figure 4. Illustration of the model for floodplain deposition.

The grain size distribution of the sediment deposited overbank is taken to be that
of the grain size distribution of that portion of the suspended sediment above the
channel bed that is also above the elevation of the floodplain. That is,

v sj
1 ⎛ (1 − ζ ) / ζ ⎞ κu * 1 / 6
c bj ∫ ⎜⎜ ⎟⎟ ζ dζ
ς f (1 − ζ ) / ζ
pl j = ⎝ b b ⎠ (22)
v sj
1 ⎛ (1 − ζ ) / ζ ⎞ κu * 1 / 6
∑j c bj ∫ςf ⎜⎜ (1 − ζ ) / ζ ⎟⎟ ζ dζ
⎝ b b ⎠

where ζf = Hb/H and Hb denotes bankfull depth. Note that Hb must be less than H, i.e.
the river must be in flood, in order to render the calculation meaningful.

In evaluating (22) the wash load is included in the calculations. That is, for range
j within the range of washload, the Rousean distribution of suspended sediment that
when integrated yields the total washload transport rate is computed, and that portion
residing above the floodplain is used as the basis for the grain size distribution of
sediment transferred to the floodplain. In this way both bed material load and wash load
contribute to floodplain deposition. Although wash load does not interact with the bed, it
is not entirely conservative due to loss to the floodplain.

The coefficient αf in (21) was obtained by calibration. For Reach 3, calibration

was performed so as to yield a set amount of sediment deposited overbank by 1992, for
which rough observations were available. For Reaches 4 and 5, calibration was
performed so as to yield a loss of total load of about 3 – 4% to overbank deposition
between Kuambit and Everill Junction under pre-mine conditions (Dietrich et al, 1999).

Spatial distribution of overbank sediment The model outlined above provides an

estimate of the rate of deposition and grain size distribution of overbank sediment, and
also provides an indication as to how this rate and distribution varies downstream. It is
not adequate, however, to determine how the thickness of the deposit varies away from
the channel. The results of the calculation were thus further analyzed assuming that the
thickness of the deposit declines exponentially away from the channel. The coefficient in
this relation was provided by W. Dietrich (personal communication). As a result, the
model provides a gross indication of the spatial variation of thickness of overbank
deposits in both the streamwise and normal coordinates.

Flood frequencies in the absence of the mine Minor adjustments in local initial
floodplain and bed elevation were made in order to calibrate the model to “observed”
pre-mine flood frequencies at several key sites down the river. In several cases the
“observed” value was based on consensus rather than hard data. There is, however,
general agreement that even under pre-mine conditions the frequency of flooding in
Reaches 3, 3a, 4 and 5, i.e. those under consideration in the model, was rather high.
Some calibrated values are presented below.

Site Percent of time the river is in flood (pre-mine)

January model May model

Konkonda 12 9
Kiunga 16 26
Kuambit 29 42
Wygerin 51 67
Manda 53 78
Obo 53 59

Treatment of Reach 3a The sediment from the Ok Tedi depositing at D’Albertis

Junction has raised the bed sufficiently so as to create substantial backwater in the
upper Fly River (Reach 3a). This backwater has not only exacerbated flooding and
dieback, but also as created a long zone of deposition of naturally derived sediment in
the upper Fly. It proved impossible to model the observed bed aggradation at Kiunga
with a model terminating at Kiunga. With this in mind an extra reach has been artificially
added upstream of Kiunga. This reach has the same bed slope and cross-sectional
shape as the reach in the vicinity of Kiunga. The length of the reach was adjusted so as
to provide reasonable results bed aggradation at Kiunga.

Model zeroing As in the case of Okgrav5, the model was zeroed so as to

produce minimal change in bed elevation, flooding frequency and rate of overbank
deposition over the 70 year calculation period in the absence of the mine. The zeroing
was performed successfully except for one feature. It proved impossible to predict the
observed pre-mine rate of fining of characteristic bed grain size between D’Albertis
Junction and Everill Junction. Although the zeroed May model does predict substantial
downstream fining between these sites, further development would have been required
to reproduce the pre-mine data.

Testing of the model The model has been calibrated for a) floodplain Manning’s
n, b) frequency of floodplain inundation under pre-mine conditions and c) rate of
overbank deposition of sediment. The model has not been calibrated for bed
aggradation rate. The values predicted by the model at Konkonda, Kuambit and Kiunga,
however, agree reasonably well with available data.

6. OUTPUT OF MODEL OkFly

 The model predicts a slug of sediment slowly migrating down the Ok Tedi-fly
River system, with increased river bed elevation, flooding and overbank deposition of
sediment persisting throughout the 70-year period (from 1985 to 2055) modeled. The
height of the leading edge of the slug of aggrading bed sediment is in excess of 3 m in
the January model, but less than 1 m in the May model. There are two main reasons for
the decrease in the height of the slug: a) the change in the cutoff size for bed material
load from 20 mµ to 43 mµ and b) the adoption of a more effective sorting routine for
suspended sediment. Both of these features tend to increase the throughput of finer
sediment from the system and reduce the adverse effects of the mine as predicted by
the May model in comparison with the January model.

The January model thus may be described as “pessimistic” and the May model
as “optimistic.” The authors of this report believe the May model to be more likely to be
closer to reality. It should be pointed out, however, that even the “optimistic” calculations
predict a continuing very strong signal of the mine in terms of bed aggradation,
increased flooding and overbank sediment deposition.

As noted earlier, all the predictions of the model are provided on, Zip discs, one
pair for the January modeling and one for the May modeling. Complete documentation
of the contents of the discs can be found in the file [Link] in the root directory of
the first disc. The following Excel files, however, which can be found under directory
General of the first disc may be of particular interest to the user

[Link] This file provides numerical information about inundation

frequencies for all schemes for the entire 70-year calculation.

[Link] This file provides numerical information on node geometry, the

assumed hydrological input and the following information for a variety of sites along the
river, for all 70 years of the simulation and for each scheme: incremental bed elevation;
incremental mean floodplain elevation, cumulative sediment delivery, percent natural
sediment in delivery, cumulative floodplain deposit mass, percent natural sediment in the
floodplain, and cumulative floodplain mass as a percent of incoming sediment to a given
reach. In most cases the sediment is fractionated according to size.

[Link] This file is a computational tool to allow the fraction of a specified

part of the floodplain which receives an accumulation depth of x mm of sediment over y
years, where x and y are user-specified.

[Link] This file presents the output in graphical form.

[Link] This file presents the incremental bed elevations at every node for
every year of the calculation for every scheme.

Some model results are summarized below in Figures 5 – 12.

 Konkonda
9
Ok Tedi 13 km up D'A.J.
Null Kiunga
8
UpFly 8 km up D'A.J.
Kuambit
7
Incremental Bed Elevation (m)
Wygerin
6 Manda
Obo
5

-1
1985 1995 2005 2015 2025 2035 2045 2055
Year

a)Figure 5a. Incremental bed elevations, May model, Null Scheme

9 Konkonda
Scheme A
Ok Tedi 13 km up D'A.J.
8 Kiunga
UpFly 8 km up D'A.J.
7
Incremental Bed Elevation (m)

Kuambit
6 Wygerin
Manda
5 Obo

-1
1985 1995 2005 2015 2025 2035 2045 2055
Year

Figure 5b. Incremental bed elevations, May model, Scheme A

 Konkonda
9
Scheme DH Ok Tedi 13 km up D'A.J.
8 Kiunga
UpFly 8 km up D'A.J.
7 Kuambit
Incremental Bed Elevation (m)

Wygerin
6 Manda
Obo
5

-1
1985 1995 2005 2015 2025 2035 2045 2055
Year

Figure 5c. Incremental bed elevations, May model, Scheme DH.

1985 1989 1990 1992

6
1995 1999 2003 2008
2013 2018 2025 2055
5
Incremental Bed Elevation, Mid Fly (m)

Null
4

-1
0 50 100 150 200 250 300 350 400
Upstream Distance From Obo (km)

Figure 6a. Bed long profiles, May model, Null Scheme.

 1985 1989 1990 1992
4
1995 1999 2003 2008
Incremental Bed Elevation, Mid Fly (m) 3.5 2013 2018 2025 2055

3
Scheme A
2.5

1.5

0.5

-0.5

-1
0 50 100 150 200 250 300 350 400
Upstream Distance From Obo (km)

Figure 6b. Bed long profiles, May model, Scheme A.

1985 1989 1990 1992

4
1995 1999 2003 2008
3.5 2013 2018 2025 2055
Incremental Bed Elevation, Mid Fly (m)

3
Scheme DH
2.5

1.5

0.5

-0.5

-1
0 50 100 150 200 250 300 350 400
Upstream Distance From Obo (km)

Figure 6c. Bed long profiles, May model, Scheme DH.

 1985 1989 1990 1992
7
1995 1999 2003 2008

6 2013 2018 2025 2055

Incremental Bed Elevation, Mid Fly (m)

Null
5

-1
0 50 100 150 200 250 300 350 400
Upstream Distance From Obo (km)

Figure 7. Bed long profiles, January model, Null Scheme.

160
< 20 micron Null, Reach 4
20 - 62.5 micron
Cumulative Floodplain Deposit (Mt)

140
62.5 - 125 micron
120 > 125 micron

100

80

60

40

20

0
1985

1989

1993

1997

2001

2005

2009

2013

2017

2021

2025

2029

2033

2037

2041

2045

2049

2053

Year

Figure 8a. Cumulative floodplain deposit in Mt, May model, Reach 4, Null Scheme.
 100
< 20 micron Scheme A, Reach 4
Cumulative Floodplain Deposit (Mt) 90 20 - 62.5 micron
80 62.5 - 125 micron
> 125 micron
70

60

50

40

30

20

10

0
1985

1989

1993

1997

2001

2005

2009

2013

2017

2021

2025

2029

2033

2037

2041

2045

2049

2053
Year

Figure 8b. Cumulative floodplain deposit in Mt, May model, Reach 4, Scheme A.

140
< 20 micron Scheme DH, Reach 4
20 - 62.5 micron
Cumulative Floodplain Deposit (Mt)

120
62.5 - 125 micron
> 125 micron
100

80

60

40

20

0
1985

1989

1993

1997

2001

2005

2009

2013

2017

2021

2025

2029

2033

2037

2041

2045

2049

2053

Year

Figure 8c. Cumulative floodplain deposit in Mt, May model, Reach 4, Scheme DH.
 Flood Frequencies, Kuambit, May Model

120

100

80
Frequency

No Mine
Null
60
A
DH
40

20

0
1980 2000 2020 2040 2060
Year

Figure 9. Flood frequencies at Kuambit, May model.

Kuambit, May predictions, incremental bed elevation relative to that

predicted for 1999 for each scheme (no average)

1.5
Incremental bed elevation, m

0.5 No mine
0 Null
A
-0.5 B
C
-1
DL
-1.5 DH

-2

-2.5
1990 2000 2010 2020 2030 2040 2050 2060
Year

Figure 10. Incremental bed elevation above the predicted 1999 level
for all schemes, May model, Kuambit.
 Predicted TSS Values Based on Sediment > 62 Microns, Kuambit, May
Model

800.0

700.0

600.0
TSS in mg/l

500.0 No Mine
Null
400.0
A
300.0 DH

200.0

100.0

0.0
1980 2000 2020 2040 2060
Year

Figure 11. Predicted total suspended sediment (TSS) concentrations based on material
finer than 62.5 microns, Kuambit, May model.

Incremental Bed Elevations, May and January Models + Marshall Data

5
Jan Model
meters

4 May Model
3 Meas 1
Meas 2
2

-1
0 100 200 300 400 500
Year

Figure 12. Comparison of measured and predicted bed profiles in terms of incremental
elevations along the Ok Tedi-Fly River system from Konkonda (far right) to Obo (far left)
in 1998. Both the predictions of the January and May models are shown.
7. REFERENCES

Brownlie, W. R. 1981 Prediction of flow depth and sediment discharge in open

channels. Report No. KH-R-43A, W. M. Keck Laboratory of Hydraulics and
Water Resources, California Institute of Technology, Pasadena, California, USA,
232 p.

Dietrich, W., Day, G. and Parker, G. 1999 The Fly River, Papua New Guinea:
Inferences about river dynamics, floodplain sedimentation and fate of sediment.
In Varieties of Fluvial Form, Miller, A. J. and Gupta, A., eds, Wiley, 345 – 376.

Parker, G., Cui, Y., Imran, J. and Dietrich, W. 1996 Flooding in the lower Ok Tedi,
Papua New Guinea, due to the disposal of mine tailings and its amelioration.
Proceedings, International Seminar on Recent Trends of Floods and Their
Preventative Measures, Hokkaido Disaster Prevention Research Center.

Parker, G. 1990 Surface-based bedload transport relation for gravel rivers. J. Hydraulic
Res., IAHR, 28(4), 417-436.

Toro-Escobar, C. M., Parker, G. and Paola, C. 1996 Transfer function for the deposition
of poorly sorted gravel in response to streambed aggradation. J. Hydraulic Res.,
IAHR, 34(1), 35-53.