A Lagrangian Method 

For Predicting Pollutant 

Dispersion in Bolinas Lagoon, 

Marin County, California 

GEOLOGICAL SURVEY PROFESSIONAL PAPER 582-B 

A Lagrangian Method 

For Predicting Pollutant 

Dispersion in Bolinas Lagoon, 

Marin County, California 

By HUGO B. FISCHER 

DISPERSION IN SURFACE WATER 

GEOLOGICAL SURVEY PROFESSIONAL PAPER 582-B 

A study of a tidal lagoon that relates observed 

dye-tracer dispersion to a numerical method for 

predicting pollutant dispersion 

UNITED STATES. GOVERNMENT PRINTING OFFICE, WASHINGTON : 1972 
UNITED STATES DEPARTMENT OF THE INTERIOR 

ROGERS C. B. MORTON, Secretary 

GEOLOGICAL SURVEY 

V. E. McKelvey, Director 

Library of Congress catalog-card No. 76-186372 

For sale by the Superintendent of Documents, U.S. Government Printing Office 

Washington, D.C. 20402 - Price 65 cents (paper cover) 

Stock Number 2401-1213 
CONTENTS 

Abstract -~------------------------------------------

Introduction ----------------------------------------

Acknowledgments -------------------------------

Description of the numerical method -----------------

General description of the computer program ____ _ 

Convective step ---------------------------------

Diffusive step ----------------------------------

Concentration-decay step ------------------------

Experiments in Bolinas Lagoon ----------------------

Description of experiments -----------·-----------

Experiment of June 18 and 19, 1968 _________ _ 

Page 

B1 

1 

2 

2 

2 

3 

5 

7 

7 

7 

7 

Experiments in Bolinas Lagoon-Continued 

Description of experiments-Continued 

Experiment of May 21 and 22, 1969 

Discussion of results ------------------------·--------

Experiment of June 18 and 19, 1968 ---·-----------

Experiment of May 21 and 22, 1969 --------------

Summary of experimental results ----------------

Numerical study of continuous discharge at the 

lagoon mouth -~----------------------------------- Conclusions _____________________________________ ..; __ _ 

Eteferences cited -------------------------------------

ILLUSTRATIONS 

Page 

B8 

8 

8 

10 

14 

14 

15 

18 

Page 

PLATE 1. Aerial photograph of Bolinas Lagoon at low tide on September 6, 1967 _________________________________ facing B4 

FIGURE 1. Diagram showing connection of one-dimensional segments for representation of Bolinas Lagoon by the numerical. program --------------------------------------------------------------------------------------- 3 

2. Diagram showing longitudinal dispersion between adjacent elements ---------------------------------------- 5 

3-9. Graphs showing: 

3. Forecast tides at mouth of Bolinas Lagoon, June 18-19, 1968 --------------------------------------- 8 

4. Observed and predicted dye concentrations at sampling sites, June 18-19, 1968 ------------------------ 9 

5. Forecast tides at mouth of Bolinas Lagoon, May 21-22, 1969 ------------------------------'-·-------- 10 

6. Flow velocity at the seaward end of segment 2, May 21-22, 1969 ------------------------------------ 11 

7. Predicted and observed dye concentrations at the seaward end of segment 2, May 21-22, 1969 --------- 11 

8. Predicted and observed longitudinal distributions of dye concentration along the axes of segments 2, 4, 

and 5 at approximately 1640 hours on May 21, 1~ho hours on May 21, and 0640 hours on May 22, 

1969 ------------------------------------------------------------------------------------------ 12 

9. Observed and predicted dye concentrations at two sampling sites, June 18-19, 1968, showing effects of 

changing the mass-transport coefficient in the numerical program ------------------------------ 13 

10. Maps showing numerically predicted concentrations resulting from a constant-rate discharge at the lagoon 

mouth, after 10 repetitions of the tidal cycle of May 21-22, 1969 ---------------------·----------------- 16. 

11. Maps showing numerically predicted concentrations resulting from an interrupted constant-rate discharge at the 

lagoon mouth, after 10 repetitions of the tidal cycle of May 21-22, 1969 ----------------------·----------- 17 

TABLES 

Page 

TABLE 1. Annotated ftow chart of the computer program _____ -------------------------------------------------------- B20 

2. Fortran listing of computer program for CDC 6400 computer ----------------------------------------------- 25 

ITI 

DISPERSION IN SURFACE WATER 

A LAGRANGIAN METHOD FOR PREDICTING 

POLLUTANT DISPERSION IN BOLIN~S LAGOON, 

MARIN COUNTY, CALIFORNIA 

By HUGO B. FISCHER 

ABSTRACT 

A numerical method is described which is capable of predicting the movement and dispersion of a pollutant in a tidal 

embayment. To use the method requires a knowledge of the 

embayment geometry and of a typical tidal cycle of watersurface elevations at various interior points. The embayment 

may be segmented into a two-dimensional pattern, although 

1low within each segment is considered to be one dimensionaL 

Each time increment includes a convective step, a diffusive 

step, and a concentration-decay step. Unlike some previous 

numerical methods, the convective step is contrived to minimize numerical dispersion. The computer program will generate the motion and dispersion of a concentrated slug of pollutant, or it will predict the concentration of pollutant at all 

points within the embayment resulting from a constant discharge at one point. 

The method is verified by predicting the dispersion of a 

slug of pollutant discharged near the mouth of Bolinas Lagoon, 

Calif. Bolinas Lagoon is typical of many tidal embayments in 

that, while it includes a number of branching eroded channels, 

most of its interior consists of mudflats which are exposed at 

low tide and shallowly covered at high tide. A field experiment is described in which a slug of Rhodamine WT dye 

tracer disperses in reasonable conformity to the numerical 

prediction. 

The numerical method is used to predict the concentration 

distribution in Bolinas Lagoon resulting from a continuous 

addition of a pollutant at the location of the existing sewage 

outfall. A study is also made of the benefit resulting from 

detention of the pollutant during a 2-hour period at low tide. 

In this case a significant concentration reduction is observed 

in those parts of the lagoon which were previously the most 

seriouely polluted. 

INTRODUCTION 

Bolinas Lagoon, Calif., is a tidal embayment approximately 3 miles long with a single opening to the 

Pacific Oeean. At low tide, water lies only in welldefined channels and pools, and most of the interior is 

exposed. At high tide, the lagoon surface presents an 

expanse of water broken only by one large and one 

small island. The tidal prism, which is approximately 

80 percent of the high-tide volume, fills and empties 

the lagoon each tidal cycle by flowing up and down 

the eroded channels. The shallow covering of mudflats at high tide and their exposure at low tide produce an ideal habitat for many kinds of bird and sea 

life. At one time Bolinas oysters were renowned in 

San Francisco restaurants, and clamming is still popular even though prohibited by health regulations during part of the year. Part of the lagoon has been purchased by the Audubon Society as a bird refuge. Of 

equal interest, though less well known to tourists, are 

a colony of seals in possession of a part of the large 

central island and a shark spawning ground in the 

northern basin. 

Bolinas Lagoon, like almost all areas of open water 

near large urban centers, has attracted the recreationalist and summer-home owner, and with him have 

come the problems of pollution. Incoming tides carry 

the pollutants which enter the lagoon at its mouth up 

the lagoon channels and distribute them throughout 

the interior. The purpose of this report is to provide 

a method for estimating pollutant concentrations within the lagoon and for determining how much pollution 

would be reduced by proposed changes. 

Previously published analytical methods for predicting pollutant dispersion were not devised to cope 

with the complicated geometry of an embayment such 

as Bolinas Lagoon. The studies of Taylor (1954) for 

flow in pipes, Elder (1959) for flow down an infinitely 

wide open channel, Fischer ( 1967, 1968) for uniform 

steady flow in natural streams, and Holley, Harleman, 

and Fischer (1970) for oscillating flow all require that 

the flow cross section be uniform along its length and 

Bl 
B2 DISPERSION IN SURFACE WATER 

constant in time. In Bolinas Lagoon the flow divides 

and redivides into ever-smaller channels as it penetrates the landward reaches, and the cross sections of 

the channels change from well defined, narrow, and 

shallow at low tide to wide and irregular at high tide. 

An analysis is required that is capable of including 

these wide ranges of geometry. 

Two numerical methods have been published to 

describe pollutant dispersion in a class of two-dimensional estuaries. Thomann and Sobel (1964) divided 

an estuary into subareas in plan, permitting transfer 

of mass between subareas according to empirical masstransfer coefficients. Orlob, Sellek, Shubinski, Walsh, 

and Stann (1968) used a similar method in which convection is included directly and is obtained from a 

numerical model of the hydraulics of the estuary. In 

both methods complete mixing is assumed within subareas during each time step. This mixing induces a 

numerical dispersion dependent on the size of subarea 

and length of time step, which is additive to the empirical diffusion coefficients between subareas. How 

much of the dispersion between subareas is included 

in the numerical dispersion of the computer program 

and how much is reflected in the mass-transfer coefficients is hard to determine. 

The first section of this report describes a numerical 

method that substantially reduces numerical dispersion. The program convects the pollutant at the mean 

flow velocity through a series of interconnected· channels and accepts changes in channel geometry caused 

by changes in tidal level. The program also includes, 

in a diffusive step, dispersion caused primarily by the 

shear effect. The shear effect is the variation in velocity between different points in the flow cross section 

and is described in detail in the analytical studies 

already referenced. To use the program the bottom 

geometry of the embayment must be known, as well 

as tidal elevations throughout a typical tidal cycle at 

a number of interior points. The program uses this 

information to predict the flow velocities at all points 

and the distribution of concentration resulting from 

either a slug or continuous injection of pollutant. 

Accuracy of the numerical program is described in 

the section "Experiments in Bolinas Lagoon." Computer-generated dye-concentration distributions were 

compared with the results of two dye-dispersion experiments in the lagoon. The observed and computerpredicted results are reasonably similar. The computer 

program was also used to predict the concentration 

distribution caused by a continuous input of mass at 

the lagoon mouth and the concentration reductions 

achievable by a short detention period during a part 

. of the tidal cycle. These results are given in the section 

"Numerical Study of Continuous Discharge at the 

Lagoon Mouth." 

Although the program was written specifically to 

describe dispersion in Bolinas Lagoon, the input data 

format is such that it could be applied without. change 

to many other embayments of similar character. It is 

hoped that this program will become a useful tool in 

pollution studies of many coastal environments. 

ACKNOWLEDGMENTS 

Many persons contributed to this report. The field 

studies were accomplished with the assistance of personnel and equipment from the Menlo Park office of 

the U.S. Geological Survey, with the particular aid of 

Mr. Loren Young, Mr. ,Joe Robles, and Mr. John 

Limerinos. In .addition, the writer wishes to thank 

Mr. Peter Ward, Mr. M. M. Das, Mr. Larry Eronini, 

