Mathematical Modeling and Analysis of Eutrophication of Water Bodies Caused by Nutrients

Abstract. In this paper a non-linear mathematical model is proposed fo r a qualitative representation of ecosystem dynamics in a eutrophied water body. The model variables are the concentration of nutrients, densities of algal popu lation, zooplankton population, detritus and the concentration of dissolved oxygen. The mod el consists of five coupled ordinary differential equations. By using the qualitative th ory of differential equations the model steady-state dynamics are studied. Simulation an alysis is also performed to see the effect of rate of input of nutrients on different variabl es participating in the proposed model.


Introduction
Eutrophication is a process in which a water body, such as a lake, becomes rich of nutrients (nitrogen, phosphorus, etc.), from domestic drainage as well as water run off from agricultural fields.These nutrients resulting in increased production of algae and macrophytes.Due to excessive growth of macrophytes in water and algae floating on the water surface, the photosynthesis process of aquatic flora decreases leading to decreased production of oxygen in the water body.This excessive growth also decreases the oxygen transfer from air to water by diffusion.When algae and zooplankton die and sink to the bottom of the water body, their decay by bacteria reduces the concentration of dissolved oxygen in the bottom of water body further to levels which may not be sufficient to support fish life, [1][2][3][4][5][6][7][8][9][10].
Several investigators have studied the effect of discharge of nutrients in water bodies such as a lake causing eutrophication, [1, 4-6, 9, 11-20].Voinov and Tonkikh [10] have presented a model for eutrophication in macrophyte lakes assuming that the nutrient is supplied only by detritus, which is formed by death of algae and macrophytes.They have not considered the input of nutrients either from domestic drainage or from water run off from agricultural fields.Truscott and Brindley [21] presented a model for evolution of phytoplankton and zooplankton in an ocean.Franke et al. [12] proposed a physicalbiological coupled model for algal dynamics in lakes by considering interactions of phytoplankton and zooplankton and their response to physical environment.Edwards and Brindley [4] have studied zooplankton mortality and the dynamical behavior of plankton population by considering nutrients, phytoplankton and zooplankton as variables.The influence of eutrophication on air -water exchange, vertical oxygen flux, etc. in the lake Ontorio has been investigated by Dachs et al. [2].The nitrification in the water column and sediment of a lake and adjoining river system has also been studied by Paur and Auer [22].Jayaweera and Asaeda [14] have studied biomanipulation in shallow eutrophic lakes by using a mathematical model involving phytoplankton, zooplankton, detritus, bacteria and fish population but they have not considered the supply of nutrients from outside.Some other ecological modeling studies involving phytoplankton, zooplankton and nutrients, relevant to our work, have also been conducted by [23][24][25][26][27], but they have not considered the concentration of dissolved oxygen in the modeling process.
Keeping in view of the above, in this paper, we model and analyze the eutrophication of a water body when nutrients are supplied to the water body from outside by water run off from agricultural fields using ecological concepts.

Mathematical model
We consider here a lake which is being eutrophied due to overgrowth of algae and other biological species caused by discharge of nutrients from domestic drainage as well as from water run off, etc. and also from nutrients formed from detritus.The bilinear interactions of variables such as the cumulative concentration of nutrients, densities of algae (phytoplankton) and zooplankton populations, density of detritus and concentration of dissolved oxygen are considered.We consider that various nutrients are supplied to the water body from domestic drainage as well as run off from agricultural fields.These nutrients may also be supplied by death of algae and zooplankton.We assume that the phytoplankton population is wholly dependent on nutrients and is being used as a food by its predator zooplankton population.It is assumed further that the level of dissolved oxygen in the water body increases by diffusion, etc. with a constant rate as well as by photosynthesis / respiration by algae.
Let n be the cumulative concentration of various nutrients, a be the density of algae, Z be the density of the zooplankton population, S be the density of detritus and C be the concentration of dissolved oxygen (DO).We assume that the cumulative rate of discharge of nutrients into the aquatic system from outside in the water body is q, a constant which is depleted with rate (α 0 n) due to natural factors.It is further assumed that the cumulative growth rate of nutrients by detritus is (π 0 δS) and depletion of nutrients by algae is proportional to both the density of algae as well as the concentration of nutrient (i.e.na).Thus, the growth rate of algae is proportional to (na) as it is assumed to be wholly dependent on the nutrients.The natural depletion rate of algae is assumed to be proportional to its density a and its depletion rate due to crowding / food utilization is proportional to a 2 in the water body.The depletion rate of algae by its predator zooplankton is assumed to be proportional to (aZ) and hence the growth rate of zooplankton is also proportional to (aZ).The natural depletion rate of zooplankton is assumed to be proportional to its density Z and its depletion rate due to crowding is proportional to Z 2 .Since some part of natural depletions of algae and zooplankton is converted into detritus, thus the growth rate of detritus is assumed to be proportional to a and Z and its natural depletion rate is assumed to be proportional to S. We consider that the rate of growth of dissolved oxygen by various sources is q c assumed to be a constant and its natural depletion rate is proportional to its concentration C. It is further assumed that the rate of growth of dissolved oxygen by algae is proportional to a and the depletion of DO caused by its use by detritus in converting itself into nutrients is proportional to its concentration S.

