A stream pollution model with intervals for the rate coefficients

A Stream Pollution Model GARY W. HARRISON* Department of Mathematics, Receiwd

with Intervals for the Rate Coefficients

Unioersiry of Georgia, Athens, Georgia 30602

20 April 1979; revised I5 October 1979

ABSTRACT A major obstacle in using mathematical models to make management decisions is the difficulty in obtaining precise values for the parameters. A classical stream pollution model for biochemical oxygen demand and dissolved oxygen is analyzed under the assumptions that the rate coefficients may vary along the stream but remain within fixed known intervals and the initial conditions are only known to lie in certain intervals. The upper and lower envelope for the set of solutions of all systems satisfying these conditions are obtained compartmental

1.

using differential models.

inequality

techniques

that were

recently

developed

for

INTRODUCTION

The classical Streeter-Phelps [l] model uses two quantities to measure the state of a polluted stream: the biochemical oxygen demand (BOD) and the dissolved oxygen (DO). Dobbins [2] argues that for most flowing streams the effect of diffusion on the steady state is negligible, in which case the steady state equations become

(1)

where x = distance downstream, C(x) = concentration of dissolved oxygen at distance x, L(x) = biochemical oxygen demand at distance x, u = average velocity of water flow (miles/day), and C, =saturation concentration for dissolved oxygen. The equations indicate that BOD and DO react to form inert compounds at a rate k,L. In addition the BOD decreases due to sedimentation or adsorption at a rate k,L and is added to at a rate I, by

*This work was partially MATHEMATICAL QElsevier

North

supported

BIOSCIENCES Holland,

Inc., 1980

by NSF Grant

No. DEB77-02942.

111

49~11 l-120(1980) 0025-5564/80/030111+

10$02.25

112

GARY

W. HARRISON

pollution sources along the stream. The DO increases through reaeration at a rate k2[ C, - C(x)] but is removed at a rate d, due to respiration by organisms in the stream. (If plant photosynthesis exceeds respiration, dB will be negative.) The use of steady state equations to describe the BOD and DO profiles along the stream implies that conditions change slowly enough over time compared to the reaction rates that the concentrations remain close to their current steady state. Since they have been used by several other authors [e.g. 2, 3, 81, I will accept this assumption. Making the change of variables t = x/u (distance downstream in units of the distance flowed in one day) and D= C, - C (oxygen deficit), the equations (1) can be put in the form

i = -k,][ h]+[ 2) [I[ -kk,k3 d

(2)

. =d/dt. These equations can be solved to show that although the where BOD level L decreases exponentially, the oxygen deficit D initially rises until most of the BOD has been removed and then falls again. The critical requirement for survival of fish and other aquatic life is that D should not rise too high (i.e., the dissolved oxygen does not become too low). A major difficulty in using this model for decision making is that the parameters in Eq. (2) are not completely known and are in fact intrinsically random. They vary over time and over distance along the stream. Even though I have assumed that the variation over time is slow enough for the steady state equations to be a valid model, variations over distance along the stream will affect the steady state. Esen and Rathbum [3] found that the average coefficient of variation for k, and k, is about 30-35%. A number of authors [2-71 have used stochastic models to find probability distributions for the BOD and D at points downstream. Most recently Padgett [8] treated (2) as a random differential equation with random constant coefficients k, and k, and random initial condition Lo and D,. Here I take a different approach, assuming that the coefficients in (2) may be piecewise continuous functions of the distance along the stream but that they satisfy known bounds

STREAM

MODEL WITH INTERVAL

and that the initial conditions

COEFFICIENTS

113

satisfy

(4) The goal will be to construct the “envelope” of the trajectories of (2) under the restrictions (3) and (4) that is, the smallest subset E(t) c R 2 such that if (L(t), D(t)) is the solution of any differential equation of the form (2) with coefficients satisfying (3) and initial conditions satisfying (4), then (QQ,D(t)) 2.

E E(t).

COMPUTATION

OF THE ENVELOPE

OF THE TRAJECTORIES

Bounds for L(t) are easy to obtain, since i depends only on L and not on D. By the comparison theorem for one dimensional systems [9], c(t) Q L(t) < L< t) if h(O) = he, L(O) = L,, and & and L satisfy

But the equation for d contains the term k,L, and although k, = k; maximizes b for any fixed L, it also results in the fastest decrease in L. Is the maximum value of D obtained by using up L quickly to achieve the maximum rate of increase b,, or by conserving L to obtain a higher fi later on? Equation (2) has the form of a linear compartmental flow model, and the following result is proved for general linear compartmental flow models in Harrison [lo]: THEOREM Let p(t)=[pl(t),p2(t)] [A

be the solution of

Pzl

A]=P1

k?(t)+

k;(t)

-k:(t)

0

1

K(t) ’

where

-1

if

Pz(t)-P,(t)

if

P*(t)-Pl(t)
if

Mt)

if

P*(t) CO,

> 0,

> 0,

P)

