STEADY STATE MODELING PROGRAM Application Manual R. P. Canale and S. Nachiappan Department of Civil Engineering College of Engineering The University of Michigan March 1972 Sea Grant Technical Report No. 27 MICHU-SG-72-207 THE UNIVERSITY OF MICHIGAN SEA GRANT PROGRAM

The University of Michigan Sea Grant Program is a part of the National Sea Grant Program, which is maintained by the National Oceanic and Atmospheric Administration of the U. S. Department of Commerce.

TABLE OF CONTENTS INTRODUCTION 1 THEORY 2 PROGRAM APPLICATION 10 EXAMPLE 19 ACKNOWLEDGMENTS 27 REFERENCES 28

FIGURES 1. Uniform Rectangular Volume 5 2. Traverse Bay Region 21 3. Model Segmentation and Assumed Circulation in Grand Traverse Bay 22 4. Total Coliform Distribution in Grand Traverse Bay in 1964 23 5. Distribution of Total Coliform Following 90% Treatment of Cherry Wastes 24 6. Distribution of Total Coliform with 90% Treatment of Cherry Waste and New Source at Segment Number 15 25 7. Total Coliforn Distribution With No Treatment With All Discharges to Segement 35 26

INTRODUCTION A general methodology for modeling biological production in aquatic systems has been described in an earlier report by Canale (1970). This work emphasized the transient behavior of relatively complex ecosystems in single homogeneous zones. The ultimate goal of the project was to suggest how the annual cycle of phytoplankton and nutrient behavior might be simulated mathematically. The report described the complex kinetics of the system as characterized by nonlinear reaction terms and time-variable rate coefficients. A relatively simpler class of problems is also of interest. These problems are concerned with the steady state distribution of species whose decay or production tendency can be described by first order kinetic formulations. The following list of water quality variables has been traditionally analyzed using such assumptions: 1) Total Dissolved Solids or Conductivity 2) Chlorides 3) Any Conservative Chemical Species 4) Total Coliform Bacteria 5) Dissolved Oxygen 6) Biological Oxygen Demand 7) Radioactive Isotopes 8) Nitrogen (Nitrification Process) This report describes a general user-oriented program capable of calculating the three-dimensional s stadytate distribution of the above water quality variables in aquatic systems.

-2 -THEORY Mathematical models which can be useful for informed management of water resources must be based on the diverse chemical, physical, and biological mechanisms active in the system. These mechanisms are recognized by appropriate terms in equations of continuity for each chemical or biological element of interest. Essentially, two distinct types of mechanisms are recognized. First, are those processes which alter the concentration of material within a closed system due to departures of the state from a chemical equilibrium state. The study of the rate of changes toward or away from this equilibrium state is called kinetics. The kinetic expressions may be dependent on species concentration, temperature, light intensity, and pH. A second process which can bring about changes in the concentration of species results from the mechanical action of the fluid circulation and the subsequent dilution of concentration gradients. The bulk behavior of the circulation is characterized by gross convective transfer, while the random small-scale fluid movement is accounted for by dispersion coefficients and transfer due to concentration gradients alone. The above ideas are summarized by Equation 1 which is a descriptive statement of the continuity law for any material. Equation 2 expresses this same law in mathematical form for a general three-dimensional system.

Rate of Change of i Rate of Input of i Rate of Output of Rate of Production Rate of Disapwith time within a _ by convection, dis- i by convection, of i by growth, pearance of i by cell persion, sedimenta- dispersion, sedi- excretion, or uptake, predation, (weight/time) tion, or migration mentation, or + dissolution within respiration, death, from adjacent cells migration from cell or precipitation (weight/time) adjacent cells (weight/time) (weight/time) (weight/time) (1) Composition Change Hydrodynamic Mechanism Reaction Mechanism __ - V(EVC)- v-(UC) + S (2) ACCUMULATION DISPERSION BULK FLOW REACTION