and 1\{r. Harry Nichandros, graduate students at the 

University of California, all of whom voluntarily 

participated in data collection. 

Funds for computer use to test and operate the 

numerical program were provided by the Water Resources Center of the University of California under 

project WRC-165, "The Management of Surface 

Water Hydrologic Systems for Water Quality Control." Some of the boats used to collect data were 

provided without cost by members of the Bolinas 

Harbor Commission. 

DESCRWTION OF DIE NUMERICAL MEDIOD 

This paper describes a numerical program which 

simulates convection and diffusion in a tidal embayment of complicated geometry. The program is 

Lagrangian in concept, in that attention is fixed on 

the motion of identified finite elements of fluid and 

diffusion is permitted between the elements. N umerical dispersion, which occurs in the convective step if 

concentrations must be assigned to fixed points on a 

spatial grid, is minimized in this program because no 

spatial grid is established. Instead, one computes the 

exact location of each finite element. 

The program will be described in three ways: (1) 

the text describes what the program does, why it does 

it, what its limitations are, and. how it might be 

extended to other embayments, (2) table 1 gives a 

step-by-step summary of the program and its subroutines, and ( 3) table 2 gives a Fortran listing of 

the program itself. 

GENERAL DESCRIPTION OF THE COMPUTER PROGRAM 

·The numerical program is intended to describe the 

convection, dispersion, and decay of a pollutant which 

has been released, either continuously from a fixed 
METHOD FOR PltEDICTING POLLUTAN'l' DISPERSION IN BOLINAS LAGOON, CALIFORNIA B3 

source or as an instantaneous slug, into a tidal embayment. Flow in the embayment has the following characteristics: (1) It is not one dimensional in the sense 

used by Harleman (1954), Pritchard (1958), Bowden 

(1963), and others, because the flow branches into a 

number of intertwining channels. Within each of the 

many channels, however, the program treats the flow 

as being one dimensional along that channel's axis. 

(2) The tidal prism forms a substantial part of the 

total volume of the embayment at high tide, so that 

the flow geometry undergoes substantial change between low and high tide. To use the program, it is 

necessary to know the bottom geometry of the embayment and the tidal stage throughout the tidal cycle in 

each channel. 

Using maps, aerial photographs, and as much personal familiarity as possible, the program user begins 

by dividing the embayment into separate channels and 

each channel into segments. The channel segment, a 

basic division required for the program, is defined as 

a reach of channel throughout which the cross-sectional 

geometry and depth of flow are at all times approximately uniform, that is a reach of prismatic channel 

flowing at uniform depth. Flow occurs forward or 

backward in the direction of the segment axis and into 

or out of the segment ends, and the total volume of 

water contained in the segment increases or decreases 

as the tide floods or ebbs. The length of the segment is 

fixed, however, and neither flow nor diffusive mass 

transport is permitted across the segment sides. Each 

segment is connected to one segment in the sea ward 

direction (towards the mouth of the embayment) and 

to as many segments as desired in the landward direction. As an example, the division of Bolinas Lagoon 

into segments is shown on plate 1 and in figure 1. 

The water contained within each segment is divided 

along the segment axis into volume elements. A volume element is a coherent volume of water having a 

uniform concentration of pollutant and filling the 

segment cross section. Its length, along the segment 

axis, equals its volume divided by the area of the cross 

section. How these elements are moved along the segment axis and between segments will be described in 

the section "Convective Step." 

Each time step of the program consists of three 

parts: a convective step, in which volume elements are 

moved within the embayment at the mean flow velocity; a diffusive step, in which mass transfer occurs 

between adjacent elements; and a decay step, in which 

a first-order decay is simulated. Before the initial step, 

the concentrations of all the initial elements are read 

as input data. To simulate dispersion of a slug of pollutant, a high concentration is read into one element 

and zero concentration into all others. Between time 

steps, an additional mass of pollutant can be added to 

the element at any desired point within the embayment to simulate a constant rate of discharge. 

CONVECTIVE STEP 

The purpose of the convective step is to move each 

finite volume of water along its segment a distance 

given by the mean flow velocity in the segment multi6000 

0 

2200 

~Dye-injection site 

Bolinas sewage outfall- 0 

~@ 

l!'IGURE !.-Connection of one-dimensional segments for representation of Bolinas Lagoon in the numerical program. Numbers 

not enclosed in circles are segment lengths, in feet. 
B4 DISPERSION IN SURFACE WATER 

plied by the duration of the time step. This can be 1 

accomplished by knowing the entering and exiting discharges and the cross-sectional area for each segment. 

The scheme isinitiated at low tide by dividing each 

segment into a desired number of volume elements 

spaced evenly along the segment axis. The discharge 

into or out of the landward end of each segment is 

computed by the change of volume of all those segments connected landward of the segment under consideration, less the fresh-water discharge, if any. Discharge is taken as being positive in the landward direction and negative in the seaward direction. 

Each segment is considered in turn, beginning with 

segment 1, the most upstream segment, and proceeding 

to the highest numbered segment at the lagoon mouth. 

The discharge through the sea ward end is found by 

the formula 

(1) 

in which qi is the volume of water transported through 

the seaward end during the time step, q0 is the volume 

transported through the landward end, Yt+1 is the volume in the segment at the end of the time step, and 

Vt is the volume in the segment at the beginning of 

the time step. Note that q;, and q0 are volumes, not discharges, so that the duration of the time step is not 

included in equation 1. The segment volumes at the 

beginning and end of each step, Vt and Vt+1 , must be 

given as input data; these data can either be measured 

in the field or be obtained from a hydrodynamic program such as the method-of-characteristics program 

given by Lai (1967). For segments which are landward ends of the embayment (segments 1, 3, 5, 7, 8 

and 10 in fig. 1), q0 is minus the fresh-water inflow, if 

any. Fresh-water inflow is not permitted into segments 

which are connected to further landward segments; 

for these, q o of the segment being considered is equal 

to the sum of qi for the landward-connected segments. 

For instance, for the system shown in figure 2, 

qo =qi +q;, +qi , and 9 6 7 8 

qi = qo + Vt+l- Vt ; (2) 

9 9 9 9 

however, qo -fresh-water inflow into landward end of 1 

segment 1. Note that the segment numbering system 

must be arranged so that the qi for the landwardconnected segments have been computed in advance of 

computing q o for the segment under consideration. 

This can be accomplished by ordering the segments, 

with lowest order (segment 1) corresponding to the 

most landward segment and segment numbers always 

increasing sea ward. If ohe wished to consider a 

braided estuary in which one landward segment was 

connected to two sea ward segments, the program 

would have to be modified. 

The convection which takes place during each time 

step is accomplished by moving volume elements between segments. At the beginning of the time step 

each segment contains NMAX volume elements identified by the subscripts m and n, where m is the number 

of the segment and n is the number of the volume element in the segment ( NMAX is different for each seg-

,ment), V m,1 is the volume of the element at the landward end of the segment, and Vm,n is the volume element at the seaward end of the segment. Consider first 

what occurs at the seaward end of a given segment. 

Two cases must be distinguished; either the discharge 

is into the seaward end (positive) or out of the seaward 

end (negative). If the discharge is into the seaward 

end, a new volume element labeled Vm,n+1 is added 

to the segment of volume qi. The number of volume elements which the segment contains is increased by one. 

If the discharge is out of the seaward end of the segcent, that is qi is negative, the seaward volume element, 

labeled Vm,n, is reduced by the amount of ihe outflow 

volume. If the endmost element, Vm,m is not large 

enough to supply the entire outflow, it is eliminated 

entirely, and the further outflow is taken from Vm,n-1 ; 

in that case, the total number of volume elements within the segment is reduced by one. If insufficient volume 

is available in the combination of Vm,n and Vm n-1, volume is taken from Vm,n- 2, and this procedure continues 

until sufficient volume is provided. The only restriction is that during a single time step no more volume 

can flow out of either end of a segment than exists in 

the segment at the start of the step. If this does happen, the program stops, and a shorter time step must 

be used. 

A similar procedure is used for flow through the 

landward end. If flow is into the landward end, that 

is qo is negative, a new element is added whose volume, 

vm,1 , is equal to the flow into the segment. The identification of each of the remaining elements in the segment is increased by one, that is Vm,1 becomes vm,2, and 

so on, and the total number of volume elements in the 

segment is increased by one. If flow is out of the landward end, that is q0 is positive, the volume which flows 

out is taken first from vm,1, and if that is not sufficient, from vm,2, and so on. If, for instance, the outflow requires all of Vm,1, and part of Vm,2, the remaining portion of Vm, 2 becomes 'l'm,1, Vm,a, becomes vm,2, 

and so on, and the total number of elements is reduced 

by 1. 

In addition to storing the volumes v1, 1 through 

Vm,m the computer stores the concentration of mass 

contained In each volume, c1,1 through Cm,n, and the 
GEOLOGICAL SURVEY 

@ 

Segment boundary and number 

used in numerical program 

4 

'Y 

Sampling location and 

number 

0 

EXPLANATION 

Line of main channel 

1000 2000 3000 FEET 

BOLINAS LAGOON 

PROFESSIONAL PAPER 582-B PLATE I 

XA 

Position of moo red boat 

May 21- 22. 1969 

Aerial view at low tide on September 6, 1967. The red spots are dye injected by helicopter (Ritter, 1969). 
METHOD FOR PREDICTING POLLUTANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B5 

number of elements in each segment. When volume elements are moved between segments, the concentrations 

attached to the elements are moved also. For instance, 

if an outflow occurs from segment 1 to segment 2, the 

concentration of the outflow becomes the concentration of the newly established element in segment 2, 

that is c2, 1• If two segments converge to provide flow 

into a third segment, for instance segments 2 and 3 

flow into 4 in figure 1, the concentrations of the flow 

coming from segments 2 and 3 would be combined in 

forming the new element in segment 4, such that the 

total mass carried by the flow was conserved. 

The concentration of the outflow of a segment must 

be calculated from the concentrations of the elements 

participating in the outflow. If all the outflow can be 

taken from the endmost element in the segment, the 

concentration of the outflow equals that of the endmost 

element. If more volume is carried out of the segment 

than is contained in the endmost element, the outflow 

concentration must be calculated by the conservationof-mass equation, 

CtVl + •• • + CqVq, 

qo 

(3) 

in which v1 through Vq_1 are the volun1es of elements 

which are carried out in their entirety, beginning with 

the outermost element and working inward, <h through 

Cq-1 are the concentrations of those elements respectively, Cq is the concentration of the next inwardmost 

volume element, and 

Vq=qo- (vt + · · · +vq-1)· 

Note that equation 3 does present one source for numerical mixing because concentrations of several volume elements are combined into one. 

This completes the description of the convective 

step. Neither the physical location of any of the volume elements nor the velocity at which they have 

moved has been computed. Nevertheless, the procedure 

just described has accomplished what was proposed at 

the outset of this section, namely the convection of 

each volume element at the mean flow velocity of the 

segment in which it is located. If no diffusion between 

volume elements were considered, the concentrations 

established in the volume elements at the outset of the 

program would be convected up and down the embayment according to the tidal velocities, the only possibility for numerical dispersion being at the connecting 

points between segments. 

DIFFUSIVE STEP 

I~ongitudinal dispersion, which occurs between elements as they move along the channels, is caused primarily by the difference in streamwise velocity between different points on the cross section. While the 

elements are convected at the mean flow velocity, some 

individual water particles move faster than the mean 

and some slower. Particles in the deeper, faster moving 

sections are convected from one element to the one 

ahead, while particles in the slower moving zones near 

the channel banks may lag behind and be moved from 

the element they are in to the one behind. The result 

is dispersion by what is referred to as the shear effect. 

Analytic studies of the shear effect, which seek 

methods to predict the rate of longitudinal dispersion, 

were referenced in the introduction. These studies give 

methods for predicting dispersion coefficients, which 

relate the rate of mass transfer through a cross-sectional surface to the local concentration gradient. The 

dispersion coefficients given by these studies are for 

asymptotically long times; that is, the flow must remain steady for a sufficiently long time before the dispersion coefficient applies. For the typical cross section 

in Bolinas Lagoon, a sufficiently long time is on the 

order of many tidal cycles; hence, the flow does not remain steady, and the cross section does not remain constant long enough for the results of the analytic studies to apply. 

The simplest model of a diffusive process assumes 

that the transfer of mass between adjacent elements is 

proportional to the difference in concentration between 

them. The purpose of this section is to provide a way 

of estimating the constant of proportionality. Consider 

two adjacent equal-length elements as they move along 

a channel, and suppose that during a time step they 

happen to move a distance exactly equal to their own 

length. Figure 2A shows the boundaries of the element 

at the beginning of the time step, and figure 2B shows 

them at the end of the time step. The assumption embodied In the numerical program is that the bounA 

---......................... -------............... ... , ----................. , ----..... -.... ..... 

~, .... \ ' )0 }0 I 0 I I __ .,/ ... .... / __ ... / 

__ ........... --..... -- _ .... 

__ ............ --- --- B 

FIGURE 2.-Longitudinal dispersion between adjacent elements. A, Position of elements at beginning of time step. 

B, Position of element boundaries at end of time step 

(solid lines) and position of interface between :fluid particles in original elements (dashed lines) . 
B6 DISPERSION IN SURFACE WATER 

daries between elements are always normal to the 

channel axis. If, however, at the beginning of the time 

step one could dye all the water particles in element 1 

blue and all the particles in element 2 red, the boundaries between the blue and the red spots at the end of 

the time step would be approximately those shown by 

the dashed lines in figure 2B. The water particles 

which began the time step in element 2, but were in 

slow moving water near the banks, will lag behind the 

element and will be in element 1 ; similarly, some of 

the faster moving particles in element 1 will have 

moved ahead into element 2. One can estimate from 

the figure that approximately 20 percent of the water 

particles will have exchanged elements. 

At the conclusion of the time step the assumption 

required by the one dimensionality of the model is 

that each element becomes fully mixed; that is, the 

particles which enter element 2 from element 1 become evenly distributed over element 2. Although this 

assumption is inaccurate, it could only be eliminated 

by establishing a two-dimensional mesh for each segment. Such a program was written as part of this 

study, but its operation for the entire lagoon would 

have cost far more than the funds allocated to the 

project. Keeping within the one-dimensional approximation and considering transport between the two 

equal-length adjacent elements shown in figure 3, the 

adjustment in concentration in each element to account 

for mass transport between the two is given by, 

c1 = c1 + E ( c2 - c1 ) and ( 4) 

t+l f t t 

02 = c2 + E ( o1 - c2 ) , ( 5) t+l t t f 

in which c is the concentration in element 1 at 

the beginning of the time step, o is the concentration in element 1 at the beginning of the time step, 

and so forth. E is the mass-transport coefficient, a constant which can be supplied as input data to the 

computer program. Equations 4 and 5 describe the 

mass transport across one face of each element. For 

an element exchanging mass with elements on both 

sides, two faces must be considered, for instance 

o2 = c2 + E ( c1 - o2 ) + E ( o3 - c2 ) • ( 6) 

t+l t t t t 

Figure 2 suggests that E should have a value of approximately 0.2. 

Normally the lengths and volumes of adjacent elements will not be equal, and no element is likely to 

move a distance exactly equal to its own length. If 

the distance traveled by an element during the time 

step is less than its length, the mass transported across 

the boundary will be reduced because the distance 

traveled by the dashed boundary in figure 3 wil~ be 

less; similarly, if the distance traveled is more than 

the element's length, the mass transport will be increased. For want of a more exact description, this 

effect has been represented in the computer program 

by a linear relationship; the mass-transport coefficient 

is adjusted as follows: 

E'=E_I_, (7) We 

B 

in which E' is the adjusted mass-transfer coefficient, 

Tis the distance traveled by the shorter element during the time step, and We-- is the length of the shorter 

element of the pair beirig considered. 

A second adjustment to the transfer coefficient, perhaps more difficult to justify analytically but essential to properly represent the physical phenomenon, 

depends on the ratio of element length to segment 

width (or element width, which is the same thing). 

If the element is much longer than wide and if it 

moves a distance equal to its own length, cross-sectional mixing may be sufficient to reduce the rate of 

longitudinal dispersion, because a water particle origina1ly in a low velocity zone of flow may have time 

to migrate into a high velocity zone, or vice versa. A 

particle originating in the low velocity zone of element 2 may at first lag behind the mean flow and 

fall into element 1, but then migrate into the higher 

velocity zone and again pass into element 2. The narrower the element, relative to its length, the more of 

this mixing will occur because the high and low velocity zones are relatively closer together. On the other 

hand, if the element width is much greater than its 

length little mixing can occur across the width of the 

cross section~ so mass transport through the boundary 

is increased. Again, for lack of a better relationship, 

the effect is considered to be linear and the masstransport coefficient is adjusted by the ratio 

E" = E' __!!_____ We 

=E BT 

w 2 ' e 

in which B is the segment width. 

(8) 

Numerical stability requires that E" always be less 

than 0.5. If equation 8 gives a value of E" greater 

than 0.5, the value is arbitrarily reduced to 0.5. 

For each segment the program considers each adjacent pair of elements in turn, first selecting the 

shorter element for use in equation 8. In effect, the 

transfer is assumed to occur between the shorter element and its equivalent length within the longer element ; then the mass transported into the longer element is distributed uniformly therein. The other pos-
METHOD FOR PREDICTING POLLUTANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B7 

sibility, basing the transfer on the longer element, is · 

not workable because shrinking the mass transported 

into the shorter element might result in a higher concentration than previously existed in either element 

separately. The equations used in the diffusive step 

are, for the shorter element, 

C8 = C8 + E" ( Cz t+l t t 

and, for the longer element, 

Cz =cz +E"(cs - Cz ) 

t+l t t t 

Cs ) 

