Eutrophication Bifurcation Analysis for Tasik Harapan Restoration

Regime shift is characterized by large, abrupt and persistent changes in ecosystem structure and dynamics. Bifurcation analysis is commonly used to identify regime shift equilibrium states and to distinguish their stability characteristics. Eutrophication in lake, a regime shift from clear-water oligotrophic state, is a stable equilibrium state that can persist for long duration. Characterized by undesirable turbid water condition, eutrophication has been known to impair valuable ecosystem services provided by lakes worldwide. The high incidence of eutrophication in Malaysian lakes (62%) mandates urgent need for lake restoration. The three-fold objectives of this paper are (1) to develop a mathematical model for analysing bifurcation criteria in regime shift, (2) to identify regime shift thresholds and (3) to propose effective ecosystem management strategy for shallow tropical lakes such as Tasik Harapan. A mathematical model consisting of four compartments: algae, phosphorus, dissolved oxygen and biochemical oxygen demand is formulated to analyse the eutrophication dynamics in the highly eutrophic Tasik Harapan (TH), a small shallow lake in Universiti Sains Malaysia (USM). Bifurcation analysis is performed by means of XPPAUT to identify the regime shift thresholds and to determine the type of lake response. Identified as irreversible, the eutrophication state of TH mandates an urgent lake restoration program to remove nutrients in the lake. Two restoration methods for reducing nutrients are assessed: (1) flushing of lake water and (2) hypothetical application of the invasive water hyacinth. Bifurcation analysis indicates that a flushing rate exceeding 0.042 day is needed to restore TH to oligotrophic state. A complicated strategy of water hyacinth application would reduce the algae concentration from 300 μg/l to 120 μg/l after 9.6 years. A severe shortfall is the complicated and uncertain process of frequent removal of water hyacinth to prevent the lake from being overwhelmed by the invasive plants. The feasibility and sustainability of these two lake restoration methods are discussed. The insights gained would be useful to the relevant authorities in determining and implementing the best remediation measures for TH.


I. INTRODUCTION
Increasingly reported in a variety of aquatic systems worldwide, regime shift is characterized by large, abrupt and persistent changes in ecosystem structure and dynamics [1], [2]. Having more than one equilibrium state, an ecosystem may undergo a regime shift that suddenly shift from one Manuscript  equilibrium state to another when a critical threshold is crossed. Examples of regime shift include kelp transition [3], coral reef degradation [4], soil salinization [5] and lake eutrophication [6]. Regime shift may occur unexpectedly and is difficult to predict. Restoration efforts on reversing the shift may be difficult or even impossible [7], [8]. In lake eutrophication, a clear water state (oligotrophic) turns into a turbid state (eutrophic) when the nutrient level exceeds a threshold, as a result of the inflows of excessive wastes derived from human activities such as industrial, agricultural and urban domestic sewage. This would cause water quality degradation, public health risks and economic losses [9].
Several regime shift models have been used to study eutrophication worldwide [10]- [14], to propose eutrophication managements [8], [15] and to provide early warning of regime shift in lakes [16], [17]. The possibility of dramatic change in ecosystem could be illustrated through regime shift model, allowing early restoration efforts. Further, effective ecosystem managements based on the regime shift theory could be proposed. Examples of lake restoration methods include reduction of external phosphorus loading, flushing of lake water and the use of macrophyte to remove nutrients and to compete with algae. Flushing can improve water quality in eutrophic lake by removing algae and nutrients from the lake water. Another way to reduce eutrophication is to introduce macrophytes, such as water hyacinth, in the lake to remove nutrients from the lake water and to compete against algae growth. Water hyacinth is an aquatic macrophyte that belongs to the family Pontederiaceae [18]. It is a highly invasive aquatic weed due to its rapid growth rate and high adaptability to extreme conditions [18]. It can survive in both tropical and temperate climatic conditions and is able to store nutrients for later stages of life cycle. Water hyacinth has been used in phytoremediation and wastewater treatment due to its high capacity for nitrogen and phosphorus absorption, and its ability to grow in heavily polluted water [18]. Because it is highly invasive, extreme care must be taken when water hyacinth is introduced in any water body.
A study by National Hydraulic Research Institute of Malaysia (NAHRIM) indicated that 56 lakes (62%) are in eutrophic state while the balance (34 lakes or 38%) is in mesotrophic state [19]. This undesirable state of lake water quality calls for urgent need for lake restoration efforts to safeguard quality water supply. Although regime shift models are widely used to study lakes' eutrophication, most of these models were developed for temperate lakes. To address water quality in tropical lake, this paper develops a mathematical model to identify regime shift thresholds and characteristics for tropical lakes and to propose effective ecosystem management by means of regime shift theory. For this purpose, a mathematical model consisting of four compartments of algae, phosphorus, dissolved oxygen and biochemical oxygen demand is formulated and used to simulate eutrophication dynamics in the highly eutrophic Tasik Harapan⎯a small lake in Universiti Sains Malaysia (USM). Bifurcation analysis is performed by means of XPPAUT to identify the regime shift thresholds and to determine the type of lake response. To improve the eutrophication state of TH, two potential lake restoration methods are assessed: (1) flushing of lake water and (2) hypothetical application of the invasive water hyacinth. The insights derived from this analysis would assist the relevant authorities in determining the best remediation measures for TH. Sustainable development of water resources would ensure a clean water supply and maintain a healthy lake ecosystem. This is in line with Malaysia's pursuit of achieving the United Nation Sustainable Development Goals (UN SDGs), namely SDG6 which calls for clean water and sanitation [20].

