AQUATIC WATER QUALITY MODELLING

1. AQUATIC WATER QUALITY MODELLING August 8, 2007 Research Professor Tom Frisk Pirkanmaa Regional Environment Centre P.O.Box 297, FIN-33101 Tampere, Finland E-mail tom.frisk@ymparisto.fi

2. 2. MASS BALANCE MODELS FOR LAKES 2.1 General Mass balance models for lakes are simple models based on considerations of mass balances of lakes, which can be used in predicting average concentrations. Most mass balance models have been formed for total phosphorus but the same approach can be in principle applied to all substances. Mass balance models are most often steady state models. Thus they are not used in predicting temporal variations of state variables but long term average values.

3. At steady state the mass balance enquation (1.6) can be • Written as follows: • I - O - S = 0 (2.1) where : I = input of the substance considered (M T-1) O = output of the substance considered (M T-1) S = the rate of decrease of the total amount of the substance considered, i.e. net retention (M T-1). Obs! Minus sign in the equation.

4. When total phosphorus is concerned, sedimentation and • release from the sediment are in practice the significant • Processes and so the term S in Eq. (2.1) is called net • sedimentation: • S = Sbr - Sr (2.2) where: Sbr = gross sedimentation (M T-1) Sr = release from the sediment (M T-1). In mass balance models, simple hydraulic descriptions are applied. If the model is applied to a lake CSTR description is normally used. In long term consideration this is justified if the lake has regular turnover periods and the average concentrations of the whole water body are calculated.

5. In elongated lakes also PFR hydraulics or the linear Combination of CSTR and PFR can be applied. According to the CSTR assumption (Eq. 1.7) Eq.(2.1) can be written in the following form: I - Q c - S = 0 (2.3) Input (I) and discharge (Q) can be measured. Accurate measurement of net sedimentation is extremely difficult. It can be calculated by means of mass balance equations. In order to form a model on the basis of Eq. (2.3) which could be used to making water quality predictions, net sedimentation must be calculated on the basis of other variables.

6. Two approaches are available: 1) Based on sedimentation coefficients 2) Based on retention coefficients. 2.2 Models Based on Sedimentation Coefficients In the approach bases on sedimentation coefficient, net sedimentation is described as proportional to phosphorus concentration (or its function) and the volume (or area) of the lake. Vollenweider’s (1969) model is based on CSTR hydraulics and description of phosphorus net sedimentation as a first order reaction.

7. The basic equation of the model is the same as presented before (Eq. 1.11): dc V---- = I - Q c – σ c V (2.4) dt where: V = volume of the lake (L3) c = phosphorus concentration in the lake (M L-3) t = time (T) Q = discharge at the outlet of the lake (L3 T-1) σ = first order net sedimentation coefficient of phosphorus (T-1)

8. The solution of Eq. (2.4) with the initial condition c=c0 when t=t0 is the following (I and Q are constant values): Q -t(σ + ---) I I V c(t) = ---------- - (----------- - c0) e (2.5) Q + σV Q + σV At steady state (dc/dt = 0) the following solution is obtained: I css = ---------- (2.6) Q + σV where: css = phosphorus concentration at steady state (M L-3).

9. Let us use the notation ρ = Q/V. Eq. (2.5) can now be • rewritten as: • -(σ + ρ)t • c(t) = css - (css - c0)e (2.7) • where: • ρ = water renewal ceofficient or flushing coefficient (T-1) = Q/V. According to Eq. (2.7) steady state is reached after an eternally long time. The transition time between to successive steady states can be described utilizing the concept of the half-life time. Let us assume that at the original steady state phosphorus concentration is= c0 and at the new steady state (wit a new loading level) it is = css.

10. Half-life time (t½) is the time which is needed to the change that the concentration has exactly the half way value between the old and new steady states: c0 + css c(t½) = ---------- (2.8) 2 On the basis of Eq. (2.7) the following formula is easily obtained: ln 2 t½ = ------- (2.9) σ + ρ