t 

8 

(9) 

(10) 

in which the subscript l refers to the longer and s to 

the shorter element. These equations are carried out 

for each pair of elements in each segment. No diffusion 

is permitted between end elements of adjacent segInents because the increased programming complication did not seem justified by the minor increase in 

accuracy. 

Operation of the program to simulate the experiments described in the next section showed that one 

further modification was necessary to correctly simulate the diffusive process. At times of low flow, particularly in the terminal segments, it often happens 

that a number of short elements are generated at the 

ends of a segment. Since the transfer coefficient between these elements must be restricted, to preserve 

numerical stability to a value of 0.5 the rate of longitudinal diffusion is incorrectly restrained. The probem is solved, albeit at the expense of generating some 

numerical mixing, by including in the program the 

possibility of combining adjacent short elements to 

form longer ones. This occurs if either of two triggering conditions are met: if the number of elements 

in the segment is greater than 10 and the average 

element length is less than one-third the segment 

width, or if the combined length of the two endmost 

elements at either end of the segment is less than 

the average element length within that segment. For 

either condition the shortest element in the segment 

is <'ombined with the shorter of the two adjacent elements to form a single new element whose volume 

equals that of the previous two and whose concentration i::; the weighted average of the previous two. The 

elemeHt combination, when required, is made after the 

convective step but before the diffusive step; after the 

combination, the program loops back and tests whether 

another element combination in the same segment is 

required. The diffusive step is not carried out until 

the triggering conditions for element combination are 

not pr•~ent. 

CONCENTRATION-DECAY STEP 

Each time step concludes by simulating a first-order 

decay of concentration. The concentration of each element is reduced by a fixed percentage determined by 

the rate constant of the process .. For instance, if one 

wishes to simulate the decay of biochemical oxygen 

dPmand (BOD) with a rate constant of 0.1 per day 

and if the time step is 15 minutes long, the concentration in each element at the end of each time step 

should be multiplied by the factor 0.99762. This factor 

is read into the program as input data and is used to 

reduce the concentration at the end of each time step. 

EXPERIMENTS IN BOLINAS LAGOON 

DESCRIPTION OF EXPERIMENTS 

The computer program can be used to study the 

dispersion of a slug of pollutant injected as an instantaneous source, or it can determine the distribution of concentration caused by a steady source. Practical concern is usually with a steady, or quasi-steady, 

source, such as a sewage outfall. Correct operation of 

the program, however, is more easily observed by study 

of the dispersion of a slug injection, because both the 

movement of the mean location of the slug and its dispersion around the mean are 1nore easily measured. 

:Moreover, if the program can be verified for the. dispersion of a slug, it can be used with confidence to 

study the distribution from a constant-rate injection, 

because a constant-rate injection is the sum of a series 

of slug injections. The two field experiments designed 

to test the computer program were both slug injections at the mouth of Bolinas Lagoon. Both injections 

were made at the point indicated on plate 1, the first 

at approximately 1400 hours on June 18, 1968, and 

the second at 1345 hours on ~fay 21, 1969. The major 

difference between the experiments was the type of 

tide, as described in the following sections. 

EXPERIMENT OF JUNE 18 AND 19, 1968 

Forecast tides at the mouth of Bolinas Lagoon for 

.June 18-19, 1968, are shown in figure 3. On June 18, 

low slack water occurred at the injection site between 

1300 and 1315 hours. The plan was to inject the tracer 

slug, 161,6 gallons of Rhodamine WT fluorescent dye, 

as rapidly as possible across the channel cross section 

at 1345 hours. The drum of dye was carried in a rowboat towed behind a larger boat intended for use later 

for sampling. The dye drum and two operators proved 

more than a safe load for the rowboat, and the hose 

prov.ided to drain the tank was too small to carry a 

sufficient flow for even an approximately instantaneous 

injection. Although the injection began on schedule 

at 1345 hours, only approximately half the tank was 
DISPERSION IN SURFACE WATER 

6~----------~-----------,,-----------. 

1-

LLJ 

LLJ 

LL. 

c:i 5 

LLI 

1-