(7b)

114

GARY

W. HARRISON

and

k:(t)=

; 3

i

if

p,(t) 2 0,

if

PI(t) <(A

(7c)

and let

if

p,(t) > 0, P,(f) co,

if

Zf the solution [L*(t),D*(t)]

Ii*I[=

Pz( t) > 0, P2(t) co.

P)

of

-k;(t)-k:(t)

d*

if if

(74

A&)][ K]+[ ::::I

k:(t)

(*)

is nonnegatioe and (L(t), D(t)) is any nonnegative solution of (2) subject to the restrictions (3) with PI(O)L(O)

+P2(o)~(o)~Pl(o)~*(o)+P2(o)~*(o)~

(9)

then for aN t > 0,

fJ,(t)qt)+P,(t)qt) Proof:

~Pl(t)L*(t)+P2(t)o*(r).

(‘0)

It is easy to show that for all

$Ipl~+~2~l~~[P,L*+p,D*l

t > 0.

APPLICATION

The theorem places no restrictions on how the initial values p(O) must be chosen, but once they are chosen the inequality (9) can be satisfied for all solutions [L(t), D(t)] of (2) that satisfy the restriction (4) by choosing [L*(O), D*(O)] according to

L*(o)=

1fo ff -0

P,(O) 2 07 If P,(O)
D*(o)=

i

(11)

2

ff

-0

If p2(0)<0.

Pz(O)~Q

STREAM MODEL WITH INTERVAL COEFFICIENTS

01

115

6

\2-4

L= Biochemical

Oxygen

Demand

FIG. 1. A collection of bounding lines of the form given by the inequality (10) gives an estimate of the envelope in L-D space at r = 3 of the trajectories of Eq. (2) under the restrictions given in (13) for Case 1. Each bounding line is found by using a different initial condition for [pi& in Eq. (6), and the heavy dots indicate the location of the corresponding solution (L*,D*) of (8). Since the exact envelope must contain all of these points, the polygon inside the bounding lines is seen to be a good approximation of the exact envelope.

The inequality (10) then indicates that all solutions of (2) subject to the restrictions (3) and (4) lie on one side of the line p,L +p2D = b in L-D space, with b=p, L* +p2D*. Furthermore, there is at least one solution of (2) satisfying (3) and (4) which actually lies on this bounding line, namely [L*(t), D*(t)]. The line shifts, however, as t changes, because p,, p2, L*, and D* change with t. By carefully choosing a number of different initial conditions p(O), one can obtain a family of bounding lines that enclose a polygon in L-D space. Figure 1 shows that if enough bounding lines are used the resulting polygon is a good outer approximation of the envelope E(t) of the trajectories. Over any interval t, < t < t, where p,(t), p*(t), and p2(t)-p,(t) do not change signs, Eq. (6) can be solved to give p,(t) =

[

PI(h)-p2(t0)$$(ebr-l)

ew + vu I

, (12)

pz(t) =p2( t&+*(‘),

where d* = kf - kr - kf and the k: are determined

by (7). Under the initial

116

GARY

W. HARRISON

condition [p,(O),p,(O)]=[l,O], (12) yields ~i(t)=e(%l+!d’, pz(t)=O. The line p,(t)L + P,(t)D = b remains vertical, and the inequality (10) becomes L(r) < L*(t) with L*(t)=z(t) from (5). Initial conditions [ - l,O] yield L(t) a&(t) with h(t) given by (5). Thus two of the bounding lines are what we already knew they should be. The most important problem, however, is to obtain a good upper bound for D(t), which requires the bounding line to be approximately horizontal, i.e. ~,(t)<~~(t). But experience shows that it is usually impossible to find a bounding line which remains horizontal, because all solutions of (6) with p*(O)#O eventually tend toward the same direction [determined by the eigenvector corresponding to the largest eigenvalue of the matrix in (6) with appropriate values of ki*], and hence all the bounding lines tend toward the same slope, which unfortunately is far from horizontal. To overcome this difficulty, pick a point t, and let [p,(ti),pz(ti)]=[O, 11, solve (6) backward over t [or use the equations (12)] to obtain [p,(0),p2(O)], choose L*(O) and D*(O) according to (1 l), and then solve (6) and (8) forward over t to obtain a bounding line (10) which is close to horizontal for points close to ti. Repeating this procedure for a sequence of points ti gives a family of bounding lines, one of which is horizontal at each ti. Lower bounds for D(r) can be obtained by the same procedure, but starting with [pl(ti),pz(Q]=

10,- 11. 3.

EXAMPLE

Padgett [8] assumed L(O), D(O), k,, and k2 to be normally distributed with means 6.8, 0.3, 0.35, and 0.75, respectively, based on the data of Thayer and Krutchkoff [5] for the Sacramento River, and with variances 1.85, 0.01, 0.015, and 0.0506, respectively, based on the coefficients of variation found by Esen and Rathburn [3]. Case I.

For purposes

of comparison

the intervals

4.13 < L(0) < 9.46, 0.10 =GD(0) Q 0.50, 0.11 <

k,

< 0.59,

0.31 <

k,

Q 1.19

(13)

are used, which are 95% confidence intervals for each variable if they are independently normally distributed with the means and variances used by Padgett. Intervals for the other parameters are used that are consistent with

STREAM MODEL WITH INTERVAL COEFFICIENTS

the data of Thayer and Krutchkoff 0.0~

[5] for the Sacramento

117 River:

k,,
O.O
The other parameters

0.15 Q k,

~0.55,

0.38 < k,

< 1.12.

are taken to be k, = d, = 1, = 0, u = 5, and c, = 9.0.

RESULTS The results in both cases are similar and will be considered together. If p,(O)> 0 and p2(0)= 1, then from (12) pz(t) remains positive but p,(t) eventually changes from positive to negative, so that the resulting bounding line shifts from a negative slope to a large positive slope. To obtain an upper bound for D at time ti, one solves (6) or (12) and (7) backwards over t from the initial condition [p,(fi),p2(ti)]= [0, l] and finds that p,(t) and p*(t) are positive for all times previous to ti, hence from (7b, c) k,*(t) = &, and kf(t)=k, for all tti-2.298 (t >ti-2.175 in Case 2) and p,(t)>p2(t) for previous points; hence kf(t)=k, for t >ti2.298 (t > t, - 2.175 in Case 2) and kf(t) = 6, before that. In other words, the worst possible (i.e. the maximum) oxygen deficit at a given point ti happens when reaeration and sedimentation take place at the minimum rate along the entire stretch, whereas the BOD-DO reaction initially occurs at the minimum rate, but shifts to the maximum rate for the last 2.298~5 miles (2.175 X 5 miles in Case 2). For distances t > ti, however, this trajectory no longer gives the maximum D(t); a trajectory in which the BOD was “conserved” a little longer by keeping k, = &, longer before switching to k, = El gives a higher value of D(t). Minimum values of D(t) can be found in a similar manner.

118

GARY

t (5t =Miles

W. HARRISON

Downstream)

FIG. 2. The outer solid lines are bounds for the oxygen deficit D in Case I and the inner solid lines are bounds for D in Case 2 of Sec. 3. The dashed lines indicate some trajectories that actually reach the upper bound in Case 2. Use of kf = 6, produces a trajectory that follows the upper bound up until f =2.175. Trajectories that touch the upper bound at points farther along the stream are obtained by switching from k: = b, to k,L = I?, at the points indicated by the heavy dots, 2.175 units before the upper bound is reached.

FIG. 3. Each polygon approximating the envelope in L-D space of the trajectories of (2) at the indicated value of I is found as in Fig. 1. The outer polygons are for Case I and the inner polygons

for Case 2 of Sec. 3.

STREAM

MODEL WITH INTERVAL

COEFFICIENTS

119

Figure 2 shows the upper and lower envelope for o(t) for both Case 1 (outside lines) and Case 2 (inner lines), plus some trajectories that touch the upper envelope in Case 2. Figure 3 shows how the estimate of the envelope of the solutions in L-D space changes over time. Each of the polygonal regions was obtained from a family of bounding lines as in Fig. 1. 4.

COMPARISON

WITH

STOCHASTIC

MODELS

Both the interval analysis techniques used in this paper and the stochastic differential equation techniques used, for example, by Padgett [8] are attempts to deal with uncertainty in our models of natural phenomena. Under the assumptions used by Padgett there is a 66% probability that L(O), D(O), kl, and k2 will simultaneously satisfy the intervals (14) given in Case 2. One might be tempted to compare the bounds for D(t) found in Case 2 (inner solid lines of Fig. 2) with 66% confidence intervals that could be computed from the probability density functions for D(t) given in Padgett’s Fig. 5. Clearly the bounds in Fig. 2 are way outside any such confidence intervals computed from Padgett’s results. But there are several reasons why such a direct comparison is not appropriate. First, Padgett uses random but constant coefficients, whereas I have allowed coefficients that are piecewise continuous functions of distance along the stream. Second, even though at a given point in time a confidence interval might contain 66% of the trajectories of systems with parameters chosen at random from the distributions used by Padgett, it is possible that most of these trajectories leave this interval for some values of t. On the other hand, the trajectories always stay within the bound given in Figs. 2 and 3 as long as the parameter functions remain in the intervals given in (14). The same reasoning shows that confidence intervals for the parameters from a distribution that holds at a fixed point along the stream need not have the same probability of containing the parameter functions along the entire stream. One would need to know something about their transition probabilities. Finally, it must be recognized that stochastic models and interval methods try to answer two different questions. Stochastic methods aim at determining what is most likely to happen, and interval methods aim at finding the most extreme (the worst or the best) that could happen. Although the trajectories of a few systems that satisfy the conditions (13) or (14) actually touch the outer envelope shown in Figs. 2 and 3 at some point in time, most of them remain closer to the mean computed for the stochastic differential equation. Bounds computed from interval techniques offer a very conservative approach to the decision maker’s problem. When intervals for the parameters can be obtained directly from the data, the interval method avoids the difficulty of hypothesizing the proper

120

GARY W. HARRISON

type of distribution function for the parameters. It also has the advantage of being computationally more straightforward, thus making it easier to apply to large complex systems. It promises to be a useful addition to the tools for dealing with uncertainty in our models, although there is still a great deal of work to be done in removing the restriction to systems with positive trajectories and extending the techniques to nonlinear systems. REFERENCES H. W. Sweeter and E. B. Phelps, A St& of the Pollution and Naiural Purijcation of Public Health Bull. 146, U.S. Public Health Service, Washington, 1925. 2 W. E. Dobbins, BOD and oxygen relationships in streams, J. Sanitary Engineering Diuision, ASCE (SA3) 90:53-78 (1964). 3 I. I. Esen and R. E. Rathbum, A stochastic model for predicting the probability distribution of the dissolved oxygen deficit in streams, Geological Survey Professional Paper 913, U.S. Government Printing Office, Washington, D.C., 1976. D. P. Loucks and W. R. Lynn, Probabilistic models for predicting stream quality, Water Resources Research 2:593-605 (1966). R. P. Thayer and R. G. Krutchkoff, A stochastic model for BOD and DO in streams, the Ohio River,

J. Sanitary Engineering Division, ASCE

(SA3) 93:59-72

(1967).

V. Kothandaraman,

Probabilistic variations in ultimate first stage BOD, J. Sanitmy Engineering Division, ASCE (SAl) %:27-34 (1970). W. J. Padgett, G. Schultz, and C. P. Tsokos, A random differential equation approach to the probability distribution of BOD and DO in streams, SIAM J. Appl. Math. 32~467-483

(1977).

8 W. J. Padgett, A stream pollution model with random deoxygenation and reaeration coefficients, Math. Biosci. 42: 137-148 (1978). W. Walter, Differential and Integral Inequalities, Springer, New York, 1970. 10 G. W. Harrison, Compartmental models with uncertain flow rates, Math. Biosci. 9

43:131-139

(1979).