II. STUDY SITE
Tasik Harapan (TH) is a small lake in USM of Penang, Malaysia, with a surface area of 6070 m 2 (1.5 acres) and a depth of 1.0 to 1.5 m [21]. It is constructed for flood mitigation in USM in the early 1980s. The lake water has turned eutrophic and green due to excessive algae growth, as a result of nutrients accumulation in the sediment over the four decades. A high algae concentration was occasionally reported in TH, i.e., 300 μg/L, indicating a highly eutrophic lake [21]. Based on Chapra [22], a lake is classified as eutrophic if the algae concentration exceeds 10 μg/L. Fig. 1 shows the location of TH in USM.

A. A-P-DO-BOD Model
A mathematical model consisting of four compartments, algae (A), phosphorus (P), dissolved oxygen (C) and biochemical oxygen demand (L) is formulated as in (1) to (4) and used to simulate the eutrophication dynamics in TH [22]- [24]. Algae concentration is a commonly-used indicator for eutrophication and phosphorus is the limiting nutrient required by algae for growth. The sources of phosphorus include external phosphorus loading (lp), excretion associated with zooplankton grazing (egA) and recycling from sediment (rP q /(P q + n q )). The loss of phosphorus is due to flushing (h1P) and uptake by algae for growth (bA(P -paA)/(ha + P -paA)). Phosphorus is a nutrient limiting factor for algae growth and is represented by the term (P -paA)/(ha + P -paA). The algae loss is due to flushing (h1A), grazing by zooplankton (gA) and mortality (saA). The sources of dissolved oxygen (DO) include the photosynthesis by algae and atmospheric reaeration. The maximum oxygen production rate from photosynthesis at saturating lighting condition is considered and represented by the term pmax = 9.6 × 1.036 (T-20) , where the water temperature T = 27C is used. Biochemical oxygen demand (BOD) measures the amount of DO consumed by microorganisms during decomposition of organic matter. In this model, DO and BOD do not directly affect the algae and phosphorus compartments. However, their inclusion in the model is needed for model calibration involving DO. Since the DO data for TH is available and the diurnal cycles of DO is algae-dependent, the algae dynamics can be estimated from the DO data. This provides estimates of important parameter values such as external phosphorus loading rate, algae growth rate and phosphorus recycling rate for bifurcation analysis and for the assessment of lake restoration methods using A-P (1)-(2) and A-P-WH (5)-(7) models. Table I displays the definition and value of the parameters in the model.

B. A-P-WH Model
The hypothetical use of water hyacinth in removing nutrients from TH is demonstrated by the A-P-WH model in (5)- (7). In this model, the maximum carrying capacity K of water hyacinth represents the maximum population of water hyacinth that can be supported in TH. Further, the growth of water hyacinth is limited by the availability of phosphorus in the lake, which is represented by the term P/(P+hw) in (7). A constant harvesting rate of water hyacinth, hW is considered in (7) to control the water hyacinth population in TH, a complicated and uncertain process. Definition and value of the parameters for water hyacinth are displayed in Table I IV. SIMULATION STUDY

A. Parameter Estimation by DO Curve Fitting
First, curve fitting of the simulation result to match the observed DO data in TH is performed by using the Particle Swarm Optimization (PSO) in Matlab [25]. In the curve fitting process, (1)-(4) are first solved by using RK4 method. Then, interpolation is performed to estimate the value of C (dissolved oxygen). PSO is used to minimize the distance between the observed DO data and simulated C subject to a range of parameters. Optimum solution, i.e. the values of the parameters b, lp, r, ka, k1 and lBOD, representing the best fit are obtained. The parameter values obtained are shown in Table I, with the reported range of parameter values, provided for reference. Fig. 2(a) displays the curve fitting for DO, C (mg/L) in TH. The curves for algae, phosphorus and BOD, as indicated in Fig. 2(b)-2(d), are then obtained using these fitted values.  Table I.

B. Bifurcation Analysis
Bifurcation analysis is performed by means of XPPAUT to identify the regime shift thresholds. XPPAUT is an open-source package that incorporates the bifurcation package AUTO. Developed by Bard Ermentrout, from Department of Mathematics, University of Pittsburgh, XPPAUT is a general numerical tool for simulating, analysing and animating dynamical systems. Analysis such as stability analysis, bifurcation analysis, nullclines and vector fields can be performed by using XPPAUT. Designed to be user-friendly, bifurcation analysis via XPPAUT only requires the equations, parameters, variables and boundary conditions of the model. Bifurcation diagrams can then be plotted within the program using various menus and buttons. The details of XPPAUT and user guideline can be referred to Ermentrout [32].
The bifurcation diagram of algae, A against external phosphorus loading, lp is plotted in Fig. 3. From Fig. 3, TH shows irreversible behaviour because the lake is still eutrophic (A = 131 μg/L) even when lp is reduced to zero (no P inflow). Therefore, other restoration method such as flushing or biomanipulation is needed in conjunction with reduction in lp to further reduce the algae concentration in TH.
Three equilibria (two stable and one unstable) exist when 0 μg/L/d < lp < 0.01595 μg/L/d. When lp increases, A increases along the stable steady state E3 until a certain value of A is reached, at the right-side inflection point of equilibrium (lp = 0.01595 μg/L/d, A = 7.353 μg/L). Further increase in lp causes the equilibrium to bifurcate or "jump" to another stable steady state E1 and remain in E1, where eutrophication occurs. If the value of A lies between E1 and E2, it will be attracted to the stable steady state E1. The trajectories of A will never approach the unstable steady state E2. An algae concentration, A, that is below E2 with lp < 0.01595 μg/L/d will approach the stable steady state E3 (A < 7.353 μg/L), indicating mesotrophic state of the lake and oligotrophic if A < 2.6 μg/L. Based on this analysis, the regime shift threshold for the eutrophic TH is lp = 0.01595 μg/L/d. Any lp value greater than this threshold will cause a sudden increase in algae concentration to eutrophic status.
A sudden large increase from A = 7.353 μg/L to more than 132 μg/L in TH is observed in Fig. 3 when the critical threshold lp = 0.01595 μg/L/d is crossed. This is due to two reasons, i.e., the phosphorus internal recycling and zero-flushing rate considered in TH. In order to identify the effects of phosphorus internal recycling on the critical threshold, Fig. 4 is plotted by considering two different values of phosphorus internal recycling rates of r = 0.1 μg/L/d ( Fig. 4(a)) and r = 0.2 μg/L/d ( Fig. 4(b)). Fig. 4 indicates that a higher value of r causes a larger sudden increase in A and a slightly lower value of the critical threshold. A higher value of r implies that more phosphorus is recycled from the sediment into the water column and the recycled phosphorus is readily absorbed by the algae. The value of r = 0.3 μg/L/d for TH as obtained from the curve fitting of DO data is within the range of r reported from literature review as indicated in Table II. Since TH is a shallow lake (lake depth of 1.0 to 1.5 m), the phosphorus released from sediment is more available to the photic zone of the lake. Wind-induced resuspension of phosphorus may be another important factor contributing to the high eutrophication in the shallow TH. In the events of high wind, phosphorus in the sediment is resuspended into the water column. This has been identified as the main mechanism driving the phosphorus internal recycling in shallow temperate lakes by [33].   It should be noted that zero-flushing rate considered in TH causes a sudden large increase in A when the critical threshold (lp = 0.01595 μg/L/d) is exceeded. Hence, the effects of two flushing rates of h1 = 0.001 day -1 and h1 = 0.002 day -1 on critical thresholds are plotted in Fig. 5, which allows comparison with the increase in A for h1 = 0 day -1 . It should be noted that the actual practical flushing is limited, indeed zero, for TH, with no additional sources of water. Therefore, the phosphorus from external loading accumulates in the lake sediments and is readily available for the uptake by algae. From Fig. 5, the sudden increase in A is smaller with a higher flushing rate. With flushing, lake water laden with algae and phosphorus is constantly flushed out of the lake, causing a reduction in the concentration and residence time of algae and phosphorus.  Flushing is considered as one of the potential lake restoration methods for TH. In order to determine the flushing rate required to restore the lake to oligotrophic state (A < 2.6 μg/L), a bifurcation diagram of algae, A against flushing rate, h1 is plotted as displayed in Fig. 6. As expected, the algae concentration decreases when flushing rate increases. From Fig. 6, three equilibria (two stable and one unstable) exist when 0.018 d -1 < h1 <0.02575 d -1 . In this range of h1, the A will be attracted to either the stable steady state E4 (A > 10 μg/L, eutrophic) or stable steady state E6 (A < 10 μg/L, oligotrophic). If A is above E5 and below E4 (with 0.018 d -1 < h1 <0.02575 d -1 ), it approaches the stable steady state E4. An algae concentration lower than E5 (with 0.018 d -1 < h1 <0.02575 d -1 ) will decrease and approach the stable steady state E6. The algae concentration will never approach the unstable steady state E5. From Fig. 6, a flushing rate of h1 > 0.02575 d -1 would restore TH to mesotrophic state (A < 9.735 μg/L), while the oligotrophic state of TH can be achieved when h1 > 0.042 d -1 (A < 2.6 μg/L).

A. Flushing
During flushing process, lake water laden with algae and phosphorus is flushed out of the lake. Based on the bifurcation diagram in Fig. 6, a flushing rate of h1 > 0.042 d -1 would restore TH to an oligotrophic state (A < 2.6 μg/L). A flushing rate of 0.042 day -1 would require the lake replaces 4.2% of its volume per day, or a residence time of 24 days. This suggests that the entire volume of TH needs to be replaced in 24 days. For this purpose, an adequate supply of water is needed to sustain the application of flushing. It should be noted that the flushing rate of 0.042 day -1 is estimated based on an external phosphorus loading, lp = 0.3 μg/L/d. Any increase in lp in the future would require more water for flushing to achieve the same result. Hence, it is not feasible to restore TH with a flushing rate of 0.042 day -1 due to huge amount of water needed every day.

B. Water Hyacinth
Next, the hypothetical use of the invasive water hyacinth in restoring TH is discussed. In this simulation, water hyacinth is introduced in the lake on 4501 th day and the maximum water hyacinth growth rate gw = 0.11 day -1 is considered. The gw is maximum at temperature 25-28C [28,39]. Fig. 7(a) displays the graph of algae, A against time with and without harvesting of water hyacinth. When the water hyacinth is first introduced in the lake, phosphorus is consumed by the water hyacinth for growth and the phosphorus uptake by algae decreases. Consequently, A decreases until A = 250 μg/L as indicated in Fig. 7(a). However, when the water hyacinth reaches its maximum carrying capacity at 70 kg/m 2 , the phosphorus cannot be further absorbed by water hyacinth. As a result, the algae continue to grow and reach 300 μg/L again within 6 years.
In another scenario, optimal harvesting of water hyacinth is considered. Here, the optimal harvesting of water hyacinth refers to the value of water hyacinth harvesting rate, h that would result in the lowest value of A. Fig. 7(b) shows the graph of algae, A against water hyacinth harvesting rate, h. Based on Fig. 7(b), an optimal value of water hyacinth harvesting h = 0.0506 d -1 would produce the lowest value of A. Since water hyacinth absorb phosphorus to grow, harvesting water hyacinth would apportion for new growth. An application of water hyacinth with optimal harvesting rate would only reduce A from 300 μg/L to 120 μg/L after 9.6 years as indicated in Fig. 7(a), indicating a highly eutrophic state in TH.
By multiplying the surface area of TH (6070 m 2 ) with maximum carrying capacity of water hyacinth (70 kg m -2 ), this results in 424,900 kg standing crop wet weight of water hyacinth. The harvesting rate h = 0.0506 d -1 denotes that 5.06% of water hyacinth needs to be removed every day. By multiplying 0.0506 with 424,900 kg, this indicates that 21,499.94 kg (~21.5 tonne) of water hyacinth need to be removed every day for the next 9.6 years so that the lowest value of A of 120 μg/L could be achieved 9.6 years later.
Assuming that one person can harvest approximately 200 kg of water hyacinth per hour [31] and each person works for 8 hours, approximately 13 workers are needed to harvest the water hyacinth every day to reduce A to 120 μg/L. Since water hyacinth consists of 90% water, it is very heavy to transport water hyacinth for disposal or utilization purposes. This restoration method is deemed not cost-effective considering the costs of controlling water hyacinth such as labour cost, equipment cost and transportation cost, not to mention the ecological effect of the highly invasive water hyacinth. Based on this result, it is not feasible to use water hyacinth to restore TH due to several reasons. First, the water hyacinth is highly invasive due to its rapid growth rate and high adaptability to extreme conditions [18]. It might be difficult to control the water hyacinth population in TH and difficult to eradicate water hyacinth once it is established [40]. Besides, the dense mats formed by water hyacinth on the water surface would reduce light penetration to lake bottom and gas exchange on the water surface [41]. This would prevent the growth of submerged plants and would reduce the DO in the lake water. A low level of DO promotes the release of phosphorus from sediment into the water column and accelerates the eutrophication process. The floating mats of water hyacinth also decrease water currents and provide breeding grounds for vectors such as mosquitoes. Mosquito-borne diseases such as dengue and malaria might occur. A whole surface area of TH covered by floating mats of water hyacinth would expose TH to be the mosquito breeding grounds. In view of water hyacinth's mild effectiveness in improving the water quality of TH and the adverse environmental and ecological impacts, it is not sustainable to use water hyacinth as a lake restoration method for TH.

VI. CONCLUSION
A mathematical model consisting of algae, phosphorus, DO and BOD is formulated to identify regime shift thresholds for TH. With the observed DO data for TH, important parameter values such as external phosphorus loading rate, algae growth rate and phosphorus internal recycling rate are estimated through curve fitting. These fitted values are then used for bifurcation analysis, the results of which indicate that the current state of TH is irreversible. The regime shift threshold for TH is lp = 0.01595 μg/L/d. Any value of lp greater than 0.01595 μg/L/d would cause a sudden increase in algae concentration exceeding 100 μg/L. Two potential lake restoration methods, namely flushing and the hypothetical application of the invasive water hyacinth are discussed. Theoretically, a flushing rate greater than 0.042 day -1 is needed to restore TH to oligotrophic state. This means the entire volume of TH needs to be replaced every 24 days. Hence, it is not feasible to restore TH with a flushing rate of 0.042 day -1 due to the huge amount of water needed every day. An application of water hyacinth with optimal harvesting would reduce the algae concentration from 300 μg/l to only 120 μg/l after 9.6 years. The use of water hyacinth in restoring TH is not cost-effective considering the various costs involved for the next 9.6 years just to get the algae concentration down to 120 μg/l, which still indicates highly eutrophic condition. More importantly, introducing the highly invasive water hyacinth into TH would pose various adverse environmental and ecological impacts. Dredging the lake sediment appears to be the only option left.