<( 

~ 

~ 4 

0 

..J 

~ 

0 

..J 

z 3 

<( 

LLJ 

~ 

LLI 

E) 2 a:l 

<( 

z 

0 

~ 1 > LLJ 

..J 

LLJ 

L__ Durat~on of ___j I expenment ~ 

OL_----------~----------_..J----------~ 

2400 1200 

JUNE 18, 1968 

2400 1200 

JUNE 19, 1968 

FIGURE 3.-Forecast tides at mouth of Bolinas Lagoon, June 

18-19, 1968 

drained at 1410 hours. The tank was then light enough 

to be upended, and the remainder of the contents was 

emptied within several minutes. 

Before the dye injection, standard concentration 

samples were prepared at the experiment site by diluting known concentrations of dye with sea water 

taken from the mouth of the lagoon. The standard 

samples had concentrations of 0.5, 1.0, 1.25, 2.5, 5.0, 

10.0, 12.5, 25, and 50 JLg/1 (micrograms per liter). The 

bottles containing the standards were handled in 

exactly the same way as all sample bottles taken from 

the lagoon during the experiment, so that any physical 

or chemical reactions which might occur to the dye in 

the collected samples would also occur in the standards 

of known concentration. After the experiment, all 

bottles containing both samples and standards were 

taken to the hydraulics laboratory of the University 

of California, Berkeley, where they were analyzed on 

June 21, 1968, for fluorescence by using a Turner 

n1odel 111 fluorometer. The calibration of the fluorometer by the standard samples was found to be linear 

in all cases. 

During the experiment, samples were obtained by 

four parties, three in boats and one on shore. A total 

of 12 sites was to be sampled approximately every half 

hour during a 12-hour period. Each boat party attempted to sample two or three sites, while the shore 

party repeatedly circled the lagoon taking samples at 

five shore sites. The sampling sites are shown on plate 

1 ; sites 1, 2, 3, 4, and 5 were marked by buoys, lighted 

at night. Site 7 was off a pier at the town of Bolinas 

and was sampled from the pier during periods of high 

tide. Site 6 was intended to be over the mudflats in 

shallow water approximately even with site 5. At 

night, the sampling party had difficulty orienting ~tself to site 5 and determining where to sample. S1te 

3 was abandoned after dark when the sampling party 

was unable to navigate between sites 2 and 3 without 

running aground repeatedly on mudflats. Thus continuous and reliable concentration records were obtained 

only at sites 1, 2, 4, 5, and 7 and the shore-sampled 

sites 8 through 12. Whenever possible, the samples 

from shore were obtained by wading into water at 

least 2 feet deep. At site 10, however, the sampler 

often sank into the mud and was forced to obtain the 

sample just at water's edge. Sampling was continuous 

in mo'3t cases until 0200 hours on June 19. At points 

1 and 7, sampling -\vas continued until 0800 hours to 

determine how much of the dye carried out to sea 

during the ebb was returned on the following flood. 

EXPERIMENT OF MAY 21 AND 22, 1969 

The experiment of May 21-22, 1969, was conducted 

in a manner similar to that of the experiment of June 

18-19. 1968, except that only two sampling boats were 

avail~ble. An improved dye injector allowed the entire 

quantity of approximately 25 gallons of Rhodamine 

'VT dye to be injected in approximately 5 minutes, 

beginning at 1345 hours on May 21. One of the two 

sampling boats was 1noored permanently at the seaward end of segment 2, at the point marked "A" on 

· plate 1 ; at 15-minute intervals samples of water from 

the surface were obtained, and water velocity was 

measured at approximately the 0.4 depth (from the 

bottom) beneath the boat. The other sampling boat 

was used to cruise up and down the lagoon, taking 

samples along the segment axes. A traverse beginning 

at the landward end of segment 2 and proceeding 

down segment 2 and segment 4 and up segment 5 to 

its landward end could be completed with a collection 

of 15 samples in approximately 40 minutes. The axes 

of segments 6, 8, 9, 10, and 11 were also traversed at 

infrequent intervals. The moving-boat procedure was 

a change from that of the previous experiment made 

for two reasons: (1) only enough personnel were available t~ form two sampling parties and (2) by moving 

one sampling boat, that party could attempt to obtain 

an approximately instantaneous view of the spatial 

distribution of dye concentration, which would compare directly with the 'computer output. 

DISCUSSION OF RESULTS 

EXPERIMENT OF JUNE 18 AND 19, 1968 

The dve concentrations measured at all sampling 

sites during the 1968 experiment are shown in figure 4. 
METHOD FOR PREDICTING POLLUTANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B9 

,,,=:;:\ .J 12,-------------------------------------~ 

: ~~~7 L ............ _~~~. ~ 

~t~_s_IT_E_3~~~~~•-•.-•-·-•-•----------------~j 

:t SITE 4 """" .. L j 20,-------~~--------------------------~ 

16 

12 

8 

ffi 4 SITE 5 I j 0 L__--.....&-------==:::..::::.... ____________________ _J 

0::: 

LLJ 

a.. 

(/) 

:::;;; 

<{ 

0::: 

(.!) 

0 

0::: 

(.) 

i 

~ 

z. 

0 

f= 

<{ 

0::: 

1-

z 

LLJ 

(.) 

z 

0 

(.) 

~~~SITE 6 ..!\~, 

~ t SITE 7 

12 

~ f SITE~ .?::=--- L 

20 

16 f- 1, 

I \ 

12 f- }' Sf4 r- SITE 9 \,_~ 0 

-

-

-

-

16·------------------------------------~ 

I~ t SITEIO " £• l ~L~ _S_IT_E __ 11~~--.-~'--~~ ~•ko----------------_J 24.-----~~----------------------------~ 

20 

16 

12 

8 

4 

EXPLANATION 

Observed 

concentration 

o---<> 

Predicted 

o~--_.~~--~--~--~L---L---~--~--~ 

1400 1600 1800 2000 2200 2400 0200 0400 0600 0800 

JUNE. 18, 1968 JUNE 19, 1968 

FIGURE 4.-0bserved and predicted dye concentrations at 

sampling sites, June 18-19, 1968. 

The flow, as it enters the mouth of the lagoon, forks 

into two channels : the northern one known as the Bolinas channel and the southern one known as the Stinson channel. The dye slug was injected only into the 

Stinson channel on an incoming tide so that none of 

the slug entered the Bolinas channel. The flow further 

divides at the east side of the lagoon into one channel 

which .flows south towards Stinson Beach (at the south 

end of the lagoon) · and one which flows north to the 

north end of the lagoon. The dye was carried up the 

Stinson channel and divided into the two channels, 

through which high concentrations reached the southernmost and northernmost ends of the lagoon. During 

tidal outflow, the dye returned from the channels 

along the same routes by which it went in. The dye 

concentration measured at station 7 showed almost no 

dye concentration at any time during the inflow and 

outflow of the first tidal cycle; thus for these tidal 

conditions essentially none of the dye carried up the 

Stinson channel crossed over and flowed out by way 

of the Bolinas channel. 

The dashed lines in figure 4 show the concentration 

at the sampling points predicted by the computer program using as a mass-transport coefficient E = 0.2. To 

use the program, the lagoon was segmented as described previously ; widths and lengths of each segment were scaled from maps and aerial photographs, 

and tidal elevations were estimated from the forecast tide at the mouth. A time step of 15 minutes 

was adopted, and a 121h hour period of inflow and 

outflow was simulated. Since tidal elevations were not 

actually measured during the experiment, the data 

used as input to the computer program are necessarily 

inexact; nevertheless, the comparison between computer prediction and measured concentrations at most 

points is quite good. 

The· best verification of the accuracy of the diffusive 

step in the computer program is the comparison of 

predicted and observed concentration at site 1. By the 

time the dye cloud reached this point during tidal outflow, the slug had been dispersing for 12 hours. The 

slope of the rising limb of the slug is a measure of 

how much the slug has dispersed. Figure 4 shows that 

the measured and predicted slopes are almost exactly 

the same; thus dispersion has been accurately modeled. 

The convective motion of the dye cloud also seems 

to have been accurately modeled, as shown particularly 

at sites 5 and 12 by the proper occurrence of the first 

passage of the cloud and the low concentration observed at high tide. The cloud passes each site at ap-

. proximately the time and with approximately the extent predicted by the computer program. The only 

significant and repeated error is that during tidal out-
B10 DISPERSION IN SURFACE WATER 

flow the concentration at several sites (particularly 

5 and 10) decreases more quickly in the prediction 

than was actually observed. It seems likely that diffilsion in the northernmost basin, represented by segment 1 in the _computer program, is not adequately 

represented. 

No predictions are shown for sites 3, 6, 7, and 8. 

At site 3 the program predicted passage of a high 

peak of short duration prior to the time the first sample was obtained and predicted outflow of the cloud 

after sampling had been discontinued; hence no comparison of predicted and observed clouds could be 

made. At site 7, the program predicted zero concentration until low tide, when the numerical program 

was stopped; the peak observed at site 7 at 0400 hours 

on J nne 19 was caused by tidal inflow carrying in dye 

which had previously flowed out to the lagoon mouth. 

The erratic concentrations measured at site 6 probably 

result from the site's location near the boundary between segments 2 and 10, one of which is fed by the 

Bolinas channel and the other by the Stinson channel. 

The sampling site was not well marked, and the sampling party probably sampled sometimes in one segment and sometimes in the other. Site 8 "is also near 

the boundary between two segments, numbers 5 and 7; 

presumably the high concentrations observed in samples taken after 2100 hours result from water in segment 5, whose concentration is better indicated by the 

measurements at site 12. 

EXPERIMENT OF MAY 21, AND 22, 1969 

The major difference between the experiments of 

1968 and 1969 was the type of tidal cycle. The 1968 

experiment was planned to coincide with a 121f2-hour 

symmetrical tide; that is, the tide returned to its 

beginning level after 12112 hours. The 1969 experiment 

was conducted during a 25-hour nonsymmetrical cycle; as shown in figure 5, the tide first rose to a low 

high, then dropped slightly to a high low, then rose 

to a high high, and finally returned to approximately 

the beginning level. This type of tide tests more severely the approximations embodied in the computer 

program, because during the long period of high water 

significant exchange is possible across the boundaries 

of supposedly nonconnect.ed segments. Also, because of 

different velocities of propogation of the tidal wave up 

the Bolinas and Stinson channels, a net circulation of 

water flowing up channels 4 and 8 and down channel 

10 seems likely between the low high and the high 

low tides. This circulation is not permitted in the 

computer program and seems to have been a source 

of significant error. 

Prior to the actual experiment, the tidal hydro6~--------~-----------.----------,-----, 

1-

I..LI 

I.LI 

IJ... 5 

ci 

I.LI 

!;: 

3: 4 

3: 

0 

..J 

~ 3 ..J 

z 

~ 

:2 2 

I.LI 

> 0 

aJ 

<( 

z 

0 

i= 

<( 

> 

~ 0 I.LI 

2400 

I Duration of I 

f-•----experiment---+1• 

1200 

MAY 21, 1969 

2400 1200 

MAY 22, 1969 

FIGURE 5.-l!"'orecast tides at mouth of Bolinas Lagoon, May 

21-22, 1969 

dynamics of the lagoon were modeled using a segmented method -of -characteristics computer program 

similar to that described by Lai ( 1967). The hydrodynamic program, which had not been developed by 

the time of the 1968 experiment, did not permit flow 

between the Bolinas and Stinson channels, however, 

and so the results are somewhat in error. The hydrodynamic program predicted velocities and tidal elevations in each segment, and the resulting elevations 

were used as input to the dispersion program. The 

predicted mean velocity in the cross section at the 

seaward end of segment 2 is shown in figure 6, which 

also shows the velocities measured beneath the moored 

boat. One should not expect the magnitudes to be the 

same, as the observed velocity does not represent a 

cross-sectional mean. Nevertheless, the general appearance of the two plots and the times of tidal reversal 

are sufficiently similar to say that the hydrodynamic 

program appears to be reasonably accurate. 

The observed and predicted dye concentrations at 

the seaward end of segment 2 are shown in figure 7. 

As for the previous experiment, all predictions shown 

for this experiment use a mass-transport coefficient of 

E = 0.2. The peak concentration arrived sooner than 

predicted and was sharper. A partial explanation may 

be that the observed values are for one point in the 

center of the channel, whereas the prediction is of 

of the cross-sectional mean concentration. The sampling party reported that after passage of the peak 

at the boat, substantial concentrations of dye could be 
METHOD FOR PREDICTING POLLUTAN'l' DISPERSION IN BOLINAS LAGOON, CALIFORNIA . Bll 

0 

z 

0 

() +1 LlJ 

(/) 

a:: 

LlJ 

0.. 

fLlJ 

LlJ 

1..1.. 

~ 

~ 

.0 

.0 

LlJ 

u -1 0 

..J 

LlJ "C 

> 0 

0 

i:i: 

EXPLANATION 

o----o 

Predicted 

...,__..... 

Observed 

-2~------~----~~------~------~------~------~------~------~----~~------~------~----~ 

1400 1600 1800 2000 2200 2400 0200 0400 0600 0800 1000 1200 1400 

MAY 21, 1969 MAY 22, 1969 

FIGURE 6.-Flow velocity at the seaward end of segment 2 (point A in fig. 1), May 21-22, 1969. Predicted values are crosssectional mean velocities. Observed values are velocities measured at 0.4 depth directly beneath the moored boat. 

a:: 

LlJ 

fTidal~ inflow ~Tidal~ outflow 

r----Tidal 

inflow 

-------+------------------ Tidal ----------------- outflow 

:J 

a:: 

LlJ 

0.. 

(/) 

:::?; 

<( 

a:: 70 (.!) 

0 

a:: 60 S:2 

:::?; 

~ 50 

z 

0 40 i= 

<( 

a:: 

fz 30 LlJ 

() 

z 20 0 

() 

10 

0 

1400 1600 

"C "C Q) Q) 

> Q; ~ (/) "C 

.0 Q) 

~ 

0 tt 

+" «< Q) 

... :2 S? Q)+J 

1\ j ~ I \ ~ .c 

I \ ~ ~ I \ (f)-

I \ 

I \ I \ 

I \ 

I \ I \ 

1800 

\ 

~ 

' ' 

2000 

MAY 21, 1969 

"C "C 

2 Q) 

> -~ Q; "C (/) Q) .0 tt 0 

2200 2400 0200 0400 0600 

MAY 22, 1969 

EXPLANATION 

o----o 

Predicted 

..,._____. 

Observed 

0800 1000 1200 1400 

FIGURE 7.-Predicted and observed dye concentrations at the seaward end of segment 2 (point A in fig. 1), May 21-22, 1969. 

seen along the banks. If enough cross-sectional concentration measurements had been obtained to define 

the cross-sectional mean, the observed and predicted 

results would have been in better agreement. (Note, 

for instance, the comparison of the centerline and 

cross-sectional mean concentrations observed in the 

Green-Duwamish River by the writer (Fischer, 1968, 

p. A20, fig. 16).) 

According to the prediction, the main body of the 

dye cloud was expected to pass point "A" during the 

initial tidal inflow. During the partial outflow between 

the low high and the high low tides, some of the cloud 

was predicted to return to the point, so that a peak 

concentration would be recorded at the low slack 

water. During inflow to the high high tide the entire 

cloud was expected to be carried landward of the 

measuring point into segment 1. Then, during the 

main tidal outflow, the entire cloud was expected to 

pass the observation point, with the peak concentration occurring at approximately 0730 hours on May 
B12 DISPERSION IN SURFACE WATER 

22. The observed values show that the prediction was 

reasonably correct; the entire cloud did pass the observation point during the first inflow, and part of the 

cloud did return so that a peak was observed at the 

low slack water around 2200 hours on May 21. During 

the high high tide, concentrations were higher than 

had been anticipated, possibly because of mixing with 

the high concentrations expected in the landward end 

of segment 8. During the main tidal outflow, a peak 

was observed at approximately the anticipated time, 

but its magnitude was much less than expected. This 

result agrees with that of the 1968 study, which showed 

that the program underestimates the rate of mixing 

in the northernmost basin (segment 1) . 

Longitudinal dye distributions along the axes of segments 2, 4, and 5 are shown in figure 8. In agreement 

with figure 7, the observed distribution shows that 

for the first traverse, at approximately 1640 hours on 

a:: 

LLI 

1-

:::::i 

a:: 

LLI 

a.. 

(/) 

::2: 

<( 

a:: 

" 0 

a:: 

~ 

::2: 

~ 

z 

0 

f= 

<( 

a:: 

1-

z 

LLI 

(.) 

z 

0 

(.) 

160 

140 

120 

100 

80 

60 

40 

20 

0 

60 

50 

40 

30 

20 

10 

0 

40 

30 

20 

10 

0 

Landward 

direction 

o--o-.0<>--o.-o.. -o.. .... '0.--o._'Q 

Segment 2 

1000 

', 

' 

2500 

Landward 

direction 

EXPLANATION 

o-----o 

Predicted at 1640 hours, May 21 

• • 

Observed 

Start 1608, end 1643 

u----o......_..._..._ ,, "'0... 

,, ',, 

"'0...., .... , 

"'0...., 

........ 'o... .d __ ""'<Y,..cr 

EXPLANATION 

o-----o 

Predicted at 1940 hours, May 21 

• • 

Observed 

Start 1854, end 1925 

EXPLANATION 

o-----o 

Predicted at 0640 hours, May 22 

• • 

Observed 

Start 0615, end 0700 

'0-----o.. ___ -o-

---o.. ___ o()._ ,' 

---o-. I --"'"'().---.& 

Segment 4 

1000 2000 3000 4000 5000 

Landward 

direction 

5000 

DISTANCE ALONG SEGMENT AXES, IN FEET FROM LANDWARD END 

4000 3000 2000 

FIGURE 8.-Predicted and observed longitudinal distributions of dye concentration along the axes of segments 2, 4, and 5 at 

approximately 1640 hours on May 21, 1910 hours on May 21, and 0640 hours on May 22, 1969. 
METHOD l!~OR PREDICTING POLLU'l'AN'l' DISPERSION IN BOLINAS LAGOON, CALIFORNIA B13 

~fay 21, has moved more landward than was predicted. In part, this may result from errors in the 

tidal hydraulics program. Also, as in figure 7, the 

observed samples were taken on the channel centerline, whereas the prediction is for the cross-sectional 

mean. Had samples been taken along the channel sides, 

a significantly different observed graph would have 

been obtained. Unfortunately, time did not permit obtaining enough samples to define the cloud completely. 

The second traverse, at approximately 1910 hours on 

May 21, shows a good agreement between prediction 

and observation in segment 5; however, the observed 

cloud has already passed from segment 2 into segment 

1, sooner than predicted. The third traverse, taken at 

approximately 0640 hours on May 22, after the beginning of the main tidal outflow, also shows reasonable agreement with the prediction in segment 5, but 

in segments 2 and 4 significantly more mixing has 

taken place than was anticipated. In part, the low 

concentrations in segments 2 and 4 may result from 

mixing with clear water from segment 10 during the 

prolonged period of high water. 

As mentioned previously, the dye-concentration predictions shown for both experiments (figs. 4, 7, and 8) 

were made using a mass-transport coefficient of E = 0.2. 

Two additional predictions were made, one using E = 

0.0, to observe the magnitude of the numerical disper~ 

sion, and one using E = 0.4, to observe the effect of 

changing the coefficient. Figure 9 shows the results for 

E = 0.0 for two sites during the experiment of June 

18-19, 1968. The figure also shows, for comparison, the 

observed data and computer results using E = 0.2. The 

figure shows that at site 5 little dispersion is observed if E = 0.0 (if the program were perfect it would 

yield no dispersion at all). During both tidal inflow 

and outflow a high peak is observed at one of the 

hourly readings, and almost zero concentration is ob-· 

served at all other times. During tidal outflow at site 

2, two peaks are observed: the first is from concentration stored at high tide in segment 5 and the second 

is from concentration from segments 1 through 4. The 

different times of return of the two peaks do not reflect numerical dispersion, but are caused by different 

rates of convection in the different segments. Thus 

figure 9 shows that if the computer program is run 

with no dispersion intended, little dispersion results; 

computer-generated numerical dispersion is not significant in this program. 

'Vhen the computer program was run with E = 0.4, 

predictions were obtained showing somewhat, but not 

significantly, greater dispersion than for the runs with 

E = 0.2. The difficulty resulted from the adjustment 

of the transport coefficient to account for element 

100 

80 

60 

40 

a:: 

LLJ 

1-

::::i 

a:: 20 LLJ 

a. 

(/) 

~ 

<( 0 a:: 

(.!) 

0 

a:: 

~ 60 

~ 

~ 

z 50 0 

i= 

<( 

a:: 

1-

z 40 LLJ 

u 

z 

0 

u 30 

20 

10 

EXPLANATION 

...__. Observed concentration 

o---o Predicted concentration, £=0.0 

o---o Predicted concentration, E =0.2 

SITE 2 

2200 

SITE 5 

2400 0200 0400 

J U N E 1 9, 1968 

i 

1\ 

II 

I\ 

II 

1\ 

I I 

I I 

I I 

I I 

1400 1800 2200 2400 0200 0400 

JUNE 18, 1968 JUNE 19, 1968 

.l!,IGURE 9.-0bserved and predicted dye concentrations at two sampling sites, June 18-19, 1968, 

showing effects of changing the mass-transport 

coefficient in the numerical program. 

travel and segment width, as given by equation 8, 

and the requirement that the adjusted coefficient E", 

not exceed 0.5. For the runs with E = 0.2 the adjusted 

coefficient did not often exceed 0.5 ; for instance, for 

the experiment of June 18-19, 1968, out of 3,034 computations of E", 716 gave a value greater than 0.5 

(24 percent). Of the remainder, 224 were greater than 

0.4, 338 were between 0.3 and 0.4, 391 between 0.2 

and 0.3, 588 between 0.1 and 0.2, and 777 less than 

0.1. If the value of E is doubled, however, the computed values of E" also double, and many more values 

of E" greater than 0.5 are obtained. Since all these 

must be reduced to 0.5 to maintain numerical stability, 

the rate of dispersion of the cloud is not greatly increased over that for E = 0.2. 
B14 DISPERSION IN SURFACE "\VATER 

SUMMARY OF EXPERIMENTAL RESULTS 

Two dye-dispersion experiments were made using 

the same procedure in the same lagoon, but they were 

for tidal eycles of significantly different character. The 

first experiment was conducted during a symmetrical 

121j2 hour cycle, during which the tidal elevation at 

the mouth rose from a low to a high and then returned 

to the same low. The results showed that circulation, 

implying a circulatory motion of water within the 

interior of the lagoon, was not an apt description of 

the flow caused by this type of tide. Water to supply 

the high-tide volume entered the lagoon by flowing 

up the several channels. At high tide there was little 

motion and little mixing. Although perhaps different 

in small-scale detail, in gross pattern the outflow appeared to be simply the reverse of the inflow. Water 

flowed out through the same channels through which 

it flowed in. 

Ritter (1969) dropped dye from a helicopter to observe motion in the lagoon during a tidal cycle similar 

to that of the first experiment. (Plate 1 is a photograph taken during his study and shows some of his 

dye patches.) The gross patterns of flow determined 

by Ritter and during the present experiment were 

similar; for instance, Ritter stated that dye deposited 

in the Bolinas channel did not cross over into areas 

apparently filled through the Stinson channel. Ritter 

remarked that the flow appeared to be more channelized on inflow than during outflow; if true this would 

not invalidate the present approach, although perhaps 

the appropriate values of E might undergo some 

change. Unfortunately, Ritter's experiments did not 

span an entire tidal cycle, and on inflow, dye spots 

were placed only in the channels; so the details of 

the flow remain open for further study. 

The second experiment was conducted during a 25-

hour tidal cycle in which the tide remained relatively 

high during the central 12 hours. While the tide remains high~ changes in tidal elevation produce surges 

within the lagoon which appear to accomplish some 

net circulation between the Bolinas and Stinson channels. The eomputer program does not allow exchange 

between the two channels, and therefore the results 

obtained for the second experiment did not match the 

predictions as well as for the first experiment. Nevertheless, the results of both experiments verified the 

numerical predictions with reasonable accuracy. The 

results show that flow in the lagoon is channelized 

and can be segmented much as was done in the numerical analysis. Thus a flow which at first sight appears 

to be too complicated for analytical study can be 

broken into a number of isolated segments, each of 

which can be analyzed and the results combined. 

NUMERICAL STUDY OF CONTINUOUS 

DISCHARGE AT THE LAGOON MOUTH 

As an example of a practical use of the numerical 

program, a study was made of a continuous discharge 

of pollutant at the lagoon mouth. Two cases were 

studied: in one the pollutant was discharged at a 

constant rate throughout the tidal cycle, while in the 

other the discharge was interrupted during the 2 hours 

of lowest water during each 121h-hour period. The 

comparison shows how much reduction in pollutant 

concentration results from interruption of the pollutant outflow at low tide. 

For each case the computer program was run continuously for ten 25-hour tidal cycles, using the same 

tidal inflow and outflow as for the slug-discharge study 

of :May 21 and 22, 1969. The concentration of all elements inside the lagoon was initially set equal to zero, 

and an addition of a constant number of mass units 

was added to the element nearest the center of segment 12 during each 15-minute time step. The injection point corresponds approximately to the location 

of the existing Bolinas sewage outfall pipe. A decay 

factor of 0.99762 per time step was used, corresponding 

to a rate constant of 0.1 per day, as usually assumed 

for BOD. The concentration units used in this study 

are strictly arbitrary and have meaning only in relation to each other; no attempt has been made to relate 

the study to the actual amount of BOD released into 

the lagoon. 

Figures lOA and B show the concentration distribution at the low high and high high tides resulting 

from the constant-injection study. The highest concentrations are in two areas near the south and north ends 

of the branches fed by the Stinson channel. The high 

eoncentrations are apparently formed at the point of 

injection at low tide, when the volume of water contained in the channel is small and only a small quantity of water is available for dilution. As the tide begins to flood, the water containing high concentrations 

is carried inland. While the incoming discharge is 

small, the concentration of the inflow remains high, 

but as the discharge increases, the input of mass at the 

injection point remaining constant, the concentration 

of the inflow is reduced. The minimum incoming concentration occurs at the point of maximum flow discharge. As the discharge decreases with approaching 

high tide, the concentration again increases, but not 

to such a high level as at low tide because more water 

is available for dilution. 

The water containing the high concentration is convected into the lagoon, and the high concentration 

zone disperses as did the slug injection. At the low 

high tide the water which originated at the lagoon 
METHOD FOR PRl<JDlCTT;'\G I'OLLL"'l'A:\'1' lllSPERSION IN BOLINAS LAGOO:\. CALIFOR:\IA Bl5 

mouth at the low low tide has been carried almost 

to the ends of all the channels. Those areas which 

show as the highest concentrations in figure lOA represent the areas to which the water at the lagoon mouth 

at low tide is convected during the inflow. The areas 

having the lowest concentration at high tide are those 

to which water entering the lagoon during maximum 

flood discharge is convected. At high tide, an area of 

fairly high concentration exists near the lagoon mouth, 

but as soon as the tide turns the water containing this 

concentration will flow out of the lagoon. 

The concentration distribution at the high high tide 

(fig. 10) is complicated by areas of high concentration 

formed at both the low low and high low tides. The 

high concentrations fonnd at the northernmost and 

southernmost ends of the lagoon arc those originally 

formed at low low tide, and the peak concentrations 

are somewhat less than at low high tide because of 

dispersion in the period between highs. A new area of 

high concentration, formed during the high low tide, 

is found in the center of the lagoon. 

The highest concentrations found in the Bolinas 

channel are significantly less than those in areas filled 

through the Stinson channel. During the experiments, 

a sandbar crossed the Bolinas channel at the point 

where the channel joins the lagoon mouth. The sandbar restricted flow into the Bolinas channel during 

periods of low tide, so that the water containing the 

highest concentration bypassed the Bolinas channel 

and flowed almost entirely up the Stinson channel. 

The Bolinas channel >vas filled with water which entered the lagoon sometime after low tide, and, hence, 

contained lower pollutant concentrations. 

Tlw low concentrations shown in figures lOA and B 

at the northernmost and southernmost extremes of the 

lagoon result from the decay factor introduced in the 

numerical study. If the pollutant did not decay according to the given decay factor, diffusion would 

raise the concentration in these extremities to the same 

level as exists in the adjacent zones of high concentration. 

The predicted concentrations shown in figure 10 

have not been verified by experiment for two reasons: 

first. the concentration distribution depends on the 

rate of decay of the pollutant, and no dye is known 

which decays at the same rate as BOD; and second, 

the distribution predicted by the computer is for a 

repetition of a series of identical 25-hour cycles of 

the tidal flow predicted for the experiment of :May 

21-22, 1969, when~as, in reality, the tidal heights vary. 

Nevrrtheless, visual and ecological observations verify 

that the distribution of pollutants in the lagoon in 

1!)68 and 1969 was similar to that shown in the figure. 

Visually, accumulations of sewage-type refuse could 

be observed along the shore at high tide in the areas 

which are shown as zones of consistently high concentration in figures lOA and B. Ecological studies 

by Gustofson (1968) suggested that these zones were 

the least habitable to higher forms of life. Thus, although the computer progmm has not been verified 

in detail for the continuous release and indeed is based 

on an approximate lagoon geometry and only one of 

a great many possible tidal flows, the concentration 

distribution predicted by the program seems to be in 

good agreement with the distribution of pollution observed in the lagoon. 

Fignres 1L1 and B sho"· the change caused by interrnpting the pollutant discharge during the 2 hours of 

lowest tide, beginning one-half hour before low slack 

watrr and continuing until 1% hours after slack water. 

Fignre llA shows the concentration distribution at the 

same time as fignre lOA , and figure llB corresponds 

to figure lOH. The high concentrations evident in fignres lOA and B have been eliminated because thc,se 

concentrations were produced by the discharge during 

the lmv tide period. The comparison of figures lOA 

and R "·ith llA and B sho\YS how much reduction in 

pollutant levels in the more Sl'rionsly polluted areas of 

tlw lagoon mig·ht be achieved by construction of a 

small detention structure to withhold discharge of 

Bolinas se,Yage for relatively short periods at low tide. 

Again. it must be emphasized that this conclusion is 

based on a numerical model using approximate lagoon 

geometry and only one of a number of possible tidal · 

flows. Tn spite of these approximations, it seems reasonabiP to conclude that provision of a detention structnre wonld significantly reduce pollution problems 

'"ithin the lagoon. 

CONCLUSIONS 

The numerical program described predicts the effect 

of tidal inflow and outflow on the distribution of a 

pollutant in a tidal embayment. Use of the program 

requires knowledge of the embayment geometry and 

the tidal elevations at various interior points throughont tlw tidal cycle. The program models the convective motion of water within the embayment and accomplishPs the convective motion of a tracer or pollutant dissoln~d in the water whilr. minimizing numerica I dispersion. Diffnsion of the tracer between 

adjacent elements of water occurs during a diffnsi ve 

stPp at a rat<~ fll'<'.dictcd by physical considerations of 

the n:wehanics of dispersion. 

The nnmrrical program has been verified by dyetracrr experiments in two types of tidal cycle in 

Bolinas Lagoon, Calif. The observed motion and dis-
B16 DISPERSION IN SURFACE WATER 

A 

PACIFIC OCEAN 

0 1000 2000 3000 4000 FEET 

N 

PACIFIC OCEAN 

0 1000 2000 3000 4000 FEET 

EXPLANATION 

CONCENTRATION. IN 

ARBITRARY UNITS 

D 

0- 10 

10- 15 

15-20 

• >20 

Low high tide 

continuous discharge 

Stinson 

Beach 

EXPLANATION 

CONCENTRATION . IN 

ARBITRARY UNITS 

D 

0-10 

10- 15 

15-20 

High high tide 

continuous discharge 

Stinson 

Beach 

FIGURE 10.-Nnmerically predicted concentrations resulting from a constant-rate discharge at tlw 

lagoon mouth, after 10 repetitions of the ticlal cycle of ::\lay 21- 22, 1969. A, at the low high tide: 

B, at the high high ticlP. 
METHOD FOR PREDIC'l'ING POLLUTANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B17 

A 

PACIFIC OCEAN 

0 1000 2000 3000 4000 FEET 

N 

PACIFIC OCEAN 

0 1000 2000 3000 4000 FEET 

EXPLANATION 

CONCEN TRATION . IN 

ARBITRARY UNITS 

0-10 

10-15 

Low high tide 

interrupted discharge 

EXPLANATION 

CONC ENTRATION, IN 

ARBITRARY UNITS 

G 

0-10 

10-15 

High high tide 

interrupted discharge 

FIGURE 11.-Numericall~· predicted concentrations resulting from an interrupted constant-rate discharge at the lagoon mouth, after 10 repetitions of the tidal cycle of May 21-22, 1969. A, at the 

low high tide, if no pollutant is discharged during a 2-hour period at the low and high tides: 

B, at the high high tide, for the same conditions as in A. 
B18 DISPERSION IN RTJRFACE WATER 

persion of the dye cloud followed the prediction of the 

numerical program. The program was also used to 

predict the distribution of a pollutant which is continuously discharged into the mouth of the lagoon 

at the location of the existing sewage outfall pipe. 

The prediction showed that at high tide, high concentrations of sewage would exist at the northernmost 

and southernmost ends of the lagoon, while lower 

concentrations would be found in the center section. 

The prediction, while not verified by a dye-tracer experiment, seemed to conform to visual and ecological 

observations of the areas of greatest pollution. 

The numerical program was used to predict the distribution of concentration if sewage were released continuously at the location of the existing discharge 

pipe but withheld during a detention period of 2 hours 

at low tide- The concentrations predicted within the 

lagoon were significantly lower than if the sewage 

were released continuously without interruption. 

The computer program developed in this study 

should be of value to pollution studies in other tidal 

embayments whose geometry has the same character 

as that of Bolinas Lagoon. The program is suited to 

embayments in which the geometry of flow changes 

significantly between high and low tide. The program 

is one dimensional in the sense that the embayment 

must be segmented into channels which are on dimensional, but these segments may be connected to form 

a two-dimensional embayment. Thus the program is 

capable of describing dispersion in a multibranched 

tidal inlet in which large expenses of mudflat are exposed at low tide but shallowly covered at high tide. 

The program should yield reasonably accurate predictions of the dispersion of a pollutant from either an 

instantaneous slug source or from a constant-rate injection. 

REFERENCES CITED 

Bowden, K. F., 1963, The mixing processes in a tidal estua ry : 

Internat . .Tour. Air and \Vater Pollution, v. 7, p. 343- 356. 

Elder, .J. W. , 1959, The dispersion of marked fluid in turbulent 

shear flow : .Jour. F lnid Mechanics, v. 5, p. 544- 560. 

F'ischer, H. R., 1967, 'l'he mechanics of dispersion in natural 

streams: Am. Soc. Civil Engineers Proc., .Jour. Hydraulics 

Div., v. 93, no. HY6, p. 187-216. 

-·- 1968, :\Iethods fo1' predicting dispersion coefficients in 

natural streams, with application to the Green and Duwamish Rivers, Washington: U.S. Geol. Survey Prof. Paper 

582-A, 27 p. 

Gustofson, J. I<'., 1968, Ecological study, Bolinas Lagoon, Marin 

County, California, report no. 2: Resources & Ecology 

Projects, P.O. Box 87, :Mill Valley, California, 45 p. 

Harleman. D. R. I<~., 1964, The significance of longitudinal 

dispersion in the analysis of pollution in estuaries: Int. 

\Vater Pollution Research Conf., 2d, Tokyo, 1964, Proc., 

v. 1, p. 279-306. 

liolley, E. A., Harleman, D. R. F ., and I<'ischer, H. B., 1970, 

Dispersion in homogeneous estuary flow: Am. Soc. Civil 

Engineers, Proc., Jour. Hydraulics Div., v. 96, no. HY8, 

p. 1691-1709. 

Lai, Chintu, 1967, Computation of transient flows in rivers and 

estuaries by the multiple-reach method of characteristics: 

U.S. Geol. Survey Prof. Paper 575-D, p. D273-D280. 

Orlob, G. T., Selleck, R, Shnbinski, R. P., Walsh, F., and 

Stann, E., 1968, Mo(leling of water quality behavior in 

an estuarial environment. Internat. Water Pollution Research Conf., 4th Prague, Proc., p, 845--855. 

Pritchard, D. w·., 1958, A study of flushing in the Delaware 

model: Tech. rept. 7, The Chesapeake Bay Inst., Johns 

Hopkins Univ. 

Ritter, .J. R., 1969, Preliminary studies of sedimentation and 

hydrolog~· in Bolinas Lagoon, Marin County, California : 

U.S. Geol. Survey open-file report, 68 p. 

Taylor, G. I. , 1954, The dispersion of matter in turbulent flow 

through a pipe : Royal Soc. [London J Proc., ser. A, v. 

223, 6. 446--468. 

'£homann, R. V., and Sobel, M. .J., 1964, Estuarine water 

quality management and forecasting: Am. Soc. Civil Engineers Proc., .J onr. Sanitary Eng. Div., v. 90, no. SA5, p. 

9--36. 
TABLES 1 AND 2 
B20 DISPERSION IN SURFACE WATER 

'.rABLE !.-Annotated flow chart of the computer program 

[Identification of Indices : IT = number of the time step within the tidal cycle; I CYCLE = number of the tidal cycle; M = number Identifying the segment; N = volume element within the segment. Exam pies : C (M,N) Is the concentration of the Nth Yolume element In the M:th 

segment B (M,IT) Is the width of the ~Ith segment during the IT th time step of any tidal cycle] 

Reference No. 

(see table 2) 

Main program 

Read input data: 

For the program as a whole: 

90 Identification of run (IDEN); 

100 Number of segments (MMAX); 

Number of time steps in one tidal cycle (!CYCLE); 

Segment into which a continuous discharge of pollutant is to flow (MINPUT); 

Amount of mass input per time step for continuous discharge (ADD); 

Basic mass-transfer coefficient for diffusive step (E) ; 

Decay factor (DIEOFF) ; 

Position of continuous discharge within segment (XADD); 

Program control information (IPRINT,.TSTOP,NPASS1,NP ASS2). 

For each segment in turn: 

110 Number of elements initially in segment (NMAX(l\1)); 

Length of segment (XL(M)); 

Segments to which given segment is connected in landward direction (MA (M), MB (M), 

MC(M)); 

120 Segment to which given segment is connected in seaward direction (MD(M)); 

Mean depth of segment at each time step throughout tidal cycle (H (M,IT)) ; 

Width of segment at each time step throughout tidal cycle (B (M,IT)) ; 

Initial concentration of each element in segment ( C (M,N)). 

l 

130 Compute initial volume of each element of each segment by dividing total segment volume by number 

of elements initially in segment. 

l 

Read volume of fresh-water discharge per time step into those segments recmvmg a fresh-water discharge, if any. 

l 

200 Begin a time step. Increment IT (time-step index) by 1. If IT becomes greater than the number of steps 

in a tidal cycle (!CYCLE), reset IT to one (begin the next tidal cycle). At the beginning of each 

tidal cycle, increment !CYCLE by 1. When !CYCLE becomes greater than .TSTOP, stop. 

l 

209 Executive subroutine MAP for each segment in turn, beginning with segment 1 and continuing through 

segment MMAX. MAP computes the volume transferred into or ont of each segment during the time 

step and the concentration of pollutant in each volnme. 
Reference No. 

(see table S) 

METHOD FOR PREDICTING POLLUTANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B21 

TABLE 1.-Annotated flow chart of the computer program-Continued 

Main program-Continued 

212 During the first tidal cycle only, print the volume transferred across the landward end of each segment 

(QO(M) ), the volume transferred across the seaward end (QI(M) ), and the average distance 

traveled by a water particle within the segment (VEL(M) ). 

1 

220 Execute subroutine BRANCH for each segment in turn. BRANCH adds or subtracts volume elements 

from the ends of each segment, according to the inflow or outflow calculated by MAP. It also calculates the concentrations of newly formed volume elements. · 

1 

250 Execute subroutine DIFUSE for each segment in turn. DIFUSE carries out mass transfer between volume elements within each segment, as described in the text in the section "Diffusive Step." It also 

accomplishes the reduction of volume elements within a segment wh~n this is required. DIFUSE 

prints the concentrations and locations of volume elements when requested by the program contra] 

information. 

1 

255 Add the pollutant mass given as ADD to the ·volume element selected by MINPUT and XADD. This 

simulates a continuous injection at a fixed location. For simulation of a slug injectio'n ADD is set 

equal to zero. 

1 

260 Reduce the concentration of each volume element by the factor given as DIEOFF to simulate a first-order 

decay process. 

1 

Return to the beginning of the time step. 

Subroutine MAP 

70 Compute the volume transported through the landward end of the segment (QO(M)) by summing the 

volume transported through the seaward ends of the segments connected on the landward side and the 

fresh-water discharge if any, by the equation 

QO(M) = QI(Ml) + QI(M2) + QI(M3); 

Ml, M2, and ~13 are the landward-connected segments, or in the case of a fresh-water discharge QI(M3) 

is a fresh-water discharge stated in the main program. 
B22 

Reference No. 

(see table ~) 

DISPERSION IN SURFACE WATER 

TABLE l.-Annotated flow chart of the computer program-Continued 

Subroutine MAP-Continued 

1 

80 Compute the volume transported through the seaward end of the segment (QI(M)) by equation 1. 

! 

90 If the volume transported through eith~r 

end of the segment is greater than the 

volume contained in the segment, the program will not operate. Print inflow, outflow, identification~ and stop. 

! QO(~f) >0 

! 

135 If flow is out of the landward end of the 

segment, compute the concentration of the 

outgoing volume by equation 4. 

l 

QI(M) <0 

! 

245 If flow is out of the seaward end of the 

segment, compute the concentration of the 

outgoing volume by equation 4. 

~ 

! 

QO(M)LO 

l 

QI(M)>O 

300 Return to the main program. The main program increments the segment index, M, by 1 and re-executes 

subroutine MAP until all segments are complete. 

Subroutine BRANCH 

Initiate for each segment in turn from the main program. 
METHOD FOR PREDICTING POLLUTANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B23 

TABLE 1.-Annotated flow qhart of the computer program-Continued 

Subroutine BRANCH -continued 

Reference No. 

(see table 2 ) 

! ! 

QO(M)<O QO(M) =0 

! 

120/230 If flow is into the landward end, create a new 

V(M,1) and compute C(M,1) by proportioning the concentrations of the inflowing 

volumes. Increase the identifying index of 

each volume element previously in the 

segment by 1, that is, V(M,1) becomes 

V ( J\1,2) , and so forth. 

! 

QI(M) >0 QI(M) =0 

! 

500/420 If flow is into the seaward end of the segment, create a new volume element, V(M,-

N + 1), whose volume equals the inflow and 

whose concentration, C (M,N + 1), is that 

calculated for the inflow by MAP. 

! 

QO(M)>O 

~ If flow is out of the landward end, reduce 

V (M,1) by the amount of the outward 

flow. If V (M,1), is reduced to zero, reduce 

V(M,2) by the remainder, and reduce the 

identifying index of each element by 1 ; 

that is, V (M,2) becomes V (M,1), and so 

forth. If V(M,2) is reduced to zero, reduce V(M,3) by the remainder, and let 

V(M,3) become V(M,1), and so forth. 

t 

! 

QI(J\1) <0 

! 

If flow is out of the seaward end, reduce 

V(M,N) by the amount of the outward 

flow. If V (M,N) is reduced to zero, reduce V(M,N-1) by the remainder, and 

reduce the number of v.olume elements in 

the segment by 1. If V (M,N -1) is reduced to zero, reduce V(M,N -2) by the 

remainder, and so forth. 

600 Return to the main program. The main program increments the segment index, M, by 1 and re-executes subroutine BRANCH until all segments are complete. 

Subroutine DIFUSE 

Initiate from the main program for each segment in turn. 

l 

100 Compute the length (XE(N) of each volume element by dividing the volume (V(M,N)) by the 

cross-sectional area of the segment (BH (M,IT)). 
B24 DISPERSION IN SURF ACE WATER 

TABLE l.-Annotated flow chart of ~he computer program-Continued 

Subroutine DIFUSE-Continued 

Reference No. 

( see table :e ) 

l 

110 Compute the distance from the landward end of the segment to the center of each volume element 

(X(N)). 

1+--------

300 If the combined length of the two endmost elements on either end of the segment is less than the 

or average length of one element, that is less than the length of the segment divided by the num400 her of elements within the segment, combine the two endmost elements into one whose concentration is the weighted average of the two being combined. Return to the beginning of this 

step. 

120 or 130 

150 

160 

170-230 

1 

Compute the mass transfer between adjacent elements, beginning with element pair (M,2), (lV£,1) 

and continuing to pair (M,N), (M,N -1). For each pair, carry out the following steps: 

1. Determine which of the two elements is shorter. 

2. Using the length of the shorter element, compute the adjusted mass-transfer coefficient 

by equation 7. 

3. Compute the mass transfer. 

4. Adjust the concentration in each element to accomplish the mass transfer. 

1 

260 If called for by the program control information, print the location and concentration of each 

volume element. 

290 If the element being considered is one into which a continuous injection is to be put, compute the 

number of the element whose concentration is to be increased. 

1 

Return to the main program. The main program increments the segment index, M, by 1 and reexecutes subroutine DIFUSE until all segme,nts are complete. 
METHOD FOR PREDICTING POLLU'l'ANT DISPERSION IN BOLINAS LAGOON, CALIFORNIA B25 

TABLE 2.-Fortran listing of computer program for CDC 61,00 computer 

PROGRAM lAGCON(lNPUT,OUTPUT) OIMENSIO~ MA(l2),M8(12l,MC(l2l,M0(12),VEL(12),10EN(l3) 

COMMON/A/NMAX(l2l,Xl(l2),B(l2,125),H(l2,125),V(l2,50),C(l2,50), 

1QI(25),Q0(20),8H(l2,125l,C0(2l),CI(25),JP,NO(l2ltNI(l2), 

2IPRINT,IT,ITPl,EFOR,EBACK,ICYCLE 

COMMON/B/NAOO,XAOO 

COMMON/C/MINPUT 

OAT.&\ NMAX ( 1) /0/ 

c-----------------------------------------------------------------------

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

THIS PROGRAM CONVECTS FINITE ELEMENTS OF FLUID WITH ACCOMPANYING 

CONCENTRATION OF A QUALITY PARAMETER, AND ALLOWS DIFFUSION BETWEEN 

ADJACENT ELEMENTS, IN A ONE-DIMENSIONAL, MULTI-BRANCHED ESTUARY 

B 

H 

v 

c 

XL 

NMAX 

M 

MMAX 

IT 

I CYCLE 

QI 

QO 

CI 

co 

MINPUT 

ADO 

= WIDTH 

= MEAN DEPTH 

= ELEMENT VOLUME 

= CONCENTRATION IN ELEMENT 

= LENGTH OF SEMGENT 

= NUMBER OF ELEM~NTS IN SEGMENT 

= SEGMENT NUMBER 

= NUMBER OF SEGMENTS IN ESTUARY 

= NUMBER Of TIME STEP 

= NUMBER OF TIME STEPS IN ONE TIDAL CYCLE 

= FLOW INTU SEAWARD END OF A SEGMENT = FLOW OUT OF LANDWARD END OF A SEGMENT 

= CONCENTRATION OF FLOW LEAVING SEAWARD END OF A SEGMENT 

= CONCENTRATION OF FLOW LEAVI~G LANDWARD END OF A SEGMENT 

= NUMBER OF SEGMENT RECEIVING A CONSTANT DISCHARGE 

= AMOUNT OF DISCHARGE ADDED PER TIME STEP 

c-----------------------------------------------------------------------

90 

100 

105 

lLO 

112 

114 

REAo 90,(10EN(l),I=l,l3) 

FORMAT(l3A6) 

READ lOO,MMAX,ICYCLE,MI~PUT,AOO,IPRINT,JSTOP,EFOR,EBACK,OIEOFF, 

lXAOO,~PASSl,NPASS2,NREAO,ITSLUC 

FORMAT(313,Fl5.0,213,2F5.0,2Fl0.0,21l,I3,15) 

PRINT 105,(1DEN(Ii,I=l,l3),0lfOFF 

FORMAT(H ,*BEGIN EXECUTION ----- *tl3A6t* DECAY CONSTANT = *• 

lfl0.8) 

IR=ICYCLE/NREAO 

IRMl=IR-1 

NREADMl=NREA0-1 

DO 121 M=l,MMAX 

READ 110 ,NMAX(M),XL(M),MA(M),MB(M),MC(M),MO(M) 

FORMAT(lOX,I2,Fl0.0,413) 

READ 120 ,(H(M,l*lR),B(M,I*IR),I=l,NREAO) 

DO 112 IT=IR,ICYCLE,IR 

BH(M,ITl=8(M,IT)*H(M,IT) 

00 114 IT=l,IRMl 

B (M,IT)=B (M,ICYCLE)+IT*(B (M,IR)-8 tM,ICYCLE))/IR 

BH(M,ITl=BH(M,ICYCLE)+IT*(BH(M,lR)-BH(M,ICYCLE))/IR 

H(M,IT)=BH(M,IT)/B(M,IT) 

CONTINUE 

DO 116 I=l,NREADMl IPl=l*l~+l 
B26 DISPERSION IN SURFACE WATER 

TABLE 2.-Fortran listing of computer program for CDC 61,.00 computer-Continued 

IEND=I*IR+IRMl 

DO llb IT=IPl,IENO _ B (M,IT)=B (M,I*IR)+(IT-I*IR)*(B (M,I*I~+IRl-B (M,I*IR))/lR 

BHCM,ITl=BH(M,I*IR)+(IT-I*IR)*(BH(M,I*IR+IR)-BH(M,I*IR))/IR 

H(M,ITJ=BHCM,IT)/B(M,IT) 

116 CONTINUE 

120 FORMAT(l0F8.0) 

121 CONTINUE 

IFCNPASSl.NE.O) GO TO 128 

00126 M=l,MMAX 

'PRINT 122,M,Xl(M),MA(M),MB(M),MC(MJ,MO(M) 

122 FORMATClHO,*BRANCH *tl2t*t LENGTH*•f12.2t*t LANDWARD BRANCHES *• 

l3(12,H,J,• SEAWARD BRANCH *,12,/J 

PRINT 124,(HCM,IT),B(M,ITJ,IT=IR,ICYCLE,IR) 

124 FORMAT(6(2X,F6.2,Fl0.1,4X)) 

126 CONTINUE 

128 CONTINUE 

DO 140 M=1,MMAX 

VOL=BHCM,ll*XL(M) 

V(M,l)=VOL/NMAX(M) 

NT=NMAX(M) 

DO 130 N=2,NT 

130 VCM,N)=V(M,l) 

140 CONTINUE 

DO 160 M=1,MMAX 

NT=NMAXCM) 

READ l20,(C(M,N),N=l,NT) 

160 CONTINUE 

JCYCLE=l 

IT=O 

DO 165 M=1,21 

COCM)=O. 

165 CI(M)=O. 

Ql(21)=0. 

198 JP=O 

c-------~---------------------------------------------------------------

c 

C BEGIN A TIME STEP 

c 

c-----------------------------------------------------------------------

200 CONTINUE 

IT=IT+l 

JP=JP+l 

IF(IT.GT.ICYCLE)JCYCLE=JCYCLE+l 

IFIIT.GT.ICYCLE) IT=1 IF(I~.EQ.ITSLUG.AND.JCYCLE.EQ.l) GO TO 300 

201 CONTINUE 

IFIJCYCLE.GT.JSTOP)STOP 

ITPl=IT+l 

lf(IT.EQ.ICYCLE) ITP1=1 

IFCJP.LT.IPRINT) GO TO 208 

PRINT 202,JCYCLE,IT 

202 FORMAT(lHl,*RESULTS AT CYCLE *tl2t*t TIME STEP *ti3,/) 

208 CONTINUE 

DO 209 M=l,MMAX 

CALL MAP (M,MA(M),MB(M),MC(MJ,MO(M)) 
METHOD FOR PREDICTING POLLU'l'AN'l' DISPERSION IN BOLINAS I_JAGOON, CALIFORNIA B27 

TABLE 2.-Fm·tran Usting of computer program for ODO 6!,00 computer-Continued 

209 CONTINUE 

lf(JCYCLE.EQ.t.AND.NPASS2.EQ.O) GO TO 210 

IF(JP.LT.IPRINT) GO TO 218 

210 CONTINUE 

00 211 M=l,MMAX 

211 VELCM)=.5*(00(M)+QI(M))/BH(M,ITP1) 

PRINT 212,IT 

212 FORMAT(H0,2(5X,*BRANCH INFLOW OUTFLOW TRAVEL* 

1,7X),*IT = *tl2,/) 

PRINT 214,(M,QI(M),Q0(M),VEL(M),M=l,MMAk) 

214 FORMAT((X,2(5X,I4,X,Fl3.0,X,Fl4.0,X,Fll.2,7X))) 

PRINT 216 

216 FORMAT(lH ) 

218 CONTINUE 

DO 220 M=l,MMAX 

CALL BRANCHIM,MA(M),MB(M),MC(M),MO(M)) 

220 CONTINUE 

DC 250 M=l,MMAX 

CALL DIFUSf(M) 

250 CONTINUE 

CCMINPUT,NAOD)=CtMINPUT,NAOO)+AOO/VtMINPUT,NADO) 

255 CONTINUE 

DO 260 M=l,MMAX 

NT=NMAX(M) 

08 260 N=l,NT 

260 C(M,N)=OIEOFF*C(M,N) 

lf(JP.LT.IPRI~T) GO TO 200 

GO TO 198 

.300 NT=NMAX ( 9) 

C(9,NT)= 500. 

GO TO 201 

END 

SUBROUTINE MAP(M,Ml,M2,M3) 

COMMON/A/NMAXI12),Xl(l2),8(12,125),H(l2,125),V(l2,50),C{l2,50), 

1QI(25J,Q0(20),RH(l2,125),C0(2l),CI(25),JP,NO(l2l,NI(l2), 

21PRINT,IT,ITPl,EFOR,EBACK,ICYCLE 

COMMON/C/MINPUT 

N T=i\jMAX ( M) 

NO(M)=O 

NI(M)=O 

QO(M)=QI(Ml)+QI(M2)+Ql(M3) 

Ql(M)=QO(M)+(6H(M,ITPll-8H(M,ITll*Xl(M) 

VOLEL=BH(M,IT)*XL(M) 

If( QO(M).GE.VOLEL) GO TO 380 

IFC-QI(M).GE.VOLEL) GO TO 390 ~0 CONTINUE 

c-----------------------------------------------------------------------

c 

C COMPUTE CONCENTRATION OF OUTFLOW FROM LANDWARD END OF SEGMENT 

c 

c-----------------------------------------------------------------------

IFlQOtMlJ 700,200,100 

100 CONTINUE 

VOL=O. 

DO 110 N~l,NT 

NTEMP=N 
B28 DISPERSION IN SURFACE WATER 

TABLE 2.-Fortran listing of computer program tor ODO 6400 computer-Continued 

VOL=V(M,N)+VOL 

lf(OO(Ml.LT.VOL) GO TO 120 

110 CONTINUE 

120 N=NTEMP 

NO(M)=N 

SUM=O. 

IF(N.EQ.l) GO TO 135 

00 130 NTEMP=2,N 

130 SUM=SUM+C(M,NTEMP-ll*V(M,NTEMP-11 

135 CO(M)=(SUM+C(H,N)*(V(M 9 N)+QO(Ml-VOL))/QO(M) 

c---------~-----------------------------~------------------------------- c 

C COMPUTE CONCENTRATION Of OUTFLOW FROM SEAWARD END Of SEGMENT 

c 

c----------------------------------------------------------------------- 200 IF(Ql(M)) 210,300,300 

210 CONTINUE 

VOL=O. 

DO 220 NINV=lrNT 

NTEMP =NINV 

N=NT+l-NINV 

VOL=V(M,Nt+VOL 

lf(-QI(MJ.LT.VOL) GO .TO 230 

220 CONTINUE 

230 NINV=NTEMP 

NI(Ml=NINV 

SUM=O. 

IFININV.EQ.l) GO TO 245 

00 240 NTEMP=2,NINV 

NINVT = NT+l-NTEMP 

240 SUM=SUM+C(M,NINVT+ll*V(M,NINVT+ll 

245 CICM)=(SUM+CCM,Nl*(V(M,N)-Ql(M)-VOL))/(-QI(M}) 

300 RETURN 

380 lf(M.EQ.MINPUT) GO TO 500 

GO TO 400 

390 IFCM.EQ.MINPUT) GO TO 600 

400 CONTINUE 

PRINT 4lO,M,QI(M),QO(M),IT 

410 FORMATClHO,*TRANSFER EXCEEDS ELEMENT VOLUME, BRANCH *tl2, 

l*t INFLOW =*tfl3.0,*t OUTFLOW =*tfl3.0t*• TIME STEP *•12) 

STOP 

500 CONTINUE 

V(M,NT+l)=QO(M)-VOLEL 

C(M,NT+l)=O. 

QI(M)=BH(M,ITPll*XL(M) 

NT=NT+l 

NMAX(Ml=NMAX(M)+l 

GO TO 90 

600 CONTINUE 

QI(M)=-VOLEL 

QO(M)=-BH(M,ITPl)*XL(M) 

GO TO 90 

700 IF(QI(Ml).GT•O.) CO(Ml=CICM2) 

IF(Ql(M2J.GT.O •• OR.QI(M3).GT.O.) CO(M)=CI(Ml) 

GO TO 200 

END 
METHOD FOR PREDIC'l'l~G POLLUTANT DISPERSION IN BOLINAS LAGOOX, CALIFORNIA B29 

'l'ABI.E 2.-Fortran listing of computer program tor ODO 6400 oomp1tter-Continued 

SUBROUTINE BRANCHCM,Ml,M2,M3,M4) 

DIMENSION VNEW(50),CN.EW(50)_ 

COMMON/A/NMAX(l2),Xl(l2),B(12,125),H(l2,125),V(l2,50),C(l2t50), 

1Q1(25),Q0(20),BH(l2,125),COC2l),Cl(25),JP,N0(12),Nitl2), 

21PRINT,IT,ITPl,EFOR,EBACK,ICYCLE 

NT=NMAX(M) 

IF(QO(M)) 100,300,200 

100 CONTINUE 

DO 110 N=l,NT 

VNEW(N+l)=V(M,N) 

110 CNEW(N+l)=CCM,N) 

NT=NT+l 

00 120 N=2,NT 

V(M,N)=VNEW(N) 

120 C(M,Nl=CNEW(N) 

VCM,ll=-QOCMI 

C(M,l)=(Ql(Ml)*CI(Ml) + QI(M2)*CI(M2) + Ql(M3)*CI(M3))/QO(M) 

IFCQI(Mll.GT.O.) C(M,l)=CI(M2) 

lf(QICM2).GT.O •• OR.QI(M3).GT.O.) C(M,l)=CI(Ml) 

GO TO 300 

200 CONTINUE 

NOM=NO(M) 

NOMMl=NOM-1 

VOL=O. 

DO 210 N=l,NT 

VOL=VOL+V(M,N) 

IFtN.GE.NO(H)) GO TO 220 

210 CONTINUE 

220 V(M,N)=VOL-QO(M) 

00 230 NTEMP=N,NT 

VCM,NTEMP-NOMMl)=V(M,NTEMP) 

230 C(M,NTEMP-NOMMl)=C(M,NTEMP) 

NT=NT-NOMMI 

300 CONTINUE 

lf(QI(M)) 400,600,500 

400 CONTINUE 

NIM=Nl(M) 

VOL=O. 

DO 410 NINV=l,NT 

N=NT+l-NINV 

VOL=VOL+V(M,N) 

IF(NINV.GE.NIM) GO TO 420 

410 CONTINUE 

420 VCM,N) = VOL+QI(M) 

NT=NT-NIM+l 

IF(V(M,NT).LE.O.) NT =NT-1 

GO TO 600 

500 VIM,NT+l)=QI(M) 

C(M,NT+l)=CO(M4) 

NT=NT+l 

600 CONTINUE 

NMAX(M)=NT 

RETURN 

END 

SUBROUTINE DIFUSEfM) 

DIMENSION 
B30 DISPERSION IN SURFACE WATER 

TABLE 2.-Fortran listing of computer program tor ODO 6400 computer-Continued 

1XE(5Q),X(50),TR(50),EF(50),TRANS(50),CP(50) 

COMMO~/A/NMAX(l2),Xl(12).8(12tl25),H(l2,125),V(l2,50),C(l2,50), 

1Ql(25),Q0(20),8H(l2,125),C0(21J,Cit2S),JP,NO(l7),NI(12), 

2IPRINT,IT,ITPl,EFOR,EHACK,ICYCLE 

COMMON/B/NAOO,XAOO 

COMMON/C/MINPUT 

:JO CONTINUE 

NT=NMAX(M) 

NTMl=NT-1 

NT2=NT*2 

c-----------------------------------------------------------------------

c 

C COMPUTE LENGTHS OF ELEMENTS (KE) AND POSITION OF MIDPUli~TS (X) 

c 

c------------------------------~---------------------------------------- 00 100 N=l,NT 

XE(N)=V(M,N)/BH(M,ITPl) 

100 IF(XE(N).fQ.O.) XE(N)=XL(M)/100000. 

X(l)=.5*XI::(l) 

IF(NT.EQ.l) GO TO 115 

00 110 N=:?.,NT 

110 X(N)=.5*(XE(N-l)+XE(N)) + X(N-1) 

IFCNT.GT.lO.AN0.1XL(M)/NT).LT.(R(M,ITP1)*.333)) GO TO 700 

IF((XL(M}/(X(2)-X(l)}).GT.NT2) GO TO 300 

IF((·XL(M)/{X("4Tl-XtNT-l))).GT.NT2) GO TO 400 

lf(NT.GE.SO) GO TO 600 

115 CONTINUE 

c-----------------------------------------------------------------------

c 

C COMPUTE AOJUSTFO TRANSFER COEFFICIENT 

c 

c-----------------------------------------------------------------------

oo 160 NTEMP=l,NTMl 

IFtXEtNTEMP+l)-XE(NTEMP))l20,130,130 

120 N=NTEMP+l 

Go ro 140 

130 N=NTEMP 

140 CONTINUE 

TR(N)=ABS((QQ(M)*(XL(M)-X(N))+Ql(M)*X(N))/(XL{M)*BH(M,ITPl))) 

TEMP=B(M,ITPl)/XE(N) 

EF(NTEMP>=EFOR 

EF(NTEMP)=EFOR*TEMP 

150 CONTINUE 

EF(NTEMP)=EFCNTEMP)*TR(N)/XE(N) 

IF(Ef(NTEMPl.GT •• S) EFtNTEMP)=.S 

TRANS(NTEMP)=EF(NTEMP)*(C(M,NTEMP+l)-C(M,NTEMP)) 

160 CONTINUE 

c--------------------------------~--------------------------------------

c 

C ADJUST ELEMENT CONCENTRATIONS TO ACCOMPLISH DIFFtJSIVE STEP 

c 

c-----------------------------------------------------------------------

oo 190 N=l ,t~TMl 

IFIXE(N+ll-XE(N)) 170,180,180 

170 CP(N)=C(M,N)+fRANS(N)*XE(N+l)/XE(N} 

GO TO 190 
180 

190 

210 

220 

230 

1\fE'l'HOD FOR PREDICTING POLLU'l'AN'l' DISPERSION IN BOLINAH LAGOON, CALI:H'ORNIA B31 

TABLE 2.-Fo.rtran listing of computer program for ODO 6,f00 computer-Continued 

CP(N)=C(M,N)+TRANS(N) 

CONTINUE 

CPINT)=CtM,NT) 

00 230 N=2,NT 

IF(XE(N-1)-XE(N)) 210,220,220 

C(M,N)=CP(N)-TRANS(N-ll*XE(N-11/XEIN) 

GO TO 23<' 

C (M, N) =CP.( Nl -TRANS ( N-1) 

CONTINUE 

C(M,l)=CP(l) 

c-----------------------------------------------------------------------

c 

c 

c 

PRINT RESULTS AT INTERVALS GIVEN BY IPRINT 

c-----------------------~-----------------------------------------------

260 

270 

275 

280 

290 

IFIJP.LT.IPRINTt GO TO 275 

PRINT 260,M,B(M,ITPl),NT,NO(M),H(M,ITPll 1 NI(M),XL(M) 

FORMATClHO,*BRANCH *tl2t*t WIDTH =* 1 Fl0.2,8Xtl2t* ELEMENTS*, 

18X,*NO = *ti2,/,14X,*OEPTH =*,Fl0.2,27X,*NI = *·,12,/, 

ll4X,*LENGTH =*,Fl0.2,/) 

PRINT 270,(X(N),C(M,Nl,N=l,NT) 

FORMAT((5X,61Fl0.2,F8.2,2X))) 

CONTINUE 

IF(M.NE.MINPUTJ RETURN 

DO 280 N=l,NT 

NTEMP =N 

IF(X(N).GE.XAOO) GO TO 290 

CONTINUE 

NADO=NTEMP 

IF((XADO-X(NTEMP-l)).LT.(X(NTEMPJ-XAOO)) NADD=NTEMP-1 

RETURN 

c-----------------------------------------------------------------------

c 

c 

c 

c 

c 

c 

RETURN TO MAIN PROGRAM. SUBSEQUENT STEPS ACCOMPLISH REDUCTION 

OF NUMBER OF ELEMENTS IN SEGMENT WHENEVER TWO SMALL ELEMENTS 

ACCUMULATE AT EITHER END OF SEGMENT, OR WHEN NUMBER OF ELEMENTS 

IN SEGMENT WILL EXCEED ALLOWABLE STORAGE. 

c----------------------------------------------------------------------- 300 VT=V(M,l)+V(M,2t 

310 

400 

500 

520 

530 

600 

C(M,l)=(C(M,l)*V(M,l)+C(M,2)*V(M,2))/VT 

V(M,l)=VT 

DO 310 N=3tNT 

V(M,N-ll=V(M,N) 

C(M,N-l)=C(M,N) 

GO TO 500 

VT=V(M,NTMl)+V(M,NT) 

CCM,NTMl)=(V(M,NTJ*CfM,NT)+V(M,NTMl)*ClM,NTMl))/VT 

V(M,NTMl)=VT 

NMAX(M)=N-T-1 

IF(NT.EQ.l) GO TO 520 

GO TO 90 

PRINT 530 

FORMAT(lHO,*ELEME\IT REDUCTION STOP*) 

STOP 

PRINT 610 
B32 DISPERSION IN SURF ACE WATER 

TABLE 2.-Fortmn listing of computer program for ODO 6400 computer-Continued 

610 FORMAT(lHO,*FORCED ELEMENT REDUCTION*) 

GO TO 300 

700 CONTINUE 

NWANT=l 

XE:MIN=XE(l) 

DO 720 N=2,NT 

IF(XE(N).LT.XEMIN) GO TO 710 

GO TO 720 

710 NWAfiT=N 

XEMIN=XE(N) 

720 CONTINUE 

IF(NWANT.EQ.l) GO TO 300 

lf(NWANT.EQ.NT) GO TO 400 

IF(XE(NWANT+l)-XE(NWANT-1)} 760,760,740 

740 Ll=NWANT-l 

LU=NWANT 

GO TO 770 

760 LL=NWANT 

LU=NWANT+l 

770 CONTINUE 

VT=V(M,LL)+V(M,LU) 

C(M,Ll)=(V(M,Lll*C{M,LL)+VlM,LU)*C(M,LU)l/VT 

V(M,Ll)=VT 

DO 780 N=LU,NTMl 

V(M,N)=V(M,N+l) 

780 C(M,Nl=C(M,N+l) 

GO TO 500 

END