11. Even though the theoretical transition time between the successive steady states is eternal according to the model, a sufficiently accurate estimate for the transition time is obtained by dividing the half-life time (Eq. 2.9) by a coefficient with a value from 3 to 5. Vollenweider’s model can be used in calculating phosphorus concentrations. The average phosphorus concentration can be calculated with Eq. (2.6) when the average value of phosphorus concentration remains constant.

12. Example 2.1. The volume of the lake is = 50 *106 m3, discharge = 5 m3s-1, phosphorus input = 5.7 t a-1 and phosphorus net sedimentation coefficient = 1.1 a-1. Calculate the average phosphorus concentration in the lake! Eq. (2.6) is applied: 5.7 t a-1 css = ----------------------------------- = 5 m3s-1 + 1.1a-1 * 50 *106 m3 5.7 t a-1 *109 mg t-1 ------------------------------------------------------------- = 5 m3s-1 *(60*60*24* 365 s a-1) + 1.1a-1 *50* 106 m3

13. 5.7 *109 mg a-1 --------------------------------------------- = 157.68 * 106 m3 a-1 + 55 * 106 m3 a-1 26.8 mg m-3 = 27 µg l-1 I, Q and V can be measured. How about σ ? It cannot be direcly measured. However, it can be calculated on the basis of mass balance data. The following formula can be solved from Eq. (2.6): I - Q c σ = --------- (2.10) c V

14. Predictions of phosphorus concentration in changing loading conditions can be made (according to the original Vollenweider’s model) by first calculating σ in the prevailing loading condition utilizing measured values of I, Q and V and Eq. (2.10). Prediction of phosphorus concentration in a new loading condition is then calculated with Eq. (2.6) using the calculated value of σ and the new value for I. Example 2.2. The (average) phosphorus concentration of the lake of example 2.1 is measured. It is on = 25.6 µg l-1. What is the real value of the sedimentation coefficient? How would phosphorus concentration change if phosphorus concentration change if input decrease to the value of 5.0 t a-1?

15. Eq. (2.10) is used: 5.7 t a-1 - 5 m3s-1 *31.536 * 106s a-1 *25.6 mg m-3 σ = ----------------------------------------------------------- 25.6 mg m-3 * 50 * 106 m3 5.7 t a-1 *109 mg t-1 – 4.052 *109 mg a-1 = ----------------------------------------------- 1.28 *109 mg = 1.29 a-1

16. Here Eq. (2.6) is applied: 5.0 * 109 mg a-1 c = ---------------------------------------------------- 157.68 * 106 m3 a-1 + 1.29 a-1 *50* 106 m3 = 22.5 mg m-3.23 µg l-1 The concentration will be about 3 µg l-1 lower. The procedure applied above is based on the assumption that σ remains the same when loading changes. In reality this is not true but F is dependent on loading and its value is not the same for different lakes.

17. The best model for calculating the sedimentation coefficient has been presented by Canfield and Bachmann (1981). According to their statistical model sedimentation coefficient can be calculated by means of the ratio of phosphorus Input and the volume of the lake: • I β • σ = σs (----) (2.11) • V • where: • σs ja β = constants.

18. Canfield and Bachmann calculated the values of the coefficients σs ja β separately for lakes and artificial lakes. They obtained the following values: • σs (a-1) β • Lakes 0.162 0.458 • Artificial lakes 0.114 0.589 • The whole data 0.129 0.549 When the parameter values presented above are used input I must be given in 1 mg a-1 and volume V in 1 m3. Sedimentation coefficient is obtained in the unit 1 a-1.

19. The model of Canfield and Bachmann can be applied to different types of lakes. Eq. (2.11) is used in calculating σ, which can be used in Vollenweider’s model (Eq. 2.6) and the model can then be applied in new loading conditions. Example 2.3. Calculate example 2.2 using Canfield- Bachmann model! • Original loading condition • 5.7 * 109 mg a-1 • σ= 0.162 a-1( ------------------- )0.458 = 1.418 a-1 • 50 * 106 m3

