AD-A269 929

A Mathematical Model for

Oil Slick Transport and Mixing

in Rivers

Hung Tao Shen, Poojitha D. Yapa, De Sheng Wang and

Xiao Qing Yang August 1993

DTIC

ELECTE

SEP.29 W9. 3

-22571

1 1_iS

t4?
AbsMOact

The growing concern over the impacts of oil spills on aquatic environments has

led to the development of many computer models for simulating the transport

and spreading of oil slicks in surface waters. Almost all of these models were

developed for coastal environments. A few river models exist. These models only

considered the movement of surface oil slicks. In this study a two-layer model,

ROSS2, is developed for simulating oil spills in rivers. This model considers the

oil in the river to consist of a surface slick and suspended oil droplets entrained

over the depth of the flow. The oil transformation processes considered in the

model include advection, mechanical spreading, turbulent diffusion and mixing, evaporation, dissolution, emulsification, shoreline deposition and sinking.

The model can be used for simulating instantaneous or continuous spills either

on or under the water surface in rivers with or without an ice cover. The model

has been implemented for the Ohio-Monongahela-Allegheny river system and

the upper St. Lawrence River. This report describes the model formulation and

implementation. A case study is presented along with detailed explanations of

the program structure and its input and output. Although it is developed for

simulating oil spills, the model can be applied to spills of other hazardous

materials.

For convewiuon of Si metric units to U.S.JBrdtsh customary units of measurement