Analysis of equilibria
The model (1) has the following three equilibria.
Case 1. E 1 (q/α 0 , 0, 0, 0, q c /α 3 ) always exists.This equilibria of model (1) explains that if the density of algae as well as zooplankton both are not participating in the system then the equilibrium level of nutrients will reach to the value q/α 0 and the equilibrium concentration of dissolved oxygen will reach to the value q c /α 3 .Here we also note that since detritus is formed due to death of algae and zooplankton, both are not participating in the system, hence the equilibrium density of detritus will become 0.
exists, provided the following conditions are satisfied: The second equilibria E 2 of model ( 1) is obtained when algae is participating in the system whereas zooplankton population is not participating in the system.In this case the equilibrium level of concentration of nutrients, density of algae, density of detritus and concentration of dissolved oxygen will reach to the values n * 2 , a * 2 , S * 2 and C * 2 respectively.These values are explicitly given by equations ( 8), ( 12), ( 10) and ( 13) respectively.Case 3. E 3 (n * , a * , Z * , S * , C * ) exists, provided the following conditions are satisfied: The third equilibria E 3 of model ( 1) explains that if algal and zooplankton both populations are present in the water body then the equilibrium values of different variables i.e nutrients, algae, zooplankton, detritus and dissolved oxygen participating in the system will be given by n * , a * , Z * , S * and C * .These equilibrium values are explicitly given by equations ( 16), ( 14), ( 15), ( 17) and ( 18) respectively.We shall show the existence of these equilibria as follows.The equilibria of model ( 1) are obtained by the following algebraic equations: The equilibrium E 1 (q/α 0 , 0, 0, 0, q c /α 3 ) exists obviously.
Using this value of a = a * 2 > 0 in equations ( 8) and ( 10), we get which are positive.Again using the values of a * 2 and S * 2 in equation ( 11), we get which is positive if the right hand side of ( 13) is positive, gives condition (3).
Existence of E 3 .For the equilibria E 3 (n * , a * , Z * , S * , C * ), a * is the positive root of the following quadratic equation provided condition (4) is satisfied.
Substituting a = a * > 0 in equation ( 9), we get which is positive if right hand side is positive, gives condition (5).
Using this value of a * in equations ( 8) and ( 10), we get and which are positive under condition (5).Again using a * and S * in equation ( 11), we get which is positive if the right hand side is positive gives condition (6).
It may be pointed out that E 2 corresponds to the situation, where zooplankton is not present.Keeping in view that the algae is predated by zooplankton, it is noted that growth rate of algae in absence of zooplankton is greater than in its presence.Therefore using a comparison theorem, Hale [29] it can be easily concluded that the population of algae is greater in absence of zooplankton than when zooplankton is present.Thus a * 2 > a * .

Remark.
To see the effect of discharge of nutrients on eutrophication as well as on the level of dissolved oxygen, we determine the rates of change of a, Z, and S with respect to q and then corresponding change of C with respect to q.
From equations ( 14), ( 15) and ( 17) we find that da * dq , dZ * dq and dS * dq are all positive.From equation ( 18), we note that which is negative if Since most of algae float on the surface of water, therefore the oxygen formed by algae during photosynthesis may go to the atmosphere and has a little chance to get dissolve into the water, thus λ 11 is very small.Hence, the above condition ( 19) would be satisfied easily.This shows that as the cumulative rate of discharge of nutrients in the water body increases, the densities of algae, zooplankton and detritus increase but the concentration of dissolved oxygen decreases.

Stability analysis
In our analysis, we assume that all the above equilibria exist.The local stability behavior of these equilibria is studied in the following theorem.