20. 5.7 * 109 mg a-1 c = ---------------------------------------------------- 157.68 * 106 m3 a-1 + 1.42 a-1 *50* 106 m3 = 24.9 µg l-1 2) New loading condition: 5.0 * 109 mg a-1 σ = 0.162 a-1(------------------- )0.458 = 1.335 a-1 50 * 106 m3 5.0 * 109 mg a-1 c = ---------------------------------------------------- = 22.3 µg l-1 157.68 *106 m3 a-1 + 1.34 a-1 *50 * 106 m3 The phosphorus concentration would be 2.6 µg l-1 (ca. 3) lower.

21. In lakes resembling rivers or otherwise elongated lakes PFR descriptions are often applied. In the first case the Additional input along the lake is assumed to be negligible. If net sedimentation of phosphorus is assumed to obey first order kinetics concentration distribution at different points of the lake can be calculated according to Eq. (1.23):

22. - σ x/u c(x) = c0 e (1.12) where: c0 = phosphorus concentration at the starting point where x=0 (M L-3) c(x) = phosphorus concentration at the distance x (M L-3) u = average velocity of flow (L T-1) The average velocity of flow can be calculated by dividing the discharge by the cross-sectional area: u = Q/Ax (1.13) where: Q = discharge(L3 T-1), Ax = average cross-sectional area(L2)

23. The average cross-sectional area can be calculated • by dividing the volume by the length of the lake: • Ax = V/L (2.14) • where: L = the length considered (L) • When PFR is applied second order kinetic have found to be better then first order kinetics for phosphorus. The • basic equation is then the following. • dc • --- = - α c2 (2.15) • dτ

24. where: α =second order net sedimentation coefficient (M-1 L3 T-1) τ = travel time (T) = x/u • The solution of Eq. (6.15) with the initial condition c=c0 when τ = 0 is the following: • c0 • c(x) = --------------- (2.16) • 1 + α c0 x/u The values of the sedimentation coefficients σ and α can be easily calculated by means of Eqs. (2.12) and (2.16), respectively, when phosphorus concentration at two different sites is known. It should be borne in mind that also here average concentrations are referred to, not single “plugs”.

25. Example 2.4. The length of an elongated lake is 12 km and volume 18 * 106 m3. The average phosphorus concentration at the uppermost site of the lake is = 22 µg l-1 and at the outlet of the lake = 16 µg l-1. The discharge through the whole lake is = 3 m3s-1 (additional discharge from the nearest drainage basin is assumed to be negligible). Calculate the first and second order net sedimentation coefficients in the lake! The average velocity of flow can be calculated using Eqs. (2.13) and (2.14): u = QL/V = 3 m3s-1 *12 * 103 m/18* 106 m3 = 2 * 10-3 ms-1

26. 1. order (Eq. 2.12): σ = ln (c0/c(x)) * u/x = ln (22/16)* 2 * 10-3 ms-1/12 * 103 m = 0.053 * 10-6 s-1 = 1.67 a-1. 2. order (Eq. 2.16): α = u(c0 - c(x))/c0c(x)x = 2 *10-3 m s-1(22-16) mg m-3/22 * 16 mg2 m-6 *12*103 m = 0.002841 *10-6 mg-1 m3 s-1 = 0.090 mg-1 m3 a-1.

27. As already mentioned, PFR can be modified so that additional discharge can be taken into account. • If net sedimentation is described as a first order reaction • Eq. (1.24) can be written in the following form: • dc SQ SI • ---- = - ---- c + ---- - σ c (2.17) • dτ Ax Ax Where: SQ = additional discharge per unit length (of the Lake or river) (L2T-1), SI = additional input per unit length (M L-1 T-1), Ax = cross-sectional area (L2).

28. The solution of Eq. (6.17) is the following: • -Kτ -Kτ • c(τ) = c0e + cs(1 - e ) (2.18) • where: • c(τ) = concentration at travel time τ(M L-3) • c0 = c(0) K = σ + SQ/Ax cs = SI/(SQ + σ Ax) = steady state concentration (at τ = ∞)