consult Standoad Practte for Use of tfe Intetnftonal System of Units (Sl, ASTM

Standord E380-89a, published by the American Society for Testing and Materlois, 1916 Race St., Philadelphia, Pa. 19103.

Thi report Is printed on poper thot contains a minimum of 50% recycled
Special Report 93-21

US Army Corps

of Engineers

Cold Regions Research &

Engineering Laboratory

A Mathematical Model for

Oil Slick Transport and Mixing

in Rivers

Hung Tao Shen, Poojitha D. Yapa, De Sheng Wang and

Xiao Qing Yang August 1993

PWParod for

UJ& ARMY COLD REG4ONS RESEARCH AND ENGINEERING LABORATORY

ST, LAWRENCE SEAWAY DEVELOPMENT CORPORATION

APru*WiorpWA~Ir~oe at4W dfto lb Is wr*Y*Wed
PREFACE

Thirs report was prepared by I long Tau Shen, FiPiesbw, Poojitha D. Yapa, Assistant

Professor; De Sheng Wang, Research Assistant; and Xiao Qing Yang, Visiting Scholar,

Department of Civil and Environmental Engineering, Clarkson University, Potsdam,

New York.

This study was supported in part by the U.S. Army Cold Regions Research and Engineering Laboratory (CRREL), Hanover, New Hampshire, under Contract No. DACA

89-SS-K0019, and the St. Lawrence Seaway Development Corporation, Department of

Transportation, under Contract No. DTSL55-89-C-C0549. Part of the work during the

final phase of this study was carried out at the Water Resources Engineering Department, LuleS University of Technology. Support provided by the Swedish Natural Science Research Council (NFR) and the LuleAt University of Technology for Dr. Shen and

Mr. Wang is appreciated.

The authors thank Steven F. Daly of CRREL, Robert Schmitt and Lewis Kwett of the

U.S. Army Engineer District, Pittsburgh, Stephen C. Hung of the St. Lawrence Seaway

Development Corporation, and Anders Sellgren and Annika Andersson of the LuleA

University of Technology for their cooperation and assistance.

This report is one of a series of reports on numerical simulation of oil slicks in inland waterways. The series coordinator is Steven F. Daly, CRREL.

iii Lii
CONTENTS

Page

Preface .....................................................................................................................................

Introduction ...........................................................................................................................

Transport and fate of spilled oil ...................................................................................... 1

O il spill m odels ................................................. 4

M odel form ulation ................................................................................................................. 6

Governing equations ....................................................................................................... 6

M odel form ulation ........................................................................................................... 8

Com puter m odel and case study ........................................................................................ 21

M odel im plem entation .................................................................................................... 22

M odel structure ................................................................................................................. 26

Case study ......................................................................................................................... 28

Input data files ....................................................................................................................... 36

D ata file creation ............................................................................................................... 37

Input adjustm ents ............................................................................................................ 51

Sam ple input data files .................................................................................................... 52

A dditional input .............................................................................................................. 57

M odel outputs ........................................................................................................................ 58

Velocity and depth distributions ................................................................................... 58

Sim ulation results ............................................................................................................ 58

Sam ple output files ............................................................................................................ 61

Conclusion .............................................................................................................................. 64

Literature cited .................................................................................... ................................. 64

A ppendix A : D efinition of FO RTRAN variables ..................................... ....................... 67

A ppendix B: List of program s ............................................................................................. 71

A bstract .................................................................................................................................... 73

ILLUSTRATIONS

Figure

1. Physical, chemical and biological processes affecting the oil slick transformation 2

2. O il slick transform ation in rivers .................................................................................. 2

3. Structure of the sim ulation m odel ............................................................................... 5

4. Method for determining the transverse distribution of the flow in a cross section 9

5. Stream-tube velocity distributions for a reach of the Ohio River ............................. 10

6. Velocity distribution in the grid system for a reach of the Ohio River .................... 11

7. Depth distribution of the channel for a reach of the Ohio River .............................. 12

8. Definition sketch for variables used to compute the slick's aspect ratio and

orientation ................................................................................................................... 16

9. Division of a nearly circular slick into pie-shaped segments .................................... 17

10. D ivision of an elongated slick into strips .................................................................... 17

11. Ohio-Monongahela-Allegherny river system ............................................................ 20

iii
Figure Page

12. Local and global coordinate systems ............................... . ....................... 22

13. Method for determining the land grid along a series of river bends. .......... 23

14. Schematization of a river reach with an island .......................................................... 24

15. Grid boxes and river boundary representation ......................................................... 24

16. Index system for locks and dams ............................................................................. 25

17. The computer model ROSS2 ...................................... 26

18. Subrou tin es ........................................................................................ ........... .. ......... 29

19. Air temperature variation during the simulation period, Pittsburgh,

Pennsylvania ............................................................................................................... 32

20. Simulation results showing the variation of concentration distribution............. 33

21. Comparison of simulated and observed leading edge positions ............... 35

22. Comparison of simulated and observed concentration peak locations ................... 35

23. Effect of time step At on leading edge positions ........................................................ 36

24. Definition sketch of the river numbering system ...................................................... 37

25. D efiring ice regions ........................................................................................................ 41

26. Simulation of flow direction in stream tubes .............................................................. 51

TABLES

Table

1. Spreading law for oil slicks ........................................................................................... 15

2. Spread rates of oil slicks and phase transition times ................................................. 16

3. Shoreline descriptor and default parameters ............................................................. 18

4. Suggested evaporation parameters for various petroleum fractions ....................... 19

5. Suggested evaporation parameters for various crude oils ...................................... 19

6. Dissolution coefficients at 250C .................................................................................... 21

7. Relationships between local and global coordinate systems .................................... 22

8. Boundary types of locks and dams ............................................................................... 25

9. Cross references of subroutines .................................................................................... 27

10. Definition of variables in the common block .............................................................. 28

11. Observed and simulated locations of the leading edge and the concentration peak 32

12. Comparison of computing time .................................................................................... 36

Accession 'or

NTIS S R'

DTIC Tj, B 0

Unarnnounced 0

JuatifIcatlon

V .str-bution/ TIC ,

RAvalabillty _odeg ""D

il and/or

blot ISpeoial

iv
A Mathematical Model for

Oil Slick Transport and Mixing in Rivers

HUNG TAO SHEN, POOJITHA D. YAPA, DE SHENG WANG AND XIAO QING YANG

INTRODUCTION

In recent years there has been a growing concern over the increasing contamination of waterways

and shoreline areas caused by oil spills. Oil spills in inland waterways can have enormous environmental and economic impacts. Oil pollution in rivers can not only cause long-term damage to the aquatic

environment for fish and wildlife but also pose threats to water supplies for areas along the river. With

heavy industrial development in areas along rivers and inland navigation activities, major rivers in the

United States are vulnerable to oil spills. For example, in January 1988, an oil storage tank at West

Elizabeth, Pennsylvania, collapsed, spilling over 17,000 barrels of diesel oil into the Monongahela River. The oil slick reached the Ohio River near Pittsburgh, which is 25 miles downstream of the spill site,

within a day and drifted farther downstream, affecting water supply intakes along the river (Miklaucic

and Saseen 1989).

Concern over the adverse impacts of oil spills has led government agencies and private industries to

develop oil spill emergency response plans (e.g. Hung 1991). An important element in these programs

is the use of computer models to predict the movement and possible impact of an oil spill. In the event

of a real oil spill, a model can be used on a real-time basis to assist the containment and recovery of the

oil, and it can also help to guide the field data collection for detailed environmental impact analysis.

These objectives can be achieved by using the model to forecast the location and distribution of the oil.

A computer model can also be used to study scenarios of possible spills to assist in developing contingency plans and assessing likely environmental impacts.

Numerous oil spill models have been developed during the last two decades, as reviewed by Stolzenbach et al. (1977), Huang and Monastero (1982) and Spaulding (1988). Almost all of these models were

developed for coastal marine environments. In a recent study Shen and Yapa (1988) developed a model, ROSS, for surface oil slick transport in rivers. Oil slick transformation processes considered in the

model include advection, mechanical spreading, horizontal turbulent diffusion, evaporation, dissolution and shoreline deposition. In this study, a two-layer model, ROSS2, is developed to take into consideration the dispersion of oil droplets into the water column, as well as related processes including

subsurface transport, mass exchanges between the surface slick and the suspended oil droplets in the

water column, and deposition of oil to the channel bottom due to sedimentation.

Transport and fate of spilled oil

In addition to the location, size and physical-chemical properties of the spill, the fate of an oil spill is

affected by complex interrelated transport and weathering processes (Mackay and McAuliffe 1988).
SpinlSpreadingd

Oil-in-WaterMirba

Dissolution Emulsion Dgdion * Degradation

Subsurface Spill

Current ;:;No Sedimentation

Sinking

. . . BiodeoagdatlonBottom Sediment

II I -, .

Hours Days Weeks Months YearsýTime From Onset of Spill

Figure 1. Physical, chemical and biological processes affecting the oil slick transformation.

Figure 1 shows a schematic representation of the fate of oil spilled in a surface water body. Oil may

enter a water body either through a surface spill or an underwater discharge. A part of the spilled oil

will form a surface slick on the water surface and can be moved about by the action of winds, currents

or waves. Some low-molecular-weight hydrocarbons can dissolve into the underlying water column,

but most of these are lost to the atmosphere through evaporation. Evaporation of volatile components

will reduce the volume of spilled oil by as much as 50% during the first few days after the spill (Butler

et al. 1976). In turbulent waters some of the oil is emulsified and dispersed into the water column as

suspended droplets. Some of the oil droplets in the suspension may become attached to suspended

particulate matter and slowly settle to the bottom; some will rise to the water surface due to buoyancy

to form a water-in-oil emulsion. Current and waves may also drive oil onto beaches or riverbanks. For

spills in rivers the transport of oil is also affected by hydraulic structures, such as locks and dams, as

shown in Figure 2. In addition, photochemical reactions and microbial biodegradation can change the

character of the oil and reduce the amount of oil present.

Evaporation

AdvectionSpreading

"tissoluttion

Oil/Water i 

2 Emulsion currnt "I"an • °*° i\ .?••=Emulsion 0

- D ispersion

SDam Sink 2

S~Mixings(oe'

Figure 2. CV/slick transfrmio in rivm.

2
Advection

Advection is a physical process that involves drifting of the surface oil slick and the subsurface oil. It

is the main mechanism that governs the location of the oil following its discharge into the river. Once

spilled into the river, the oil will be transported downstream by the flow, with the influence of wind.

The advection of surface oil is caused by the combined effects of surface current and wind drag. The

advection of subsurface oil is the movement of suspended oil droplets entrained in the flow due to the

subsurface current.

Spreading

Spreading is the physical process that involves drifting of the surface oil slick and the subsurface oil

as it is being transported in the river. The spreading phenomena include mechanical spreading and

horizontal diffusion. Mechanical spreading is the horizontal spreading of the surface oil slick due to the

balancing forces of inertia, gravity, viscosity and surface tension. This mechanism terminates when the

thickness of the slick is reduced and the thin oil slick ruptures into small patches. Horizontal diffusion

is the spreading of oil due to the turbulent fluctuations of wind and current velocities. This is the main

mechanism that causes the spreading of both the surface and subsurface oils. Since the spreading of oil

enhances other weather processes, such as evaporation, dissolution and emulsification, it is one of the

most important processes affecting the fate of the spilled oil.

Evaporation

Evaporation occurs immediately after the spill. As the surface slick spreads, more of the hydrocarbons are exposed to the atmosphere, causing the evaporation rate to increase. The amount and rate of

evaporation depend on the percentage of light, or volatile, components in oil. Evaporation is the most

significant physical-chemical process causing the reduction in oil volume. Highly refined oil can lose

75% or more of its volume through evaporation within a matter of days.

Dissolution

Some of the same oil components that are subjected to evaporation can also dissolve into the water

column from a surface slick. Only the low-molecular-weight hydrocarbons have an appreciable solubility in water. The fraction of oil dissolved is very small compared to the oil that is evaporated. However, this extraction process can be important because of the toxicity of the dissolved fraction.

Emulsification

Dispersion, or oil-in-water emulsification, is the mixing of surface oil into the water column. Breaking waves and other surface turbulence can mix the surface oil into the water column by forming many

small globules of oil that can be rapidly dispersed over the depth of the river flow and transported by

the current. Some of the fine oil droplets in the water column may dissolve in the water colunur, at•ch

to solid particles and biodegrade. A significant portion of suspended oil droplets will rise to the water

surface when the buoyancy force is large enough to overcome the vertical mixing. This resurfacing

mechanism can increase the area of oil on the water surface and the duration of oil passage at any

particular site along the river.

Water-in-oil emulsions can also be formed, particularly with heavy crudes and residual oil. The resulting emulsions contain a large percentage of water but have a semisolid texture, often referred to as

"mousse" because of their appearance. The formation of mousse is a process of dispersion of emulsified water into the oil. This process will increase the viscosity and volume of the slick. Stable emulsions

can have a water content of up to 809/9. The mechanism of water-in-oil emulsification is not dear, alOtough it is believed that this is due to the presence of asphaltenes, waxes and surfactants and the

turblence associated with dispersion and resurfacing processes.

3
Photo-oxidation

In the presence of atmospheric oxygen, natural sunlight has sufficient energy to change the composition of the oil. This photo-oxidation process is very slow and usually becomes important a few days

after the spill. The extent and rate of photo-oxidation depends primarily on the chemical composition

and optical density of the oil. The photochemical reactions can change the interfacial properties of the

oil and may dissolve toxic organic species into the water column. However, its effect on the oil slick

transport process is not very important.

Shoreline deposition and sedimentation

Oil may be brought to the riverbanks and deposited along the shoreline, to be later re-entrained into

the river current. This process can significantly affect the distribution of oil and should be modeled

Some of the suspended oil droplets may also sink to the riverbed. This sinking or sedimentation process occurs due to an increase in density of the oil, resulting from either the evaporation and dissolution of lighter fractions of the oil, or adherence of the oil droplets onto suspended sediment. The oil

deposited on the channel bottom may be moved laterally or resuspended, or it may undergo further

biological or physical-chemical reaction. Little is known about the ultimate fate of the sedimented oil.

Biodegradation

All of the processes just described, except possibly photo-oxidation, can only redistribute the oil.

They cannot remove the hydrocarbon from the environment. Real degradation takes place only through

biochemical oxidation. This biodegradation process, which may continue for years after a spill, is the

principal long-term means of removing the spilled oil from the environment.

Oil spill models

Most of the existing oil spill models simulate only the advection and spreading of a surface slick.

Other models deal extensively with physical-chemical processes but lack the component for simulating the movement of the slick. Only in recent models have the incorporation of both transport and

weathering processes been attempted (Huang and Monastero 1982, Spaulding 1988). Since there is a

significant lack of data for a reliable analytical formulation to be established for many of the weathering

processes, it is impractical to include all of them into an oil spill simulation model. It would be more

useful to include the most significant processes, i.e. those accounting for the bulk of the oil, while omitting others so that uncertainty in the outcome can be reduced. Almost all of the existing oil slick models

were developed for coastal marine environments.

Only a few models were developed for inland waters (Huang and Monastero 1982). Tsahalis (1979b)

developed a simulation model for predicting the transport, spreading and associated shoreline contamination of a surface oil slick in a river. In his model the current velocity distribution in the river is

calculated by empirical relationships determined from field data with some modifications for the secondary current in river bends. The oil slick is assumed to remain circular, w'ith its radius calculi ted

according to Fay's spreading laws (Fay 1969,1971). The drift velocity of the slick center is determined

by formulas derived by Tsahalis (1979a) from laboratory experiments. Fingas and Sydor (1980) developed a two-dimensional model for surface oil slicks in a short river reach. The current distribution was

determined using the finite-difference scheme of Leendertse (1970). The entire oil slick volume is represented by a large number of individual parcels. The drift velocities of these parcels are determined by

the wind factor approach. A random fluctuation component is included to represent horizontal diffusion. The area of each oil parcel is calculated by Fay's spreading laws for circular slicks.

Recently Shen and Yapa (1988) developed a computer model, ROSS, for surface slicks in rivers. The

model considered the effect of ice covers. In this model the oil slick is considered to be composed of a

large number of discrete parcels with equal volume. These parcels are tracked for their positions and

aibjected to an equal mass loss rate during the simulation. The two-dimensional velocity distribution

of the river water is computed using a separate hydrodynamic model. The advection of each parcel in

4
INPUT

"* Model parameters: Evaporation. ditsolution, dispersion, resurfacing. sedimeitation

"* Oil spill information: Spill k,' -- on and duration, volume and piperties of oil

"* River and weather darn "-,.ometry. wind velocity, discharge and water levels, air

temperature, ice 'ýon, on, river coordinate system transformation

DETERMINE DEPTH AND VELOCITY DISTRIBUTIONS

"• Two-dimensional depth distribution "* Two-dimensional depth-averaged velocity distribution

ADVECTION AND HORIZONTAL DIFFUSION

"• Surface layer

" Suspended layer

EMULSIFICATION, RESURFACING AND SEDIMENTATION

"• Some surface oil is emulsified into water "• Some suspended oil rises to surface and some sinks to the bottom

EVAPORATION AND DISSOLUTION

(only for surface layer)

SHORELINE DEPOSITION

(only for surface layer)

1,

MECHANICAL SPREADING

(only for surface layer)

OUTPUT

"* Volume of oil parcels and distributions in surface and suspended layers

"* Amount and location of deposition on shorelines "* Amounts of evaporation, dissolution, emulsification, resurfacing and sedimentation

Figure 3. Structure of the simulation model

the slick is determined from the wind velocity and the computed current velocity, using the wind factor

approach. The spreading of the slick is simulated by considering both mechanical spreading and horizontal turbulent diffusion. The model developed by Fay (1969, 1971) is used to simulate mechanical

spreading, while a random-walk simulation is used to simulate spreading due to surface turbulence.

The loss of oil due to the shoreline deposition is calculated according to the oil-retention capability of

the shorelines. The losses of oil due to evaporation and dissolution are calculated based on empirical

formulations that consider the effects of slick area, wind velocity, temperature and oil properties. None

of the above models considered the motion of oil in suspension and the exchange between the sirface

oil and the suspended oiL Sedimentatiom of oil to the bottom was also not considered.

In the present study the model ROS is generalized by considering the movement of oil both on and

under the water surface. The present model, ROSS2, is primarily designed for the transport of oil spilled

in rivers; however, it can also be used for spills of other hazardous materials. A brief outline of the

structure of the model is presented in Figure 3.

5
In ROSS2 the spilled oil is divided into a surface layer and a suspended layer, with continuous

exchanges of oil mass between the two layers. In the simulation, oil in both layers is discretized into a

large number of parcels that are tracked for their positions and volumes at each time step during the

simulation. Besides the two-dimensional velocity distribution, the distribution of flow depth is also

determined for simulating transport in the suspended layer. Because of the exchange between the surface layer and the suspended layer, the number of oil parcels, as well as the volume of each parcel, will

change from time to time. The advection of oil parcels in the suspended layer is determined by the

depth-averaged current velocity, and those on the surface are determined by the wind factor approach.

The model developed by Fay (1969, 1971) is used to simulate mechanical spreading, and a randomwalk simulation is used to account for horizontal spreading due to turbulent diffusion. The loss of oil

in the surface layer due to evaporation and dissolution is calculated based on analytical and empirical

formulations, which consider the effects of slick area, wind velocity, temperature and oil properties.

The loss of oil due to shoreline deposition is calculated according to the oil-retention capability of the

shorelines. The amounts of vertical dispersion into the water column, resurfacing and deposition to the

channel bottom are calculated based on the depth-averaged oil concentration in the suspended layer.

Weathering processes that occur long after the onset of the spill are not well understood and less significant, so they are not considered in this model. This is also justified from the operational point of view,

since the main objective of the model is for short-term tactical forecasts to assist in designing response

measures in the event of a major spill.

The present simulation model has been applied to the upper St. Lawrence River and the OhioMonongahela-Allegheny river system. This report will use the Ohio-Monongahela-Allegheny river

system as an example to illustrate the model implementation in detail. The Ohio-Monongahela-Allegheny river system is more complex and difficult to implement, so it serves as a better example.

MODEL FORMULATION

Governing equations

In the present model, oil in the river is considered to consist of a surface slick coating on the water

surface and a mixed layer containing suspended oil droplets extending over the depth of the flow, with

continuous exchange betw een the two layers. The thickness of the surface layer is assumed to be negligible in comparison to the suspended layer. The equation of motion for oil in the surface layer can be

written as:

ac + a (UC)+ a ( aS (D a1) + a tD ay

at ax 5;TX ax OY a

+ ,E1 Vb Cv] z =o -'Y Cs - Ca SE + Ms (x, y) - Ds (x, y) (1)

where x, y, t = space and time variables

z = vertical coordinate measured downward from the water surface

C, = local volumetric oil concentration in the surface layer per unit surface area

c, = local volumetric concentration of oil in the suspended layer per unit volume of water

u, v, = components of surface drift velocity in x and y directions, respectively

Dx, DY, = diffusion coefficients in the x and y directions, respectively

cy, = coefficient representing the probability of deposition of an oil droplet reaching the

water surface

Vb = buoyant velocity of suspended oil globules

y = coefficient describing the rate at which the surface oil is dispersed into the water

column

Ca = area concentration of oil in the surface layer

6
SE = rate of evaporation and dissolution per unit area of the surface slick

M, = effect on the distribution of surface oil by mechanical spreading

Ds = effect on the distribution of surface oil by shoreline deposition.

The concentration distribution in the suspended layer can be described by:

a( h + 5 (u C ,h) + d (vC ,.h)

aýxhD,"''ý) ý (hiDyajV- n1ý Vbcll.I = + 7tC_" -PfICA:=,, Q)

where C, = depth-averaged volumetric concentration of oil in the suspended laver

h = flow depth

u, v = components of depth-averaged river current in x and t! directions, respectivekl

P•1 = coefficient defining the rate of net oil deposition on the riverbed per unit area.

Equation 2 can be simplified by introducing the mass conservation equation of water:

a3h + I (uh ) + a (vh )0=. (3)

at ax Wy

In Lagrangian form the simplified equation for the suspended layer becomes:

Dt hi Lax k ax ay~ 'CT l OCv~ - 1--[ •)+• (h Dx•cYh ay•a(h yfClh_ •Wbv+ C -s•vC (4)

where QD= L+u a +v ý . Dt at ax ay

zCVC = Clcv(• =0;

PC, (=3c~Z=h)IA.

The movement of spilled oil in the river as described by eq I and 4 is mainly governed by advection

and diffusion. These equations can be solved when the flow condition is known. Finite-difference and

finite-element methods for solving advection-diffusion equations exist. These methods, however, often suffer from numerical diffusion problems and aemand excessive computational time and memory

when applied to long river reaches. In this study a modified Lagrangian discrete-parcel method is used.

The method is a generalization of the method used in the ROSS model (Shen et al. 1990).

Lagrangian discrete-parcel algorithm

In the Lagrangian discrete-parcel algorithm, both the oil on the water surface and the oil in the suspension are represented as ensembles of a large number of small parcels. Each parcel has a set of timedependent spatial coordinates and a mass associated with it. The movement of each parcel in fle river

is affected by the wind, the water current and the concentration of surrounding parcels. During each

time step, all the oil parcels are first displaced according to the current velocity and a fluctuation component at their respective locations. The turbulent fluctuation component, which represents the effect

of horizontal diffusion, is related to the diffusion coefficient based on the ranidom-walk analogy. After

all the oil parcels are displaced according to advection and diffusion, further modifications are introduced to account for all other transport and weathering processes. Due to these processes, especially

"the exchanges between the two layers, the mass of each oil parcel will be different and will change from

7
time to time. The number of parcels will also increase with time during the simulation. If a large number of parcels are released in the river, and their discrete paths and masses are followed and recorded

as functions of time relative to a reference grid system in space, then the density distribution of t

particles in the water body can be interpreted to give the concentration ot the oil.

The approach requires an efficient bookkeeping procedure rather than the solution of a large matrix

associated with a conventional Eulerian finite-difference or finite-element method. The algorithm is

inherently stable with respect to time, although the time step should be compatible with the grid size

and velocity for numerical accuracy. Since the movement of each parcel in the nver depends on the

distribution of the entire ensemble on the water surface and in the suspension, all parcels must be

traced at each time level before proceeding to the next.

The detailed structure and implementation of the present numerical model will be discussed and

presented later. In the following sections the analytical formulations used for each component of the

model will be discussed.

Model formulation

Current velocity

Since the water current affects advection, spreading and the exchange of oil between the two layers,

the distribution of the magnitude and direction of the current and the flow depth must be determined

first. Numerous numerical methods exist in the literature for determining the two-dimensional flow

distribution in shallow waters (Leendretse 1970, Hamilton et al. 1982). However, this type of method is

time consuming and impractical for long rivers. A quasi-two-dimensional stream-tube method is used

in the present study. In the present approach the time-dependent discharge distribution Q(xt) along

the river is first obtained through the use of a one-dimensional hydraulic transient model, which was

developed based on the St. Venant equations:

aQ +aA =0 (5)

ax at

S( Q) 2aA + 2Q a +gAaH -gA(So-Sf)=O (6)

at ( ax A ax (a)

where x longitudinal distance along the river

A = flow cross-sectional area

H = water level

S,= bed slope

Sf = friction slope.

The friction slope can be calculated by Manning's equation:

sf = (7)

2.21 A 2R 4/3

in which R is the hydraulic radius and nb is the Manning's roughness coefficient of the bed. For an icecovered reach, the composite Manning's coefficient, which accounts for the resistance of ice cover and

the riverbed, should be used instead of n. The composite Manning's coefficient can be calculated by

the Belokon--Sabaneev formula (Larsen 1969):

[ný +nb (8) L 2 J

The hydraulic radius can be assumed to be half of the flow depth for an ice-covered reach in wide

rivers.

8
Once the one-dimensional solution is obtained, the discharge Q can then be distributed across the

width of the river using a stream-tube model (Shen and Ackermann 1980) to give the two-dimensional

velocity distribution at selected cross sections.

For a channel cross section, as shown in Figure 4, the transverse distribution of the flow can be determined using the stream-tube method. The cross section is first discretized into trapezoidal elements. By

applying Manning's equation, the ratio of discharges between the entire cross section and a partial

cross section can be written as:

Qp _ =1 Ap Rp/3 (9)

Q N A R 2/3

where P = number of trapezoids in the partial cross section

N = total number of trapezoids describing the entire cross section

Qp = cumulative discharge up to and including the Pth trapezoid

Q = total discharge through the entire cross section

A;, Rp= area and hydraulic radius of the pth trapezoid, respectively.

The cumulative discharge Qp can be computed by first rewriting eq 9 as:

P

Qp FQ Y, ApR 2/3 (10)

p=1

with

F Q= (11) ,N 2/3

Yp=IA"pR'p

and

QP - QP-i = FQ ApRp2/ý (12)

Based on the computed distribution of Qp, stream-tube boundaries within the cross section can be

determined by simple interpolation. Once the stream-tube boundaries are located, the flow through

each stream tube is then divided by the cross-sectional area of the stream tube to obtain the depthaveraged velocity, and the area of the stream tube is divided by its width to obtain the width-averaged

depth. The velocity and depth are then assigned to the center of the stream tube. By applying this

procedure to successive cross sections along the river, a two-dimensional depth-averaged velocity and

Op-i:

a. Cumulative discharges auw a b. Area and discharge for the Pth trapezoid.

partial cross section.

Figure 4. Method for determining the transverse distribution oftheflow in a cross section.

9
Ohio River 0BonI '

Figure 5. Stream-tube velocity distributions for a reach of the Ohio River.

width-averaged depth distributions can be obtained. The direction of a velocity vector is the same as

the vector connecting center points of a stream tube at successive cross sections in a relatively straight

reach, or perpendicular to the cross section in a river bend. As an example, the simulated depth-averaged velocity distribution for a reach of the Ohio River is shown in Figure 5.

Once the velocity and depth distributions along stream tubes are established, the distribution of

velocity for all the points in a predefined grid system can be obtained through linear interpolations.

Samples of calculated velocity and depth distributions are shown in Figures 6 and 7. The validity of the

stream-tube method has been examined by Shen and Ackermann (1980).

Advection

During each time step At of the oil spill simulation, the displacement AS of an oil parcel is calculated by numerically integrating the drift velocity Vt over the time period At. When the At value is

large, subintervals 8tk are needed for numerical accuracy. In this case the displacement during the time

interval At is

- K.

AS Y, Vk 8tk (13)

k=1

where Vk = drift velocity of the oil parcel during the time interval Stk

AS = displacement during the time interval At

K

I 8tk = At.

k=l

The values of Stk should satisfy the condition (Roache 1972, Cheng et al. 1984)

8t < u + ._. -k (14)

where Uk and zk are the x and y components of the velocity.

Open uater condition. In the present model the advective velocity of each oil parcel is computed as:

Vt ft + (15)

10
OHIO RIVER

Velocity (It/s)

EJ1.05 -1.37

~J1.38 -1.69 I

S170 -2.01

M2.02 -2.33

IM 2.34 - 2.64

GrOfn 0 Weirton

Figure 6. Velocity distribution in the grid system for a reach of the Ohio River.

where Vt is the drift velocity of an oil parcel and V and V' are the mean and turbulent fluctuation

components of the drift velocity. The component V'is included to simulate the horizontal diffusion of

the oi parcels. The formulation for V' will be discuslsed later. In the suspended layer the mean drift

velocity is equal to the depth-averaged current velocity Vc, determined from the stream-tube simulationt. On the surface layer the mean drift velocity is conisidered to be a vector sumn of a wind-induced

and water-current-induced drift (Stoizenbach et A. 1977). The mean drift velocity V. on the surface can

be calculated by the formula:

VS. = 40ewVW+ qCVC (16)

where Vw, wind velocity at 10 m above the water surface
59

Depth (It) r.. 0.00- 5-00

0~~ Weito 501 - 10.00

10 2,01 - 25.00

m 25.01 - 27.00

27.01 - 30.00

M 30.01 - 34.88

Figure 7. Depth distribution of the channel for a reach of the Ohio River.

Vc = depth-averaged current velocity

(xw = wind drift factor to account for the drift of the surface slick due to wind

ac = correction factor for the current on the surface.

A value of 0.03 has been commonly used for ua (IHuang and Monastero 1982). The wind speed on the

river is usually different from that at a nearby meteorological station (Weeks and Dingman 1972). A

correction factor should be used when the meteorological station is situated away from the river.

Assuming a logarithmic velocity profile, the surface velocity of the water current can be related to

the depth-averaged current velocity by the relationship:

ac == 1 j+ui (17) Vc KVc

where V0 = surface velocity

u. = shear velocity

ic = Karman constant (0.4).

The ratio ;ý usually vares between 1.1 and 12.

ke-covera d cndition. Oil spill movement under ice is a topic that has received little theoretical or

laboratory treatmemt. For the suspended layer the transport of oil droplets can be treated by the advective diffuskm formulation, as described in eq 4. The advection velocity is equal to the depth-averaged

current velocity. The advection of surface oil on the underside of the ice cover is, however, more compfimbid. As s e by the wexrimental study of Cox and Schultz (1981), an ice cover may be classified into euue categories for slick advection. The ice cover may be smooth, contain small roughness

eenumts, or contwn large roughness elements. The following discussion summarizes the results of

ths lmt ed experlnental studies and serves as the basis of the present simulation.

12
Under smooth ice covers with no current, the oil will rest at an equilibrium thickness that is described by the empirical equation

Seq = 1.67- 8.5(Apw) (18)

where 8eq = static equilibrium slick thickness (cm)

A = (pw - p0 )/pw, the relative density difference between oil and water

pw and p0 = water and oil densities, respectively (g/cm3

).

An ice cover is considered to be smooth when the height of the ice roughness is smaller than the equilibrium thickness of the oil, 8q.

The current velocity at which the oil begins to move along an ice cover is called threshold velocity,

Uth. For a smooth cover the value of Uth was empirically determined to be a function of the oil viscosity

pO and is given as

U th = 305.79/(88.68 -,o) (19)

where Uth has the units of cm/s and A. has units of g/cm-s. Viscosities for crude and fuel oils fall in the

range of 5-50 centipoises (1 cp = 1.0 x 10-2 g/cm-s = 2.4 lb/ft-hr).

A rough ice cover has the ability to retain the oil between the roughness elements. As the current

velocity is increased, the oil will creep along the upstream face of the roughness element until it spills

over the element and moves downstream. The threshold current velocity at which the oil will move

downstream under a rough ice cover is called the failure velocity, Ufi:

UH = 1.5(2 (P'Pw) [,o/wg (p. + p.)31/2)1/2 (20)

where %Yw is the oil-water interfacial tension. The failure velocity Ufl is the velocity above which no oil

can be contained upstream of a large roughness element.

If the depth-averaged current velocity is larger than the threshold velocity, Uth or Ufl, the relationship between the current velocity and the slick velocity is given as:

= K (21) SVCJ 0.11.5 F2e +1.10.5

with

F8 = (22)

(A g~e)l /

2

where V = mean oil drift velocity

c = u- ent speed

K = friction amplification factor

F5 = slick densimetric Froude number

g = gravitational acceleration.

The value of K is a function of the roughness height of the cover. With the limited data available, K is

assumed to vary linearly between 1.0 for a smooth cover and 2.6 for an ice cover with Manning's roughries coefficient n1 = 0.055.

13
Horizontal diffusion

The term V' in eq 15 accounts for the horizontal turbulent diffusion due to the turbulent fluctuation

of the drift velocity. Based on the random-walk analysis (Fischer et al. 1979):

V"= (4Dw/St )1/2 (23)

and

V'= V"Rnei' (24)

where 5t = time step

DT = diffusion coefficient

Rn = normally distributed random number with a mean value of 0 and a standard deviation

of 1.

The directional angle 0' is assumed to be a uniformly distributed random angle ranging between 0 and cr.

In rivers the diffusion coefficient is affected by the shear velocity u. and the depth of flow h (Fischer

et al. 1979):

DT= 0.6h U-. (25)

Flume experiments by Sayre and Chang (1969) indicated that DT for surface dispersants can also be

calculated by eq 25. Using Manning's equation, we can rewrite eq 25 as:

DT = 0.40nbVgl/ 2h5 /

6 (26)

for free surface flow (ft2 /s) and

DT = 0.64 iiVg12 h5 /

6 (27)

for ice-covered flows (ft2 /s).

Mechanical spreading

The mechanical spreading process only affects the spreading of surface slicks. Formulations used in

the model will be discussed in this section.

Spreading in open water. A number of theories have been proposed for the process of mechanical

spreading in open waters (Blokker 1964, Fay 1969,1971, Hoult and Suchon 1970, Mackay et al. 1980b).

In this study Fay's spreading theory (1971) is used because this theory is based on a rather comprehensive description of the spreading mechanism and has been verified by laboratory experiments (Fay

1971, Hoult and Suchon 1970) and other analytical solutions (Fannelop and Waldman 1972).

Fay's spreading theory is derived for single-component, constant-volume slicks with idealized configurations in quiescent water. This theory considers the spreading of oil as a result of two driving

forces--gravity and surface tension-counterbalanced by inertia and viscous forces. The spreading of

an oil slick is considered to pass through three phases. In the beginning phase, only gravity and inertia

forces are important. In the intermediate phase, gravity and viscous forces dominate. The final phase is

governed by the balance between surface tension and viscous forces.

Formulas for both one-dimensional spreading and radial spreading at different stages are summarized in Table 1. The spreading rate during each phase can be obtained by taking time derivatives of the

formulas given in Table 1. The equations for spreading rates are summarized in Table 2. The time rate

of chand of the oil volume in these equations represents the change due to weathering, including

14
Table 1. Spreading law for oil slicks. (After Fay 1971, Hoult 1979, Waldman et al. 1973.)

Spreading Width Radius

phase L, R,

Gravity-inertia 1.39 (AgAt 2

)1/3 1.14 (AgVt 2)1/4

Gravity-viscous 1.39 (AgA 2 t

3 /2 V-1/2)1/4 0.98 (AgV 2 t 3/2 V-1/2)1/6

Surface tension-viscous 1.43 (o;2 t

3 p2 v-1)1/ 4 1.60 (c;2 t

3 p-2 v-1)1,' 4

A 0.5 volume of the oil per unit length of oil slick.

v kinematic viscosity of water.

evaporation and vertical mixing, and the changes in oil volume distribution in various parts of the

slick. hi addition to the rates of spreading at different phases, the times at which phase transitios occtr

also need to be determined. These transition times can be obtained by letting equations in the appropriate phases equal each other and solving for the time. Equations for the transition times are summarized

in Table 2.

Fay (1969, 1971) observed that the changes in slick properties caused by weathering may result in

the eventual cessation of mechanical spreading. Based on a number of fields observations, Fay suggested that

Af =105 V31 4 (28)

where Af is the final slick area in m2 and V is the total volume of the slick in m3. In this study the cessation of mechanical spreading is considered to occur when the slick thickness reduces to 10-5 V"4 m.

The formulas presented in Tables 1 and 2 were derived for simple slick geometries that exist under

idealized conditions. In the simulation model the radial spreading formulas are used when the slick is

nearly circular, and the one-dimensional formulas are used when the slick area is elongated. A slick is

considered to be elongated when the aspect ratio of the slick area is greater than 3. The aspect ratio

refers to the length-to-width ratio of the slick and the orientation refers to the angle 0 between the major

axis of the slick and the x axis, as shown in Figure 8.

The orientation of the slick is computed using the moments of inertia and the product of inertia of

the slick. The angle e can be calculated using:

tan (20) = • (29)

Ix -ly

where

IX =LP y= XxýPY= TY(30)

where L4 and ly are the moments of inertia of the oil slick with respect to the x and y axes, respectively,

and P, is the product of inertia. Once the orientation is known, the x and y coordinates of all the particles in the slick are transformed into coordinates in an x 'and y' coordinate system in which the x 'axis

is the major axis of the slick. The aspect ratio is computed by the following equation, with x' and y'

coordinates of all oil particles:

Aspect ratio = (31)

15
Table 2. Spread rates of oil slicks and phase transition times. (After Shen

et al. 1990.)

One-dimensional slicks, dLeldt

Spreading phase Spreading rate

Gravity-inertia 2 kio (Ag)/ 3 (L-1/ 3 dL t2/3 + L2

/

3 t-1/3)

3 ~121// dt-58d 3

Gravity-viscous k2o (Agv1/) (t I U-5/8 + dLt 3/8)

Surface tension-viscous k3o ((p2)"2 (tVY 11 4

4

Circular slicks, dR~dt

Spreading phase Spreading rate

Gravity-inertia /A (t2 3) dV

Gravity-viscous k2v (Agv12)l'6 1 dl 4 (3d t)V/3

Surface tension-viscous k2, (ap;1) (tvV"

Times of phase transition

Transition One-dimensional spreading Radial spreading

Gravity-inertia to (L0v0 [Ag32)'fl( 2Ž I V(. 1/3

gravity-viscous (k2i) tvWAg l

Gravity-viscous to (2&/3 (AgL4 d2 2P v/2I)2/3 (k2v)9 (_•_) (V2VAg)1/3

Surface tension-viscous W30 (k2t) (r

Note:

L2 = 0.5 volume of oil per unit length along the major axis of the slick (L is a characteristic

length); kIo = 139, k2 = 139, k30 = 1.43.

V = Total volume of oil slick; k2 = 1.14, k2 = 0.98, k2t = 1.60.

-- Figure 8. Definition sketch for variables used to compute the slick's aspect ratio and orientation.

16
Centroid of - -

Figure 9. Division of a nearly circular slick into pie - - S... Sick

shaped segments.

Oil Slick Strips

Strip rT "r 1 --.

Centroids

Figure 10. Division ofan elongated slick into

strips. 7,7 7 7; 77,

Direction of Spread J

This aspect ratio is used to determine whether the spreading of the oil slick is radial or one-dimensional. The ube of 3 as the aspect ratio for transition from radial spreading to one-dimensional spreading is subjective but gives reasonable results. For a nearly circular slick the slick area is divided into

eight pie-shaped sectors, as shown in Figure 9. For an elongated slick the entire slick is segmentized

into a series of strips, as shown in Figure 10. For both cases the spreading rate of each segment is calculated using Fay's formula shown in Table 2, independent of other segments in the slick. This simplification assumes that concentration gradients at the boundaries between neighboring segments are negligible for mechanical spreading (Shen and Yapa 1988). The segmentization also allows for different spreading rates in different regions of the slick, thus providing a more realistic description of the field situation. For a nearly circular slick the rate of outward movement of an oil parcel along the radial direction

in a particular pie-shaped sector located at a distance r from the centroid of the slick can be calculated

from the spreading rate of the mean radius r of the sector.

dr - r d. (32)

dt r dt

Parcels scattered at distances far away from the main slick will be excluded from this process since

small isolated patches of oil will not be subjected to mechanical spreading. In this study, parcels that

account for the outer 5% of the total slick volume are excluded. This is equivalent to excluding the

parcels located at a radial distance greater than 2.27 from the centroid of the slick, as determined from

numerical experiments for the spreading of circular slicks.

For an elongated slick the rate of outward movement in the width direction of an oil parcel in a

segment with mean length 1, located at distance I from the centroid of the segment, can be calculated

from the rate of spreading of the mean length of the segment as:'

A = L df (33) dt [ dt

Parnels for the outer 5% of the slick volume are not subjected to mechanical spreading.

Spreading under an ice cover. Hoult et al. (1975) suggested that appreciable mechanical spreading will

only occur during a continuous spill, because the oil thickness is stabilized when the equilibrium thick17
ness for the flow condition is reached. There is no pressure gradient or surface tension force to cause

further spreading of the oil. The ail reaches an equilibrium state when cavities formed by the ice roughness contain a volume of oL tha: can decrease only with a significant increase in the current speed. Since

a continuous spill will repeatedly add oil to 111 the cavities, the excess oil will effectively spread to the

empty neighboring cavities and establish an equilibrium state there.

The formula used to model mechanical spreading for a continuous spill under ice is:

r = 0.25(agQ02)l/ 6t2/3 (34)

where r = slick radius

A = (p - po)/pw, the relative density difference

g = gravitational acceleration

Qo = average volume flow rate from the beginning of the spill

h' = half of the root-mean-square roughness height of the ice cover.

This equation is a result of balancing the frictional drag from the ice cover with the pressure drop that

occurs as the oil flows into open cavities. In the simulation models it is assumed that no mechanical

spreading will occur for an instantaneous spill or after the continuous oil discharge stops. If the oil discharge is in progress and the slick is nearly circular, the mechanical spreading will be calculated by eq

34.

Shoreline boundary conditions

When oil reaches a shoreline, it will deposit along it. The deposited oil may be re-entrained into the

river by the current. The rate of re-entrainment can be formulated according to the "vulnerability" index of the shoreline (Gundlach and Hayes 1978). This index is designed to reflect the environmental

sensitivity of the shoreline to oil pollution. For beaches of different vulnerability indices, Torgrimson

(1984) suggested "half-life" values to describe the ability of the shore to retain the oil. Half-life is a parameter that describes the "absorbency" of the shoreline by describing the rate of re-entrainment of oil

after it has landed at a shoreline location. Table 3 presents half-lives for different types of shorelines

along with their vulnerability indices. According to the half-life formulation, the fraction of the beached

oil re-entrained into the water body during each time step dVb is calculated from:

dVb = 1- e-kAt = 1 - 0O.5At/ (35)

Vb

where Vb is the volume of oil on the beach. The decay Table 3. Shoreline descriptor and default

constant k can be expressed in terms of the half-life X as: parameters.

Half-life Vulnerability

[-In (11/2)] /X (36) Shoredine descriptor (hr) index

Instantaneous rejection* 0 0

For the purpose of calculating evaporation rate, the Exposed headland I 1

beached oil is assumed to occupy a predefined width of Wave-cut platform 1 2

each computational grid along the shoreline. This width Pocket beach 24 3

is determined by the user according to the field condi- and gbeach 24 4

Sand and gravel beach 24 5

tiom. Sand and cobble beach 8760 6

Exposed tide flats 1 7

Eapotion and dissolution Sheltered rock shore 8760 8

Sheltered tide flats 8760 9

Ezvqvetion. Evaporation is the main process that caus- Sheltered marsh 8760 10

es the low in oil volume during weathering. Mackay et Land 8760 11

al. (1980a) developed a formula for calculating the rate This option is added in the present computer model.

18
of evaporation. The volume fraction of oil evaporated is determined as

F = UE(1/C)[InPo + ln(CKEt + I/Po) (37)

where c•E = modification coefficient introduced in the present study to account for the effect of the

change in oil properties on the evaporation rate due to emulsification and other weathering processes

E = KEt = evaporative exposure term, which varies with time and environmental conditions

KE = KMAV/(RTVo)

KM = mass transfer coefficient (mi/s) (0.0025 LU 2d,

A = spill area (m2)

v = molar volume (m3/mole)

R = gas constant [82 x 10-6 atm n3/ (mole-K)]

T = surface temperature of the oil (K), which is generally close to the ambient air temperature TE

Vo = volume (i 3).

The initial vapor pressure P0 in atm at the temperature TE is determined from

In P0 = 10.6[1 - To/TE] (38)

where To is the initial boiling point (K). The constant C can be determined by the relationship TEC =

constant. Values of C for TE = 283 K are given in Table 4. For crude oils of different API index values, C

at TE = 283 K are given in Table 5 along with T0 .This table can be replaced by the following functional

relationships obtained through curve fitting (Shen et al. 1990):

C = 1158.9API -1.1435 (39)

and

To = 542.6 - 30.275API + 1.565API 2 + 34.39API 3 + 0.0002604API 4 (40)

The API index and the specific gravity of the oil are related by

Specific gravity = 141.5/(API + 1315). (41)

Table 4. Suggested evaporation parameters for Table 5. Suggested evaporation

various petroleum fractions (TE = 283 K). parameters for various crude oils

(After Mackay et al. 1980b.) (TE = 283 K). (After Mackay et al.

1980.)

Oil type T, (K) C P, (atm)

Motor gasoline API (g cnr 3

) C T, P,

Summer 314 5.99 0.313

Winter 308 6.23 0.39 10 1.0 89.2 366 0.044

Aviation gasoline 341 2.81 0.12 12 0.986 69.4 348 0.088

Diesel fuel 496 5.57 3A x 1"-4 15 0.966 52.1 339 0.13

Jet fuel 418 5.06 6.0 x 10-3 20 0.934 34.7 329 0.18

No. 2 furnace oil 465 7.88 1.1 x 10V 25 0.904 27.2 330 0.17

Lube (heavy and light) 583 8.61 1.32 x 10-3 30 0.876 22.33 325 0.21

Heavy gas oil 633 8.99 2.0 x 10-6 35 0.850 1935 314 0.31

Residuals 783 3.37 735 x 10-p 40 0.825 17.9 304 0.45

Ligt gas oil .473 6.37 8.1 x 10-4 45 0.802 16.4 283 1,004

19
"

'' ~o-s I Z I"

~j~D i

S.. .. : .I%-- ,r

- 0S _ _ _ _ _ _ _ _ _

2 - j
The molar volume of oil is required in eq 37. The value, ' Aie molar volume can vary between 150 x

10-6 and 600 x 10-I m3 /mole. Typical values of molecular weight of various oil components are given

by Moore et al. (1973). The compositions of several different oils are given along with the concentration

and molecular weight Mi of each component. Based on these the molecular volume can be computed

(Shen et aL. 1990). The volume V0 for each time step is determined by

V0 V0 + AVm (42)

where Vo = effective initial volume of surface oil

V0 = initial volume of oil spilled on the water surface

AVm = volume increase in the surface layer due to the exchange between the two layers.

Dissolution. Cohen et al. (1980) developed a method to calculate the total dissolution rate N:

N = KASS (43)

where N = total dissolution rate of the slick (g/hr)

K = dissolution mass transfer coefficient, assumed to be 0.01 mr/hr

A, = slick area (M2

)

S = oil solubility in water.

Huang and Monastero (1982) suggested that for a typical oil the solubility can be calculated by

S = Soe-0.1t (44)

where So is the solubility for fresh oil and t is time in hours. Huang and Monastero suggested a typical

value of 30 g/m 3 for So. In the present study the formula of Lu and Polak (1973) is used, which provides more information on solubility. In this method the rate of dissolution is calculated as

rd = cde -o.t (45)

where rd is the rate of dissolution (mg/m2-day). For three oil samples tested the coefficient c and d are

given in Table 6.

Table 6. Dissolution coefficients at 25°C.

(After Lu and Polak 1973.)

Oil c d KS, - cd

type API (mg/m2

) (/day) (11g2 hr)

No. 2 fuel oil 35.5 1043 0.423 0.0184

Crude oil 38.6 8915 2.380 0.884

Bunker C oil 14.8 459 0.503 0.0104

Mixing, emulsification and sedimentation

Terms representing mass exchanges in the surface and the suspended layers in eq 1 and 4 are calculated at the end of each time step. The volume of oil mixing into the suspended layer per time step in a

grid is calculated as yCsAxAyAt. Similarly the volume of suspended oil resurfaced is calculated as

aCrVVbA~xAyAt. The volume of oil settled to the channel bed is calculated as rCVhAxAyAt.

COMPUTER MODEL AND CASE STUDY

Based on the analytical formulation presented in the previous chapter, we developed a computer

model for the Ohio-Monongahela-Atlegheny River (OMA) system (Fig. 11). This chapter describes the

21
model along with a case study. It should be not- Table 7. Relationships between local and global

ed that the following discussions deal with coordinate systems.

many programming details. For a user to ap- Reion Orlgipii of local Rofatiii, a•g'lh

ply the model to the OMA system, a menu-driv- no. coordinates (X,,y, deV,',', Per

en model (Yapa et al. 1990a) provides a simple 1 0.0,0.0 0.0 Ohio

means for spill simulations. 2 348200.0, 138800.0 -90.0 Ohlo

3 3754000, -42200.0 -90.0 Ohio

Model implementation 4 392000.0, -84800.0 -180.0 .lononghela

5 315800.0, -10460.0 -180.0 Monionglhcla

The g-rid sstern 6 383200K0, -36800.0 360.0 Allegheny

7 592000.0, -768000.0 90.0 AIh'gheov

In the computer simulation a two-dimensional reference grid is used to discretize the river. This reference grid is required to identify the or x

locations and velocities of oil parcels. The grid

boxes in this model are squares of equal size.

The iocation of each grid box is identified by its :

x and y ndices.

The reference grid described is adequate for

rivers with only one orientation, e.g. from east

to west, or with only small variations in their .

orientation. Unfortunately the OMA river sys- -:

tem has a complicated configuration with ;U•

changing orientations. Seven sets of local coor- g4'k

dinates are introduced, as shown in Figure 12, • °0

for different regions in the system to accommodate the changing orientations. During the sim- 3

ulation, transformation between a local coordinate system anct the global system is required. 4i>.

The local x-axis is set to be in the main flow di- •';:-

rection toward upstream. Relationships be- ' ,,'i'.

tween the global and local coordinate systems r

are given in Table 7. x, M land'z Monongahela River

Besides the fixed-grid systems a temporary V 3 and ": Allegheny River

moving grid system traveling with the slick is 07

used in the model to save memory space. This

moving grid system will just be long enough to ...

cover the length of the oil slick and the travel

distance in one time step of the Bow computa- Figure 12. Local and global coordinate sl Jems.

tion. In a long simulation the slick may become

too long to be covered by the moving grid because of bank deposition and rejection. For this case there

are two options, one of which is to make the length of the moving grid long enough to cover the longest

slick during the spill or the total length of the river. However, this option is limited by the n'emory

space of microcomputers. Hence, the second option is used. In this option the model selects a reasonable length of moving grid, both for computer memory and for the length of the oil slick, for example,

30 miles long. The moving grid will follow the slick's leading edge and cover a 30-mile-long segment of

the river upstream from the leading edge. The oil outside the moving grid, if it exists, will generally be

inpia~nt.

A series of successive river bends can cause difficulty in defining the river boundary. In Figure 13

there are four continuous iver bends. There are more than three sets of boundary points that intersect

with a local j grid line. In this model all land boundary points, except for the two extreme end points,

22
Local x or i

Grid Directon x

Local y or 1

Grid Direction

0l Bank Grid Point

a• Combination

Flow

Figure 13. Method for determining the land grid along a series of river bends.

are treated as island boundary points. The boundary grids along a series of bends can be divided into

three types, as shown in Figure 13. They are island grids, bank grids and their combinations.

Release of oil

The oil released into the river is represented by a user-specified number of particles. Each particle

can be considered as a parcel of oil. The user can select a particle number up to 1000 for each layer at the

beginning of the simulation. During the simulation the number may increase to 3000 or even more than

5000 due to the exchanges between the surface and the suspended layers. For simulation accuracy,

approximately 500 particles in each layer is suggested to be used at the beginning of the simulation.

However, in an emergency situation the user may reduce the particle number to 50 or less to reduce the

computing time. The volume of oil released into each layer is determined by the user's input. Parcels in

the same layer have the same volume at the beginning. The total number of parcels in each layer is the

same at the beginning, so that the oil is released at a constant rate in each layer. For example, consider a

spill of 10,000 gallons, with 7000 gallons in the surface layer and 3000 gallons in the suspended layer. If

the user chooses 1000 parcels in each layer, then each surface oil parcel will contain 7 gallons, and a

suspended one will contain 3 gallons at the beginning.

An oil spill will be treated either as an instantaneous spill or a continuous spill. When the spill duration is zero, the spill will be treated as instantaneous. Otherwise, it will be treated as a continuous spill.

Velocity and depth distributions

Analysis of the transport of oil spill parcels in a river requires well-defined water velocity and depth

distributions. To calculate these distributions, various information is needed, including the total discharge, river stage, cross-section geometry, shoreline characteristics, ice condition and wind condition.

The velocity and depth distributions in the river are calculated in two stages. First, a one-dimensional

unsteady flow model is used to generate discharge and river stages at specified cross sections referred

to as nodes. For the OMA model the Flowsed model is used for this purpose (Johnson 1982). The oil

spill model uses this information to compute two-dimensional velocity and depth distributions.

In the oil spill model, branches are used to describe the channel configuration. The beginning and

end of these branches are generally selected to coincide with nodes of the one-dimensional unsteady

flow model. In cases where it is necessary to have more branches in the oil spill model than in the

23
Flow •

S• " l i Flow

33

a. Nodes. b. Branches.

Figure 14. Schematization of a river reach with an island.

unsteady flow model, water levels at intermediate points are obtained through interpolation. When an

island is encountered, there will be branches on each side of the island (Fig. 14). Each of these branches

includes a portion of the stream tubes. The information available is the magnitude and direction, i.e.

the x and y components, of velocities at stream-tube centers at all cross sections.

The assignment of velocity and depth through the entire river requires the establishment of a grid

system, as shown in Figure 15. Velocities and depths obtained from the stream-tube analysis are interpolated to all grid boxes. River boundaries for this grid system are defined by boundary grid boxes. For

each grid along the x-direction, two corresponding y-grid boxes are used to define the location of the

riverbanks. In places where there are islands, more y-grid boxes are used to define the upper and lower

boundaries of islandcs. The model can handle any number of islands in the river. However, if more

islands are intersected by a vertical line along which the x coordinate is constant, more space will be

needed for the input data, and simulation speed will be slowed down. Defining islands in the grid

system is independent of defining islands in a branch configuration. This allows for a greater flexibility

in the simulation of oil slick transformation. For example, a small island that covers only three or four

grid boxes, which may be too small to be defined in the branch configuration, can be easily defined in

the grid box system. All islands must be defined in the grid box system.

Extra

igriub(i) Grid Layer River

/ Bank

19

17

16 W 1-l-1"9

1 1 13 31 ý1 13 13

12

12 12 12

10 i 2, i)

S9 9%'•d , I- igdb(2, j) • 8 a a ••

7 2 3 46 78 77

- - igrub(1, r)

4 4 4 ',, '6 5 igdb(i, i)

2 2:IX

1 2 3 4 5 8 7 9 9 10 11 1213 14 15 1817IS19 20 2122 2324 25

i - grids

Figure 15. Grid boxes and river boundary representation.

24
The velocity and depth within a box are assumed to be uniform. Only boxes found within the defined river boundaries are assigned velocities and depths. Generally spacing between cross sections in

the stream-tube computation is much larger than the grid size Ax. The following procedure is adopted

for assigning velocity and depth to grid boxes based on the velocity and depth at cross sections. In the

following discussion the term "velocity/depth point" is used to represent a point at which a computed

velocity/depth is located. At first all velocity and depth points computed from the stream-tube analysis will be assigned to the boxes within which they lie. If more than one velocity and depth point fall

into a box, the average value of velocity/depth is assigned to that box. In the next step, more velocity/

depth points will be obtained through interpolation between two successive cross sections. The number of interpolation points between two successive cross sections can be changed through user-defined

input. These velocity/depth points obtained through interpolation will now be assigned to grid boxes

using the same method as before. After this, if there is any box without an assigned velocity/depth, an

average velocity/depth calculated based on neighboring boxes will be assigned. Upon the completion

of these steps, the entire river is scanned to ensure that velocity/depth points are assigned to all boxes.

Oil slick orientation

The surface slick can spread either axisymmetrically or one-dimensionally. The determination of

which type of spreading prevails depends on the shape of the slick, i.e. the aspect ratio. Therefore, after

the oil particles moved to their new positions, the centroid, aspect ratio and orientation of the oil slick

are calculated. These calculations will consider only oil particles in the water grids.

Table ary"ye 8. of Bound- locks : ,.

and dams.

Code Rejection

no. rate

13 2

14 3

15 4

16 5

16

Figure 16. Indexsystem for locks and dams.

• =- Lock and Dam

13 = Particle moves opposite local x-coordinate.

14 . Particle moves along local y-axis.

15 . Particle moves opposite local y-axis,

16 = Particle moves along local x-axis.

Boundary conditions and locks and dams

Twelve half-life values were defined in Table 3 for shoreline boundaries. Oil slicks passing a lock

and dam will be completely mixed into the water column. The program is directed to this condition

through the use of shoreline boundary codes. Four new boundary condition codes are added. As defined in Table 8, code numbers 13-16 refer to lock and dams with four different orientations (Fig. 16). A

set of special rejection rate codes are given for these conditions to direct the program to well-mixed

conditions. Oil parcels deposited on the banks at a lock and dam are instantly rejected into the water.

This is achieved by using a set of shoreline codes 0 and 6 as in Table 8. This set of codes can also be used

at other shoreline locations for instantaneous rejection if needed.

25
CONTINUED

Choose river combination:

1. Ohio-Monongahela

2. Ohio- Allegheny Call advsur and advsus to advect particles both in surface and

suspended layer

erunning mode: Call edam to mix surface particle into water body, if surface

Determine rparticle passing dams and locks

0simulate from t = 0 Call orientto calculate spill center on surface and in

1 to start from previously stopped point suspension, and establish if axisymmetrical or onedimensional spreading prevails

Call pievolto calculate oil volume on surface in each pie area Read: if axisymmetrical spreading prevails

Oil spill constants Call strvol to compute oil volume on surface in each strip if

River division data one-dimensional spreading prevails

Spill site, volume, initial particle no.

Determine the type of oil spill: Instantaneous or continuous Call emulfn to calculate emulsification

Determine initial moving Detemineiniial grds:iO, grids: i0, lend endCall ovig Call redepn evaporto to calculate compute evaporation resurfacing and sedimentation

Check if oil spill site is correct Call disolto compute evaporation

If it is a continuation run, read back the previous results Cal disolu to compute dissolution

For uf= 1, no. of unsteady flow step Call lumperto compute particle number, total volume, position

index for each grid in order to save time for calculating

exchanges between two layers If Jut> 1, check oil particles out moving grids Call exchan to do exchange between surface layer and

Reset moving grid for next time step: il, lend suspended layer

Call bndrej if particles have hit the shore and determine if

particle will be retained or rejected

For the reach covered by moving grids: Call sprdax to spread oil particles in axisymmetrical case

Read cross section locations and connection information Call sprdldto spread particles in one-dimensional case

Read cross section geometries Call rearrg to lump small particles and divide large particles

Read boundary grids for river shore, islands and lock and dam Call ploter to calculate one-dimensional concentration

Read oil retention parameters for boundary grids distribution

Read ice region data

.... - Output information for generating numerical plot of oil spill

Call vetdep to calculate 2-D velocity and simulation and one concentration distribution

depth distribution within moving grids

For i=1. no. of oil spill step

ff it is continuous oil spill, and the spill is still progressing, call

presur and prsus to release oil particles at equal Intervals

Figure 17. The computer model ROSS2.

Model strcue

The computer code is written in standard Fortran 77. No compiler-specific extensions have been

used. The code has been tested only on 32-bit machines with Fortran 77 compilers that support 32-bit

computations. The program can be run on other machines (e.g. 24-bit) with minor modifications. Appendix B gives a listing of all the subroutines used. The model is named ROSS2, and the model structure it dsown in Figure 17. All variables are dimrensioned to have sufficient storage to run for each

segment of the river system. The velocity and depth distributions will be kept only for a reach long

26
enough for the oil slick to travel one un- Table 9. Cross references of subroutines.

steady flow calculation time step. The Program Cal(I Input fib's

model can be run for a new spill simula- ross2 input, ckout, positn, veldep, ohndc, River division data,

tion or continue from a previously ter- ohvad, advsar, presus, advsur, ad'us, spill data,

minated run by changing weather and bndloc, check], ckout, con, rejiloc, shoreline condition,

emdain, orient, pievol, strvol, emufii, ice regions oil characteristics. redepn, areacl, evapor, disolu, lumiper,

exchan, bndrej, shift, sprdax, prilif,

Subroutine cross reference rearrg, ploter

The model consists of the main pro- advsur gauss, randu, bndrej, refloc None

advsus gauss, randu, bndrej, rejloc None

gram ROSS2 and 37 subroutines, as sum- areacl None None

marized in Table 9. An alphabetically ax- bndloc None None

ranged listing of the subroutines is giv- cndral rcfc None

calcij None None en in Appendix B. calcIm None None

The common block includes arrays check1 calcij None

and their definitions shared by subrou- ckout calcij,calctm None

con conchk Previous simulation

tines. Sizes or arrays can be adjusted in results

this block according to the moving grid conchk calcij Previous simulation

results size, maximum particle number, river disolu None None

segment number, maximum island emdam None None

number and stream-tube number. Table emulfn None None

evapor None None 10 gives definitions of important van- exchan None None

ables included in the common block, input None Spill constants, river divison data

When the model is applied to other riv- None river N one

ers, dimensions should be adjusted. Oth- lumper1 None None

er variables used in the program are pre- ohndc None Flow data

ohvad None None sented in Appendix A. pievol None None

ploter None None

Main program posit2 None None

The main program is used to estab- presur advsur, gauss, randu, bndrej, rejloc None

.is te presus advsus, gauss, randu, bndrej, rejloc None

rearrg lumperl None

sides interactions going through the relca None None

common block. The major functions of reloca None None

reloxy None None the main program can be summarized reixy None None

as follows: shift None None

*Read oil spill information in the data sprdax None None sprdld None None

file xxx.spl. strvol None None

oRead boundary grid box numbers veldep ohndc, ohvad Geometry data,

and set boundary conditions from flow data

xxx.zgr and xxx.bnd.

'Call con to read back the simulation results if it continues a job from a previously terminated run.

'Call ckout to adjust the moving grid position based on the oil slick position.

*Call veldep to calculate velocity and depth distributions of the river reach superimposed by moving

grids during every unsteady flow time step.

oCall presus and presur to release oil particles at equal intervals for the oil spill duration.

*Call advsus and advsur to advect oil particles in both the suspended layer and the surface layer.

'Call emdam to mix oil particles into the water column when passing locks and dams.

OCall orient to calculate the oil slick orientation, determining whether the slick is quasi-one-dimensional or quasi-axisymmetric; call pievol to calculate the volume of each pie segment for an axisymmetric slick; or call struol to calculate the volume in each small strip for a one-dimensional slick. All

of these are for mechanical spreading calculations.

*Call emufin and redepn to compute oil volumes into the water column, resurfaced to the water sur27
Table 10. Definition of variables in the common block (dimension parameters).

Variable Value Definition

numbm 37 Maximum number of branches of river segments

numdep 66 Maximum number of sounding depth locations at all cross sections

numice 10 Maximum number of ice regions

numpar 15000 Maximum number of oil particles

numpie 8 Number of pies into which the axisymmetrical slick is divided

numvl0 2 Number of positions where the velocities are assigned zero

numpos 26000 Maximum number of moving grids

numsec 330 Maximum number of cross sections among river segments

numstb 11 Number of stream tubes

nwth 10 Number of shorelines, including islands, locks and dams

nudxx 1100 Maximum grid number in x direction of moving grids

nurmxgd 4600 Total number of grids in the x direction

nseg 5 Maximum number of river segments

nisl 4 Maximum number of islands a cross section can meet

nout 800 Maximum number of particles out of moving grids that can be plotted

face and deposited onto the riverbed during the time step, as well as the cumulative values.

*Call aread to calculate the slick area in water, under ice and on land for the calculation of evaporation and dissolution.

e Call evapor and disolu to calculate evaporation and dissolution.

*Call lumper and exchan to calculate oil exchanges between the suspended layer and the surface

layer.

*Call sprdax for an axisymmetrical oil slick, or call sprdld for a one-dimensional slick, to calculate

mechanical spreading.

*Call bndrej to reject oil particles back into the river according to the bank retarding condition (halflife).

*Call rearrg to adjust particle sizes, ie. lump small particles togptber and divide large particles to

satisfy the requirements of the random-walk simulation and the computer memory space.

*Call ploter to calculate the one-dimensional cross-sectional average oil concentration distribution

along the river.

-*Save the simulation results into corresponding output files at a designated output frequency.

Subroutines

Since all subroutines and related variables are listed in Appendices A and B, and some simple subroutines can be easily understood with comments in the listing, no further explanation will be given.

Block diagrams for some main subroutines are given in Figure 18 to assist in understanding the programCase study

The model is applied to the oil spill in the Monongahela River on 2 January 1988 (Miklaucic and

Saseen 1989). Simulation results for the leading edge and concentration distribution of the oil plume

will be compared with the field data.

Model parameters

A number of empirical parameters, which can affect the simulation results, are to be selected by the

user. In this section the effects of these parameters will be discussed to assist the user in making appropriate selections.

One of the parameters Aat will affect the movement of the slick is the shoreline boundary condition.

Under an instantaneous-rejection shoreline boundary condition, oil will not deposit on riverbanks. Partice moved to the shore will be rejected to water instantly. This will give a shorter slick, and its leading

edge will advance at a faster rate.

Thoetcaly the simuLition sbmuld not be affected by the spill time step. However, numerically a

28
START

Set river region for velocity and depth

distribution calculation (veldep)

Input data:

a Geometry data of river cross sections (veldep)

e Discharge and water level at nodes of unsteady

flow model (ohndc) * Interpolate discharge and water level to stream- a. Velocity and depth computations.

tube cross sections (ohvad)

Stream-tube velocity and depth:

* Calculate stream-tube velocity and depth (ohvad)

* Determine velocity direction (ohvad)

Grid velocity and depth:

L Unearly interpolate stream-tube velocity and depth

to grids

* Search zero-velocity grids and interpolate velocity

and depth with the four surrounding grids (ohvad)

INext Parcels

Save the results for oil advection and

mixing calculation (ohvad) within oil Yes Release oil parcels by spill duration the same time interval

Next small time step No

EDDetermine advection velocity:

"" If under ice cover, check moving criteria

"* If in open water, consider wind effect

I Advect oil parcels by a small time step within which

oil parcels cannot travel more than one grid size

S Check boundary conditions:

9 It oil passes a lock or dam, mix it into water column

9 If oil hits land, mark It as land parcel

b. & "da aid adWcdion ofi parcels.

Figasm I8. Subroutine.

29
ISTART

"* Calculate emulsification in each grid (emulfn) 1

"* Calculate resurfacing in each grid (redepn) I

"* Calculate evaporation in each grid (evapor)

c. Mixing computation.

* Add one particle in each grid suspended layer with

the total volume of emulsification and dissolution

* Add one particle in each grid in surface layer with

a total volume of resurfacing

"* Adjust particle volume in surface layer

" Adjust particle volume in suspended layer S

"* Adjust volume of particles on land

END Calculate orientation of oil slick

Axisymmetric No Apc ai e One-dimensional

spreading spreading

angled wedges and calculate and calculate volume change

volume change rate rate for spreading

Calculate average spreading Calculate average spreading

rate of a pie segment rate of a strip

Sprea oil particles proportional to the Spread oil particles proportional to the1

ddistancep from the center of the slick distance from the center line of the strip1

Check oil particles; passng Check oil particles passing

dam or lock and hitting land dam or lock and hitting land

Next wedge I F Next strip

d. MedatW spreading computatikm.

Figum 18 (cont'd). Submtutines.

30
TSTART

Lump small particles:

"* Check and lump small particles in each grid of

surface layer

"* Check and lump small particles in each grid of

suspended layer

"* Check and lump small particles in land grids

Remove oil particles with zero volume in surface

layer, suspended layer and on land

Divide large particles:

"* Check and divide large particles in surface layer

and redistribute them randomly in each grid

"* Check and divide large particles in suspended layer

and redistribute them randomly in each grid

"* Check and divide large particles on land and

redistribute them randomly in each grid

f Recount the number of particles on land and in water

and renumber them for next time step

N e. Particle size adjustment.

Figure 18 (cont'd).

large time step can lead to inaccurate results, especially for the non-instantaneous-rejection condition.

In the model, boundary rejection, exchange between layers, mechanical spreading, evaporation and

dissolution are done once every spill time step to save computation time and memory. The spill time

step can affect the slick length through all of these simulation components. If the time step is too large,

discontinuous distributions of the portions of oil deposited on shoreline, on the water surface and in

the suspension will appear. When the instantaneous-rejection condition is used, the effect of the spill

time step is small

The random-walk theory requires a large number of particles with equal mass. Numerical experiments show that the equal mass requirement can be relaxed if the particle number is large enough. In

the present model the particle number increases when the particle mass becomes more uneven due to

the exchange between the surface and suspended layers The model also keeps at least one particle in

each grid that contains oil.

To demonstrate the effects of the shoreline condition, simulations for the time step At and the initial

particle number No are carried out for the 1988 oil spill in the Monongahela River for a duration of 10

days. The model parameters used in these simulations are:

" Buoyant velocity: v = 0.01 ft/s;

" Resurfacing coeffiient: a = 1.0;

• m iation coefficient P = 10-1/s;

* Riverbed deposition coefficient: = 10/s;

31
"* Initial particle number No = 40, or 400 each layer;

"* Spill time step: At = 5, 15 or 30 minutes; and

"S Evaporation modification coefficient E = 0.3.

Ice covers can affect the movement of oil in the river. According to the model the advection volocity

of the oil slick under a smooth ice cover will only be reduced by 10-15%. The ice cover has very little

effect in the advection of suspended oil. For the 1988 spill the air temperature was below freezing, as

shown in Figure 19. Ice may have existed in the Ohio River during the later part of the simulation

period. However, since the ice covers formed in the Ohio River are generally smooth, the existence of

an ice cover will be neglected in the sample simulations for simplicity.

30

S20

E 10

0-

-10

0 2 4 6 8 10 12 14

Days from January 1, 1988

Figure 19. Air temperature variation during the simulation period,

Pittsburgh, Pennsylvania.

Figure 20 shows simulation results with the idealized instantaneous-rejection shoreline condition

and a partial-rejection shoreline condition determined for the OMA system. Table 11 and Figures 21

and 22 compare simulated and observed locations of the leading edge of the slick and the concentration peak. These results show that the slick locations simulated using the instantaneous-rejection shoreline condition do not agree with the observed locations, although the errors are relatively small. Simulations using the partial-rejection shoreline condition compared well with the field observation. SimuTable 11. Observed and simulated locations of the leading edge and the concentration peak.*

Leading edge Concentration peak Leading edge C.;.,n.rati,;;;

(rimer mile) (river mile) (river mile) (river mile)

Date Observed Simulated Observd Simulated Date Observed Simulated Obseroed Simulated

Instantaneous shoreline rejection Partial shoreline rejection

2 Jan -25.0 -25.0 -25.0 -25.0 2 Jan -25.0 -25.0 -25.0 -25.0

3 Jan 0.0 -1.21 -4.5 -2.54 3 Jan 0.0 0.41 -4.5 -2.11

4 Jan 25.0 22.24 no 19.04 4 Jan 25.0 21.06 no 14.06

5 Jan 45.0 52.14 no 49.90 5 Jan 45.0 43.11 no 32.34

6 Jan 65+ 74.32 no 72.58 6 Jan 65t 58.82 60.0 52.43

78" 74.32 60.0 72.58 780* 58.82 60.0 52.43

7 Jan 73t 98.89 65.0 %.07 7 Jan 73* 77.62 65.0 64.06

88" 9869 77.0 96.07 88** 77.62 77.0 64.06

8 Jan 89.0* 115.8 77.0 112.9 8 Jan 89.0* 90.43 77.0 78.41

"99.0" 115.8 84.0 112.9 99.0" 90.43 84.0 78.41

9 Jan 94. 125.1 no 122.6 9 Jan 94.0* 101.5 86.0 8539

108" 125.1 86.0 122.6 108" 101.5 86.0 85.39

10 Jan 112.0 1252+ 96.0 125.1+ 10 Jan 112.0 109.7 96.0 93.51

0 Initial partide numbe in each layer is 400, spill time step is 900 s.

f Vory an Al.s-Kedidom (1989).

-Rai"yet aL (199).

32
2.0 x 103 I I I I I I I I I I i I I I

Monongahela Ohio

1.6

.0

a. Instantaneous shore- ' 1.2.

line rejection, N, = 400, 0

At= 5 min. oo

(D

E 0.8 -.. U C M

0.4

2.0x 10 3

Monongahela Ohio

1.6-

.0

b. Instantaneous shore- • 1.2 -

line rejection, N =400,

At= 15 min. (3

"E0.8 .-

-)--l i i -

2.0 x103

1

Monongahela Ohio

1.6

C

c. Instantaneous shore- B 1.2 -

line rejection, N, = 40,

At = 30 min. ( -

0.8 -I Vo o N 00 M

0.4 I

Figure 20. Simulation results -

showing the variation of con- .. L

centratoin dsibution. -so so 1o0 150

River Mile

33
2.0x 103 I I I I I I

Monongahela Ohio

1.6

0

-t S1.2 d. Partial shoreline rejection, N,= 400, At C

0 = 5 ln.

S 0.8

0

C14 0q r"t• 4 0 C•oO a)

0.4

2.0x 103 I I I

Monongahela Ohio

1.6 I

CB

c0

im 1.2 e. Partial shoreline re- C

jection, NIV= 400, At =

o 15 min.

E 0.8 -

o., •1 , •. .' -' • ' • "" 0. JIc A-c I, CO 0

0.4

2.0x 103

Monongahela Ohio

1.6

"1.2 f. Partial shoreline reC jection, N, = 400, At S- = 30 min.

0.8

WC I'- CO C -" I -, - -• •-~- -•

0.4 - Figure 20 (cont'd). I Simulation results

showing the variation

0 0 100 150 of concentration disRiver Mile tribution.

34
U1/5 E. Liverpool

New Cumberland

116 'preirton et Beaver 117 ".Wefalon 17 ,. Steubenville V;8 O/04•• IA *- talso

101/70

1/9 ;.,0' Warwood 116

1t10 q• 1/9 ," 18Wheeling

11,0 Moundsville Esol

-Hannibal StationEmon

1/3 N1 Pittsburgh

L;0 2

e Observed (Vicory and Ashles-Kedziora 1989) North Braddock

O Observed (Berkey et al. 1989) 1/2

* Simulated a Elizabeth

N. Charierol Spill Site

Maxwell

."oongahe\]

Figure 21. Comparison of simulated and observed leading edge positions; N, = 400, At 300s.

O\.io

1/6 E. Liverpool

New Cumberland &

1/5

W/7eirton er Beaver1/9 Steubenville falls

1/10 Warwood 114U

.0 Wheeling Moundsville

Hannibal Station Emswch

* Pittsburgh

1/3 1 - ,,-2

* Observed (Vicory and Ashles -Kedziora 1989) North Braddock

o Obseved (Berkey et al. 1989) 1/2

8 Simulated Elizabeth

N. Charierol Spill Site

Maxwell

%'0,7ngaht'

Figure 22. Comparison of simulated and observed concentration peak locations; N, = 400, At 300 s.

lations using small particle numbers generally resulted in shorter slicks with higher concentration. The

effect of time t(ep At is shown in Figure 23. This figure shows that the slick slowed down with the

increase in At. However, this effect is small.

A sumnmary of computing time used for all of the sample simulations run on a SUN SPAC stationl+

is givenin Table 12. This t, "show6 that the computing time decreases with a reduction in particle

35
120

• __ t = 8 days

100

•" a Simulated

w

0 Observed

80

0-

m 0

0 10 20 30

At, Spill lime Step (min.)

Figure 23. Effect of time step At on leading edge positions.

Table 12. Comparison of computing time.

Parameters CPU time (s) Parameters CPU time (s)

Instantaneous rejection i'mtial rejection

N,=400;At=300s 23,855 N, = 400; A, =300s 67,354

No = 400; At =900 s 10,797 No = 400; At = 900 s 32,248

No = 400; At =1800 s 3,612 N0 = 400; At = 1800 s 19,258

number and an increase in At. The use of the instantaneous-rejection shoreline condition reduces the

computing time by approximately 60%. It is recommended that, for field applications in theOMA river

system, the partial-rejection shoreline condition with At = 900 s and NO = 400 be used, considering the

simulation accuracy and computing time. The ratio between the real time and the computing time is

about 1:12. This ratio is smaller at the early stages of simulation. For emergency applications an instantaneous-rejection shoreline condition along with a smaller No and a larger At should be used. The ratio

between the real time and the computing time is about 1:240 when At = 30 minutes and No = 40 are

used.

Applications to non-oil spills

As mentioned earlier the model can be applied to non-oil spills. Non-oil spills can be simulated by

inputting appropriate values for the spill parameters described earlier. For neutrally buoyant miscible

materials, the buoyant velocity should be set to zero.

INPBQ DATA FILES

There are two categories of inputs that are required to run the model. The first category can be conidered to be fixed data for a given river readc. This is the information required to describe the shoreline and cros&-sectnal geometry of the river. Normally there is no need to adjust these data. The secnd catety includes various parameters that may be adjusted according to the flow condition and

Vill a rstica. To completely understand the setup of the data files, it is helpful to go through the

s -te procedre that follows. Sample input data files are given later. If the user is only interested

m, ad**% pammeters, haging the spill locatioi or establishing new stage and discharge conditions

ota pmvlouy modeled river, no new data files need be created.

36
Data file creation

For demonstration purposes the Ohio-Monongahela river combination will be used as an example.

This river shows the complexity of a river system that L Ite model can handle. Where the data fies for

the Ohio-Monongahela river combination are not sufficient to describe the detail, other river data will

be presented. Eight data files exist for inputting information into the computer model:

xxx.geo River geometry, cross sections and branch geometry

xxx.igr Boundary grid box numbers

xxx.ice Ice regions, ice roughness and oil viscosity

xxx.flo Water level and discharge at nodes determined from the unsteady flow model

xxx.bnd Index for half-life value assignments to shore grids

xxx.spl Oil parameters, spill location, wind condition and ambient temperature.

xxx.ndc Information to transform computed unsteady flow results to what are needed in the oil spill

model flow computation

xxx.msc River segment definitions, parameters for exchange between two layers and information needed to coordinate transformation data.

The first step in preparing the data file is to make a sketch similar to Figure 24 to describe the branch

and cross-section numbering system used in the data file. The procedure for obtaining this sketch is as

follows:

1. Determine the number of river segments. In this model the entire Ohio-Monongahela-Allegheny

river system is divided into seven segments, as described earlier. The numbering of segments

depends on the river combination. The Ohio-Monongahela has five river segments numbered

from downstream to upstream. Each segment has its own local coordinate system, which can be

transformed to the global coordinate system using the transformation data in data file ohmon.int.

All geometry data willbe prepared segment by segment with independent numbering. To ensure

that all water grids are covered, cross-section numbers will have overlaps in the vicinity of the

division line between river segments.

2. Determine the number of branches. Each branch must contain at least two unsteady flow nodes

and one stream-tube cross section.

3. Number the cross sections in consecutive order from the downstream end to the upstream end of

the segment. Around islands, first number the bottom side up to the last cross section prior to the

confluence into a single channel, then go back to the top side and continue the numbering sequence.

29

25

* Unsteady Flow Nodes 21 22

Stream-tube Cross Sections 20 0010

18

175 15

16 12

Figure 24. Definition setch of the rkvr numbering system.

37
4. Determine the number of stream tubes in each branch. The total number of stream tubes always

remains the same. The stream tubes must be divided according to the ratio of the flow split around

the island.

Once the branches and cross sections for the river are established, scaled maps of the river and cross

sections are used to:

1. Establish the global x-y Cartesian coordinates, to be superimposed over the river.

2. Establish the reference coordinates on the left bank (looking downstream).

3. Digitize the entire river boundary and cross-section locations simultaneously.

4. Digitize the cross-section geometries by measuring a transverse distance from the reference point

on the left bank to a corresponding sounding point.

5. Determine the way to assign flow direction to each cross section. For the cross sections along river

bends, use 1 as the direction index, which implies that the flow directions will be perpendicular to

the cross section. Along straight reaches, especially in reaches with islands, 0 is chosen as the flow

direction index. This implies that the flow direction will be along the line connecting corresponding stream-tube centers in two neighboring cross sections.

6. Schematize boundary boxes, i.e. for every x-grid there exists corresponding upper and lower riverbanks and island boundary boxes.

The files are generally broken up into blocks and cards. A block covers a group of data that may

contain one or more card types. A card type is one line of specific data, which is sometimes repeated.

(For example, Block 4 in xxx.geo has Cards 1 and 2, where Card 2 is repeated as many times as needed.)

By inspecting the example of a card and comparing it to the complete sample data set, it is easy to see

how the entire file comes together.

Most of the data read into the model are in list-directed I/O (free format). If column numbers are

shown, the data must be formatted accordingly. Otherwise it is necessary to put only one space or

comma between each number in a card.

In the following sections, data files of the Ohio-Monongahela river combination will be used as

examples to explain the setup of data cards.

Geometry data- ohmon.geo

Ohmon.geo consists of five segments of geometry data for the five corresponding river segments.

Each segment has four blocks. None of the data in this file needs to be adjusted from one spill to the

next. The user may choose to add additional cross sections, change the number of branches, etc. However, the user is cautioned to consult this section before making such changes. Since all segments have the same

format, only one will be explained here.

ohmon.geo Block 1, branch and grid information

Example:

24 16

Variable Type and Column

npe length number Dscription

nbench Inte - Number of branches

Idntm Inte - Number of velocity interpolations between cross sections

in a stream tube

38
Card 2 (one number for each branch)

Example:

12 22 32 44 45 64 89 104 105 121 131 138 139 141 150 208 210 218 228 233

Variable Type and Column

name length number Description

lcstsq(i) Integer - Last cross section in each branch. In the last branch, use the

second to last cross section. There must be a total of nbranch

data. If the line is not long enough, continue to the next line

ohmon.geo: Block 2, cross section and connection information

Card 1 (one card for each cross section)

Example:

1 16840.66 7741.73 -0.44720 0.89443 11 11 2 0 0.0 0

Variable Type and Column

name length number Description

Integer - Cross-section number

x Real - x coordinate (in global coordinate system) of reference point

on left bank

y Real - y coordinate (in global coordinate system) of reference point

on left bank

cst Real - Cosine value that the angle cross section makes with positive x-axis (global)

snt Real - Sine value that the angle cross section makes with positive

x-axis (global)

nstube(i) Integer - Number of stream tubes at current cross section

numcon(i) Integer - If all stream tubes continue to next cross section undivided,

= 11; if stream tubes divided into two channels from a main

channel, = 12; if stream tubes from this channel and another

channel connect to next section that is in the main channel, =

21

nfirco(i) Integer Next cross section connecting to the current cross section.

For a divided channel around an island, this represents the

first cross section connected to in the lower division from

the main channel cross section

nseco(i) Integer For a divided channel around an island, this represents the

first cross section connected to the upper division from the

main channel cross section (if no island, = 0). If this is the

first cross section in the upper branch, = 888. If this is the last

in the upper branch, = 999.

Atub(i) Real If a branch has small islands, = discharge fraction on left side

of the island; if no island, = 0.0

lbend(i) Integer Flow direction index; for cross sections in a river bend, =1;

otherwise =0

39
ohmon.geo: Block 3, cross-section geometry

Card 1

Example:

6 39

Variable Type and Column

name length number Description

Integer - Cross-section number (for checking)

nslsct(j) Integer - Number of sounding depths (points along the bed in the

transverse direction) used to describe the cross-section geometry

Card 2 (as many cards as required to input all sets of ywid.z)

Example:

0.0 660.0 540.0 660.0 540.0 596.0 ...

Variable Type and Column

name length number Description

ywid(ij) Real - Transverse distance from the reference shore to the]h sounding depth in ith cross section

z(i,j) Real - Bed elevation at the jth transverse point (sounding depth)

for the jth transverse point (sounding depth) for the jth cross

section

Note: Block 3 must be repeated lcstsq(nbmch) times (Le. number of cross sections defined).

ohmon.geo- Block 4, defines all grid boxes that are to be assigned zero velocity

Card1

Example:

0

Variable Type and Column

name length number Description

nzrvb Integer Number of boxes to assign zero velocities

Card 2

Variable Type and Column

name length number Description

izrbx Integer - x-grid number of ith box to have zero velocity

izrby Integer - y-grid number of AtI box to have zero velocity

Note There are must be nzrvb pairs of izrbx(i) and izrby(i). Data can be continued to as many lines

asneeded.

For the other river segment, just repeat block 1 to block 4. This is the end of data file ohmon.geo.

a c•datan A04 oiMojee

The rhmm file contaifs information identifying ice regions. An ice region is defined by a range of

gWd boes that are eWovd; using the example in Figure 25 an ice region may be identified as extenh* fton rid (14,8) to grid (18,7). The ice region then covers every grid from (14,8) to (18,7) as

"shown in the figure. This tnfomatiorn is used when determining if spreading and advection takes place

under an ic cover or on open water.

40
17

161

m is, o _ Ice Region (14 8) to (18. 7) 14, ,

131

00512.5e eio

-- • c e R g i on2 • "( 2 3 , 7 ) t o ( 2 3 , 1 1 )

10

7

6

e20 25 Figure 25. Defining ice regions.

ohmon.ice: Block 1

Card 2

Example:

0235 12.5

Variable Type and Column

name length number Description

anice amiuo Real Real - Manning's n for ice roughness - Viscosity of oil (gm-cm/s) (poise)

Example:

2

Variable Type and Column

name length number Description

nicerg Integer - Total number of ice regions

Example:

14 8 18 7

Variabl Typ and Column

name ltengt number oecretion

nicexl(i) niceyl(i) Integer Intege - x grid at the beginning of ice region - y grid at the beginning of ice region

ncey2i) ntWy'(i Integer - x grid at the end of ice region Integer - y grid at the end of ice region

This is tie end of data file ohmon.ice.

Rmowdatad hn=Ao

The Ohnf file contains the unsteady flow time step, the water level and the discharge at each

nkode kt ft river as defined by the onedimaesional unsteady flow model (Yang et al. 1990). Also in41
cluded are the ice conditions for each cross section in the river. These data are separated from the ice

region data in ohmon.ice. The oil spill simulation model converts this information into boundary conditions for each river branch.

This file consists of three blocks of information. All blocks are listed below with descriptions and

corresponding components. Blocks 2 and 3 must be repeated every time when the velocities and depths

are updated in the model, i.e. every time step of the one-dimensional flow model. Therefore, the data in

this fie need to be adjusted on a more regular basis.

If the format or column numbers are shown, the data must be formatted accordingly. Otherwise it '•

necessary to have only one space or comma between the data.

For the Ohio-Monongahela-Allegheny river system, the difference in branch configuration between

the oil spill model and the unsteady flow model is built into the oil spill model through subroutine

ohndc. The number of branches in this data file is the same as the number of branches in the unsteady

flow model. The discharges and water levels at the nodes of the branches of this data file should be

interpolated to the cross sections of the oil spill model. If the two models have the same number of

branches, then subroutine ohndc should be changed to the following three lines, which are already in

the program. These lines are presently given as comments, which can easily be activated:

do 230 i=l, nsect

read(7,*) wlsct(i), qsct(i)

230 continue

ohmon.flo: Block 1, time step in one-dimensional model

Qr 1d

Example:

6.0

Variable Type and Column

name length number Description

ufdt Real - Unsteady flow time step (hr)

Qtd 2_(Comments line)

Example:

branch qu qd hu hd region 11/1/88 hour=0.0

These are the variable name, date and time of the data to be given in the following block.

ohmoflo. Block 2, discharge and water level

Q& I (one card for one branch)

Example:

1 52810. 52810. 625.19 625.00

2 51761. 52810. 625.25 625.19

Variable Type and Column

name lenth number Description

i n - Branch number for checking

qut) Real - Discharge at upstream end of each branch (cfs)

qd(M Real - Discharge at downstream end of each anch (ce)

wu(l) Reel - Water level at upstream end of each branch (ft)

wd(i) Rel - Water level at upstream end of each branch (ft)

Noe. Card 1 needs to be repeated as many times as the number of brandwb of this river segment.

42
ohmon.flo: Block 3, ice conditions (thickness information)

Card 1

Example:

1

Variable Type and Column

name length number Description

icinfo Integer - Nnmber of cros, sections with ice-cover-d condition. If no

ice-covered section exist, set icinfo =1 and then in card 2, define the section to be open

Card 2

Example:

2 open

Variable Type and Column

name length number Description

is i4 ID number of the cross section with ice-covered conditions

word a4 Cross-section ice cover condition, "full" = fully covered, "part" = partially covered, "open" = open water. If word =

"full," then the card has only one value and that is the ice

thickness across the river for that x-section. If word = "part,"

then card 3 must have an ice thickness defined at each vertical line defining the x-section = nslsct(is). If all numbers don't

fit into one card, as many cards as necessary may be u--ed.

Card 3 (for fully covered cross section)

Card 3 must follow every card 2 if word = "full" or "part." If word = "open," card 3 is not needed.

Example.

0.6

Variable Type and Column

name length number Description

fulti Real - Ice thickness (ft) of fully covered cross section. Only one value is read as input, and it will be assigned to the entire cross

section

Card 3 (for partially covered cross section)

Example:

0.3 OA 0-5 0.5 0.4 0.2

Variable Type and Column

name length number Description

tice(ij) Real - Ice thickness (ft) of partially covered cross section. There must

be one value for each sounding depth of the cross section

with one additional value for the extreme left of the cross

section where depth sounding is not input through data because it is always assumed to be zero

Note: Block 2 followed by 3 must be repeated for every one-dimensional unsteady flow time step.

"This is the end of data file ohmon.flo.

43
Boundary grid data: ohmon.igr

ohmon.igr Block 1, boundary grid boxes

Card 1

Example:

343 121 130 125 125

Variable Type and Column

name length number Description

i i4 - x-grid box i umbwr

igrilb(i) i5 - y-direction grid box number of lower river boundary for ith

x-grid (water-side grid)

igriub(i) i5 y-direction grid box number of upper river boundary for ith

x-grid (water-side grid)

igrlbj,i) i5 y-direction grid box number of lower boundary of jth island

for t1h x-grid (land-side grid)

igrub(j,i) i5 y-direction grid box number of upper boundary of jth island

for it*h x-g•id (land-side grid)

Note: Block I must repeat ntogridx (no. of grids in x-direction) times.

This is the end of data file ohmon.igr.

Shoreline boundary condition data: ohmon.bnd

The ohmon.bnd file consists of three blocks of information AX blocks are listed below with components and description. Most of the data read into the model are in free format. If the format or column

numbers are shown, the data must be formatted accordingly. Otherwise it is necessary to have only

one space or comma between the data. All fies with the bnd extension follow the same format.

ohmon.bnd: Block 1, half-life index.

CArd I (one card for each range of grid boxes)

Example:

1 1 81 1

Variable Type and Column

name length number Description

k Integer Shore number. Lower y (river) = 1; upper y (river) = 2; lower

y (first island) = 3; upper y (first island) = 4; etc.

Ifrom Integer Beginning limit (grid box number) for half-life designation

to shore

Ito Integer Ending limit (grid box number) for half-life designation to

shore

icode Integer Integer identifying which of the 16 half-life values are to be

assigned to a grid; for determining the value for lock and

dams, refer to Figures 15 and 16

ard lkt (must be included)

0 0 0 0

This is the end of data file ohmon.bnd.

O0l spill dabi, c&uwnopl

TMw dman Mfle cosmigs of two blocks of infomation with a varying or constant number of cards

kt eh block& AU blocks are td below with the components and description. Most of the data read

44
into the model are in free format. If the format or column numbers are shown, the data must be formatted accordingly. Otherwise it is necessary to have only one space or comma between the data. All files

with the spl extension follow the same format.

ohmon.spl: Block 1, Oil spill and simulation parameters.

Card 1 (river name-for identification only)

Example:

Ohio Monongahela

Variable Type and Column

name length number Description

text a20 - Any text to identify the river system

Card a (Type of oil, for identification only)

Example:

Fuel No. 2

Variable Type and Column

name length number Description

fueltp a20 - Text for identifying the oil type

Card 3

Example:

24.0 4 1 1 1 0 1800.0 -1.0

Variable Type and Column

name length number Description

totime Real - Time of oil spill simulation (hr)

ievery Integer - Frequency of obtaining output from simulation results from

ohmon and plotter. For example, if ievery = 4, output is written every four time steps

ipoti Integer - Two options: one (1) results in output of fixed data, such as

cross-section geometry data and shoreline conditions; zero

(0) results in no output

ipot2 Integer - Two options: one (1) results in output of computed velocities to a data file for plotting; zero (0) results in no printout

ipot3 Integer - Two options: one (1) results in output of particle location to

a data file for plotting; zero (0) results in no printout

ipot4 Integer - Two options: one (1) results in output of number plot of particle distribution to a data file for plotting; zero (0) results in

no printout

spitim Real - Duration of oil spill in seconds. If = 0, the spill will be treated as instantaneous; otherwise it will be continuous

diffur Real - Horizontal diffusion coefficient (ft2

/s) for river. If the default

formulation is derired, set this value to -1.0

45
Card 4

Example:

400 705000. 900. .84 1.411e-5 7.55e-4 1.14 1.1 1.6 1.4 1.4 1.4

Variable Type and Column

name length number Description

ntotal Integer - Number of parcels defined in each layer at the beginning of the simulation (currently 400)

spvol Real - Total volume of oil spilled (gal.)

spgoil Real - Specific gravity of oil

aniu Real - Kinematic viscosity of water (ft2/s)

sigma Real - Surface tension of oil (lb/ft)

ak2i Real - Fay's gravity-inertia phase spreading coefficient (axisymmetrical)

ak2v Real - Fay's gravity-viscous phase spreading coefficient (axisymmetrical)

ak2t Real - Fay's surface tension-viscous spreading coefficient (axisymmetrical)

akcl0 Real - Fay's or Waldman's gravity-inertia spreading phase coefficient (onedimensional)

akc20 Real - Gravity-viscous phase spreading coefficient (one-dimensional)

akc3O Real - Surface tension-viscous phase spreading coefficient (one-dimensional)

Card5

Example:

342650. -97800. .7063e-2 .1873e-27.88 465.0

Variable Type and Column

name length number Description

spx Real - x-coordinate of spill site (ft)

spy Real - y-coordinate of spill site (ft)

vmuni Real - Molar volume of oil (ft3/mol)

soluni Real - Solubility of fresh oil (lb/ft3)

cevp Real - Coefficient C of evaporation characteristics of oil

toevp Real - Boiling point temperature of oil (K). If you define a value of less

than 1.0 for toevp, the program will define the evaporation characteristics, using fitted curves. Therefore, the input vaiues of cevp and

toevp have no influence on computations, although they are read

ohmon.spl: Block 2, wind and temperature data.

Qs 1 (1 card for each time step

Example0.2 11.80 21.50 24.00

Variable Type and Column

name length number Description

vwc Real Empirical coefficient for modifying the wind speed on land to that

on river

vwmag Real - Wind speed (ft/s)

tdeta Real - Wind direction (angle measured clockwise from the north in de- gree)

tetwf Real Air temperature (OF)

NoWe Cad 5 must be repeated for every time step of the simulation.

This is the end of data file ohmon.spl.

4'
Bridging data between the unsteady flow model and the oil spill model: ohmon.ndc

The ohmon.ndc file consists of two blocks of information with a varying or nonvarying number of

cards in each block. All blocks are listed below with the components and description. Most of the data

read into the model are in free format. If the format or column numbers are shown, the data must be

formatted accordingly. Otherwise it is necessary to have only one space or comma between the data.

All files with the ndc extension follow the same format.

ohmon.ndc: Block 1, stream-tube cross-section distribution

Card 1 (Comments)

Example:

Region 1: part 1

Variable Type and Column

name length number Description

dummy a4 - Dummy variable for skipping a line in the data file

Card 2

Example:

5 1 45

Variable Type and Column

name length number Description

Integer - Branch number in this river segment

ns Integer - Number of stream-tube cross sections in this branch

jar(j) Integer - List of numbers of all stream-tube cross sections in the branch

Note: Card 2 must repeat for all branches in this segment.

ohmon.ndc: Block 2, location of stream-tube cross sections

ard 1

Example:

2 2.11 1 0.09

Variable Type and Column

name length number Description

j Integer - Cross-section number

ihup Integer - Number of upstream unsteady flow node

xup Real - Distance between cross section and up.tream node

ihdp Integer - Number of downstream unsteady flow node

xdp Real - Distance between cross section and downstream node

Note: Card 2 must repeat as many times as the number of cross sections in this river segment.

Repeat blocks I and 2 as many times as the number of river segments. This is the end of data file

ohmon.ndc.

Mis/elkneous data- ohmon.msc

"The ohmon.msc file consists of two blocks of information with a varying or nonvarying number of

car& in each block All blocks are listed below with the components and description. Most of the data

read fto the model are in free format If the tormat or column numbers are shown, the data must be

format[ aordingly. Odtewise it is necessary to have only one space or comma between the data.

All fies with the mwc extuen follow the same format

47
ohmon.msc: Block 1, some miscellaneous data.

Card.1 (Bank rejection options)

Example:

inst

Variable Type and Column

name length number Description

reJp a4 - Bank rejection parameter; for instantaneous rejection, =

"inst;" if = "full," then the oil deposited on the shoreline during each time step will be rejected completely at the end of

the time step; for other rejection rates, = "part"

Card 2 (Output options)

Example:

0

Variable Type and Column

name length number Description

ipropt Integer Two options: larger than 1, output grid locity, stream-tube

velocity to data files for plotting; less than 1, no output to

data files

Card 3

Example:

1.0

Variable Type and Column

name length number Description

*ha Real - Resurfacing coefficient, a; if all particles that rise to the surface layer stay in surface layer, = 1; otherwise, less than 1.0

Card 4

Example:

Le-5

Variable Type and Column

name length number Description

beta Real - Deposition coefficient

Card 5

Example:.

0.001

Variable Type and Column

name length number Description

gama Real - Emulsification coefficient (1/s)

Card6

Example:

0.028

Variable Type and Column

e ngth number Description

amnwn Real - Mann s coefficient of riverbed

48
Card 7

Example:

32.2

Variable Type and Column

name length number Description

gravac Real - Gravitational acceleration (ft/s 2)

Card 8

Example:

0.01

Variable Type and Column

name length number Description

vbuoy Real - Buoyant velocity (ft/s)

Card 9

Example:

1.1

Variable Type and Column

name length number Description

alpc Real - Coefficient for converting depth-averaged velocity to surface

velocity, a,

Card 10

Example:

0.03

Variable Type and Column

name length number Description

alpw Real - Wind effect coefficient on water surface velocity, a,

Card 11

Example:

05 05

Variable Type and Column

name length number Description

ratsur Real - Fraction of oil spilled on water surface

ratsus Real - Fraction of oil spilled in the suspended layer

Card 12

Example:

0.3

Variable Type and Column

name length number Description

alevp Real - Coefficient for modifying evaporation due to emulsification,

aE

449
Card 13

Example:

1000 200.

Variable Type and Column

name length number Description

ntxgrid Integer - Total x-grid of moving grids

dx Real - Grid size (ft)

Card 14

Example:

5

Variable Type and Column

name length number Description

nr Integer Total number of river segments

ohmon.msc: Block 2, coordinate transformation information.

Card 1

Example:

129.

Variable Type and Column

name length number Description

sO Real - River mileage of the downstream boundary of the study

reach

Card 2

Example:

0.0 0.0 348200.0 144000.

Variable Type and Column

name length number Description

xmin(i) Real - Minimum x-coordinate of ith river segment

ymin(i) Real - MIindmum y-coordinate of ith river segment

xmax(i) Real - Maximum x-coordinate of ith river segment

ymax(i) Real - Maximum y-coordinate of ith river segment

Cold I

Example:

0.0 0.0 0.0

Variable Type and Column

name length number Description

x0(i) Real - x-origin (coordinate) of the iAh local coordinate system

yO(i) Real - y-origin (coordinate) of the ith local coordinate system

beta0(i) Real - Rotation angle of the ith local coordinate x positive with global x positive; "-" for clockwise and "+" for counterclockwise

50
Card 4

Example:

1741 46.

Variable Type and Column

name length number Description

ngrid(i) Integer - Maximum x-grid number in the ith river segment

sbm(i) Real - River mileage at the end of the ith river segment

Note: Repeat cards 3 and 4 times as many times as the number of river segments.

This is the end of data file olhmon.msc.

Input adjustments

For a river that already has the necessary input files, very little modification in input is needed to run

the model with different field conditions, such as new river discharge, different oil properties, new spill

location, etc. Input data that are most likely to require modification are cited below with some guidelines and suggested values. Readers should refer to the preceding section for the format of the data.

Velocity

Interpolating stream-tube velocities to boxes depends on the variable kintm (ohmon.geo, Block 1, Card

1) and ibend(i) (ohmon.geo, Block 2, Card 1). As kintm increases, more velocity interpolations are made

between cross sections. If the cross sections are spaced far apart, kintm should be in the range of 5-8 or

more; otherwise kintm can be 5 or less. The problem with using too small a value for kintm is that

more grid boxes will end up without assigned

velocities after the interpolation stage. They will

be assigned a velocity based on their neighboring boxes. This is less accurate than assigning a

velocity based on interpolation. One should try

to assign as many boxes as possible in the interpolation stage. On the other hand, if the value of

kintm is too large, it will take more computation

time. Besides the value of velocities, appropriate

simulation of the direction of velocity is also important. Ibend(i) supplies two options for determining velocity directions. When ibend(i) = 0, the

velocity will be along the line connecting two corresponding stream-tube centers in neighboring -lo

cross sections. This usually gives a fairly accurate

flow direction. However, for cross sections in a

river bend, as shown in Figure 26, ibend(i) should

be 1. This makes the velocity direction perpendicular to the cross section.

Another factor that affects assigning velocities

to grid boxes i the number of stream tubes se- itend. 1

lected. The choice of the number of stream tubes --- ib . 0

(nstube, ohmon.geo, Block 2, Card 1) will affect the

computer time and the accuracy. For a river with

highly irregular cross sections or when a high

degree of accuracy is required, more stream tubes

should be used. Figure 26. Simulation offlow direction in stream tubes.

51
Stage-discharge data

Stage and discharges (ohmon.flo, Block 1, Card 2) in river branches are unlikely to remain constant

over time. Therefore, this information will require updates as the need arises. The one-dimensional

unsteady flow model can be used to obtain this information.

Oil spill parameters

The oil spill parameters (ohmon.spl, Block 1, Cards 2,3,4,5) will require the most modification from

spill to spill. Each spill will have its own simulation duration, simulation time step, printing options,

number of particles, volume, physical properties and location, so it will be necessary to adjust these

data each time the simulation is done for a new spin. The suggested guidelines for these paranmeters

follow.

The initial total particle number, ntotal, is about 500; a maximum of 1000 is possible in the current

version. During the simulation the number may increase to several thousand because of the exchange

between the two layers. A larger particle number will give a better result at the expense of a longer

execution time. The time step suggested for spildt is 900 s (15 minutes).

For other spill parameters the following values may be used:

spgoil = 0.7-0.98

aniu = 1.411 x 10-5 ft2 /s (at 50 0F)

sigma = 206 x 10-5 (Ib/ft)

ak2i= 1.14

ak2v 0.98

ak2t= 1.6

akcl0 = 1.39

akc20 = 1.39

akc30 = 1.43

vmnui= 0.7063 x 10- 2 (ft3 /mol)

soluni = 0.1873 x 10-2 (b/f,).

Wind data

The computer model reads the wind data (ohmon.spl, Block 2, Card I) at each time step. This allows

for changing the wind speed and direction in the simulation. A correction is needed because the wind

speed on land and that on the river are different. If no better information is available, a value of 0.2 may

be used for vwc (ohmon.spl, Block 2, Card 1).

Sample input data files

All file names to be opened for input and output will be in the file ohmon.fnm. Sample input data files

included in this section are:

1. ohmon.geo

River geometry, cross sections and branch geometry of the Ohio-Monongahela rivers.

2. ohmon.igr

Boundary grid data for the Ohio-Monongahela rivers.

3. ohmon.ice

Ice parameters, ice regions and ice thickness in the Ohio-Monongahela rivers. This file is needed

even when the model is run for open-water conditions. Samples will be given for the ice-covered

conditio and the open-water condition.

4. ohmoPnflo

Water level and discharge at nodes from the unsteady flow model for the Ohio-Monongahela

rivemr

5. o•nmo.bnd

Half-life index for shoreline grids of the Ohio-Monongahela rivers.

52
6. ohmon.spl

Oil parameters, spin location, wind condition and ambient air temperature.

7. ohmon.ndc

Information for converting unsteady flow computation results to what are needed in oil spill model

flow computation for the Ohio-Monongahela rivers.

8. ohmon.msc

River division data, exchange parameters for the two layers and information for coordinating transformation data for the Ohio-Monongahela rivers.

ohmon.fnmThu Jan 24 08:59:15 1991 1

ohmon. fnm:

lohxnon geo

2oilprt.out

7ohmon. flo

8ohmon. bnd

9ohmon. ndc

llohmon. sur

13ohmon. sus

12ohmon. spl

14ohmon. ice

20ohmon .msc

22concen.dat

29ohmon. igr

28vel.pos

30dep.pos

ohmon.geo:

24 16 region 1

12 22 32 44 45 64 89 104 108 121 131 138 139 141 150 165 177 194 204 208

210 218 228 233

1 16840.66 7741.73 -0.44720 0.89443 11 11 2 0 0.0 0

2 18925.39 7018.47 -0.44720 0.89443 11 11 1 0 0.0 0

3 20144.21 8923.44 -0.56714 0.82362 11 11 2 0 0.0 0

1 10

0.0 660.0 540.0 660.0 540.0 596.0 600.0 596.0 650.0 596.0

700.0596.0 800.0 596.0 900.0 596.0 900.0 660.0 1274.0 660.0

2 20

0.0 665.0 1200.0 665.0 1280.0 660.0 1380.0 625.0 1659.0 625.0

1800.0616.0 1900.0 594.0 2200.0 594.0 2420.0 580.0 2770.0 581.0

2940.0583.0 3080.0 590.0 3140.0 590.0 3240.0 606.0 3300.0 606.0

3380.0633.0 3830.0 635.0 3870.0 640.0 4030.0 645.0 4500.0 665.0

0

23 5 region 2

53
ohmon.igr Thu Jan 24 09:01:241991 1

ohmon.igr

1 15 19

2 15 19

3 15 19

4 15 19

5 14 19

6 14 19

7 14 19

8 14 19

9 14 19

10 14 19

11 14 19

12 14 19

13 14 19

14 14 19

15 14 19

16 14 20

17 14 20

18 14 20

19 14 20

20 14 20

21 14 21

22 15 22

23 15 22

24 15 23

25 15 23

26 16 23

ohmon. ice:

0.035 0.84

0

0.035 0.84

0

0.03S 0.84

0

0.035 0.84

0

0.035 0.84

0

0.035 0.84

0

0.035 0.84

0

0.035 0.84

0

0.035 0.84

0

0.035 0.84

0

0.035 0.84

54
0

0.035 0.84

0

0.035 0.84

0

ohmon.floThuJan 24 09:03:01 1991

ohmon. flo:

6.0 flow time step

branch qu qd hu hd region 1 1/ 1/88 hour-0.

1 52810. 52810. 625.19 625.00

2 51761. 52810. 625.25 625.19

3 50721. 51761. 625.37 625.25

23 48876. 48878. 669.09 668.89

24 48877. 48876. 669.30 669.09

1

2 open

branch qu qd hu hd region 2 l/ 1/88 hour = 0.

1 48877. 48876. 669.30 669.09

2 48860. 48877. 669.45 669.30

3 48850. 48860. 670.50 669.45

4 48460. 48850. 670. 87 670.60

ohmon .bnd:

1 1 81 1

2 1 81 2

1 82 84 2

2 82 84 2

3 82 84 2

4 82 84 2

1 85 34 2

2 85 34 2

1 3433582

2 3433582

3 343 358 2

4 343 358 2

1 359 530 1

2 359 530 1

1 531 532 2

2 531 532 2

3 531 532 2

4 531 532 2

1 5335361

2 5335361

1 5375382

2 5375382

3 5375362

4 5375382

1 5397401

2' 5397401

55
ohmon.splThuJan 24 09:04:27 1991 1

ohmon.spl:

Ohio Monongahela

No. 2 Fuel oil

240.0 288 1 1 1 03600. -1.0 totime, ievery,

400 705000.0 300. 84 1.411E-5 7.55E-4 .14 1.1 1.6 1.4 1.4 1.4

342650.-97800..7063E-02.1873E-02 7.88 465.0

0.2 11.80 210.50 24.00

0.2 4.21 100.50 20.00

0.2 0.00 0.00 17.50

0.2 2.53 40.50 15.50

0.2 2.53 40.50 14.50

0.2 0.00 00.00 18.00

0.2 0.00 00.00 25.00

0.2 3.37 60.50 29.50

0.2 3.37 60.50 30.50

0.2 2.53 90.00 29.50

0.2 7.59 200.50 28.50

0.2 10.96 230.00 28.00

0.2 14.33 230.SO 27.50

0.2 28.66 240.50 27.00

0.2 36.25 250.50 27.50

0.2 37.93 280.00 25.50

0.2 33.72 280.$0 19.50

0.2 26.13 270.50 13.50

0.2 25.29 280.00 9.00

ohmon.ndc:

Region 1:part 1

1 12 1 2 3 4 5 6 7 8 9 10 11 12

2 10 13 14 15 16 17 18 19 20 21 22

3 10 23 24 25 26 27 28 29 30 31 32

4 12 33 34 35 36 37 38 39 40 41 42 43 44

5 1 45

22 8 211 212 213 214 215 216 217 218

23 10 219 220 221 222 223 224 225 226 227 228

24 6 229 230 231 232 233 234

Next we have located x-section

1 2 2.11 1 0.09

2 2 2.00 1 0.20

231 25 2.50 24 1.00

232 25 2.04 24 1.46

233 25 1.50 24 2.00

234 25 1.02 24 2.48

Region 2: part 1

1 8 1 2 3 4 5 6 7 8

56
ohmon.msc Sun Jan 20 16:42:57 1991 1

part full for rejrat = 1 ,inst for inst. rej., part for calculated;

0 =ipropt, >1 output vel. and depth in grid, opposite in inverse.

1. =alpha, resurfacing

1.e-5 =beta, deposition

0.00001 =gamma, emulsification

0.028 =amngn

32.2 =gravac

0.01 =vbuoy

1.1 =alpc

0.03 =alpw

0.8 0.2 =ratsur, ratsus

0.2 =alevp, evap. modified coeff.

1000 200. =ntxgrid, dx the x-grid number of working length

5 =nr

129. sO

0. 0. 348200. 144000. xmin,ymin, xma.x,ymax

0. 0. 0. xO,yO,betaO

1741 46. ngrid(1), sbrn(1)

348200. -36800. 432000. 139200. xmin,ymin,xmax,yymax

348200. 139200. -90. xO,yO,beta0

2621 0.0 ngrid(2), sbrn(2)

378400. -84800. 400000. -36800. xmin,ymin,xmax,ymax

378400. -36800. -90. xO,yO,betaO

2861 -11.5 ngrid(3), sbrn(3)

315800. -132000. 392400. -84800. xmin, ymin, Xmax, ymax

392400. -84800. -180. xO,yO,beta0

3244 -40.5 ngrid(4), sbrn(4)

52000. -220000. 315;00. -104600. xmin,ymin,xmax,ymax

315800. -104600. -180. xO,yO,betaO

4563 -132. ngrid(5), sbrn(5)

2288 3288

Additional input

Besides the input data files mentioned above, there are two additional sets of data to be input on the

screen after the start of the simulation.

Rier combinations

The model is designed to simulate spills that occur either in the Monongahela River or in the Allegheny River. After the user starts running the oil spill model, the screen will show:

*********** ***************************

"* select river combination: *

"* for Ohio-Monongahela River input: 1 *

"* for Ohio-Allegheny River input: 2 *

Choose 1 for an oil spill in the Monongahela River or 2 for a spill in the Allegheny River. If the spill is

in the Ohio River, choose either 1 or 2. However, 2 will make take less computing time, since the data

files on the geometry and shoreline boundary conditions are smaller.

57
Continuation runs

The model can be used either for a simulation from the beginnirig of a spill or for continuing a previously stopped simulation. The user will be required to choose from these two options after choosing

the river combination. The screen will show the following:

* The model can be run from t = 0 or *

* restarted after a stop *

* Enter 0 or 1 and hit return

* 0 for a run starting from t =0 *

* 1 to start from a previously stopped point *

MODEL OUTPUTS

Output from the model can be directed to different devices, such as a console, a printer or files.

Velocity and depth distributions

Flow velocity and depth, both along the stream tubes and in the grids, are calculated in subroutine

veldep. If needed the user can choose the printout option (ipropt> 1, ohmon.int, Block I, Card 2) to save

all simulation results into corresponding files.

"* strm.veh: Save coordinate (xy) and velocity (v., vy) at stream-tube locations.

"* strn.dep: Save coordinate (xy) and depth (h) at stream-tube locations.

"* grid.velh Save coordinate (xy) and velocity (vx, vy) at grid box locations.

"* grid.dep: Save coordinate (xy) and depth (h) at grid box locations.

All of these data are used in advection and mixing calculations. They cannot be saved in the microcomputer version because of the limitation on the number of files to be opened. The arrangement of

data files is shown as follows:

Data card of velocity for grids and stream tubes:

x y v. vY

Data card of depth for grids and stream tubes:

x y h

where coordinates and depth are in feet and velocity is in feet per second.

The cards for velocity and depth are repeated as many times as the distribution points. In this model, only the velocities and depths within the moving grid at each flow time step will be saved into files

to save memory and computer time.

Simulation results

Oil spill simulation results will be output to files on the selected output frequency (ievery, ohmon.spl,

Block 1, Card 2). These results include the following files:

*oilprt.out: Spill site and oil parameters. A summary of distributions of the oil, including updated

total volume (for checking purposes), volume in the surface layer, volume in the suspended layer,

total volume evaporated and oil evaporated in the current time step, total volume dissolved and

volume disolved during this time step, volume on land, volume deposited onto the riverbed, and

volume emulified and resurfaced.

*o•mon.sur: Locations and volumes of surface oil particles, including updated total particle number

and the center of slick

f 58
ohmon.sur: Block 1

Card 1

Example:

Ohio Monongahela Surface No. 2 Fuel Oil

Variable Type and Column

name length number Description

text a20 11 River name

sur a14 31 Indicate results of surface slick

fueltp a20 45 Name of spilled oil

Card 2

Example:

881.25

Variable Type and Column

name length number Description

vave Average volume of particles both on surface and in suspension at the beginning for graphics purpose (gal.)

Card 3

Example:

400 705000. 900. 0.81 0.141E--04 0.755E-03 0.706E-02 0.187E-02 0.OOOE+00

564000.000 141000.000 342650 -97800

Variable Type and Column

name length number Description

ntotal Integer - Total number of particles defined in each layer (e.g. 400)

spvol Real - Total volume of oil spilled

spgoil Real - Specific gravity of oil

aniu Real - Kinematic viscosity of water (ft2/s)

sigma Real - Surface tension of oil (lb/ft)

vmuni - Molar volume of oil (ft3/mol)

soluni - Solubility of fresh oil (1b5

/f-t)

amuni - Absolute viscosity of oil (lb/ft2)

- Total volume of oil on water surface at the beginning of spill

(gal.) - Total volume of oil in suspension at the beginning of spill

(gal.)

59
ohmon.sur: Block 2

Card 1

Example:

342650 -97800

Variable Type and Column

name length number Description

ipartx(1) i8 - x-coordinate of the spill site (ft)

iparty(1) i8 9 y-coordinate of the spill site (ft)

Card2

Example:

2300 391145 -77512 24.0000 -1.6 -0.5

Variable Type and Column

name length number Description

nwtl Integer - Total number of surface oil particles at the present time

ipartx(1) i8 - x-coordinate of surface slick center (ft)

iparty(1) i8 9 y-coordinate of surface slick center (ft)

ttt Real - Completed simulation time (hr)

vwx Real - Present x component of wind speed (mph)

vwy Real - Present y component of wind speed (mph)

Card3

Example:

391163 -85810 538.35 391160 -85399 150.81 391185 -84711 163.11391370

-85965 132.45

This card has four sets of values [x,y, vopsur(i)] to indicate the coordinates and volume of each particle;

x and y in are in feet, and the volume of the ith surface particle vopsur(i) is in gallons. This card will

repeat for all particles.

Block 2 will be repeated as many times as desired output times.

*ohmon.sus: Locations and volumes of suspended oil particles, including the updated total particle

number and the spill center. The output format is the same as ohmon.sur.

econcendat: One-dimensional cross-section averaged concentration distribution (including the surface layer and the suspended layer), indicating the volume concentration at different river miles.

concendat: Block 1

Example:

& time, JO, disO= 24.0000 205933.9545

where & is for a graphics purpose; time is the total length of the completed simulation (hr); iO is the

intl x-grid number of the moving grids; disO is the river mile of the downstream edge of the moving

g"id (mile).

60
Card 2

Example:

2.827 0.5278E-02

Varýible Type and Column

name length number Description

dis Real River mile at which cross-section-averaged volume concentration is calculated

conc2 Real Volume concentration at river mile distance (10-3)

Block I will be repeated as many times as the designed output number.

concen.sur: One-dimensional surface concentration distribution, indicating the volume concentration

at different river miles along the river. This file will not be available in the output for a microcomputer. The format of this file is the same as concen.dat.

The graphics output that corresponds to these outputs can be obtained by using the graphics program.

Sample output files

oilprt.out WedJan 30 11:53:42 1991 1

oilprt.out:

1 Ohio Monongahela

**** * **** * *** **

* continuous spill*

*at *

*342650.,-97800. *

*for 60. min. *

simulation period = 240.0 hrs

characteristics of spill

no. of particles : 400

oil spilled : 705000. gals of No. 2 Fuel Oil

dt for spill simulation : 300. secs.

specific gravity of oil : 0.84 (api index = 37.0)

kinematic visco. of water: 0.1411E-04 sq ft/sec

surface tension : 0.7550E-03 lbs/ft

spreading coefficients

k2i k2v k2t clo c20 c30

1.14 1.10 1.60 1.40 1.40 1.40

molar volume : 0.7063E-02 cu ft/mol

solubility of ZzasL oil : 0.1873E-02 lbs/cu ft

viscosity of oil : 0.00 gms/cm-sec

Manning's roughness of ice:0.000

61
api option is not selected-evap. constants are c= 7.88t0= 465.0

surface diffusion-default formulation is used

time, iO, jend, nhitb=

8.00 2064 3064 1346

slick condition at the end of this time step

total volume (cu. ft) = 94244.5

total volume on land (cu. ft) = 44775.

total on surface = 46411. total in suspension = 714.53

evaporation (cu. ft) this step = 18.694total evap.(cu. ft) 2222.5

dissolved this step = 0.58905 total = 29.087

total resurfacing (cu. ft) = .12147E+06

total deposition (cu. ft) = 121.39

total emulsification (cu. ft) = 13582.

total out of range on surface (cu. ft = .66263E-02

total out of range in suspended (cu. ft) = .54000E-02

:oncen.dat ThuJan 24 09:12:37 1991 1

concen.dat:

6 time, iO, disO= 24.0000 2059 31.5455

0.4180 0.5278E-02

0.3801 0.1285E-03

0.3378 0.1198E-01

0.2152 0.6928E-04

0.1616 0.4378E-01

0.1193 0.1687E-01

0.7693E-01 0.3473E-01

0.2336E-01 0.1634

-0.1899E-01 0.1171

-0.7256E-01 0.3773

-0.1261 0.1599

-0.1944 0.9520E-01

-0.2627 0.7579E-01

-0.8132 0.7076E-01

-0.85C6 0.7076E-01

-0.9576 1.052

-1.060 0.9058E-01

-1.128 0.1246

-1.213 0.2245

-1.281 3.160

-1.869 1.859

-1.937 2.506

-2.006 1.563

-2.059 2.579

-2.113 4.183

-2.166 0.9436E-01

-2.274 0.4299E-03

-2.327 2.119

-2.3R1 0.1234E-04

-2.419 1.403

-2.472 0.1794E-01

62
-2.510 0.6395

-2.564 0.166-

-2.602 0.3221

-2.639 0.2070

-2.677 0.1014

-2.715 0.5485E-01

-2.753 0.6229E-01

-2.791 0.7658E-01

-2.833 0.1138

-2.871 0.4878E-01

-2.909 0.1256E-01

-2.947 0.2762E-01

-3.000 0.5685E-01

ohmon.surThu Jan 24 09:10:02 1991 1

ohmon.sur:

Ohio Monongahela Surface No. 2 Fuel Oil

881.250

400 705000. 300. 0.84 0.141E-04 0.75'E-030.706E-02 0.187E-020.000E+00

564000.000 141000.00

342650 -97800

3228 383560 -63631 24.0000 0.3 -1.7

385603 -55130 127.61 385022 -55737 151.35 384943 -55816 122.80 385406 -55514 137.42

3R5453 -55484 90.32 385005 -55790 130.41 384217 -56667 167.51 384269 -56384 102.60

384204 -56198 145.88 384541 -S5888 99.56 385486 -55416 225.97 384341 -56130 230.82

384681 -55980 115.19 384548 -56119 218.39 384121 -56382 91.33 384184 -56360 162.97

384164 -56405 92.04 384390 -56438 167.40 384514 -S6237 216.95 384947 -55902 149.75

384212 -56492 115.06 383893 -56521 114 22 384098 -56850 92.86 384082 -56556 170.92

383831 -56706 102.75 384108 -56482 108.71 384059 -56695 209.84 383949 -56903 500.66

383958 -56767 121.98 384161 -56739 100.85 383714 -57187 114.55 384053 -56653 131.79

383848 -57152 91.06 383700 -57115 100.12 383734 -56833 118.72 384327 -56450 119.88

384393 -56345 132.88 383628 -57093 152.93 383296 -57726 301.19 383217 -57443 321.61

385530 -55327 90.76 383568 -57324 304.25 384432 -56297 96.78 383300 -57414 192.12

ohmon.sus:

Ohio Monongahela Suspended No. 2 Fuel Oil

881.250

400 705000. 300. 0.84 0.141E-04 0.755E-030.706E-02 0.187E-020.OOOE+00

564000.000 141000.00

342650 -97800

324 383388 -63317 24.0000 0.3 -1.7

381100 -60059 107.54 380927 -62207 103.10 381289 -62940 118.05 381496 -63364 88.44

381612 -63445 93.64 383286 -64942 102.53 385900 -C'ýP0 1.47 385700 -55100 1.57

385500 -55300 1.23 385500 -55500 4.27 385100 -55700 6.67 385300 -55700 0.57

384500 -55900 1.30 384700 -55900 9.48 384900 -55900 1.15 384300 -56100 6.66

384500 -56100 9.24 383900 -56300 0.08 384100 -56300 11.60 384300 -56300 5.42

384500 -56300 1.51 383700 -56500 0.05 383900 --56500 3.79 384100 -56500 3.25

384300 -56500 1.67 383500 -56700 0.06 383700 -56700 5.42 383900 -56700 21.05

384100 -56700 3.49 384300 -56700 0.52 383500 -56900 9.22 383700 -56900 13.31

383900 -56900 5.42 384100 -56900 0.58 383300 -57100 20.31 383500 -57100 13.78

383700 -57100 10.03 383900 -57100 0.55 383300 -57300 21.70 383500 -57300 7.66

383700 -57300 3.99 383100 -57500 25.84 383300 -57500 18.08 383500 -57500 1.74

63
CONCLUSION

This report describes the model formulation and implementation of a two-layer model, ROSS2, for

simulating oil spills in rivers. The model considers the oil in the river as consisting of both the surface

slick and suspended oil droplets mixed over the depth of the flow. The oil transformation processes

considered in the model include adv'ection, mechanical spreading, turbulent diffusion and mixing, evaporation, dissolution, shoreline deposition and sinking.

The model can be used for instantaneous or continuous spills, either on or under the water surface

in rivers with or without an ice cover. Although it is developed for simulating oil spills, the model can

be applied to spills of other hazardous materials.

The model has been implemented for the Ohio-Monongahela-Allegheny river system and the upper St. Lawrence River. A case study for the spill on the Monongahela River in January 1988 is presented along with detailed explanations of the program structure and its input and output. The simulation

results compared well with field observations.

LITERATURE CITED

Berkey, E., ILL. Price, S.M. Creeger and W.M. Johnston (1989) Assessment of environmental effects from

the January 2, 1988 diesel oil spill into the Monongahela River. In Proceedings of Pittsburgh's Oil Spill,

Pittsburgh, Pennsylvania.

Blokker, P.C. (1964) Spreading and evaporation of petroleum products on water. In Proceedings, 4th

International Harbour Congress, Antwerp, p. 911-919.

Butler, J., J. Morris and T. Slecter (1976) The fate of petroleum in the open ocean. In Sources, Effects and

Sinks of Hydrocarbons in the Aquatic Environment. The American Institute of Biological Sciences, p. 287-

297.

Cheng, R.T., V. Casulli and S.N. Milford (1984) Eulerian-Lagrangian solution of the convectiondispersion equation in natural coordinates. Water Resources Research, 20(7): 944-952.

Coher, Y., D. Mackay and W.Y. Shiu (1980) Mass transfer rates between oil slick and water. The Canadian Journal of Chemical Engineering, 58(October).

Cox, J.C., and LA. Schultz (1981) The containment of oil spilled under rough ice. In Proceedings of the

1931 Oil Spill Conference, American Petroleum Institute, Washington, D.C., p. 203-298.

Fannelop, T.K., and G.D. Waldman (1972) Dynamics of oil slicks. American Institute of Aeronautics journal, 10- 506-510.

Fay, J.A. (1969) The spread of oil slicks on a calm sea. In Oil on the Sea (D. Hoult, Ed.). New York.

Plenum Press, p. 53-64.

Fay, J.A (1971) Physical processes in the spread of oil on a water surface. In Proceedings of the Joint Conference on Prevention and Control of Oil Spills, Washington, D.C., American Petroleum Institute, p. 463-

467.

Fingas, M., and M. Sydor (1980) Development of an oil spill model for the St. Lawrence River. Inland

Waters Directorate, Water Management and Planning Branch, Technical Bulletin No. 116.

Fischer, H.B., E.J. List, R.C.Y. Koh, J. Imberger and N.H. Brooks (1979) Mixing in Inland and Coastal Waters. New York: Academic Press.

Gundlach, E.R., and M.D. Hayes (1978) Vulnerability of coastal environments to oil spill impacts. Marine Technological Society Journal, 12(4): 18-27.

Hamilton, G.D., CW. Soileau and A.D. Stroud (1982) Numerical modeling study of Lake Pontchautrain. Journal of the Waterways, Ports, Coastal and Ocean Division, ASCE, 108: 49-64.

Hoult, D.P. (1972) Oil spreading on the sea. Annual Review of Fluid Mechanics, p. 341-367.

Hoult, D.P., and W. Suchon (1970) The spread of oil in a channel. Department of Mechanical Engineerin& Massachusetts Institute of Technology, Cambridge, Massachusetts.

Hoult, DY., S. Wolfe, S. O'Dea and J.P. Patureau (1975) Oil in the Arctic. Massachusetts Institute of

Techtology, Cambridge, Massachusetts, Report No. CG-D-96-75.

64
Huang, J.C. (1983) A review of the state-of-the-art of oil spill fate/behavior models. In Proceedings of the

1983 Oil Spill Conference, American Petroleum Institute, Washington, D.C., p. 313-322.

Huang, J.C., and F.C. Monastero (1982) Review of the state-of-the-art of oil spill simulation models. Final

Report submitted to the American Petroleum Institute, June.

Hung, S. C. (1991) Seaway emergency response plan. In Proceedings, 1991 Oil Spill Conference, San Diego.

Johnson, B. H. (1982) Development of a numerical modeling capability for the computation of unsteady

flow on the Ohio River and its major tributaries. USA Waterways Experimental Station, Technical Report HL-82-20.

Larsen, PA. (1969) Head losses caused by an ice cover on open channels. Journal of Boston Society of Civil

Engineers, p. 45-67.

Leendertse, J.J. (1970) A water-quality simulation model for well-mixed estuaries and coastal seas. Volume I, Principles of computation. The Rand Corporation, Santa Monica, California, Report RM-6230-

RC.

Lu, C.Y., and J. Polak (1973) A study of solubility of oil in water. Environmental Protection Service,

Canada, Report EP5-4-EC-75-1.

Mackay, D., and C.D. McAuliffe 0988) Fate of hydrocarbons discharged at sea. Oil and Chemical Pollution, 5:1-20.

Mackay, D., S. Paterson and S. Nadeau (1980a) Calculation of the evaporation rate of volatile liquids.

Proceedings of the 1980 National Conference on Control of Hazardous Material Spills, Louisville, Kentucky, p.

361-368.

Mackay, D., S. Paterson and K. Trudel (1980b) A mathematical model of oil spill behavior. Environmental Protection Service, Fisheries and Environment Canada.

Miklaucic, F.A., and J. Saseen (1989) Oil spill, Floreffe, Pennsylvania--Case history and response evaluation. Proceedings, Oil Spill Conference, San Antonio, Texas.

Moore, S. F., R.L. Dwyer and AM. Katz (1973) A preliminary assessment of the environmental vulnerability of Machias Bay, Maine, to oil supertankers. Parson Laboratory, Massachusetts Institute of Technology, Report No. 162.

Roache, R.J. (1972) Computational Fluid Dynamics. Albuquerque, New Mexico: Hermosa Publishers.

Sayre, W.W., and FM. Chang (1969) A laboratory investigation of open-channel dispersion processes for

dissolved, suspended, and floating dispersants US. Geological Survey, Professional Paper 433-E, p. 71.

Shen, H.T., P.D. Yapa and M.D. Petrold (1990) Simulation of oil slick transport in Great Lakes connecting channels. USA Cold Regions Research and Engineering Laboratory, CRREL Report 90-1.

Shen, H.T., and P.D. Yapa (1988) Oil slick transport in rivers. Journal of Hydraulic Engineering, ASCE,

114(5): 529-543.

Shen, H.T., and N.L. Ackemnnan (1980) Wintertime flow distribution in river channels. Journal of the

Hydraulics Division, ASCE, 106(HY5): 805-817.

Stolzenbach, K.D., O.S. Madsen, E.E. Adams, AM. Pollack and C.K. Cooper (1977) A review and evaluation of basic techniques for predicting the behavior of surface oil slicks. Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, Report No. 22.

Spaulding, M.L. (1988) A state-of the-art review of oil spill trajectory and fate modeling. Oil and Chemical

Pollution, 4: 39-55.

Tsahalis, D.T. (1979a) Theoretical and experimental study of wind and wave induced drift. Journal of

Physia O graphy, 9: 1243-1257.

Tsadalis, D.T. (1979b) Contingency planning for oil spills: R1VERSPIL-A river simulation model. In

Proceedings of the 1979 Oil Spill Conference, Washington, D.C. American Petroleum Institute.

TorSrimscn, G.M. (1984) The on-scene-spill model: A user's guide. NOAA Hazardous Materials Response Branch, Technical Report (unpublished).

Vlcmy, AMH, and L Ahles-Kedzora (1989) Interstate in-stream tracking of the January 1988 Ashland oil

qpill to the Motongahela and Ohio rivers. In Proceedings of Pittsburgh's Oil Spill, Pittsburgh, Pennsylvania.

Waidman, GA.. RLA. Johnson and P.C. Smith (1973) The spreading and transport of oil slicks on the

65
open ocean in the presence of wind, waves, and currents. Coast Guard Report No. CG-D-1 7-73, AD763 926, July.

Weeks, W.F., and S.L. Dingman (1972) Thermal modification of river ice covers: Progress and problems.

Proceedings,The Role of Snow and Ice in Hydrology, Banff, p. 1427-1435.

Yang, X.Q., P.D. Yapa, H.T. Shen and D.S. Wang (1990) Unsteady flow simulation for the Ohio River

system. Department of Civil and Environmental Engineering, Clarkson University, Report 90-10.

Yapa, P.D., D.S. Wang, H.T. Shen, X.Q. Yang and J.B. Perry (1990a) An integrated oil spill model for

Ohio-Monongahela-Allegheny river system (User's guide). Department of Civil and Environmental

Engineering, Clarkson University, Report 90-11.

Yapa, P.D., H.T. Shen, X.Q. Yang, D.S. Wang and J.B. Perry (1990b) User's manual for an integrated oil

spill model for St. Lawrence River. Department of Civil and Environmental Engineering, Clarkson University, Report 90-7.

66
APPENDIX A: DEFINITION OF FORTRAN VARIABLES

The program was written in FORTRAN 77 and was tested on the f77 compiler. In this appendix,

important FORTRAN variables are given. For more detailed instructions on the usage of the program,

refer to the main text. Descriptions of the main program, subroutines and preparation of data files are

also given in the text. The program listing is provided in Appendix B.

Variable Symbol Type Description

akcl0 c10 real Fay's coefficient for gravity-inertial phase one-dimensional

spreading

akc20 c20 real Fay's coefficient for gravity-viscous phase one-dimensional

spreading

akc30 c30 real Fay's coefficient for surface tension- viscous phase one-dimensional spreading

ak2i k• real Fay's coefficient for gravity-inertial phase axisymmetric spreading

ak2t k2t real Fay's coefficient for surface tension-viscous phase axisymmetric

spreading

ak2v k2v real Fay's coefficient for gravity-viscous phase axisymmetric spreading

alevp SE real Coefficient for modifying evaporation due to emulsification (- 0.2-

0.5)

alpc sC real Weighting factor for the rivei current velocity to compute the surface drift velocity (- 1.1)

alpha a real Probability factor for re-entrainment of oil globules from the suspended layer to the surface layer

alpw aw real Weighting factor for the wind velocity to compute the surface drift

velocity (- 0.03)

amngn ii real Manning's roughness coefficient for the riverbed

angle real Acute angle of the major axis of the slick with the x-axis (radians)

anice ni real Manning's roughness coefficient for ice at the ice/water interface

amiuo 9 real Absolute viscosity of oil (gm-cm/s)

amuni t real Absolute viscosity of oil (Qbf-s/ft2)

aniu nw real Kinematic viscosity of water (ft2/s)

api API real American Petroleum Institute index for oils

apitem v-eal Temporary storage for API values

beta real Fractional empirical factor for bottom deposition of oil

cevp C real Evaporation coefficient C

cordlb(m) complex Complex variable giving x and y coordinate locating cross section on reference shore

cordv(m,n) complex Coordinates at which vstrm(mn) is acting

ddd DT real Horizontal diffusion coefficient

depth(n) h real Depth of the river at the nth grid box center (ft)

dstrm(m,n) real Flow depth at mth cross section and nth stream tube (ft)

difx DX real Diffusion coefficient in the x direction

dify D real Diffusion coefficient in the y direction

dvdt dv/dt real Rate of change of oil volume

dx real Width of the grid in x or y direction (since square grids are used,

the widths are equal) (ft)

67
fevpl real Fraction of oil evaporated at the end of previous time step

fevp2 real Fraction of oil evaporated at the end of present time step

fueltp character Type of oil

full character Character variable

gama Y real Empirical fraction for oil exchange from surface to suspended

layer

gravac g real Gravitational acceleration (ft/s 2)

hlife(17) X real Half-life of the shoreline for oil retention (hr)

ibend(n) integer Index for the cross section; if 0, flow direction will be along the

connection line of corresponding stream tubes in two neighboring cross sections; if 1, flow direction (shore or island)

idam(n) integer Particle numbers that have hit locks and dams; stored in the array where n is from 1 to ndam

idum(n) integer Temporary storage

ievery integer Frequency of obtaining output; e.g. value of 2 gives output

every other time step

igrilb(i) integer Water-side j grid number (y direction) of the lower boundary of

the river for the tih grid in the x direction

igriub(i) integer Water-side j grid number (y direction) of the upper boundary of

the river for the ith grid in the x direction

igrlb(k,i) integer Land-side j grid number (y direction) of lower boundary of kth

island (if any) for the zjh grid in the x direction

igrub(k,i) integer Land-side j grid number (y direction) of the upper boundary of

kdh island (if any) for the jth grid in the x direction

ihitb(n) integer Particle numbers that have hit land; stored in the array where n

is from 1 to nhitb

imovin(n) integer Particle numbers that are in the water; stored in the array where

n is from 1 to nmovin

indpm integer Index for printout; if 0, no printout; if 1, printout

indxld integer Index for spreading; if 0, axisymmetric spreading; if 1, one-dimensional spreading

ipos integer Grid box number

ipoll(n) integer Beginning grid box number of the nth ice region

ipos2(n) integer Ending grid box number of the nth ice region

izrbx(n) integer x grid number of the nth grid box to have zero velocity

izrby(n) integer y grid number of the nth grid box to have zero velocity

jdum(11) integer Temporary storage

kintm integer Number of interpolations between two velocity points of two

consecutive cross sections in the same stream tube

lcstsq(m) integer Last cross-section number of branch m

Ulmax integer Maximum 1th grid of the present moving grid

11min integer Minimum t1h grid of the present moving grid

Imax integer Maximum -fli grid of the present oil slick

1min integer Minimum 11h grid of the present oil slick

nbmch integer Number of branches used in oil spill model

ndam integer Total number of particles that have hit locks and dams

nfizrc(n) integer Next cross section connecting to cross section in question. For a

divided channel 'around island, this represents the first cross section cmuected to in the lower division from the main channel

CroL sectimn

NO& itg Total number of grids in the x direction

inteer Total number of particles that have hit land

nker ite Total number of ice regions

68

54,9 : '• -'
nicexl(n) integer x grid box of the starting point of ice region n

nicex2(n) integer x grid box of the ending point of ice region n

niceyl(n) integer y grid box of the starting point of ice region n

nicey2(n) integer y grid box of the ending point of ice region n

nmovin integer Total number of moving particles

nseco(n) integer For a divided channel around an island, this represents the first

cross section connected to in the upper division from the main

channel cross section (if no island, = 0; if lower division complete

and returning back to upper division, = 888; if both divisions are

complete and resuming main channel, = 999)

nsisct(n) integer Number of sounding depths used to describe the channel geometry at the nth cross section

nstube integer Total number of stream tubes at the nth cross section

ntpar integer Total number of particles

numcon(n) integer Condition number of cross section n. If all stream tubes continue

to next cross section undivided, = 11; if stream tubes divide into

two channels from main channel, = 12; if stream tubes from divided channels connect back to main channel, = 21

nzrvb integer Total number of boxes to assign zero velocities

parsur(n) complex x and y coordinate location of the nth particle in the surface layer

(ft, ft)

parsus(n) complex x and y coordinate location of the nth particle in the suspended

layer (ft, ft)

preldt real Time step for releasing new particles in case of continuous spill

(in one spildt several preldt occur) (s)

qsct(n) real Discharge at the nth cross section (ft3/s)

radius(n) real Distance of the nth particle from the spill centroid in case of axisymmetric spreadingor the distance of the nth particle from the

strip centroid in case of one-dimensional spreading

sctang(n) real Angle the nth cross section makes with positive x axis (radians)

sigma a real Surface tension (Ibf/ft)

slickr(n) real Slick radius of the nth pie segment (ft)

solbit real Solubility of oil (gm/m 3)

soluni real Solubility of oil (lb/ft3 )

spaice real Surface area of spill under ice (ft2)

sparea real Open surface area of spill (ft2)

spcen complex x and y coordinates of the spill centroid (ft, ft)

spcenO complex x and y coordinates of the spill centroid at the beginning of simulation (ft, ft)

spgoil real Specific gravity of oil

spildt At real Oil spill time step (s)

spIrat real Rate of oil spill in case of continuous spill (ft3/s)

spltim real Duration of oil spill in case of continuous spill (s)

spvol real Total volume of oil spilled (U.S. gal.)

text(n) character Storage array for character string

tice(mn) real Thickness of ice at the mth cross section an n

t

h sounding depth

location (ft)

timet real Total time from the beginning of simulation to the end of present

time step (s)

totme real Total time for which the simulation has to be performed (hr)

to(evp real Total volume of oil evaporated until the end of present time step

(ft3)

"tb (4i) real Oil rejection rate from the shore. The dimension 4 is for the four

69
shores (lower: 1, upper: 2, lower island: 3, upper island: 4). The

dimension i represents the ith grid in the x direction

tOevp real Initial boiling point of oil (K)

tOuni real Initial boiling point of oil ('R)

vbuoy Vb real Buoyant velocity of the oil globules in the suspended layer (ft/s)

vcar(n) complex x and y components of velocity in the nth grid box (ft/s, ft/s)

vmuni real Molar volume of oil (m3 /mol)

vmol real Molar volume of oil (m3 /mol)

volevp real Volume of oil evaporated during the present time step (ft 3)

volpie(n) real Volume of oil in the nth pie multiplied by 8 (ft3 )

volspi real Total volume of oil spilled (ft3)

vopsur(n) real Oil particle volume of nth surface particle

vopsus(n) real Oil particle volume of nm suspended particle

volspl real Imaginary original volume of oil spilled (m3

)

vpsur( real Initial particle volume in surface layer (ft3)

vtsurO real Initial total volume in surface layer (ft)

vtsusO real Initial total volume in suspended layer (ftW)

vstrm(m,n) complex x and y components of the stream velocity at the rnth cross section

and nth stream ube (ft/s, ft/s)

vwind complex x and y components of wind velocity (ft/s, ft/s)

vwmag real Magnitude of the wind velocity (ft/s)

vzero real Volume of the spill (ftW)

wsct(n) real Water level at the nth cross section (ft)

word character Character variable

yshift(n) real Distance the nth particle is shifted in the negative y direction if

there is an island (ft)

ywid(m,n) real Distance along the mth cross section from the reference bank to

the nt sounding depth location (ft)

z(m,n) real nth sounding depth for the mth cross section (ft)

zd(n) real Reference datum for nth section from which the sounding depth

is evaluated(ft)

70
APPENDIX B: LIST OF PROGRAMS

The computer program is written in FORTRAN 77 and has been tested with several Fortran

compilers. The program usage is discussed in the main text and Appendix A. Model output (results

files) are described in the text. The main program, ROSS2, with all of the following subroutines will

be saved on a floppy disk with input and output sample data for the users' reference.

ross2

advsur advsus areacl bndloc bndrej calcij

check checki chif ckout con conchk

disolu emdam emulfn evapor exchan input

lumper ohndc orient pievol ploter positn

presur presus random rearrg redepn rejloc

reloca shift sprdld sprdax strvol ohvad

veldep

71
RO Form Approved REPORT DOCUMENTATION PAGE OMB No. 07010788

Pubic reporting Ibrden for this collection of information is estimated to average 1 hour per response, includrng the time for reviewing instructions, searching existing data sources, gathering and

maintaining the data needed, and compateing and reviewing the collection of information. Send comments regarding this burden e-timate or any other aspect of th:s collection of informalon. including suggestion for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports. 12t5 Jefferson Davws H-gnway, Suae 1204 Arritglon.

VA 22202-4302, and to the Office of Management and Budget. Papeoerok Reduction Project (0704-0188). Washington, DC 20503.

1. AGENCY USE ONLY (Leave biank) 12. REPORT DATE 3.REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

A Mathematical Model for Oil Slick Transport and Mixing in Rivers

Contract No. DACA89-SS-K0019

6. AUTHO-S

Contract No. DTSL55-89-C-C0549

Hung Tao Shen, Poojitha D. Yapa, De Sheng Wang and Xiao Qing Yang

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

REPORT NUMBER

Department of Civil and Water Resources Engineering Department

Environmental Engineering LuleA University of Technology Special Report 93-21

Clarkson University LuleA, Sweden

Potsdam, New York

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

U.S. Army Cold Regions Research St. Lawrence Seaway

and Engineering Laboratory Development Corporation

72 Lyme Road Department of Transportation

Hanover, New Hampshire 03755-1290 New York 13772-0520

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

Approved for public release; distribution is unlimited.

Available from NTIS, Springfield, Virginia 22161

13. ABSTRACT (Maximum 200 words)

The growing concern over the impacts of oil spills on aquatic environments has led to the development of

many computer models for simulating the transport and spreading of oil slicks in surface waters. Almost all

of these models were developed for coastal environments. A few river models exist. These models only considered the movement of surface oil slicks. In this study a two-layer model, ROSS2, is developed for simulating oil spills in rivers. This model considers the oil in the river to consist of a surface slick and suspended oil

droplets entrained over the depth of the flow. The oil transformation processes considered in the model indude advection, mechanical spreading, turbulent diffusion and mixing, evaporation, dissolution, emulsification, shoreline deposition and sinking. The model can be used for simulating instantaneous or continuous

spills either on or under the water surface in rivers with or without an ice cover. The model has been implemented for the Ohio-Monongahela-Allegheny river system and the upper St. Lawrence River. This report

describes the model formulation and implementation. A case study is presented along with detailed explanations of the program structure and its input and output. Although it is developed for simulating oil spills,

the model can be applied to spills of other hazardous materials.

14. SI•BECT TERMS 15. NUMBER OF PAGES

com~ter ()l splls78 odes Oil Rivers 16. PRICE CODE

17. SECURITY CLAFSIFICAT1ON 118. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT

OF REPORT OF THIS PAGE OF ABSTRACT

UNCLASSI D UNCLASSIFIED UNCLASSIFIED UL

NON7504140S0M0 Standard Form 298 (Rev. 2-89)

Piewcitbe byANt SIW Z39-18

29flI02