-4 -A direct solution of Equation 2 for natural systems is not possible. Therefore in practice it is necessary to use approximations which are equivalent to considering a continuous body of water as a series of finite interconnected segments as shown in Figure 1. In this case, the steady-state continuity equation with first order kinetics reduces to the following: dCk Vk dt = [-qkj (k Ck+kj C )+Ekj (-Ck) -VkKkCk+Wk dt kj kj kj j j 3 k jk (3) where: C = concentration of water quality variable in segment k, (mg/l) Vk = volume of segment k, (cft) Q -= net flow from segment k to segment j (positive outward) (cubic feet per second) akj = finite difference weight given by ratio of flow to dispersion, 0<a<l, dimensionless kj = kj Ekj = dispersion mixing coefficient between segments k and K j (cft/day) = E A tj/L kj kj kj Kk= first order reaction coefficient in segment k for water quality variable C (day-l) E = dispersion coefficient between segments k and j (square miles/day) Ak. = cross-sectional area between segments k and j (square feet) Lk = average of characteristic lengths of segment k and k j (feet) Wk = source (or sink) of variable C in segment k (pounds/ (day) C. = concentration of water quality variable in adjacent segments (mg/l) The above units for the parameters are those commonly used in practice and may be changed according to user preference.

Iv -( C k 3 F Ad Ck2 - C 9 ~~~~ ~~~~~k4 V~ ~k Ck Figure 1. Uniform rectangular volume, Ck6 ~~~~~~~~Ck ~~~~~

-6 -If all terms involving the variable C are grouped on the left hand sides, Equation 4 is obtained: akkCk+Eakj C Wk kk k kjj k( j 0 J (4) Where: akk = (Qkja kj+Ek +V k1k akj kj kj kj If n segments are used, then a series of n equations can be written for the n unknowns, C, C2, C3 ---C. Boundary conditions apply at the interface of the n segments and outer sections of the waterbody. The flow between the boundary and the section is designated Qkk. When Qkk is positive (leaving the section) akk= (Qkjakj+Ek )+kKk+kQkkkkkk (5) and Wk= Wk+(Ekk-Qkkkk)CB: (6) where C is the boundary concentration of the variable C and must be B known. If Qkk is negative then, akk = (kjakjj +Ek +VkK+Qkk kk+Ekk (7) and Wk = Wk+(Ekk-QkkOkk)CG' (8)

-7 -The n equations would then be given by: allC1 + a12C2 + I... alnC = W1 a21C1 + a22C2 + 2 + a2nC - W2 a C + a C +...... + a C = W nl 1 n2 2 nn n n where the boundary conditions are incorporated into Wk's. There are a number of numerical procedures for solving such a set of n simultaneous equations in n unknowns. The SSIP uses an MTS subroutine called SLE1 which uses the Gaussian elimination technique. Equation 3 is suitable for single dependent variables that are not forced by outputs from other quality systems. Examples of such a type are chloride, coliform, and BOD. Other water quality variables such as dissolved oxygen are coupled to other systems. The utilization of oxygen depends on the distribution of BOD, and therefore it is necessary to first obtain the distribution of BOD and then use these results in a continuity equation for D.O. Thus, the mass balance equation for D.O. is: dC dk 0 S Qk (a kjCk+kjCj +Ekj(Cj-Ck) ]VkKak(sk- k k dkLkWk (10) dt j where Ck is the saturation value of D.O., Kk is the reaeration coeffisk ak cient in segment k, Kdk is the deoxygenation coefficient, L is the biochemical oxygen demand and +W is now interpreted as sources and sinks of D.O. such as benthal demands and photosynthetic production or respiration. With Lk known from previous calculations, the final solution of Equation 10 is similar to solution of Equation 3.

-8 -As spatial approximations to derivatives have been used in Equations 3 and 10, some errors are introduced into the analysis. One of the errors is "psuedo or numerical dispersion." It appears due to the assumption of completely mixed finite volumes. Numerical dispersion is defined by Equation 11. numkj kj kj (k (11) When akj = 1/2, E is zero. kJ 'numk. On the other hand, it can be shown that for a positive solution the terms off the main diagonal in the left-hand side of Equation 9 should be non-positive. This condition is satisfied if l kj>l kj Akj / / kj I kj| (12) Writing Equation 12 in another way, it is seen that Lkj must be chosen such that, <E kjAkj 1 kj JQkjj (1lkj) (13) For the case of zero numerical dispersion kj = 1/2 and 2E.A - < k... kj- IQkj (14) If akj is set equal to 1/2 and Lkj is chosen in such a way as to satisfy Equation 14, then it may be necessary to handle many segments.

-9 -However, if akj differs very much from 1/2, then numerical dispersion would be high. Further, making akj = 1/2 does not imply the best solution for the case of unequal-sized segments. In SSMP akj is first set equal to L. kj LkLj (15) In other words, the segment whose center is nearer to the interface would have more weightage in determining the concentration at the interface in Equation 3. The value of ckj is then checked against Equation 12. If it is not satisfied, then Ekj ak. is made equal to 1 - QT Choosing proper spatial grids for approximations to the differential equations is still very much an art. The more numerous the segments in a model, the more accurate the resulting solution. However, in such cases the computer costs may be very high, so a compromise is necessary. Considerations of computer size, nature of problems, degree of accuracy, simplicity of the resulting finite difference equations, and availability of verifying field data all influence the choice. For additional details concerning these questions the reader is referred to Thomann (1971).

-10 -PROGRAM APPLICATION The subroutine SSMP (Steady State Modeling Program) is on permanent file under the computer center user number SBEH. The following MTS cards are required to access and execute SSMP. column 1 on card $SIGNONACCNO (your user number) PASSWORD (your password) $RUNASBEH:SSMP (SSMP DATA) $ENDFILE $SIGNOFF Note: A refers to a blank space. No other blanks should be included between characters.

Data Columns Variable Mode Description Card on Card Name No. ~1 1-80 Title A A maximum of 80 columns is used to describe the model (Literal) or run. 2 This card gives the user options with print out information. All the variables on this card have a value of zero (0) or one (1). 5 KEY(l) I 0 - A value of zero for this variable gives (Integer) a print out of the number of segments in the ( simulation run and the scale factors KEY(1) employed. 1 - A value of one for this variable suppresses the above output. 10 KEY(2) I 0 - A value of zero for this variable gives a print out of the values of the interfacial areas, dispersion factors, flows, and the KEY(\) boundary segments for each segment in the jhLL {,/) \. simulation run. 1 - A value of one for this variable suppresses the above print out. 15 KEY(3) I 0 - A value of zero for this variable gives a print out of the values for the character( istic lengths and volumes for all the segKEY(3) ments in the simulation run. 1 - A value of one for this variable suppresses the above output. 20 KEY(4) I 0 - A value of zero for this variable gives a print out of the boundary conditions used KEY(4)in the simulation run. KEY (4) 1 - A value of one for this variable suppresses the above output.

Data Columns Variable Mode Description Card on Card Name No. 2 25 KEY(S) I 0 - A value of zero for this variable suppresses the following output. KEY(5) 1 - A value of one for this variable gives a print out of the values of alpha internally computed for all the interfaces in the simulation run. 35 KEY(7) I 0 - A value of zero for this variable gives a print out of the values of the load inputs KEY(7) used in the run. 1 - A value of one for this variable suppresses the above output. 40 KEY(8) I 0 - A value of zero for this variable gives a print out of the values of the reaction coefficients used in all segments. If an oxygen vKEY (8) deficit run is required the values of ) reaeration coefficients are also displayed. 1 - A value of one for this variable suppresses the above output. 45 KEY(9) I 0 - A value of zero for this variable suppresses the following print out. K 9 ^ 1 - A value of one for this variable gives a KEY9)print out of the part of the A matrix requested by next card. In a deficit run, the same part of the second matrix is also displayed.

Data Columns Variable Mode Description Card on Card Name No. 2 50 KEY(10) I 0 - A value of zero for this variable suppresses the following print out. KEY(10-) 1 - A value of one for this variable gives a I print out of the right hand side vector, W in [A]C = W. If a deficit run is required, a second right hand side vector is also displayed. 80 KEY(16) I 0 - A value of zero for this variable means Y( no deficit run is requested. KEY(16) 1 - A value of one for this variable means deficit run is requested. 3 If KEY(9) is zero skip this card. If KEY(9) is one the user has the option of choosing some part of the A matrix for display. 1-5 IMIN* I IMIN = row of A matrix after which output is requested. 6-10 IMAX* I IMAX = row of A matrix before which output is requested. Example: If A is a 100 by 100 matrix then: (a) A value of IMIN = 1, and IMAX = 10, displays the first 10 rows of A matrix. (b) A value of IMIN = 1 and IMAX = 100 prints out the entire matrix.

Data Columns Variable Mode Description Card on Card Name!No. 4 1-3 N* I N = number of segments used in the model (maximum of 100). 6-15 Area Scale G Scale factor which multiplies all areas. (Real E or F type) 16-25 E Scale G Scale factor which multiplies all dispersion coefficients. 26-35 Q Scale G Scale factor which multiplies all flows. 36-45 Length G Scale factor which multiplies all characteristic lengths Scale of the segments. 46-55 Volume G Scale factor which multiplies all volumes. Scale 56-65 BC Scale G Scale factor which multiplies all boundary concentrations. 5 1-5 Area F Each segment has a maximum of 6 interfaces. Since there (F-type are 3 interfaces per card, 2 cards per section must be real) specified even if less than 4 interfaces are used (the second card would then be a blank). The segments must be 6-10 E F inputted in ascending order. 11-20 Q F Area = Interface areas between segments. 21-25 CORSEG* I E = Dispersion coefficient between segments. 26-30 Area F Q = Flow between segments. 31-35 E F CORSEG = Segment number that shares interface with present segment. 36-45 Q F The convention for the sign of flow Q.. between segments i and j is as follows:

Data Columns Variable Mode Description Card on Card Name No. 5 46-50 CORSEG* I If the flow is going from i to j Q.. is positive. 12 51-55 Area F If the flow is going from j to i Q..ij is negative. 56-60 E F If all the E values are same then enter the values in columns 6 to 10 of the first card (i.e., card #5) and 61-70 Q F leave the spaces for other values' blank. An identical double subscript indicates a boundary interface. To 71-75 CORSEG* I determine the sign of the flow, assume j to be outside the model. Thus, Qll is negative if the flow is from outside the model, and positive if the flow leaves segment 1. NOTE: Only the j part of subscript is specified by CORSEG. The i part of the subscript is assumed to be the segment number as determined by the card on which Q is found. Example: Suppose segment number 50 has a boundary interface. Then all the boundary values (area, E and Q) can be indicated either on 99th or 100th card. For example, boundary interface area on columns 51-55, E on 56-60, Q on 61-70 and CORSEG on 71-75 on the 99th card. CORSEG in this case is 50. IA 1-10 Length F This is the first card after all the areas, E's, Q's, and CORSEG's have been specified. The lengths and volume of 11-20 Volume F each segment are placed in the same order as the segments (i.e., ascending order). There are 4 lengths and volumes per card. 61-70 Length F 71-80 Volume F

Data Columns Variable Mode Description Card on Card Name No. 1B This is the first card after all lengths and volumes have been recorded. 1-3 NUMBC* I NUMBC = Number of segments which have boundary concentrations. 1C 1-7 BC F If NUMBC on card number 1B is zero, then skip this card. 8-10 IBCSEG* I If NUMBC on card number 1B is not zero, then input the boundary concentrations and corresponding segment numbers. 71-77 BC F Remember NUMBC on card lB specifies the number of boundary concentrations, and each should be accounted for on card 78-80 IBCSEG* I 1C. Eight boundary concentrations and segments are allowed per card. BC = boundary concentration IBCSEG = corresponding segment number 1D 1-10 Load G After the A matrix has been identified the source vector is formed. Loads are inputted with their corresponding 11-13 ISEG* I segment. Use a blank card for either no loads or after the last load card. In other words, a blank card is required to stop the loading input. One load per card.

Data Columns Variable Mode Description Card on Card Name No. 1E 1-5 K = First order reaction coefficients. Sixteen coefficients per card are permitted. Segments are in consecu6-10 K F tive order. Note: If all the K's are same, enter that value in column 1-5 of the first card and leave the remaining spaces blank. (All segments must be accounted for by using the correct number of blank cards.) 76-80 FOLLOWIING CARDS ARE NECESSARY ONLY IF OXYGEN DEFICIT RUN IS REQWESTED 2A 1-10 CSK F Saturated dissolved oxygen concentration at the tempera- t ture used in the run. 2B 1-5 KAK = First order reaeration coefficient. 6-10 Description is the same as for card number 1E. KAK F 76-80 2C 1-10 Load G Oxygen loads for different segments are recorded. 11-13 ISEG* I Description same as in card number 1D.

Data Columns Variable Mode Description Card on Card Name No. 2D 1-3 NUMBC* I Number of segments which have boundary concentrations for the deficit run. 2E 1-7 BC F If NUMBC on card 2D is zero, then skip this card. 8-10 IBCSEG* I If NUMBC is not zero, then input the boundary concentrations for the deficit run and the corresponding segment numbers. 71-77 BC F Description same as in card number 1C. BC must be expressed in the same units as in card number 1C. 78-80 IBCSEG* I *Variables denoted in this manner must be right justified within the specified field.

-19 -EXAMPLE Figure 2 is a map of Grand Traverse Bay which has been modeled employing the techniques described in this report. Shown is the location of the major discharges of pollution from two cherry processing plants and the Boardman River. Figure 3 shows the map of the lower part of the west arm of the Bay with the completely mixed segmentation scheme superimposed upon it. It is noted that forty-eight sections were defined and that segment sizes are smallest in the vicinity of point load discharges. The segment size increases as the boundary with the upper Bay is approached. Figure 3 also illustrates the circulation, river, and waste flow routing as estimated from a preliminary hydrodynamic model by Green (1971). For the calculations that follow the first order decay of total coliform was characterized using a reaction rate coefficient of.76 day1 and the dispersive mixing using E = 1 miles /day. Figure 4 shows the calculated total coliform distribution for 1964 in the Bay using surveyed coliform loading at the points shown in Figure 2. It is noted that shore-line concentrations are well above values recommended for total body contact recreational use. An example of the distribution of coliform that would result from 90 percent treatment of the existing cherry waste is shown in Figure 5. Such reduction

-20 -accompanied by the simultaneous introduction of a new coliform source at Cedar Creek is shown in Figure 6. The effect of a simple diversion of all wastes to deeper water by a 2000 ft. outfall pipe is illustrated in Figure 7. The above calculations are based on preliminary information and are intended only to serve as an example illustrating the utility of the computer program described herein. The successful application of these techniques is dependent upon "hard" verification of both circulation and water quality models with observed field data. In addition, the value of coefficients such as K and E must be consistent with user experience and laboratory testing. Thus caution must be coupled with sound professional engineering judgment to avoid erroneous or premature conclusions regarding the state of water quality in a given body of water.

-21 -Lake Michigan Grand Traverse Bay East Arm West Arm L,:.~.-4- J Mi les 0 5 Be a r dm an River TFrayerse~ City: Figure2 - Traverse Bay Region

OIi lid~~~~~~~~~~~~~~~~~~~~~v U~.) CD~~~~~~( PC) r\) CD ~~~C 'iCD~~~~f.A). ~ ~.0~~ L C+- N) r P. Itn ~~~~~1 ~ ~ ~ ~ ~ ~ = it j Et Co VD~~~~c cn ~~~~~~~~~~~~3C )\J A b~d co Qb ~q ID or r ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J — J~~~~~~~~~~~~~~~~~~~~~-so" r ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~~~~~~~~~~~~~~~r J. d- I-,!-, C- 0 -Clb O c+-~,0 a +3~~~~~~ — mt~j ~

~- *t96 tt *9 -TiD. ut - uot.nqTcJ szT.(I m;JOJTToo3 t113o0 - at an3.lF o 0 0 n I o o o o o 0 % 0 ^0 0 0 0 0 oo o o o 0 o o vn o.0 0 0 00 0 0 ) h r-i ~U ---?r ----e --— i~yp~_ C) C) O OO O ^^^/ - OO -i "- b *~~~, A,-.00* I o' T 0 0 1~/T/a!~ZL /c,,, ~ 4,,"~~~~~~~~~r I~00T./TT~~~~~~~~~~~~~~or

'SUKOTTO9;'ao 00" 0n 6 O O~ 0 C) 0 OO N~~~~~~N 0 T ~~/ 007:~~~~~~~~~~~~~~00 -* - ^^ IT^ - ^~~~~ ^~^ ^ ^.E WAS Tr wo-*~"^ OT JL " ~~~~~~~~~~~~"t~ twm "01

ST' a qnrn 3 a B 'i; 3;z~toS'Ei pti-. 194saiMa l3J 0 0~o low* OT~~~~~~~~~~~ Nl~~~~il 007: urrr~~~~~~~~~~~~~~~~~o, wrrr I~~~~~~~~~~~40 " ^-1^^ ^.I dop~~~~~~~~~~~~~~~~~~~~ oo^^^^' \ '^ ^ ^\ ^^ ^ ^o SHOP Or~~~~~~~~~~~~~~~~~~~O II "^~P ~I f ~9sl ^^ ^\^^O /~e- 41x ^ I~ ^ 1 / ^ ^ ^^ ^^J^~~~~OT. 7 ^ ^^f ^ "<.?.^~

(DST/ j / OO ( ^ I ^ -N /^ \ v. t q r00 1CI~~~~~lr, ~ angrP ^^OOT/TT~~~~~~~~~o ~[1^^ \~~~~~OO -98 -crrr~~~~~p -' N~ / oorg~~~~~~~~~~~~~~~~~~~~o 0 1 ~ 9 -9 —~~- /rr D~O -VaP, -V / TtUdoT/TTT A-PO, ~ ~ ~ O 'TMOOT/T~~~~~~~~~~~~~a;)~~*

-27 -ACKNOWLEDGMENTS The writers acknowledge the Sea Grant Program which has supported the effort required to write the computer program and prepare this manual. The computer program itself was largely the work of Mr. S. Nachiappan.

-28 -REFERENCES Canale, R. P., "A Methodology for Mathematical Modeling of Biological Production," Report to the University of Michigan Sea Grant Program (1970). Green, A. W., Jr., Unpublished progress report to the University of Michigan Sea Grant Program (1971). Thomann, R. V., "Time Variable Water Quality Models —Estuaries, Harbors, and Off-Shore Waters," from Advanced Topics in Mathematical Modeling of Natural Systems, 16th Summer Institute in Water Pollution Control.(1971).

UNIVERSITY OF MICHIGAN 3 9015 02651 8616