current research
Mathematical Model for Fate of Pollutants in Near-Shore Waters Walter J. Wnek and Edward G. Fochtman' IIT Research Institute, 10 W. 35th St., Chicago, IL 60616
Mathematical concepts necessary to calculate the fate of pollutants in the near-shore waters of a lake are developed and an example calculation involving three separate sources calculated. The model considers variable diffusion coefficients and variable velocities. Use of transforms and Green's function avoids the finite differences approach, permits rapid calculation, and results in a stable answer for long distances from the source.
The equation based on this concept, which governs distribution of a pollutant, is given by the partial differential equation : dc dt
+ Dz-aa cx +
dC by
cy-
+ Cz-dCb z
=
subject to the initial condition
L
ake Michigan has long been used as a means for disposal of sewage and industrial waste. There is a clear need to reverse the trend of deteriorating water quality and to preserve this valuable resource. Recent efforts to abate pollution have focused attention on the complex nature of the problem and the lack of reliable information on the fate and persistence of pollutants introduced into the lake. Such information is needed to direct antipollution efforts rationally so that those steps can be taken that will result in the greatest effect for their cost. While many scientists have studied this problem on both an experimental and theoretical basis, there is still a great need for more definitive predictive techniques. In this work, available mathematical techniques were modified to develop a program which requires little computer time and is more amenable to prediction of conditions at greater distances or times from the source. Effects of multiple sources, boundaries, and reaction or decay rates were also considered.
C(X,Y,Z,O)
= 0
and the boundary condition of no mass flux across the bounding surface S
where c is concentration; x , y , z are coordinates measured respectively in the direction of the shore, across the lake, and vertically downward (Figure 1); t is time; n is normal to S; L ; ~cy, , and uz are velocities for the x,y,z directions, respectively, functions of x , y , z ; K,, K,, K, are turbulent dispersion coefficients for x,y,z directions respectively, functions of x,y,z; and R is rate of reaction of pollutant with its environment. This equation is arrived at by performing a material balance on a differential element of volume anywhere in the lake.
Qualitative Description of Model
At the beginning of the study, possible dispersion mechanisms were selected, and consequences of these evaluatede.g., molecular diffusion, velocity gradients, and turbulence. The mechanism which provided the best explanation of the dispersal process was one based on the concept of dispersion as an eddy diffusion process. Although phenomena of the diffusion type have been observed in many other fields, differences in these fields from pollution of Lake Michigan were great enough so that developments in them could not directly be taken and applied to the process of dispersion in a lake.
Surface
r-Y
Bottom
To whom correspondence should be addressed.
Figure 1. Diagram of coordinate system for Lake Michigan Volume 6, Number 4, April 1972 331
First term represents local time rate of change of concentration, the next three terms define changes in concentration due to advection, the next three terms account for dispersive action of the pollutant due to turbulence, and the last term allows for depletion of the pollutant due to chemical reaction with its environment. Although several analytical solutions to Equation 1 are available in the literature (Crank, 1957; Carslaw and Jaeger, 1959), they could not be used owing to the simplifying assumptions made, such as constant dispersion coefficients and velocities, no reaction term, no boundaries, and a twodimensional coordinate system. Another method used to solve Equation 1 is to use a numerical technique such as finite differences (Water Resource Engineers, 1966; 1968a,b). Unfortunately, this method of solution requires tremendous computer time and tends to be unstable for predictions far in the future. It was possible to remove these stringent restrictions and reduce the computation time to a few seconds by developing an analytical solution by particular transformations and the concept of Green’s functions. The model of Lake Michigan can calculate the amount of pollution due to any number of sources on the lake at any time and point in the shore waters of the lake. The dispersion coefficients, velocities, and reaction term can be arbitrary functions. The model also allows for boundaries such as the lake’s bottom and shoreline, thermoclines, and thermobars. An expanded description of the model and its development along with an example problem are presented. Quantitatice Description of Model
Correlations for Eddy Diffusion Coefficients. Eddy diffusion coefficients are characteristic parameters of a given lake and are empirically determined. Although few experimental values have been determined (FWPCA,1967) for the Great Lakes, enough data have been taken on other bodies of water such as rivers (Glover, 1964), oceans (Hill, 1966; Koczy, 1958; Ozmidov, 1958; Foxworthy and Kneeling, 1969) and other lakes (Hutchinson, 1957) so that order-of-magnitude estimates can be made from the following correlations. Longitudinal dispersion means dispersion in the direction of flow. It includes the combined action of turbulence and variation of velocity with position. Glover (1964) gives data expressed as eddy diffusion coefficients for rivers which can be used to account for turbulence due to shear flow. Diffusivity in longitudinal direction due to flow in the x direction is
RH (ft) H (ft) L (ft) N,
= = = =
+
HL/2 ( H L ) (hydraulic radius) depth of lake or thermocline width of lake or thermobar Manning’s constant ~ 0 . 0 2 5
Similar formulas apply for flow in the y and z directions. To account for the variation of these eddy diffusion coefficients with distance due to transition from small to large eddy motions inherent in homogeneous turbulence, the following correlations from the oceanography literature (Hill, 1966) may be used. The longitudinal and transverse coefficients decrease with depth according to: -~
A K iAz r Y
~xlsurrace
Kvlsuriace
AKJAz
O
s 4.6 x
a
10-3 ft-1
Longitudinal and transverse coefficients on the surface increase with distance from the source up to the length scale of large eddy motion (Hill, 1966; Ozmidov, 1958)
where r (ft) source and
=
d x 2
+ y 2 is the absolute distance from the 5 cs 5 2.5
1X
ft2I3/sec
X
For Lake Michigan, this length scale is about 800 ft so that for distances greater than 800 ft, diffusivities do not vary (FWPCA,1967). Combining these relationships, the expressions for the eddy diffusion coefficients as functions of the three space variables are
’
+ csr4I3(1- az) K y = K,o + cSr4l3(1 - az) K,
=
Kz
= CKZHL;*
Kx0
For values of r greater than 800 ft, r should be set equal to 800. Although the model can include flows in more than one direction, these coefficients include only velocity in the x direction, along the shore. The effect of shear is incorporated in K, so that K, should be larger than Ky and K,, since the latter include only the effect of turbulence. To get a feeling of the magnitude of these diffusivities, let us calculate them under the following conditions : z = 20 ft
r > 800 ft H = 27ft uz = 0.072 ft/sec
A similar formula for transverse dispersion is K,o (ft*/sec) = cK~RHu*
The value of the Chezy coefficient and the shear velocity are
K , (ft2/sec) = cKZHrj*
Nc
where cKY,cKZ = empirical constants for a given body of water whose values are approximately 500, 0.72, and 0.067, respectively [see Glover (1964) for further details] u* (ftisec) = l/s u,/N,, (shear velocity) g = gravitational constant uT (ftisec) = inshore velocity in x direction N , = 1.49 (Chezy coefficient)
(Hill, 1966; FWPCA,1967)
-a
where
Kzo (ft2/sec) = cK~RHv*
and for vertical dispersion
=
=
1.49(27)O. 0.025
=
104
cK,,
332 Environmental Science & Technology
so that K,o = 500(27)(3.82 X
=
51.5
=
0.074
K,o
=
0.72(27)(3.82 X
Kz
=
0.067(20)(3.82 X loW3)= 5.1 X
Values of the diffusivities are then
+ + +
51.5 2.5 X 10-S(800)4'a(l- 4.6 X 51.5 16.6 = 68.1 ft2/sec K,, = 0.074 16.6 = 16.7 ft2/sec Kz = 5.1 X ft2/sec Kz
= =
X 20)
One may compare these values with radial diffusivity measured for Lakes Michigan and Erie which is a combination of Kz and K y (Table I). Correlations for K , and Ku are seen to be in agreement with the data. Although there is no equivalent data for K,, its correlation has been confirmed by numerous measurements on rivers (Glover, 1964) and is in agreement with ocean data (Koczy, 1958). Modification and Solution of the General Equation for Lake OF RATEEQUATION. Experimentally, it has Michigan. FORM been found that for several nonconservative materials introduced in the lake such as BOD, R(c) is best represented by a first-order reaction term as R(c) = kc
(2)
where k is a constant having units l/time and is on the order of O.l/day and that the velocities and dispersion coefficients are only functions of the coordinates and not concentration. Using these facts, it is seen that Equation 1 is simplified to a second-order linear partial differential equation for which an analytical solution appears feasible. By rearranging Equation 1, making two changes of variables, and using separation of variables, Equation 1 can be further simplified so that an analytical solution can be obtained by means of the method of Green's functions. REDUCTION TO THE CONSTANT COEFFICIENT EQUATION WITH A MODIFIED SOURCE TERM.Let mean velocities and diffusivities be defined as:
Table I. Experimental Values of the'Radial Diffusivity K, for Lakes Michigan and Erie (FWPCA, Lake Michigan Currents, 1967) Depth, ft z, ft u, ft/sec K,, ftZ/sec 27 20 0.098 31.2 27 5 0.111 59.1 27 20 0.072 39.8 25 5 0,101 44.2 42 20 0.066 34.4 42 20 0.098 35.5
i
= z
- uzt
Then bc bc =bXi b i i
i = 1,2,3
and (7) Also b_c at
-- bc - + at
bc b_i t_= bc i = l a;;., b t --
~
- bC
(8)
vi -
bii
Substituting Equations 5-8 into Equation 4 (9) ANALYSIS OF DIFFUSION TERMS.Make the following change of variables:
xi = " dFi
i = 1,2,3
Then where i = 1,2,3 represents x,y,z. The subscripts 0 and P denote the source and the point where concentration is desired, respectively. u p is the volume element equal to ( x p - XO)( y p
-
yo)
(zp
-zo).
Introducing these quantities into Equation 1 yields the constant coefficient equation
1 __ b2c --b 2 C = Ki bXiz Substituting Equation 10 into Equation 9
SIMPLIFYING TRANSFORMATION IN TERMSOF CONSTANT. Let
where
c =
e-kg;
THE
RATE (12)
Then, substituting Equation 12 into Equation 11 and is to be treated as a source term. The rationale behind this is that the solution for an arbitrary source term can be readily written in terms of the solution for an instantaneous point source which is designated as the Green's function. It should be noted that S is a measure of the deviation of the true solution from that for the mean coefficient case. For this reason, S can be considered as perturbing the mean coefficient solution. ANALYSISOF CONVECTION TERMS.Make the following change of variables. A
x =x
.. y
=
y
- vzt - v,t
(5a) (5b)
This equation is now in the form of the simple diffusion equation, and the methods available for solving it can be used. The solution for an instantaneous point source will be found from which the solution for an arbitrary source term can be generated. SEPARATION OF VARIABLES. Assuming a solution of the form
2
= fi(x,t)fi(y,t)h(z,t>
(14)
and substituting it into Equation 13 without the source term Volume 6,Number 4, April 1972 333
c =
Me-'
'
8 ( ~ t ) ~dI?,&Rz /*
exp [ - ( x - 0,t)2/4R,t -
0, - OUt)'/4Rut - (Z - Ozt)'/4 Rzt] (23) The Green's function is the solution for a unit instantaneous point source released at x',y',z',t' and no decay. Thus, G(x,y,z,t; x',y',z',t') =
exP TRANSIENT SOLUTION FOR Equation 15 is f.
A
{-
X
8[r(t - t ' ) ] 3 / ' d ~ , ~ u R z
[(x - x') - Oz(t - t'N2 4 R , ( t - t')
- 1-0, - y ' ) - O,(t - t')p 4 R,(t - t')
[ ( z - z') -
POINTSOURCE.A solution of
- O,(t - t')
?')I2
4&(t
A - t"2 exp [ - ( X i 2 / 4 tll
(16)
where A is an arbitrary constant. Substituting Equation 16 into Equation 14 c* = -A -e x p [ - ( X 2 / 4 t )
tal 2
- ( Y 2 / 4 t )- ( Z 2 / 4 t ) ]
(17)
Substituting Equation 17 into Equation 12 Ae-'' c = - p i 2 exp [ - ( X 2 / 4 t )
- ( Y 2 / 4t) - (Z2/4t)]
R,Rzx2
(24)
+ R,R,y2 + R,R,zz 4
b = tiz2/4R, (18)
}
TRANSIENT SOLUTION FOR A CONTINUOUS SOURCE.The solution for an arbitrary source is obtained from the solution for the instantaneous source by integrating it over instantaneous sources released every dt. However, it will be convenient first to rearrange Equation 23. Let a=
Transforming back to the 18 becomes
1
R,R,R,
+ Ou2/4I?, + iiz2/4Rz+ k
Then
:,>,; coordinate system, Equation
dg
- [(&i)2/4t]
c = Azexp{-[(")'/4t] taiz
exp
-
(-
- b t ) (25)
Let the pollutant be introduced at rate cp(t)into the lake continuously at (O,O,O) where [ is the time of introduction. Then the effect of a quantity of pollutant p d t on concentration at the point x,y,z at time t i s given by
or c = Ae-' -- exp [-(22/4 R,t) - 6 2 / 4 Rut>- (22/4 R,t)I t3/2
(20)
c=
Transforming back to the x,y,z,t coordinate system, Equation 20 becomes c=-
Ae-" pi2 exp 1-(x - gzt)*/4R,t - ( y - ii,t)2/4 Rut (z
- OZt)2/4R,tl (21)
If the rate then
cp
The constant A depends on the boundary conditions. To obtain the Green's function, the boundary condition is selected as an instantaneous point source of amount M at (O,O,O). Later, provision is made for sources not on the shore. With the reaction term k = 0, the following relation holds
M =
J-mm
J-mm
J-mm
c/k =
o dxdydz
M
. I
,.l
(22)
which states that the amount of material released must be contained somewhere in the system if no reaction occurs. Substituting Equation 21 into Equation 22, integrating, and solving for A (Crank, 1957) A =
is taken as a constant or a suitable mean value,
p
=
Then du
ab =
u3/ad[ and Equation 27 becomes
~_
8 ( i ~ t ) ~dR&k /2
Therefore, the solution for an instantaneous source is 334 Environmental Science & Technology
The integral in Equation 28 has been evaluated (Abramowitz and Stegun, 1965) so that Equation 28 becomes
y*21
y.41
Figure 2. Pictorial representation of the method of images
where
For a continuous source released at t = t oat x = XO,y = yo, z = zg, Equation 29 should be modified by the substitution x=x-x0,
y=y-yo,
z=z-z0,
t=t-ro
STEADY-STATE SOLUTION FOR A CONTINUOUS SOURCE.To obtain the steady-state solution, take the limit of Equation 28 as t + to obtain
and can be satisfied if the concentration profile for an infinite system is considered to be reflected first at z = h, again at z = 0, then at x = I1 and so on, then at y = 12, again at y = 0, then at y = 12, and so on with the result of each successive reflection being superimposed on the original profile. This will produce a double infinite series for the concentration profile in a finite system. An equivalent approach to this problem is to use the method of images in which an infinite system of sources is placed outside of the finite system such that their fluxes cancel at the boundaries. In addition, the velocity field for the images is replaced by the negative of that for the source for those components not parallel to the boundary. To show that this method does indeed give zero flux at the boundaries, let us evaluate the flux at y = 0 for the simple case of an instantaneous point source at x = 0, y = yo,z = 0 and a reflecting boundary at y = 0. Only one image source at y = -yo with a negative velocity in they direction is required. Using Equation 23, the concentration is 2 Me-"
c =
m, [ -:;;
(X
- ozt)z
8 (d)a'2 {exp CORRECTION DUE TO THE NONCONSTANT COEFFICIENTS. If the term en'sin Equation 13 is considered as a source distributed over the space and time domains, then a first-order correction may be applied to the solution for the case of mean coefficients by adding the term AC
= e-'
l'dr'
LmLw J-mw
dY
=0
- O,t)Z
4 Ret
4 R,t
so that the flux is D,c
- dc - KU-
dY
(X
=
2 Me-" 8 (rt)a'2 d
K Z ,
- &t)2 - ( Z - v,t)2
X
2 (Y - Yo -
m
x
where c is the solution for the case of mean coefficients. In general, this integral must be evaluated numerically. However, since the velocities and diffusivities are not strong functions of position and the solutions are not extremely sensitive with respect to them, the corrections would not be expected to be a great amount. FINITE BOUNDARIES. To account for finite boundaries (such as the surface and bottom of the lake, a thermocline, the shorelines, or a thermobar), the procedure is as follows: If the pollutant is of the type that deposits on the bottom such as suspended solids, so that the boundaries are absorbing ones, the previous solutions may be used. However, if it is reflected at the boundaries as dissolved solids are, then the condition that there is no flow of pollutant across these boundaries must be satisfied. These conditions for rectangular boundaries are expressible as
- oc - K,-
(Z
-(Y
en"S(c,x ',t ')G(x,t ;x',t ')dx 'dy 'dz '
D,C
--I.
a t y = Oandy = 12
4 R,t
4 Rut
Substituting y = 0
('" - K$)l,:o
=
2 Me-" 8 ( ~ t ) ~ dl?&,R, / 2 X
(Yo
+ v,t)*
Before applying these methods to the problem at hand, let us further demonstrate the method for the problem of a source at y = 0 in a one-dimensional medium with boundaries at y = 0 and y = 1. Iff(y,t) is the solution for the semiinfinite case, then it can be seen from Figure 2 that the solution for the finite case is Volume 6, Number 4, April 1972 335
lmaqe 5
Image 2
flux due to image 2 across the bottom, etc. As shown in Carslaw and Jaeger, 1959, if 0 = r / m where m is a positive integer, the 2 mth image will coincide with the source so that a closed system is obtained. MULTIPLE SOURCES. To treat multiple sources, one calculates the contribution for each source as given by Equation 32 and adds them to obtain the combined effect. This is justified because the general partial differential equation is linear. Example Problem
Image 6
\
\
Image I
Bottom
It is desired to know the individual and combined effects of three outfalls situated on the shoreline of Lake Michigan on the concentration of BOD along the shoreline and 50 ft out from shore. Locations, concentrations of BOD, and flow rates of the three outfalls are:
Figure 3. Location of images for a sloping bottom
5
f(z
n=m
- 2 no
where the velocity field for the images is the negative of that for the source, and the summation starts at n = 0 and progresses 1 , -1, 2, -2 . . . . +n, -n. Thus, the solution for the finite case satisfying the boundary conditions expressed by Equation 31 for a source located at (xo,y o , z o )and released at to may be written down as C(X,Y,Z,t) =
5 5
x
m = - - m n=--m
+ (2 mIz - z + (2 nI1 - zo), t - to] + y + (2 mlz - Y O ) ,z + (2 nh + t - to1 + Y + (2 mlz + YO),z + (2 nh t - to1 + Y + (2 mlz + Y O ) ,z + (2 nh + t - to]} (32)
{ c [ x - xo, Y
yo),
- XO, c[x - X O , c[x - X O , c[x
Outfall 1: x = 0,20 ppm, 10,000 gpm Outfall 2: x = 2000 ft, 50 ppm, 50,000 gpm Outfall 3: x = 5000 ft, 30 ppm, 200,000 gpm
ZO), ZO),
ZO),
where the negative of is to be used for images in the y direction and the negative of or for images in the z direction. C can be either Equation 23 for the instantaneous source case, Equation 29 for the transient continuous source case, or Equation 30 for the steady-state continuous source case. Another example of the use of image sources is the case of a sloping bottom as shown in Figure 3. The reasoning behind the system of images shown is as follows: Image 1 cancels the flux due to the source across the bottom; image 2 cancels the flux due to image 1 across the surface; image 3 cancels the
During the time of year to be considered, velocity along the shore of the lake is 1 ft/sec, mean depth is 30 ft, and there is a thermobar 500 ft offshore. The rate constant for the decomposition of BOD is 0.1 day. Information to be obtained are: (1) the steady-state contribution of each outfall to the BOD and their total for y = 0 and y = 50 ft from x = -5000 ft to x = 20,000 at 1000 ft increments, and (2) concentration history due to outfalls at the shoreline 5000 ft downstream for outfall 3. Solution: The requested information is given in graphical form in Figures 4-6. They were obtained for the conditions: (1) reflecting boundaries at y = 0 and 500 ft and at z = 0 and z = 30 ft; (2) uz = 1 ft/sec, uy = ur = 0; and (3) mean diffusivities using the same parameters as in the sample calculation. Nomenclature = concentration
C
c ~ ~ , c K= ~ ,empirical c ~ ~ = empirical C.
constants, dimensionless constant, ft 2' */set g = gravitational constant H = depth, ft k = first-order rate constant, sec-' K,,K,,K, = turbulent dispersion coefficients for the x,y,z directions, respectively, ft */set Kzo,KVo = turbulent dispersion coefficients due to shear for the x , y directions, respectively, ft*/sec
VI
Source 2
!/--
Source 3
-
Figure 4. Concentration of BOD along shoreline for y = z = 0 at steady state
0
f 6 g 5 $ 4
3 2 I
O -2000
-
2000
4000
6000 8000 10,oOO 12pOO 14,000 16,000 l8,OOO 2
336 Environmental Science & Technology
00
6
5
-5
4
Figure 5. Concentration of BOD along shoreline for y ft and z = 0 at steady state
s
‘;3
e
=
50
0
8 2
I
0
Figure 6. Transient mcentration of BOD at shoreline due to source 3 SO00 ft downstream from source 3
2
4
6
8
IO
12
14
16
18
20
Time ( h r )
width, ft normal to S = Chezy coefficient = Manning’s constant = ft = rate of reaction of pollutant with its environment = hydraulic radius, ft = bounding surface = time, sec = shear velocity, ft/sec = velocities for the x,y,z directions, respectively, ft/sec = coordinate measured in direction of shore, ft = coordinate measured in direction across lake, ft = coordinate measured in direction vertically downward, ft = empirical constant, ft-1 =
=
-\/2+.,
Literature Cited Abramowitz, M., Stegun, I. A,, “Handbook of Mathematical Functions,’’ pp 304, Dover, New York, NY (1965). Carslaw, H. S., Jaeger, J. C., “Conduction of Heat in Solids,” Oxford Univ. Press, New York, NY (1.959). Crank, J., “The Mathematics of Diffusion,” Oxford Univ. Press, New York, NY (1957). Foxworthy, J. E., Kneeling, H. R., “Eddy Diffusion and
Bacterial Reduction in Waste Fields in the Ocean,” Research Grant WP00931 FWPCA, conducted through Allan Hancock Foundation by University of Southern California, Report 69-1, September 1969. FWPCA, Lake Michigan Currents, Chapt. 8, November 1967. Glover, R . E., “Dispersion of Dissolved or Suspended Materials in Flowing Streams.” Geol. Survey Professional Paper, 433-B (1964)Hill, M. N., “The Sea,” Vol. 1, Section VI, Interscience, New York. NY. 1966. Hutchinson, ’G. E., “A Treatise on Limnology,” Vol. I, Wiley, New York, NY, 1957. Koczy, F. F., “Natural Radium as a Tracer in the Ocean,” Proc. UN Int. Conf. Peaceful Uses A t . Energy 2nd, 18, 351 (1958). Ozmidov, R. V., “Computation of Horizontal Eddy Diffusivity of Dye Tracers in the Sea,” Dokl. Akad. Nauk. SSSR, 120,761 (1958). Water Resources Engineers, Inc., “Applications of Mathematical Models for Prediction of the Thermal and Quality Behavior of Lake Washington, A Report Prepared for the Washington Pollution Control Commission,’’ Olympia, Wash., November 1968a. Water Resources Engineers, Inc., “A Hydraulic Water Quality Model of Suisun and San Pablo Bays, A Report of an Investigation Conducted for FWPCA,” March 1966. Water Resources Engineers, Inc., “Prediction of Thermal Energy Distribution in Streams and Reservoirs,” prepared for the Dept. of Fish and Game, California, August 1968b. Receivedfor reoiew March 3,1971. AcceptedSeptember 3,1971.
Volume 6 , Number 4, April 1972 337