29. The first order net sedimentation coefficient can be calculated on The basis of Eq. (2.18): 1 c0 – cs SQ σ = --- ln ( ---------- ) - ----- (2.18b) τ c(τ) – cs Ax The solution (2.18) is based on the assumption that the cross-sectional area is the same throughout the whole lake. If there is also additional input (SQ) average velocity of flow (u) cannot remain the same but it must increase along the lake.

30. It can be calculated as follows: • Q0 SQ • u(x) = ---- + ---- x (2.19) • Ax Ax • where: • u(x) = average velocity of flow at the distance x • from the starting point (L T-1) Q0= discharge at the starting point(L3 T-1)

31. Travel time (τ) to the distance x is calculated in the following way: x dx Ax SQ τ(x) = ∫ ----- = --- ln(1 + ---- x) (2.20) 0 u(x) SQ Q0 where: τ(x) = travel time (detention time) to the distance x from the starting point (T).

32. If net sedimentation is described as a second order Reaction the basic equation is the following: dc SQ SI --- = - --- c + --- - α c2 (2.21) dτ Ax Ax where: α = second order net sedimentation coefficient (M-1 L3 T-1) For solving Eq. (2.21) the derivative of the left hand side is set as zero and the second degree equation is solved.

33. The solutions are: • r1 = -(SQ/Ax + ((SQ/Ax)2 + 4α SI/Ax)½)/2α (2.22) • and • r2 = -(SQ/Ax - ((SQ/Ax)2 + 4α SI/Ax)½)/2α (2.23) • where : • r1 and r2 = the zero points of the second degree polynom • on the left hand side of Eq. (2.21)

34. Of these zero points r1 is always negative but r2 is always positive and it represents the steady state concentration which would be reached if the travel time were eternal positive (similar to cs in Eq. 2.18). Solution of Eq. (2.21) has the form: -(r1 - r2)ατ r1(c0 - r2) - r2(c0 - r1)e c(τ) = ------------------------------------------- (2.24) -(r1 - r2)ατ (c0 - r2) - (c0 - r1)e If c0 = r2 the concentration remains the same (c0) in the whole system.

35. Travel times τ for different distances x can be again calculated using the following equation: x dx Ax SQ τ(x) = ∫ ----- = --- ln(1 + ---- x) (2.25) 0 u(x) SQ Q0

36. Fig. 2.1 (Frisk 1989)

37. 2.3 Models Based on Retention Coefficients In models based on sedimentation coefficients, net sedimentation is described as proportional to phosphorus concentration, whereas in models based on retention coefficients net sedimentation is described as proportional to input: • S = R I (2.26) • where: • R = retention coefficient