Theorem 1. The equilibrium E
The proof is given in Appendix II.
In the following we prove that E 3 is nonlinearly stable.For this we need the following lemma, following, Freedman and So [30], Shukla and Dubey [8,31].Lemma 1.The set Ω is a region of attraction for all solutions initiating in the positive octant. where Proof.Adding the first four equations of model ( 1), we get By using comparison theorem, Hale [29] it can be easily concluded From the last equation of model ( 1), we have Using the maximum value of a, we get Again using comparison theorem, we get It is noted here that if this condition is satisfied then (20) would always be satisfied.It is further noted that this condition is feasible for very large α 0 and small π 0 .
The proof is given in Appendix III.
The above theorems imply that under certain condition, the system variables would attain their equilibrium values and the densities of algae, zooplankton, detritus increase as the cumulative rate of discharge of nutrients increases but the concentration of dissolved oxygen decreases.

Numerical example
To check the feasibility of our analysis regarding the existence of E 3 and corresponding stability conditions, we conduct some numerical computation by choosing the following values of the parameters in model (1).
It is found that under the above set of parameters, conditions for the existence of interior equilibrium E 3 (n * , a * , Z * , S * , C * ) are satisfied and E 3 is given by n * = 3.093580, a * = 1.137432,Z * = 0.318716, S * = 0.655267, C * = 1.096836.
Thus, it has three real and two complex eigenvalues, which are either negative or have negative real parts.Hence E 3 is locally stable.
It is pointed out here that for the above set of parameters, the conditions for local stability (20) and nonlinear stability ( 22) are also satisfied.
Further for the above set of parameters, a computer generated graph of n verses C is shown in Fig. 1, which indicates the global stability of (n * , C * ) in n-C plane.In Figs.2-5, we see the effect of q on a, Z, S and C, by keeping other parameters constants (as given in (23).We observe here that as the rate of input of nutrients q increases, the density of algae, zooplankton and detritus increases whereas the concentration of dissolved oxygen C decreases.In these figures we are also seeing that if   the rate of introduction of nutrients q = 0, then the density of algae, zooplankton and detritus are tends to zero, whereas the concentration of dissolved oxygen reaches to its maximum level.This result is clear because of the fact that the nutrients then formed from detritus will not be sufficient for the growth of algal population, zooplankton population and detritus.

Conclusion
In this paper, a non-linear mathematical model is proposed and analyzed to study the eutrophication of a lake due to overgrowth of algae, zooplankton and other biological species caused by excessive supply of nutrients from water run off from agricultural fields and other sources.The modeling analysis shows that as the supply of nutrients in the water body increases, the densities of algae and other biological species increases initiating the eutrophication process.It has also been shown that the density of detritus increases correspondingly leading to a decrease in the concentration of dissolved oxygen (as oxygen, formed due to photosynthesis by floating algae, goes to the atmosphere and does not affect the concentration of dissolved oxygen).By simulation analysis it has been shown that if the rate of input of nutrients from outside is high then the concentration of dissolved oxygen may become negligible in a water body.

Appendix I: Stability of E
The general variational matrix for the model ( 1) is given as follows: We shall use eigenvalues theory to study the local stability behavior of equilibria E i (i = 1, 2).Let M 1 be the matrix obtained from M after substituting for E i (i = 1, 2).
For the stability of E 1 (q/α 0 , 0, 0, 0, q c /α 3 ), we note that one of the eigenvalues of M 1 is θ1β1q−α0α1 α0 , which is clearly positive whenever E 2 exists (condition (2)).Hence E 1 is unstable whenever E 2 exists.For the stability of E 2 (n * 2 , a * 2 , 0, S * 2 , C * 2 ), we note that one of the eigenvalues of M 2 is θ 2 β 2 a * 2 − α 2 , which is positive whenever E 3 exists (as a * 2 > a * and condition (5)).This shows that E 2 is unstable whenever E 3 exists.As we cannot say much about the behavior of E 3 from the corresponding variational matrix M 3 , we study the behavior of this equilibrium by using Liapunov's method.
From conditions (A.3)-(A.7),we can choose a positive m 3 and hence m 4 if the following condition is satisfied:

Fig. 2 .
Fig. 2. Variation of a w.r.t.t for different q.Fig.3. Variation of Z w.r.t.t for different q.

Fig. 3 .
Fig. 2. Variation of a w.r.t.t for different q.Fig.3. Variation of Z w.r.t.t for different q.

Fig. 4 .
Fig. 4. Variation of S w.r.t.t for different q.Fig. 5. Variation of C w.r.t.t for different q.
β 1 a * < From conditions (A.10)-(A.12),we can choose a positive m 3 and hence m 4 if the following condition is satisfied: