Developing a simple dynamic simulation model of phosphorus in Lake
Ontario
Hamid Abdolabadi, sMojtabe: Ardestani®, Amin Sarang °, John C. Little*
* Ph.D Candidate in Envi duate Faculty of Envi t, University of Tehran.
’ Graduate Faculty of Envi , University of Tehran, Tehran, Iran,
© Dey of Civil and Envi 1 Engineering, Virginia Tech, Blacksburg, USA.
*Corresponding Author Tel: +98-915-188-1801; Fax: +98-21-44731988
E-mail address: h.abdolabadi@ut.ac.ir
Abstract
Food cycle creates a mechanism through which nutrients and other components of life are
ilable for living or i Improper fi ion of one cony of the cycle may lead to
emerging disruption in the life cycle of organisms. In this paper, the dynamic behavior of
phosphorus was studied by considering three state variables: productive organisms (organic
phosphorus), dead or is (organic phosphorus) and inorganic phosphorus. The system
dynamics model intelligibly provides information about changing state variables levels with
respect to all interactions and feedbacks. The model was applied to Lake Ontario and run for a
year with daily time steps. The model demonstrates acceptable performance in estimating the
variables’ concentration. Stratif ion and mixing dition have a significant effect on the
variables’ concentration during two periods (Day 100 to 158 and day 315 to 335) leading toa
decrease in the soluble reactive phosphorus concentration of 80% in the hypolimnion
(compared to the epilimnion) as well as an increase of phytoplankton concentration from 0.2 to
0.4 (mg/) in the epilimnion.
Keywords: Phosphorus cycle, System Dynamics, Phytoplankton, Lake Ontario.
I. Introduction
The food cycle is an ecological process in which nutrients change from one form to another.
The existence of the ecosystem depends on this cycle. Without it, productive organisms would
not be able to acquire the nutrients necessary for their survival (Foéllmi, 1996).
Several processes including solar radiation, prey-predator dynamics, and hydrological and
climatic parameters directly affect the food cycle. However, solar energy has the most impact on
the natural process of growth and death in nature. The sun supplies the required energy for the
plants (producers) and they absorb light to photosynthesize (Deaton & Winebrake, 2000).
During photosynthesis, the plants convert light, carbon dioxide, and mineral nutrients to
chemical energy so as to produce more complex compounds of carbon, which may then be used
as a fuel source for growth. The consumers eat these producers and break down the complex
compounds within them to obtain energy and nutrients. After the consumers die, they will be
decomposed, and returned to the ecosystem, in the form of a primary nutrient, which can again
be used by the producers.
The primary cycle of biomass in the environment is more complex. To simply describe the
dynamics of the phosphorus cycle, we should consider at least three storages including the
storage of productive organisms (eutilents could be found within them in their living form),
dead or i and d or i (bacteria turn the dead organism into primary
mineral nutrients) (Féllmi, 1996).
Studying the phosphorus cycle is of p ignij for ing water quality. Over
the past decades eutrophication has been one of ‘the most common problems affecting water
quality of lakes (Andersen et al., 2002; Chapra, 1997; Jong et al., 2002). Eutrophication is a
complex natural process occurring gradually due to nutrient enrichment in water bodies which
can irreversibly affect the ec Phosphorous has been recognized as the main It
factor for algal growth (Gilbert et al., 2010). In most cases, as the conc ion of phosphorous
declines, the growth rate will also decline due to lack of available nutrients. Therefore, studying
the phosphorous cycle in water bodies could be beneficial to clearly understand a simple food
web process.
Over the last years, making use of object-oriented models has become very common in analyzing
complex phenomena. Such models create a flexible and user-friendly framework to develop up-
scaled models for analyzing complex systems (Ahmad & Simonovic, 2000). System dynamics is
a feedback based method, which usually does not require advanced mathematical descriptions.
Some adi of such simulatic hods are appropriate under ding of a ph
high speed numerical calculation, modeling reliability, and the possibility of expanding and
easily changing model structure (Loucks et al., 1981; Simonovic & Fahmy, 1999). System
dynamics approach has been used for eutroy models, predicting long-term water quality
changes, analysis of the Quality Control policies in river basins, water allocation, and the
operation of multi-purpose dams )Mitra & Flynn, 2010; Geene, 1996; Vezjak et al., 1998). In
this study, the ation of phosphorus was simulated with daily time step in three forms
ie. productive or is (organic phosphorus), dead or i (organic phosphorus) and
inorganic phosphorus. We also used system dj ics to simulate the behavior of phosphorous
in lakes during the annual stratification and mixing cycles.
2. Material and Methods
2.1. Model description
7 he system dynamics approach is a powerful tool for developing object-oriented models to
over time that involve feedback effects (Sterman, 2000;
Elshorbagy & Ormsbee, 2006). In this study, system dynamics was used to simulate the
phosphorus cycle. It is a powerful and simple approach to model complex systems. This
approach works based on feedbacks among variables and regulates itself with respect to the past
behavior of the system (Forrester, 2007). In these models, an initial causal loop diagram oft he
desired problem is drawn and a graphical diagram of the model is then plotted using stocks,
flows, arrows, and converters (Ford, 1999).
The causal loop diagram of the phosphorus cycle model consists of reinforcing and balancing
loops. The combination of positive and negative feedback loops allows the system to reach steady
state. Figure I shows the causal loop diagrams of the main state variables in the system and
their feedback loops.
3. Theory: subsystems of the phosphorus cycle simulation model
As the phytoplankton concentration increases, the concentration of soluble reactive phosphorus
will decrease due to consumption. Increasing the phytoplankton concentration leads to greater
nutrient uptake. As a result, reducing nutrient concentration slows the growth of phytoplankton
(balancing loop). On the other hand, d posit of phytoplank increases the
concentration of organic comy ds. Additionally, if soluble phosphorus centration is
on the increase, the nutrient concentration will increase and lly the phytoplankton
concentration will rise (reinforcing loop). All the paramount parameters, their relationships,
and their impacts on the state variables were determined and the structure of the model was
developed for each state variable based on its control equations.
In this study, the soluble reactive phosphorus, soluble reactive phosphorus and
h lank were idered as the key state variables in the epilimnion and hypolimnion.
T he model was applied to Lake Ontario. To demonstrate the validation and accuracy of the
model, we used a range of factors and coefficients suggested by Chapra (1999). This model
S I temperature ch The stratification divides the lake into two layers (the
epilimnion and hypolimnion), which were idered to be pletely mixed vol
Turbulent diffusion connects the two layers and the inputs and outputs flow from the
epilimnion. Table 1 defines all symbols used in the following the equations.
+
Soluble reactive phosphorus Phytoplankton
B (organic)
(inorganic phosphorus)
Productive organism
(organic)
‘Dead organism
Fig. 1. The causal loop diagram of phosphorus cycle.
cycle model
Table 1. Values of all the parameters used in the p
Parameters Symbol Value | Parameters Symbol Value
Concentration of phytoplankton, Concentration of SRP (mgP m-
(mgChla m3) 4 3) p
Phytoplankton growth rate @ € r 2 A —
720°C (a) aa) Thermocline area (m?) , 185*10
Phytoplankton losses due to k a . V —
vebpiration and excretion (a) be 0.025 | Epilimnion volume (m') ‘ 254*10
ppreplankion settling velocity (m 0.2 | Hypolimnion volume (m?) V, 14*19"!
Optimal light level (ly €") I, 350 | Out flow (m? yr!) Q 212*10°
Temprature factor 0 1.066. | Thermocline diffusion: ¥,
Attenuation of growth due to light, Summer-stratified (cm? 1) 9.13
Attenuation of growth due to 4 cmcepnril 7 2B
phasphorsiis d, Winter-stratified (cm? s!)
Light extinction due to factors K a2 | Start of summer stratification iii
other than phytoplankton (m!) e ‘ @ ss
Phosphorous half-saturation k 3 Time of establish stratification, sg
constant (ugP L"!) 2 (a) tes
Photoperiod (sunlight fraction of Time | Onset of end of stratification
a J ; Loos 315
day) Series | (d) oe
Extinction coefficient k, End of stratification (d) t 20
Variable load (mg a!) Ww Epilimnion thickness (m) H, 20
Concentration of NSRP (mgP m3) C,.,, Hypolimnion thickness (m) H, 82
NSRP settling velocity (m a") Vn 0.2 _ | The stoichiometric coefficient
for the conversion of a 1
Decomposition rate for NSRP (dx 0.1 | Phosphorous to phytoplankton“
y) t " (mgP/mg Cht!)
Equations 1 and 2 are mass balances for the epilimnion and hypolimnion (Chapra, 1997). The
subscripts e and h designate epilimnion and hypolimnion for state variables of the phosphorus
cycle model.
Ve =W— Qe, +¥A le C)ES, o
it
y, fo =0.4 (0-0) £5, Q)
1 UAE“ O ES;
3.1. Lake volume
The volumes of the epilimnion and hypolimnion were calculated based on the thermal
stratification of the lake. In the model, a state variable was used to simulate the lake volume.
3.2. Phosphorus cycle
3.2.1. Phytoplankton
The growth of phytoplankton is a function of temperature, light, nutrients, and algae. The
concentration of algae was repres d by the ation of chlorophyll a (Asmala, 2011;
Flynn, 2010).
Equation 3 describes the model of phytoplankton growth. The concentration of phytoplank
in the hypolimnion was calculated based on the settling process from the epilimnion and the
effective diffusion between the two layers (Chapra, 1997).
de
ic.
V, 2 = Kae (Ts gy E) Vee Kine¥ Cue + A, Ca Cad) Va NC 3)
egg (Ts Cpe 1) = Ky 291-066" $, dp 4
The Michaelis-Menten equation was used to account for nutrient limitations, as follows:
¢
op an (3)
sp tsp
The following equation considers the effect of light limitation on the ph lank
growth rate.
4 _2.718f ap ala goth aeep ata hats o
“RA | OL “UT,
k, =k! +0.0088c, + 0.054c,°°
”
igure 2 shows the stock and flow dit for in us cycle model.
Fi 2 shows the stock and. diag phytoplank phosph Te del.
<@
2. . ez 2
Aap — oe Epean eye Epi Out Phyto ——krae
iciotien Epilimnion krae
<Time
hy Temperature», HyPoliminion
temperature in lake
Hy In Phyto Quah] Hy Ont Phyt ee
Fig. 2. Stock and flow diagram for phytoplankton.
3.2.3. Non-soluble reactive phosphorous (particulate phosphorus)
Phosphorous is classified into two groups of reactive and soluble reactive phosphorous.
This classification is due to the method of measuring the phosphorous concentration in the field.
In most cases, the phosphorous concentration is measured as reactive soluble phosphorous and
total phosphorous. The ation of soluble reactive phosphorous is equal to the
difference between these stocks. The particulate phosphorus increases through the
decomposition of the phytoplank and its ation is controlled through the
decomposition, settling, and decay processes. Eventually, the particulate phosphorus transforms
into reactive soluble phosphorous. Equations 9 and 10 describe the non-soluble reactive
phosphorous model.
Mnrsrpe
Ve Wage Oye eg K pV Cag KV Cog: + A (Cog Corse) eA Cap 7)
Ay
Vy a = Wan + ya K na Can Kn Cars +24 (Carpe Cnr) +¥,ACune Vp Cane 20)
3.2.3. Soluble reactive phosphorus
Soluble reactive phosphorus is one of the main limiting factors that affect the food chain.
Organic phosphorus is converted into soluble reactive phosphorus, which is used by
phytoplankton (Chapra, 1997). Figure 3 shows the stock and flow diagram for soluble and non-
soluble reactive phosphorus in the model.
dC,
V, a = Wore — Dope — K geV Cope + KreVCnsrpe + UA: (Cop ~ Corpo) dy)
Co,
Yi, ap = oh + Kei V Cop + UA, (Cope = Capt (12)
K te,
ioe
ee aaa
Epitimnion __Kre look up
a
o—_e Sf —
Hy Setting
velocity
Hypolimixion SRP
Hy In SRP eee By Out sRP
Fig. 3. Stock and flow diagram for phosphorus.
4. Study site
Lake Ontario is one of the five Great Lakes of North America and it is the 14th largest lake in
the world. It is bounded on the north and south by the Canadian ince of Ontario and
on the south by the American state of New York. The Great Lakes watershed is a region of high
biodiversity and Lake Ontario is important because of its diversity of birds, fish, reptiles,
amphibians, and plants. The lake's primary source is the Niagara River, which drains Lake
Erie, while the St Lawrence River serves as the outlet. The drainage basin covers 64,030 km?
and 49% of it is forested, 39% is agricul l, while the ining 13% is urban (Agency
USEPA, 1998). The lake has an important freshwater fishery, although it has been negatively
affected by water pollution (Christie, 1974). Figure 4 shows Lake Ontario.
o_ 30 60, Paeeasraae:
ee Kilometers
Fig. 4. Lake Ontario.
5. Calibration and validation
To analyze the model behavior and characterize the impact of the parameters on the model’s
output, a sensitivity analysis of the model was conducted. In this regard, in each model run, the
effect of each parameter’s ch was q ified, while k 1g the other parameters constant.
To determine the most important parameters, a relative ty factor was calculated using
quation 13 (Shirmoh di et al., 2006).
00, P
F. 13)
sa Sp *O (13)
In this equation, “O” is the model’s output, and “P” is the model input parameter. The model’s
parameters were calibrated manually making use of a trial and error method. Table 2 indicates
the model’s sensitivity analysis for an average ation of ph lank The
coefficient reveals that the factors most affecting the lake’s phytoplankton concentration were
the algae growth rate, phytoplankton decay rate, light extinction, and the de rate of
non- soluble reactive phosphorus (NSRP). The range of parameters changes on account of the
trial and error, and the closeness of the model’s average output and the observed data.
Table 2. The sensitivity coefficient of the phytoy ion in the lake.
Parameters Symbol Variation range FE
Phytoplankton growth rate Kga,20 0.8 — 1.6 0.623
Phytoplankton losses due to ka 0.02 — 0.04 0.572
respiration and excretion
Light extinction due to factors ki 0.1-.0.4 0.532
other than phytoplankton e
Decomposition rate for NSRP k, 0.05 — 0.2 0.323
After analyzing the sensitivity, the model’s error was measured using the Nash-Sutcliff
efficiency coefficient, the Pierson correlation coefficient, and the standard error. Nash-Sutcliff
efficiency is defined in equation 14:
Yo-¢,7
E,,=1-———_ d4)
YC, = Cone)?
Where C, is observed variable and Cn is simulated variable. Nash-Sutcliffe efficiencies can
range from —« to 1. An efficiency of 1 (E = 1) corresponds to a perfect match of estimated
outcomes to the observed data. An efficiency of 0 (E = 0) indicates that the model predictions
are as precise as the mean of the observed data (Nash & Sutcliff, 1970). According to Table 3
the appropriate accordance of the model’s output with the observed data is indicated by a value
of the Nash-Sutcliff coefficient (Ens) close to 1, the data correlation coefficient (R), and also a
low standard error for the model’s calibration and verification periods. The result indicates a
reliable model definition. Figure 5 shows observed data and simulation results for soluble
reactive phosphorus: in the epilimnion and hypolimnion of Lake Ontario. It can be seen that the
correl stil d and observed data for phosphorus is very close in both layers.
Table 3. The results of Nash-Sutcliff coeffici ion coefficient, and the standard error.
Layer .
Criterion ’ Se Ens
Epilimnion 0.966 15 0.84
Hypolimnion 0.86 1.07 0.68
18
16
Soluble reactive phosphorus (mg Picubic m)
oO
0 30 cc) 90 120 150 180 210 240 270 300 330 360
Time (Day)
——Fpilimnion soluble reactive phosphorus. wwe Hypolimnion soluble reactive phosphorus
> Epilimnion SRP (observed) © Hypolimnion SRP (observed)
Fig. 5. The observed and simulated concentrations of the soluble reactive phosphorus.
Specific tests can be utilized to assess the accuracy of a model, based on the system dynamics
approach, as described below.
5.1. The dimensional analysis test
The unit consistency of each state variable in the model was analyzed by the Units Check order.
If there is an inconsistency in units, a unit error will emerge. Although a dimensional error does
not affect the numerical calculation and the model’s outcomes, it can lead to misunderstanding
for complex systems. In this study, there is no unit error in the model.
5.2. Extreme conditions test
Every model must cope with extreme circumstances during the simulation period. For this
reason, the input load and initial concentration of SRP, NSRP, and phytoplankton load were
considered as 0. Accordingly, the phytopl Co ation decreases from the initial
amount of 1 to 0.4 mg Chla per cubic meter (Figure 6). In addition, the SRP and NSRP
concentrations increase during the simulation owing to d position of the phytoplank
Consequently, the model works accurately under extreme conditions.
12 1.2
Phytoplankton (mg chl a/cubie m)
Phosphorus (mg P/cubic m)
Time (Day)
—Phytoplankton = = soluble reactive ees N luble reactive NSRP)
Fig 6. The results of state variable concentration under extreme Conditions.
6. Results and Discussion
After developing the model structure for the epilimnion and hypolimnion of the lake based on
their equations, the concentration of soluble reactive phosphorous, non-soluble reactive
hosphorous, total phosphorous, and algae (phytoplankton) stocks were lyzed from January
to December. To assess the simulation results, we first analyze the behavior of the state variables
separately for each layer. Figure 7, shows the simulation results in the epilimnion.
, i
: 3
Bg &
3 et
Fy
g 4 g
z a
23 A
£ =
£2
1
0 0
0 30 60 90 120 150 180 210 240 270 300 330 360
Time (Day)
— Phytoplankton ~~ soluble reactive phosphorus(SRP)
Non-soluble reactive phosphorus(NSRP) = = ~ Total phosphorus(TP)
Fig 7. Simulation results for state variables in the epilimnion,
According to Fig. 7, the phytoplankton growth rate increases slightly during the first two months
of the year. This might be due to the gradual increase in light and the photoperiodic ratio. These
factors approximately neutralize and control the growth of phytoplankton.
Th he phytoplankton concentration increases abruptly due to the abundance of soluble reactive
phorous and a di rise in temperature and solar radiati The phytoplank
‘concentration reaches a peak on about Day 160 when the stratific n develops, and
afterwards, until the end of the thermal stratification, has a gradual decreasing trend. As pointed
out in the modeling section, the phytoplank decay rate p (K,), which includes
respiration, decomposition and excretion processes, is the most important factor controlling the
ion. This p increases when the temperature increases based
on the Michael-Manten equation. Therefore, the reason for the gradual decrease in the
h lankton concentration during the could be due to the increase of this parameter.
Eventually, the interaction of both feedback loops could be recognized as the main reason for
the maximum phytoplankton concentration at 9 (mg Chla/cubic m) on about Day 160. The
soluble reactive phosphorous is a primary nutrient. for algae growth. Its concentration declines
rapidly as the algae begins to grow. It hes the ion of 2 (mg P/cubic m)
on about Day 160 when the ph lank ation peaks. The trade-off between
phytoplankton growth and death causes SRP to gradually reduce in the epilimnion throughout
the summer up to about Day 315.
Since the level of soluble reactive phosphorous increases with phytoplankton decay, its
concentration increases slowly because of the slow growth of algae and the continued entry of
input loads. After Day 160, the soluble reactive phosphorous concentration decreases in
the epilimnion, alongside the reduction in the phytoplankton ation, and slowly
continues through Day 315.
Strong vertical mixing occurs during the late when the ] plank ‘ation
declines and the SRP and NSRP concentrations increase. Figure 8, shows the simulation results
for state variables in the hypolimnion.
2
Phosphorus (mg P/cubic m)
3
g
3
3
7
&
5
?
a
a4
Es
ar 30 60 90 ey nr rT es ee ne a ee
Time (Day)
— Phytoplankton ~~ soluble reactive phosphorus(SRP)
Non-soluble reactive phosphorus(NSRP) = = =Total phosphorus(TP)
Fig 8. The results of phosphorus si ion in the hypolit
The phytoplankton concentration reaches a maximum of 5.5 (mg Chla/cubic m) on about day
150 in the hypolimnion, and after that, its concentration signif ly decreases b of light
limitation (caused by the algal bloom in the epilimnion) and decomposition and excretion
processes. During the summer stratification, temperature remains almost unchanged in this
layer. As a result, the ph ation in the hy; is smaller than in the
epilimnion. Soluble reactive phosphorous increases in the hypolimnion due to the reduction of
phosphorous consumption, the settling of algae and soluble reactive phosph Ss from the
epilimnion, and the decomposition process. The other important point is the fluctuation in the
phosphorous concentration in the epilimnion and hypolimnion on Day 315. As the water is
getting cold in the upper layer in the late summer, it will be denser which causes destratification
of the lake.
The soluble reactive phosphorous concentration in the hypolimnion is higher than the
epilimnion owing to settling of the soluble reactive phosphorous from the epilimnion, and
phytoplankton decay in the hypolimnion. The soluble reactive phosphorous concentration
gradually decreases during the stratification period based upon the reduction in phytoplankton
, and lly on Day 315 it equals the soluble reactive phosphorous
concentration in the epilimnion at around 9 (mg P/cubic m).
The total phosphorous concentration in the epilimnion (accumulation of the soluble and non-
soluble reactive phosphorous) decreases during the thermal stratification period due to the
phytoplankton consumption (Fig. 7). Moreover, it declines in the hypolimnion in two time
periods (Fig. 8). The first of these is the time of establishing stratification in the lake between
Day 100 and Day158, and the second period is the time of turning over between Day 315 and
Day335.
7. Conclusion
In this study, the phosphorous cycle in a lake was simulated using a system dynamics approach.
The results show that the model has the capability to simulate the behavior of the phosphorous
cycle. In this model, the overall trends and fl i of the phytoplankton (productive
organisms), non-soluble reactive phosphorus (dead or i and soluble reactive
phosphorus (inorganic phosphorous) were studied. Despite the similar overall behavior of the
variables in the upper and lower layer, the results revealed that the peak concentration of the
phytoplankton was smaller in the hypolimnion. It was also shown that the concentration of state
variables were affected by the strati ion and the complete mixing
The model could be used to de ate the ii ions and feedbacks of iy of the
phosphorous cycle. It also can provide students, decision-makers and opera's with an
interactive learning envir to eval the effecti of s. Finally,
adding a more comprehensive carbon cycle and including the release of ‘phosphorous ‘from the
sediments would improve the model reliability, although this will increase the model complexity
and the need for more accurate data.
References
Ahmad, S., Simonovic, S.P., 2000. Dy i deling of flood policies. In
Proceedings of the 18th International Conference of The System Dynamics Society:
Sustainability in the Third Millennium, Bergen, Norway, 6-10.
Anderson, D.M., Glibert, P.M., Burkholder, J.M., 2002. Harmful algal blooms and
eutrophication: nutrient sources, c ition, and s. Estuaries 25, 704-726.
Asmala, E., Saikku, L., Vienonen, S., 201 1. Import-export balance of nitrogen and phosphorus
in food, fodder and fertilizers in the Baltic Sea drainage area. Science of the Total Environment
409(23), 4917-4922.
Chapra, S.D., 1997. Surface water quality modeling. Ith Edition, Mc Graw-Hill, Inc.
Deaton, M. L., Winebrake, J. J., 2000. Dy i deling of envir I Springer
Verlag.
Elshorbagy, A., Ormsbee, L., 2006. Object-oriented modeling approach to surface water quality
management. Environmental Modeling Software 21(5), 689-698.
Flynn, K.J., 2010. Ecological modelling in a sea of variable stoichi ry; dysfi ionality and
the legacy of Redfield and Monod. Progress in Oceanography 84, 52-1 65;
Féllmi, K., 1996. The phosphorus cycle, phosphog is and marine phosph ich der
Earth-Science Reviews, 40(1): 55-124.
Ford, A., 1999. Modeling the envi : an introduction to system d, i deling of
environmental systems, Island press.
Forrester, J.W., 2007. System dynamics a personal view of the first fifty years. System Dynamics
Review 23(2-3), 345-358.
Geene, B., 1996. Eutrophication of Droodzen reservoir in Iran. Disseratition for the degree of
doctor of science (tecnology), Wageningen.
Glibert, P.M., Allen, JL, Bouwman, L., Brown, C., Flynn, K.J., Lewitus, A. and Madden, C.
2010. Modeling of HABs and eutr i status, adi hall Journal of Marine
Systems, 83, 262-275.
Jonge, V. N. de, Elliott, M. , Orive, E., 2002. Causes, historical development, effects and future
hallenges of a envir | problem: eutr . Hydrobiologia 475/476, 1-19.
Loucks, D. P., Stedinger, J.R., Haith, D.A., 1981. Water resources systems planning and
analysis. Ith Ed., Prentice Hall, Englewood Cliffs, New York.
Mitra, A., Flynn, K.J., 2010. Modelling mixotrophy in harmful algal blooms: More or less the
sum of the parts? Journal of Marine Systems 83, 158-169.
Shirmohammadi, A., Chaubey, I., Harmel, R. D., Bosch, D. D., Muiioz-Carpena, R., Sexton, A.,
Arabi, M., Wolfe, M. L., Frankenberger, J., 2006. Uncertainty in TMDL models, Transactions
of the ASAE 49(4), 1033-1049.
Simonovic, S.P., Fahmy, H., 1999. A new modeling approach for water resources policy
analysis. Water Resources Reseatch 35 (1), 295-304.
Sterman, J. D., 2000. Busi ics; system thinking and modeling for a complex world.
Irwin McGraw-Hill Publication.
Vezjak. M., Savsek. T., Stuhler. E. A., 1998. System dynamics of eutrophication process in
Lakes. European Journal of Operational Research 109, 442-451.