38. When Eqs. (2.269 and (2.3) are applied together we • receive the formula for calculating concentration: • c = (1 - R) I/Q (2.27) • Retention coefficient can be calculated on the basis o • mass balance data using Eqs. (2.1) and (2.26). Steady • state is a requirement. • R = (I - O)/I (2.28)

39. There are many models for calculating the retention coefficient using the basic variables (I,Q,V,A). In the models developed in the USA and Canada its most often assumed that input (I) has no effect on the value of retention coefficient. Eg. Chapra (1975) has presented the following model: ν R = -------- (2.29) ν+ qs where: qs = areal hydraulic loading (L T-1) = Q/A, A = area of the lake (L2), ν= apparent settling velocity (L T-1)

40. According to Chapra (1975) ν is constant and at least almost • the same for the different lakes (16 m a-1). Eq. (2.29) is in fact • based on a kinetic description where net sedimentation is • described as proportional to phosphorus concentration • and the area (not volume) of the lake: • S = ν c A (2.30) Thus apparent settling velocity (ν) and sedimentation coefficient according to Vollenweider’s (1969) model (σ) have the following relationship: • σ = ν/z (2.31) • where: • z = average depth of the lake (L) = V/A.

41. However, apparent settling velocity is not constant but its has been found to vary between 10 - 20 m a-1. Lappalainen (1974) presented a model for calculating retention coefficient in which input has effect on the value of R. In 1977 Lappalainen rewrote his model in the following form: • (cI - 6)T • R = 0.9 ------------------ (2.32) • 200 + (cI - 6)T • where: cI = mixing concentration (or initial concentration of phosphorus = I/Q, must be given in 1 mg m-3

42. T = theoretical detention time or residence time = V/Q, must be given in months (in calculations 1 month = 30 d and 1 a = 12 months). Eq. (2.32) can also be simplified to the following form (Frisk 1978): • cI T • R = 0.9 -------------- (2.33) • 280 + cI T

43. Fig, 2.2 Chapras’s (1975) model: 1: ν = 16 m a-1, 2: ν = 10 m a-1 (Frisk 1978) R = f( (cI – 6)T) according to Lappalainen (1974)

44. Phosphorus concentration in the lake can be calculated by means of Eq. (2.27) utilizing Eq. (2.32.), (2.33) or (2.29). For calculating retention coefficients there are also many other models. Example 2.5. Calculate the estimate of the impact of the reduction of phosphorus loading in Examples 2.1 – 2.3 using models based on retention coefficient! The average depth of the lake z = 10 m. The basic equation is (2.27). Mixing concentrations (cI=I/Q) In the original (cI1) and new (cI2) loading situation are first calculated:

45. 5.7 t a-1 *109 mg t-1 cI1 = ----------------------------------- = 36.149 mg m-3 5 m3 s-1 *365* 86400 s a-1 5,0 t a-1 *109 mg t-1 cI2 = ----------------------------------- = 31.710 mg m-3 5 m3 s-1 *365* 86400 s a-1 1) Chapra’s model, Eq. (2.29): 50*106 m3 A = V/z = --------------- = 5* 106 m2 10 m

46. 5 m3 s-1 *365* 86400 s a-1 qs = Q/A = ----------------------------------- = 31.536 m a-1 5 106 m2 16 m a-1 R = ------------------------------- = 0.3366 16 m a-1 + 31.536 m a-1 Original loading situation: c1 = (1 – 0.3366)* 36.149 mg m-3 = 23.98 mg m-3 According to the model R does not change when I changes c2 = (1 – 0.3366)* 31.710 mg m-3 = 21.04 mg m-3

47. The concentration will be 2.9 mg m-3 (ca. 3 µg l-1) lower. 2) Lappalainen’s model, Eq. (2.32): 50* 106 m3 T = V/Q = ----------------------------------- = 3.858 month 5 m3 s-1 * 86400 * 30 s month-1 Original loading situation: (cI1 - 6)T = (36.149 - 6)mg m-3 *3.858 month = 116.31 mg m-3 month

48. 116,31 R1 = 0.9 ------------------ = 0.3309 200 + 116.31 c1 = (1 – 0.3309)*36.149 mg m-3 = 24.19 mg m-3 New loading situation: (cI2 - 6)T = (31.710 - 6)mg m-3 *3.858 month = 99.189 mg m-3 month 99.189 R2 = 0.9 ----------------- = 0.2984 200 + 99.189

49. c2 = (1 – 0.2984)* 31.710 mg m-3 = 22.25 mg m-3 The concentration will be 1.9 mg m-3 (ca. 2 µg l-1) lower. 3) Lappalainen’s model, Eq. (2.33): Original loading situation: cI1T = 36.149 mg m-3 *3.858 month = 139.46 mg m-3 month 139.46 R1 = 0.9 ------------------ = 0.2992 280 + 139.46 c1 = (1 – 0.2992)* 36.149 mg m-3 = 25.33 mg m-3

50. New loading situation: cI2T= 31.710 g m-3 *3.858 month = 122.38 mg m-3 month 122.38 R2 = 0.9 ------------------ = 0.2737 280 + 122.38 c2 = (1 – 0.27347* 31.710 mg m-3 = 23.03 mg m-3 The concentration will be 2.3mg m-3 (ca. 2 µg l-1) lower.