Simulating Nutrients and Phytoplankton Dynamics in Lakes: Model Development and Applications

: Due to eutrophication, many lakes require periodic management and restoration, which becomes unpredictable due to internal nutrient loading. To provide better lake management and restoration strategies, a deterministic, one-dimensional water quality model MINLAKE2020 was modiﬁed from daily MINLAKE2012 by incorporating chlorophyll-a, nutrients, and biochemical oxygen demand models into the regional year-around temperature and dissolved oxygen (DO) model. MINLAKE2020 was applied to six lakes (varying depth and trophic status) in Minnesota focusing on studying the internal nutrient dynamics. The average root-mean-square errors (RMSEs) of simulated water temperature and DO in six lakes are 1.51 ◦ C and 2.33 mg/L, respectively, when compared with proﬁle data over 2–4 years. The average RMSE of DO simulation decreased by 24.2% when compared to the MINLAKE2012 model. The internal nutrient dynamics was studied by analyzing time series of phosphorus, chlorophyll-a, and DO over several years and by perform-ing a sensitivity analysis of model parameters. A long-term simulation (20 years) of Lake Elmo shows that the simulated phosphorus release from sediment under the anoxic condition results in surface phosphorus increase, which matches with the observed trends. An average internal phosphorus loading increase of 92.3 kg/year increased the average daily phosphorus concentration by 0.0087 mg/L.


Introduction
Eutrophication has been a threat to waterbodies since the beginning of the twentieth century in industrialized countries [1][2][3][4]. A large proportion of the anthropogenic increase in nitrogen and phosphorus flux due to industrialization is delivered to ground or surface waters through direct runoff, human and animal wastes, and atmospheric deposition. Over time, excess nutrients are transported to waterbodies [5,6]. When a waterbody undergoes any human-influenced ecosystem changes such as nutrient loading, extreme weather events, or invasive organisms; algal species (cyanobacteria) can form dense overgrowths known as algal blooms. Since these blooms can produce toxins that are harmful to people and wildlife they are often referred to as harmful algal blooms (HAB). HABs cause undesirable changes in aquatic resources such as reduced water clarity, hypoxia, fish kills, loss of biodiversity, and an increase in nuisance species [7,8]. Oxygen is consumed by both living and dead algae which results in low oxygen concentration in lakes typically known as hypoxia (below 2-4 mg/L of dissolved oxygen (DO)). Eutrophication also has a detrimental effect on human health through increased exposure to cyanobacteria toxins [9,10], nitrites, and nitrates [7,8] in drinking water. Since most cities use surface water as the drinking water source, HABs can cause serious problems of off-flavor odor and taste (sometimes described as earthy or musty). In some cases, drinking water no longer remains safe to drink and complete remediation is needed. For example, the state of Ohio committed to In MINLAKE2020, the chlorophyll-a model was modified to overcome many of the limitations of previous MINLAKE models. The model can simulate up to three algal groups (diatoms, green algae, and blue-green algae) and the shift from diatoms early in the season to green and then blue-green algae during the summer. The algal groups are distinguished by different rates of photosynthesis, respiration, settling, zooplankton grazing, and different nutrient requirements. A schematic diagram of the phytoplankton cycle applicable to all algal groups is presented in Figure 1. Phytoplankton growth depends on the maximum growth rate of the algae (Gmax), half-saturation coefficients for nutrients, water temperature, solar irradiance, external nutrient concentrations, and the current Chla concentration. The maximum growth rate of algae varies for different classes of algae. The algal growth limitation by nutrients is modeled using a Michaelis-Menten equation [50]: where f(S) is dimensionless, S is the concentration of the nutrient (P, N, or Si) in water (mg/L), and KS is the half-saturation constant for the nutrient (mg/L). Algal growth dependence on water temperature is modeled by equations given by Lehman et al. [51]: T is the water temperature (°C). The maximum growth occurs at an optimal temperature, Topt (°C), and the growth rate decreases both above and below Topt. Tmin (°C) is a low density (#/m 3 ), ΔTg (<1 d) is time that zooplankton spends in a layer during the night to graze (d), V(IZ)/V(I) is the day depth (layer IZ) and the layer I volume ratio, and CF is the unit conversion from L to m 3 (=0.001).
Zooplankton grazing is only simulated to represent the dynamics of algae. Zooplankton grazing rates vary with different classes of algae. For example, zooplankton is more likely to feed on green algae than blue-green algae [23]. ) (10) where ΔTg(I) is the time of grazing for layer I (d), TD is the photoperiod (h) from sunrise to sunset, Chla(k, I) is the chlorophyll-a concentration of phytoplankton group k in layer I, and IZ is the day-depth layer, where DO ≥ 0.5 mg/L. Zooplankton grazing of phytoplankton (Chla) occurs during nocturnal migration at the day depth. The nocturnal grazing rate is calculated for each layer (I = 1, …, IZ) between the day depth and the surface using the volume day depth/layer ratio Vr = V(IZ)/V(I) [43].

Zooplankton Simulation
MINLAKE2020 includes a zooplankton model to simulate (1) Chla lost by zooplankton grazing and (2) DO consumed by zooplankton respiration. A single class of zooplankton is simulated in the lake environment each day. During the day, zooplankton retreats to deeper water seeking refuge from visual predators. They begin to rise to the surface at dusk while grazing and return to deeper layers at dawn (Figure 2). Zooplankton activity of these two periods is treated separately in the model. Zooplankton is assumed to have a constant reproduction rate and a time-varying predation rate for determining the zooplankton population ZP(t) as a function of time (t, day). The day depth, light level at the day depth, and predation on zooplankton are calculated as the first step in zooplankton simulation. The day depth of zooplankton is identified as the deepest layer in which the DO concentration is greater than 0.5 mg/L and therefore changes with time depending on DO vertical depth distribution. Below the day depth, the grazing is zero. In MINLAKE2020, dual effects of seasonal predation and light limitation are included to simulate biomanipulation techniques related to methods of increasing the zooplankton population [54]. The dominant zooplankton predators are assumed to be visual predators and zooplankton predation only occurs in the daytime. The light limitation assumes a linear variation of predation between two light levels [55]. phosphorus is released from the detrital mass through decay. Though diffusion of phosphorus occurs between layers, phosphorus is also transported indirectly between layers by phytoplankton and detritus settling. Figure 3 graphically represents/summarizes the SRP fate and transport modeled by MINLAKE2020. Phosphorus, accumulated from the detrital biomass, sediment release (zero-order kinetics), and respiration, are used by the phytoplankton, in the presence of sunlight, for growth. Algae need both nitrogen and phosphorus for growth. However, phosphorus is particularly important for algal growth as it is usually in short supply compared to other nutrients. If it is assumed that nitrogen is in abundant supply, phosphorus becomes the only limiting nutrient for algal growth (green and blue-green algae), which is modeled in the application of MINLAKE2020. Uptake depends on the maximum growth rate, light limitation, nutrient limitation, Chla concentration at that time, and the yield ratio of P to Chla as shown in Equation (15). The differential equation representing SRP fate and transport in a layer is given as Equation (15) In MINLAKE2020 phytoplankton growth is simulated by external nutrient limitation and the model does not allow for nitrogen fixation (storage of excess nitrogen for later use) as MINLAKE88 did [49]. MINLAKE2020 simulates the release of ammonia from dead algae indirectly through detrital decay. In the normal range of ammonia concentration, when KNH(n) is small, the Michaelis-Menten ratio for ammonia in Equation (16) will be close to 1.0, and in Equation (17) one minus the Michaelis-Menten ratio will be small; therefore, it has preferential uptake of ammonia over nitrate.

BOD Simulation
BOD is an important parameter in the DO, phosphorus, and nitrogen cycles. The microbial decay or decomposition of organic matter, which is detritus from the mortality of phytoplankton in MINLAKE2020, consumes oxygen, and therefore, the amount of organic matter is represented as BOD, an oxygen equivalent. However, DO directly affects biological decay processes and phosphorus release under anoxic conditions. In the regional DO and the model does not allow for nitrogen fixation (storage of excess nitrogen for later use) as MINLAKE88 did [49]. MINLAKE2020 simulates the release of ammonia from dead algae indirectly through detrital decay. In the normal range of ammonia concentration, when KNH(n) is small, the Michaelis-Menten ratio for ammonia in Equation (16) will be close to 1.0, and in Equation (17) one minus the Michaelis-Menten ratio will be small; therefore, it has preferential uptake of ammonia over nitrate.

BOD Simulation
BOD is an important parameter in the DO, phosphorus, and nitrogen cycles. The microbial decay or decomposition of organic matter, which is detritus from the mortality of phytoplankton in MINLAKE2020, consumes oxygen, and therefore, the amount of organic matter is represented as BOD, an oxygen equivalent. However, DO directly affects biological decay processes and phosphorus release under anoxic conditions. In the regional DO model [22], a constant rate for BOD was used for each simulation lake depending on the lake's trophic state. For the year-round lake water quality model [46], different constant rates of BOD for the open water seasons and winter ice cover periods were used over multiple years. However, BOD is an important parameter for nutrient cycles and DO cycle and is affected by mortality, organic decay, diffusion, and advection ( Figure 6) and simulated separately in MINLAKE2020. The differential equation representing BOD in a layer is given as Equation (18): diffusion advection

Interaction and Connection among Modeling Variables
This study is focused on the internal nutrient dynamics and its interaction/connection to the phytoplankton, BOD and DO dynamics in the lakes; therefore, the inflow and outflow sub-model was disabled. The MINLAKE2020 DO model is different from that of MINLAKE2012, where daily Chla and BOD were specified as model input based on lake trophic status. For MINLAKE2020, several modifications were made to simulate phosphorus, nitrogen, Chla, and BOD concentrations and zooplankton activity in the DO simulation. The photosynthetic oxygen production calculation was modified, and the nitrification and zooplankton respiration were added to the DO model while the rest of the terms were the same as in MINLAKE2012. In MINLAKE 2020, photosynthetic oxygen production is dependent on the simulated daily Chla concentration rather than the data driven Chla pattern and annual mean concentration used in MINLAKE2012.
The diffusion of DO occurs between the water layers in the metalimnion and hypolimnion. The spring and fall overturn periods for temperate lakes completely mix water columns and make all modeling variables distributed uniformly along with the depth. Since zooplankton usually spend the largest amount of their time in the day depth layer, zooplankton respiration is simulated at the day depth only. BOD removes oxygen from the water layer through the decay of detritus which is accumulated from phytoplankton mortality and removed by settling ( Figure 6). Nitrification removes oxygen from the water layer through the conversion of ammonia to nitrite then to nitrate. Nitrification is applied to the DO model only if nitrogen is simulated.
From Figures 1-7 and from Equations (9)- (19), one can see the complex interactions and connections among modeling variables. For example, when the settling velocity of the phytoplankton group is increased, it reduces Chla concentration and in turn, affects detritus (BOD) and phosphorus, but phosphorus goes back to affect phytoplankton growth. Additionally, DO concentrations in different water layers connect/integrate changes in water temperature, nutrients, and phytoplankton dynamics ( Figure 1). The sediment oxygen demand is an important factor for DO mass balance and results in the anoxic condition in the hypolimnion that leads to phosphorus release from sediment. Therefore, exploring/understanding the internal cycles/dynamics of nutrients/ phytoplankton and their interactions in different lakes helps us to gain insights of lake ecosystem and then develop appropriate restoration strategies.

MINLAKE Overview
MINLAKE models use the basic one-dimensional advection-diffusion equation to simulate the dynamic variations of state variables in horizontal layers of a lake.
where C is the concentration of a state variable, v is the vertical settling velocity of the particulate form of some of the state variables (v = 0 for dissolved); z is the vertical coordinate measured positively downward; K z is the vertical turbulent diffusion coefficient; and A horizontal area of the control volume. The one-dimensional vertical advectiondiffusion equation is solved using a series of layers characterized by depth from the water surface, thickness, layer volume, and horizontal areas. The MINLAKE model developed by Riley and Stefan [43] (called MINLAKE88 in this paper) was capable of simulating water temperature and DO profiles with three levels of complexity for phytoplankton. The first and simplest form is a single algal group model with productivity distributed in the mixed layer [43]. The second level of complexity has the algal growth term unique to each layer using Michaelis-Menten equation for a single algal group. The level three considers up to three groups of algae, phosphorus, nitrogen, BOD for lakes during the open water season. MINLAKE88 is not widely used because of its inability to perform multiple-year simulations. Due to a lack of available data, the nutrient model could not be developed further. Gu and Stefan [21] included an ice-cover period simulation in MINLAKE88 to model snow thickness, ice cover thickness, and water temperature but no other state variables. In 1991, Hondzo and Stefan introduced a more general water temperature simulation model for MINLAKE which can be applied to a wide variety of lakes and regions [44]. An important modification of MINLAKE was accomplished in 1994 when Fang [45] developed the regional dissolved oxygen model and combined it with MINLAKE to study the impact of global climate warming on lake water quality and fish habitat in Minnesota lakes. The regional lake model was developed to model different types (categorizing by stratification strength and eutrophication) of lakes by maximum depths (shallow, medium-depth, and deep), surface area (small, medium area, and large), and trophic status (eutrophic, mesotrophic, and oligotrophic). The regional DO model is a simplified version of the DO model that did not simulate daily nutrients and Chla concentrations but used the annual mean Chlorophyll-a concentration with seasonal variation patterns [22] to specify daily Chla for DO simulation. Separate sub-models for winter conditions were developed and integrated with MINLAKE96 to simulate water temperature and DO year-round over many years as long as daily weather data are available [46]. The MINLAKE96 model was further modified and refined as it was used in a study in 2010, to simulate water quality conditions in cisco lakes, which are typically deep mesotrophic or oligotrophic lakes [47]. The MINLAKE96 model was further modified by West-Mack to simulate phosphorus, nitrogen, and chlorophyll-a simulation in 1998 [23,48] and the updated model is called MINLAKE98 here. The mass balance equations for chlorophyll-a, phosphorus, nitrogen, and DO were modified from those of MINLAKE88. The model was tested on three lakes and produced satisfactory results but due to lack of observed data and inability of the model to run for multiple years, the nutrient model was not further developed. The source code of MINLAKE98 was also lost and no longer available for further improvement. MINLAKE96 was further modified in 2018 to calculate the hourly water temperature and DO using hourly weather conditions [24,49]. MINLAKE2012 is the most recent version of the daily MINLAKE model and can simulate year-round water temperature and DO in various lakes of different regions, using an Excel spreadsheet as the user interface. For all MINLAKE model variants, water temperature is simulated first by solving the following heat transport equation.
where T w (z, t) is the water temperature in • C, which is a function of depth (z in m) and time (t in d); A(z) (m 2 ) is the horizontal area for each layer of water as a function of the depth; K z (m 2 /day) is the vertical turbulent heat diffusion coefficient which is a function of depth and time; ρC p (J/m 3 -• C) represents the heat capacity of water per unit volume; H w (J/m 3 -day) is the heat source and/or sink term per unit volume of water. Determination of the turbulent diffusion coefficient is discussed in detail by Fang [45]. In the regional daily MINLAKE model, the vertical heat diffusion coefficient K z for epilimnion and hypolimnion is calculated using the following equation: where K z is the vertical diffusion coefficient in cm 2 /s (1 cm 2 /s = 8.64 m 2 /day = 0.36 m 2 /h), A s is the surface area of the lake (km 2 ) and N 2 is the Brunt-Vaisala stability frequency of the stratification (s −2 ). In the epilimnion, N 2 was set at a minimum value of 0.000075 [42]. Equations (1) and (2) are solved numerically using an implicit finite difference scheme and a Gaussian elimination method with time steps of one day.

Phytoplankton Simulation
In MINLAKE2020, the chlorophyll-a model was modified to overcome many of the limitations of previous MINLAKE models. The model can simulate up to three algal groups (diatoms, green algae, and blue-green algae) and the shift from diatoms early in the season to green and then blue-green algae during the summer. The algal groups are distinguished by different rates of photosynthesis, respiration, settling, zooplankton grazing, and different nutrient requirements. A schematic diagram of the phytoplankton cycle applicable to all algal groups is presented in Figure 1.
Phytoplankton growth depends on the maximum growth rate of the algae (G max ), half-saturation coefficients for nutrients, water temperature, solar irradiance, external nutrient concentrations, and the current Chla concentration. The maximum growth rate of algae varies for different classes of algae. The algal growth limitation by nutrients is modeled using a Michaelis-Menten equation [50]: where f (S) is dimensionless, S is the concentration of the nutrient (P, N, or Si) in water (mg/L), and K S is the half-saturation constant for the nutrient (mg/L). Algal growth dependence on water temperature is modeled by equations given by Lehman et al. [51]: T is the water temperature ( • C). The maximum growth occurs at an optimal temperature, T opt ( • C), and the growth rate decreases both above and below T opt . T min ( • C) is a low temperature at which phytoplankton growth is reduced to 90% from the optimum. T max ( • C) is the high temperature at which growth is reduced by 90%. Phytoplankton growth is usually limited by available light, which is a function of the depth and the light attenuation coefficients due to water (K w ) and algae (K c ). The Haldane equation [52] is used in MINLAKE2020 to calculate the light limitation for algal growth, which was also used in MINLAKE88, MINLAKE98, and the regional DO model: where f (L) is the light limitation coefficient (dimensionless), I(z) is the photosynthetically active radiation (PAR) as a function of depth, TD is the photoperiod (h) as a function of Julian day, RAD is daily solar radiation, and K 1 and K 2 are the light limitation and inhibition coefficients [53], respectively. The units for I(z), K 1 , and K 2 are in µE/m 2 /s. In MINLAKE88 and MINLAKE98, light limitation and inhibition coefficients were specified by the user. In MINLAKE2020, light limitation and inhibition coefficients are calculated using the same equations as in the regional DO model [23]. Phytoplankton populations are removed from a water column by four processes: respiration, mortality, settling, and zooplankton grazing ( Figure 1). Each phytoplankton population is assigned a fixed or calibrated respiration rate, mortality rate, and settling rate. Respiration affects the available phosphorus and DO immediately whereas mortality contributes with a time lag through detrital decay. In MINLAKE2020, a single class of zooplankton is simulated. To simulate the Chla lost by grazing of zooplankton, the zooplankton population (ZP(t), #/m 3 ) is simulated in a separate subroutine on daily basis. Grazing is assumed to take place in the evening when zooplankton rises to the upper layers and is dependent on the temperature and Chla concentration. The Michaelis-Menten equation is used to simulate the effect of Chla concentration on grazing. It is assumed that no grazing occurs below a threshold Chla concentration (Chla min in Equation (9)): rate. Respiration affects the available phosphorus and DO immediately whereas mortality contributes with a time lag through detrital decay. In MINLAKE2020, a single class of zooplankton is simulated. To simulate the Chla lost by grazing of zooplankton, the zooplankton population (ZP(t), #/m 3 ) is simulated in a separate subroutine on daily basis. Grazing is assumed to take place in the evening when zooplankton rises to the upper layers and is dependent on the temperature and Chla concentration. The Michaelis-Menten equation is used to simulate the effect of Chla concentration on grazing. It is assumed that no grazing occurs below a threshold Chla concentration (Chlamin in Equation (9)): where Chla is the chlorophyll-a concentration (mg/L), vc is the phytoplankton settling velocity (m/d), Km, Kr, θm, and θr are mortality and respiration rates (d −1 ) and corresponding temperature adjustment coefficients (dimensionless), respectively. Gmax is the maximum growth rate of phytoplankton (d −1 ), GRmax is maximum grazing rate (mg Chla/individual zooplankton per d), Chlamin is minimum chlorophyll concentration for grazing to occur grazing by zooplankton growth diffusion mortality settling respiration (9) where Chla is the chlorophyll-a concentration (mg/L), v c is the phytoplankton settling velocity (m/d), K m , K r , θ m , and θ r are mortality and respiration rates (d −1 ) and corresponding temperature adjustment coefficients (dimensionless), respectively. G max is the maximum growth rate of phytoplankton (d −1 ), GR max is maximum grazing rate (mg Chla/individual zooplankton per d), Chla min is minimum chlorophyll concentration for grazing to occur (mg/L), K gchla is half-saturation constant for grazing (mg/L), K P and K N are the half-saturation constants for phosphorus and nitrogen (mg/L), respectively; P and N are the available concentration of phosphorus and nitrogen (mg/L), respectively; ZP is the zooplankton density (#/m 3 ), ∆T g (<1 d) is time that zooplankton spends in a layer during the night to graze (d), V(IZ)/V(I) is the day depth (layer IZ) and the layer I volume ratio, and CF is the unit conversion from L to m 3 (=0.001). Zooplankton grazing is only simulated to represent the dynamics of algae. Zooplankton grazing rates vary with different classes of algae. For example, zooplankton is more likely to feed on green algae than blue-green algae [23].
where ∆T g (I) is the time of grazing for layer I (d), TD is the photoperiod (h) from sunrise to sunset, Chla(k, I) is the chlorophyll-a concentration of phytoplankton group k in layer I, and IZ is the day-depth layer, where DO ≥ 0.5 mg/L. Zooplankton grazing of phytoplankton (Chla) occurs during nocturnal migration at the day depth. The nocturnal grazing rate is calculated for each layer (I = 1, . . . , IZ) between the day depth and the surface using the volume day depth/layer ratio V r = V(IZ)/V(I) [43].

Zooplankton Simulation
MINLAKE2020 includes a zooplankton model to simulate (1) Chla lost by zooplankton grazing and (2) DO consumed by zooplankton respiration. A single class of zooplankton is simulated in the lake environment each day. During the day, zooplankton retreats to deeper water seeking refuge from visual predators. They begin to rise to the surface at dusk while grazing and return to deeper layers at dawn (Figure 2). Zooplankton activity of these two periods is treated separately in the model.
Zooplankton is assumed to have a constant reproduction rate and a time-varying predation rate for determining the zooplankton population ZP(t) as a function of time (t, day). The day depth, light level at the day depth, and predation on zooplankton are calculated as the first step in zooplankton simulation. The day depth of zooplankton is identified as the deepest layer in which the DO concentration is greater than 0.5 mg/L and therefore changes with time depending on DO vertical depth distribution. Below the day depth, the grazing is zero. In MINLAKE2020, dual effects of seasonal predation and light limitation are included to simulate biomanipulation techniques related to methods of increasing the zooplankton population [54]. The dominant zooplankton predators are assumed to be visual predators and zooplankton predation only occurs in the daytime. The light limitation assumes a linear variation of predation between two light levels [55].
where PD d is the daytime predation rate (d −1 ) and P d is the daily predation rate (d −1 ) calculated using Equation (12). XI is the light intensity at the day depth (µE/m 2 /s), XI min is the light intensity at which no predation occurs (µE/m 2 /s) and XI max is light intensity above which predation is not light inhibited (µE/m 2 /s). When the light intensity XI is less than XI min or larger than XI max , the ratio in Equation (11) is reset to zero or one, respectively. Both daytime and nocturnal predation are calculated in MINLAKE2020. Daytime predation combines light limitation with a time-varying maximum predation rate. A linear function is used to calculate the time-varying daytime predation rate P d given in Equation (12) when Julian day DY is between DY min and DY max .
where P min and P max are minimum and maximum predation rates (d −1 ), respectively. DY min is Julian day of the last day of minimum predation rate and DY max is Julian day of the beginning of maximum predation rate: P d = P min when DY < DY min and P d = P max when DY > DY max . For example, for Elmo Lake, West and Stefan [23] set P min , P max , DY min , and DY max as 0.05 d −1 , 0.7 d −1 , 110 (20 April), and 140 (20 May); it means P min occur on and before 20 April and P max occur on and after 20 May. Zooplankton density in the daytime is determined using Equation (13) including first-order reproduction and daytime predation: where Repro is the reproduction rate (dimensionless) and ZP(t − 1) is the zooplankton density in the previous day. Nocturnal predation occurs during nocturnal migration at the day depth. The nocturnal predation rate is calculated for each layer between the day depth and the surface.
where PD(t, I) is the nocturnal predation rate in layer I during migration (#/m 3 ) and PD n is the nocturnal predation rate (d −1 ) as a constant input parameter. ZP(t) is the zooplankton population in the day depth layer calculated using Equation (13), and ZP min is the minimum concentration of zooplankton for predation to occur (#/m 3 ). V(IZ) and V(I) are the volumes of the day depth layer and layer I (m 3 ), respectively. ∆T g (I) is the time that zooplankton spent in layer I during the night (d). The daytime ZP(t) minus PD(t, I) gives the zooplankton population for the next day. The second part of zooplankton simulation is the simulation of phytoplankton grazing by zooplankton which begins with a vertical rise in the evening. Temperature and Chla concentration affect grazing in the water layer (see Equation (9)). A Michaelis-Menten ratio is used to express the effect of Chla concentration on grazing with a refugium effect for no grazing below a threshold Chla concentration (Chla min in Equation (9)).

Phosphorus Simulation
In most cases, phosphorus is known to be the primary nutrient controlling the trophic state of lakes in the Upper Midwest USA and Canada [56]. Phytoplankton can only use the soluble reactive phosphorus (SRP) which is composed of orthophosphate and polyphosphate ions. The model only simulates the readily accessible phosphorus and indirectly models organic phosphorus as detritus. Phytoplankton growth removes SRP from the water. Respiration releases phosphorus into the water column. Mortality does not directly release phosphorus to the water column but contributes to the detrital mass (BOD); phosphorus is released from the detrital mass through decay. Though diffusion of phosphorus occurs between layers, phosphorus is also transported indirectly between layers by phytoplankton and detritus settling. Figure 3 graphically represents/summarizes the SRP fate and transport modeled by MINLAKE2020.
Phosphorus, accumulated from the detrital biomass, sediment release (zero-order kinetics), and respiration, are used by the phytoplankton, in the presence of sunlight, for growth. Algae need both nitrogen and phosphorus for growth. However, phosphorus is particularly important for algal growth as it is usually in short supply compared to other nutrients. If it is assumed that nitrogen is in abundant supply, phosphorus becomes the only limiting nutrient for algal growth (green and blue-green algae), which is modeled in the application of MINLAKE2020. Uptake depends on the maximum growth rate, light limitation, nutrient limitation, Chla concentration at that time, and the yield ratio of P to Chla as shown in Equation (15).
The differential equation representing SRP fate and transport in a layer is given as Equation (15): phosphorus is released from the detrital mass through decay. Though diffusion of phosphorus occurs between layers, phosphorus is also transported indirectly between layers by phytoplankton and detritus settling. Figure 3 graphically represents/summarizes the SRP fate and transport modeled by MINLAKE2020.
Phosphorus, accumulated from the detrital biomass, sediment release (zero-order kinetics), and respiration, are used by the phytoplankton, in the presence of sunlight, for growth. Algae need both nitrogen and phosphorus for growth. However, phosphorus is particularly important for algal growth as it is usually in short supply compared to other nutrients. If it is assumed that nitrogen is in abundant supply, phosphorus becomes the only limiting nutrient for algal growth (green and blue-green algae), which is modeled in the application of MINLAKE2020. Uptake depends on the maximum growth rate, light limitation, nutrient limitation, Chla concentration at that time, and the yield ratio of P to Chla as shown in Equation (15). The differential equation representing SRP fate and transport in a layer is given as Equation (15): where KBOD and θr are decay rate (d −1 ) and corresponding temperature adjustment coefficient, respectively; SP is the rate of phosphorus released at the water-sediment interface (g P/m 2 /d) at anoxic condition and it is calibrated against available phosphorus/Chla/DO detrital decay respiration sediment release diffusion growth of phytoplankton (15) where K BOD and θ r are decay rate (d −1 ) and corresponding temperature adjustment coefficient, respectively; S P is the rate of phosphorus released at the water-sediment interface (g P/m 2 /d) at anoxic condition and it is calibrated against available phosphorus/Chla/DO profiles, Y PChla is mass yield ratio of phosphorus to chlorophyll, and Y PBOD is mass yield ratio of phosphorus to BOD. A phosphorus/chlorophyll yield coefficient (Y PChla ) is used to determine the amount of phosphorus consumed during photosynthesis as well as the amount of phosphorus released during algal respiration. In MINLAKE98, the value of Y PChla was derived from the mass yield coefficient of phosphorus to BOD divided by the mass yield coefficient of chlorophyll to BOD [23]. This value is 1.1 mg P/mg Chla for Y PChla , which was assumed to be constant in MINLAKE98, and is close to that presented by Thomann and Mueller [57] of 1.0 µg P/µg Chla. In MINLAKE2020, Y Pchla is set to 1.1 mg P/mg Chla for all simulated lakes. However, the phytoplankton biomass does not depend on phosphorus solely, it also depends on the nitrogen concentration.
There are three source terms for phosphorus: detrital decay (death of phytoplankton), sediment release, and phytoplankton respiration. For many lakes which have a history of progressive eutrophication, the lake sediments have now become the primary source of phosphorus to the water. If the sediment-water interface is anoxic, phosphate ions go to the water at an increased rate, depending upon the concentration difference between porewaters and the overlying water [58]. MINLAKE2020 simulates SRP release back to water when the DO concentration becomes zero.
The daily zero-order internal phosphorus release (same as Equation (15)) was included in MINLAKE88 [20] and MINLAKE98 [46] before. CE-QUAL-W2 [59] includes the zeroorder and first-order phosphorus release from sediment. EFDC model [60] has a governing equation for total phosphate including a sediment-water exchange flux of phosphate (g P/m 2 /d) for the control volume at the bottom. ELCOM-CAEDYM [29] simulates phosphorus release from sediment as a function of temperature, DO and pH. There are some empirical models that quantified long-term internal phosphorus release for the whole lake (not for mass balance equation for a water layer or control volume), for example, Nurnberg [61] determined the internal phosphorus release per year after analyzing the data from various lakes. Stigebrandt et al. [62] and Stigebrandt and Andersson [63]

Nitrogen Simulation
When phosphorus is in excess (e.g., P >> K p ), nitrogen can become the limiting nutrient for algae growth, which is not a usual scenario in many lakes. Nitrogen is available in two forms (ammonium NH 4 , and nitrite plus nitrate, represented as NO 2-3 in this paper). MINLAKE2020 models both ammonia and NO 2-3 separately. Schematic diagrams of the NH 4 and NO 2-3 sub models are given in Figures 4 and 5, respectively. Nitrogen N in Equations (9), (15), and (19) is the sum of NH 4 and NO 2-3 concentrations (in mg N/L).
Diffusion of ammonia and NO 2-3 occurs between layers. Nitrification is the biological oxidation process (assuming DO) of ammonia to nitrite followed by the faster oxidation of the nitrite to nitrate so that nitrite and nitrate are often modeled as one variable. The uptake of ammonia and nitrate due to phytoplankton growth is calculated as a zero-order sink term with a preference in the uptake of ammonium over the uptake of nitrate. It directly links to the growth of phytoplankton. Respiration of phytoplankton releases ammonia (Equation (16)) but not nitrate. Mortality of phytoplankton does not directly release ammonia or nitrate to the water but contributes to the detrital mass (BOD), and then ammonia is released through the decay of the detrital mass, as shown in Equation (16). Therefore, ammonia is also transported indirectly between layers by phytoplankton (Equation (9)) and detritus settling (Equation (18)). Ammonia and nitrite releases from the sediment are also modeled when DO is low (depending on the DO concentration of the overlying water). The model does not consider the atmospheric deposition of ammonia or NO 2-3 and denitrification (the reduction of nitrate to nitrogen gas, N 2 ). monia is released through the decay of the detrital mass, as shown in Equation (16). Therefore, ammonia is also transported indirectly between layers by phytoplankton (Equation (9)) and detritus settling (Equation (18)). Ammonia and nitrite releases from the sediment are also modeled when DO is low (depending on the DO concentration of the overlying water). The model does not consider the atmospheric deposition of ammonia or NO2-3 and denitrification (the reduction of nitrate to nitrogen gas, N2).
Equations (16) and (17) provide the governing equations to model NH4 and NO2-3 in each water layer. KNI and θNI are the first-order nitrification rate (d −1 ) and the temperature adjustment coefficient, respectively; KNH(n) and KTN(n) are the half-saturation constants for preferential uptake of ammonia over nitrate and nitrogen uptake for each algae class, respectively; SNH, SNO, θNH, and θNO are ammonia and nitrite release rates from sediment (g N/m 2 /d) and corresponding temperature adjustment coefficients, respectively; YNHBOD, YN-  (16) monia is released through the decay of the detrital mass, as shown in Equation (16). There-fore, ammonia is also transported indirectly between layers by phytoplankton (Equation (9)) and detritus settling (Equation (18)). Ammonia and nitrite releases from the sediment are also modeled when DO is low (depending on the DO concentration of the overlying water). The model does not consider the atmospheric deposition of ammonia or NO2-3 and denitrification (the reduction of nitrate to nitrogen gas, N2).
Equations (16) and (17) provide the governing equations to model NH4 and NO2-3 in each water layer. KNI and θNI are the first-order nitrification rate (d −1 ) and the temperature adjustment coefficient, respectively; KNH(n) and KTN(n) are the half-saturation constants for preferential uptake of ammonia over nitrate and nitrogen uptake for each algae class, respectively; SNH, SNO, θNH, and θNO are ammonia and nitrite release rates from sediment (g N/m 2 /d) and corresponding temperature adjustment coefficients, respectively; YNHBOD, YN-HChla, and YNOChla are the mass yield ratios of ammonia to BOD, ammonia to chlorophyll, Equations (16) and (17) provide the governing equations to model NH 4 and NO 2-3 in each water layer. K NI and θ NI are the first-order nitrification rate (d −1 ) and the temperature adjustment coefficient, respectively; K NH(n) and K TN(n) are the half-saturation constants for preferential uptake of ammonia over nitrate and nitrogen uptake for each algae class, respectively; S NH , S NO , θ NH , and θ NO are ammonia and nitrite release rates from sediment (g N/m 2 /d) and corresponding temperature adjustment coefficients, respectively; Y NHBOD , Y NHChla , and Y NOChla are the mass yield ratios of ammonia to BOD, ammonia to chlorophyll, and nitrate to chlorophyll, respectively; and K BOD and θ BOD are the first-order BOD or detritus decay rate (d −1 ) and the temperature adjustment coefficient, respectively.
In MINLAKE2020 phytoplankton growth is simulated by external nutrient limitation and the model does not allow for nitrogen fixation (storage of excess nitrogen for later use) as MINLAKE88 did [49]. MINLAKE2020 simulates the release of ammonia from dead algae indirectly through detrital decay. In the normal range of ammonia concentration, when K NH(n) is small, the Michaelis-Menten ratio for ammonia in Equation (16) will be close to 1.0, and in Equation (17) one minus the Michaelis-Menten ratio will be small; therefore, it has preferential uptake of ammonia over nitrate.

BOD Simulation
BOD is an important parameter in the DO, phosphorus, and nitrogen cycles. The microbial decay or decomposition of organic matter, which is detritus from the mortality of phytoplankton in MINLAKE2020, consumes oxygen, and therefore, the amount of organic matter is represented as BOD, an oxygen equivalent. However, DO directly affects biological decay processes and phosphorus release under anoxic conditions. In the regional DO model [22], a constant rate for BOD was used for each simulation lake depending on the lake's trophic state. For the year-round lake water quality model [46], different constant rates of BOD for the open water seasons and winter ice cover periods were used over multiple years. However, BOD is an important parameter for nutrient cycles and DO cycle and is affected by mortality, organic decay, diffusion, and advection ( Figure 6) and simulated separately in MINLAKE2020.
The differential equation representing BOD in a layer is given as Equation (18): model [22], a constant rate for BOD was used for each simulation lake depending on the lake's trophic state. For the year-round lake water quality model [46], different constant rates of BOD for the open water seasons and winter ice cover periods were used over multiple years. However, BOD is an important parameter for nutrient cycles and DO cycle and is affected by mortality, organic decay, diffusion, and advection ( Figure 6) and simulated separately in MINLAKE2020. The differential equation representing BOD in a layer is given as Equation (18): where vB and YCHBOD represent detritus settling velocity (m/d) and the mass yield ratio of Chla to BOD, respectively. In the model, BOD is increased from two sources ( Figure 6). First, detritus travels to adjacent layers via diffusion. Secondly, the mortality of the phytoplankton adds to the detritus concentration. BOD has two sink terms ( Figure 6): advection (settling in the vertical direction) and organic decay. Detritus falls from the concerning layers and goes to another layer or the sediment. The microbial decay of organic matter is a function of the detrital mass expressed in oxygen equivalents (BOD). The mortality of cells and a fraction of the grazed phytoplankton are converted from Chla concentrations to oxygen equivalents using the constant carbon/Chla ratio and stoichiometric relationships. The result is a one-to-one correspondence between detrital decay and the utilization of oxygen. This cycle is very important for nutrient calculation as it directly adds to the nutrient load through detrital decay.
diffusion advection phytoplankton mortality detrital decay (18) where v B and Y CHBOD represent detritus settling velocity (m/d) and the mass yield ratio of Chla to BOD, respectively. In the model, BOD is increased from two sources ( Figure 6). First, detritus travels to adjacent layers via diffusion. Secondly, the mortality of the phytoplankton adds to the detritus concentration. BOD has two sink terms ( Figure 6): advection (settling in the vertical direction) and organic decay. Detritus falls from the concerning layers and goes to another layer or the sediment. The microbial decay of organic matter is a function of the detrital mass expressed in oxygen equivalents (BOD). The mortality of cells and a fraction of the grazed phytoplankton are converted from Chla concentrations to oxygen equivalents using the constant carbon/Chla ratio and stoichiometric relationships. The result is a one-to-one correspondence between detrital decay and the utilization of oxygen. This cycle is very important for nutrient calculation as it directly adds to the nutrient load through detrital decay.

DO Simulation
DO is one of the vital parameters of lake water quality simulation. Aquatic organisms and fish depend on the availability of DO in the waterbody [64]. A schematic diagram of the processes contributing to the DO concentration is given in Figure 7. It shows that DO is added to a water layer through diffusion and photosynthesis; and is removed by respiration of algae and zooplankton, detrital decay (BOD), sediment oxygen demand (SOD), and nitrification. The surface reaeration can add or remove DO depending on whether surface DO is less or greater than saturated DO (a function of surface temperature and lake elevation). Phytoplankton (modeled as Chla) growth can add DO to the water layer through photosynthesis to the point where water could be supersaturated with DO in some cases. These dynamic processes can happen over time scales of less than one day (the time step of the MINLAKE2020 simulation). Therefore, the model outputs DO profiles as an integration of different physical (e.g., mixing), chemical, and biological processes over the day. DO removal from the water layer through phytoplankton respiration is simulated to occur at a constant rate throughout the day while photosynthesis occurs only during the hours with solar radiation. In MINLAKE2020, the adjustments for low DO levels on SOD, BOD, and algal respiration follow Edwards and Owens's [65] formula. SOD is calculated for each layer, and it is treated as a sink term in the one-dimensional (vertical) transport equation [22,57]. Oxygen uptake of the sediment depends on the area and composition of bottom materials in contact with the water [22,66].
The differential equation representing DO dynamics in a layer is given as Equation (19): surface DO is less or greater than saturated DO (a function of surface temperature and lake elevation). Phytoplankton (modeled as Chla) growth can add DO to the water layer through photosynthesis to the point where water could be supersaturated with DO in some cases. These dynamic processes can happen over time scales of less than one day (the time step of the MINLAKE2020 simulation). Therefore, the model outputs DO profiles as an integration of different physical (e.g., mixing), chemical, and biological processes over the day. DO removal from the water layer through phytoplankton respiration is simulated to occur at a constant rate throughout the day while photosynthesis occurs only during the hours with solar radiation. In MINLAKE2020, the adjustments for low DO levels on SOD, BOD, and algal respiration follow Edwards and Owens's [65] formula. SOD is calculated for each layer, and it is treated as a sink term in the one-dimensional (vertical) transport equation [22,57]. Oxygen uptake of the sediment depends on the area and composition of bottom materials in contact with the water [22,66]. The differential equation representing DO dynamics in a layer is given as Equation (19): where Kzr and θzr are the zooplankton respiration rate (d −1 ) and the temperature adjustment coefficient, respectively; KNH(n) and KTN(n) are the half-saturation constant for preferential uptake of ammonia over nitrate and for nitrogen uptake for each algae class, respectively; Sb and θSOD are sediment oxygen demand of sediment (g O/m 2 /d) and corresponding temperature adjustment coefficient, respectively; ke is the surface oxygen transfer coefficient (m/d), A(1) and V(1) are horizontal area (m 2 ) and lake volume (m 3 ) for the first or surface layer; YNHO2 and YCHO2 are the mass yield ratios of ammonia to oxygen and chlorophyll to oxygen, respectively. Calibration of the sediment oxygen demand is very important for simulating DO in the hypolimnion.
diffusion BOD sediment oxygen demand respiration and photosynthesis nitrification reaeration (surface only) zooplankton respiration (19) where K zr and θ zr are the zooplankton respiration rate (d −1 ) and the temperature adjustment coefficient, respectively; K NH(n) and K TN(n) are the half-saturation constant for preferential uptake of ammonia over nitrate and for nitrogen uptake for each algae class, respectively; S b and θ SOD are sediment oxygen demand of sediment (g O/m 2 /d) and corresponding temperature adjustment coefficient, respectively; k e is the surface oxygen transfer coefficient (m/d), A(1) and V(1) are horizontal area (m 2 ) and lake volume (m 3 ) for the first or surface layer; Y NHO2 and Y CHO2 are the mass yield ratios of ammonia to oxygen and chlorophyll to oxygen, respectively. Calibration of the sediment oxygen demand is very important for simulating DO in the hypolimnion.

Interaction and Connection among Modeling Variables
This study is focused on the internal nutrient dynamics and its interaction/connection to the phytoplankton, BOD and DO dynamics in the lakes; therefore, the inflow and outflow sub-model was disabled. The MINLAKE2020 DO model is different from that of MINLAKE2012, where daily Chla and BOD were specified as model input based on lake trophic status. For MINLAKE2020, several modifications were made to simulate phosphorus, nitrogen, Chla, and BOD concentrations and zooplankton activity in the DO simulation. The photosynthetic oxygen production calculation was modified, and the nitrification and zooplankton respiration were added to the DO model while the rest of the terms were the same as in MINLAKE2012. In MINLAKE 2020, photosynthetic oxygen production is dependent on the simulated daily Chla concentration rather than the data driven Chla pattern and annual mean concentration used in MINLAKE2012.
The diffusion of DO occurs between the water layers in the metalimnion and hypolimnion. The spring and fall overturn periods for temperate lakes completely mix water columns and make all modeling variables distributed uniformly along with the depth. Since zooplankton usually spend the largest amount of their time in the day depth layer, zooplankton respiration is simulated at the day depth only. BOD removes oxygen from the water layer through the decay of detritus which is accumulated from phytoplankton mortality and removed by settling ( Figure 6). Nitrification removes oxygen from the water layer through the conversion of ammonia to nitrite then to nitrate. Nitrification is applied to the DO model only if nitrogen is simulated.
From Figures 1-7 and from Equations (9)- (19), one can see the complex interactions and connections among modeling variables. For example, when the settling velocity of the phytoplankton group is increased, it reduces Chla concentration and in turn, affects detritus (BOD) and phosphorus, but phosphorus goes back to affect phytoplankton growth. Additionally, DO concentrations in different water layers connect/integrate changes in water temperature, nutrients, and phytoplankton dynamics (Figure 1). The sediment oxygen demand is an important factor for DO mass balance and results in the anoxic condition in the hypolimnion that leads to phosphorus release from sediment. Therefore, exploring/understanding the internal cycles/dynamics of nutrients/phytoplankton and their interactions in different lakes helps us to gain insights of lake ecosystem and then develop appropriate restoration strategies.

Model Coefficients and Parameters
MINLAKE2020 model was designed to simulate small lakes (A s < 25 km 2 ) provided that the user specifies the input data and calibration parameters accordingly. Variation in lake characteristics is reflected in model input data/parameters. Lake bathymetry and weather data (depending on lake geographic location) need to be supplied to the model. For many lakes in Minnesota, depth or elevation contour lines can be downloaded from the Minnesota Department of Natural Resources (MN DNR) LakeFinder website (https://www.dnr.state.mn.us/lakefind/index.html, accessed on 5 June 2019), from which horizontal areas at different elevations/depths can be determined. The weather data include daily air temperature ( • F), dew point temperature ( • F), wind speed (mph), solar radiation (Langley), sunshine percentage, and precipitation including rainfall (cm) and snowfall (mm). To facilitate the comparison of DO simulation results, in MINLAKE2020, the user can run the DO simulation in two ways: (1) using MINLAKE 2012 regional DO model, i.e., nutrients and Chla are not simulated (called the RegDO model); (2) using MIN-LAKE2020 (called the NCDO model). This model first simulates nutrients and Chla and then DO [67]. Table 1 lists nutrient, Chla, and DO calibration parameters in MINLAKE2020 and includes descriptions and effects on specific model results. The number of algal classes and the light attenuation coefficient are important input parameters for the simulation. Moreover, the snow and ice model require various coefficients (e.g., snow and ice density, specific heat, thermal conductivity, etc.) as input parameters, which has been well tested in the previous studies [68]. The temperature adjustment coefficients for BOD, photosynthesis, respiration, and sediment oxygen demand were set to 1.047 [69], 1.066 [57], 1.047 [20], 1.065 [70], respectively. Mass ratios of Chla to oxygen, phosphorus to oxygen, and phosphorus to chlorophyll-a are 0.0083, 0.0091, and 1.1, respectively. The inclusion of phytoplankton and zooplankton simulation in MINLAKE2020 calls for many additional input parameters or model coefficients. Most of the zooplankton-related coefficients were taken from West and Stefan's study [23]. Zooplankton respiration rate was set to 0.002 d −1 and the reproduction rate was 0.02 d −1 . Minimum light intensity for zooplankton predation, light intensity for maximum zooplankton predation, Julian day for the end of low predation period, and Julian day for the beginning of maximum predation were set to 0 µE (m 2 s) −1 , 0.1 µE (m 2 s) −1 , 110 (20 April), and 140 (20 May), respectively. Minimum seasonal day time predation rate, maximum seasonal day time predation rate, and overnight predation rate were set to 0.05 d −1 , 0.7 d −1 , and 0.03 d −1 , respectively. For any lake, these parameters were kept constant whereas some of the model parameters were updated/calibrated based on the comparison of simulation results with observed data. The first four parameters were used in the regional DO model, but the respiration rate K r was not calibrated (constant).

Lakes Simulated
In this study, six lakes ( Table 2) were selected for the model calibration, the sensitivity analysis of model parameters/coefficients, and understanding the interaction/connection among seasonal variations of nutrients and chlorophyll-a: two shallow lakes (Pearl and Carrie), two medium-depth lakes (Riley and Thrush), and two deep lakes (Carlos and Elmo). The maximum depths range from 5.6 to 50.0 m. All six lakes are located in northeastern Minnesota since they have the necessary data for the study. The nearest weather station to the lake was selected for providing weather data: St. Cloud Regional Airport for Pearl, Carrie, and Carlos lakes; Duluth International Airport for Thrush and Riley lakes; and Minneapolis/St. Paul International Airport for Elmo Lake. The geometry ratio (GR = A s 0.25 /H max , A s in m 2 , and H max in m being the surface area and the maximum depth of the lake) is a very important characteristic parameter of a lake that is related to stratification, lake habitat, etc. [71]. The lake geometry ratio is between 0.75 and 7.53. The lower the geometry ratio, the stronger the lake stratification. Two medium-depth (6 m < H max ≤ 20 m, [72]) and two deep lakes (H max > 20 m, [72]) selected for the study are strongly stratified (geometry ratio less than 3), one medium-depth and one shallow lake are weakly stratified (geometry ratio between 3 and 10). Based on observed Chla concentration, Pearl, Carrie, and Riley lakes are eutrophic (mean Chla > 10 µg/L [73]), Elmo and Carlos are mesotrophic lakes (mean Chla between 4 and 10 µg/L [73]), and Thrush is an oligotrophic lake. The nutrient model was calibrated and validated based on available measured water temperature, chlorophyll-a, phosphorus and DO profile data on particular days, downloaded from the LakeFinder website.

Model Calibration
For MINLAKE 2020 model application to the six study lakes, the temperature model was calibrated first and then the nutrient model was calibrated. Temperature model calibration ensured that thermal and mixing dynamics were modeled accurately because water temperature and mixing dynamics directly affect nutrients, Chla, and zooplankton processes (Equations (9)- (19)). The wind sheltering coefficient and the multiplier for diffusion coefficient in metalimnion are the main calibration parameters for temperature modeling. Although MINLAKE2020 has an integrated nitrogen model, for this study, only phosphorus was simulated since phosphorus is the limiting nutrient in these six lakes. Green algae and blue-green algae were simulated separately and then combined to represent the total chlorophyll-a concentration. The MINLAKE2020 development also included the inflow and outflow subroutines from MINLAKE88, which were tested/verified to ensure they function properly; however, the inflow/outflow function was not activated for the simulation of the six study lakes ( Table 2). Certain inflow and outflow were reported for Carrie, Pearl, and Carlos Lake [74][75][76] whereas inflows in the other three lakes were minor, and the inflow quality (nutrients and phytoplankton) data were scarce. These approximations are appropriate since the study objective is to examine/understand the internal dynamics and cycles of nutrients over multiple years in six lakes with different stratification and trophic characteristics. Figure 8 shows an example of the calibration results of water temperature and DO time series at two depths (1 and 7 m) at Lake Carrie including measured data. During 2008-2010, Lake Carrie had measured water temperature and DO profile data for 36 days or 342 data points in total. MINLAKE2020 simulated water temperature and DO with a root mean square error of 1.75 • C and 1.95 mg/L, respectively. Corresponding regression coefficients of measured versus simulated (R 2 ) are 0.99 and 0.93, respectively. The statistical results summarized in Table 3 show that for the six lakes, MINLAKE2020 model performed better than MINLAKE2012, especially for DO simulations, when simulated profiles were compared with observed profiles. The main reason for this improvement is the simulation of Chla concentration on daily time step rather than using the specified pattern of observed data. Table 3 shows the model performance improved significantly with the NCDO model in Carlos and Thrush lakes. The average root-mean-square error (RMSE) of DO simulations in six lakes from MINLAKE2020 decreased by 24.2%, and average Nash-Sutcliffe efficiency (NSE) [77] also increased with respect to MINLAKE2012. Chlorophyll concentration affects the solar radiation attenuation in the water column (Equation (8)) and then affects water temperature simulation as shown in Table 3 even though the average RMSE, NSE, and R 2 for temperature (regression coefficient of measured versus simulated) from the two models are almost the same.
For Chla simulation of six lakes, RMSE ranges from 0.0006 to 0.0276 mg/L. Figure 9 shows an example comparison between simulated Chla of MINLAKE2020 (NCDO model) and specified Chla from MINLAKE2012 (RegDO model) for Lake Elmo. From 1980 to 2018, there were 74 days in 10 years with measured temperature and DO profiles (1506 data points) for model calibration but no profile data in 2007 and 2009. The average chlorophyll concentration was 0.0075 mg/L in 74 days but 0.0036 mg/L over 10 days in 2008. Even though there were no profile data in 2009, we identified some Chla data in 2009 as shown in Figure 9. MINLAKE2012 uses the annual mean Chla concentration and seasonal variation patterns [22] (depending on trophic status) to specify daily Chla for DO simulation. Therefore, the RegDO model had higher Chla in 2007 and 2009 due to the lack of available profile data in these two years whereas the NCDO model predicts the Chla reasonably well in 2008 and 2009 when comparing with data. NCDO model in Carlos and Thrush lakes. The average root-mean-square error (RMSE) of DO simulations in six lakes from MINLAKE2020 decreased by 24.2%, and average Nash-Sutcliffe efficiency (NSE) [77] also increased with respect to MINLAKE2012. Chlorophyll concentration affects the solar radiation attenuation in the water column (Equation (8)) and then affects water temperature simulation as shown in Table 3 even though the average RMSE, NSE, and R 2 for temperature (regression coefficient of measured versus simulated) from the two models are almost the same.   Some lakes, such as Lake Elmo and Carrie Lake (Table 3), do not exhibit a noticeable change in simulated DO concentrations based on the model used for simulation. Some lakes, e.g., Pearl Lake, exhibit a noticeable change in simulated DO concentration depending on the model (RegDO or NCDO model). Figure 10a When BOD is simulated, the winter DO decreases, predicted by MINLAKE2020, are smaller than those predicted by MINLAKE2012 when BOD is specified as a part of the model inputs. There are very limited data for P comparison between simulated and observed, and RMSE ranges from 0.005 to 0.036 mg/L. DO simulation is an overall model performance indication (Table 3).

Chla and Phosphorus Profiles
Lake Elmo was extensively monitored in 1988 by the Metropolitan Council [65] and had measured phosphorus concentration data at five depths (0, 8, 16, 24, and 32 m) and DO data at 31 depths for open water season. The comparison between simulated and observed concentrations for Chla, phosphorus, and DO on three days in 1988 is presented in Figure 11. Since Elmo is a deep lake and solar radiation cannot penetrate below the euphotic zone, the Chla concentration becomes zero in the deep layers. On 11 April 1988, the lake is more or less well mixed and phosphorus concentration did not vary much throughout the depth but DO concentration gradually declined along with depth due to the con-

Chla and Phosphorus Profiles
Lake Elmo was extensively monitored in 1988 by the Metropolitan Council [65] and had measured phosphorus concentration data at five depths (0, 8, 16, 24, and 32 m) and DO data at 31 depths for open water season. The comparison between simulated and observed concentrations for Chla, phosphorus, and DO on three days in 1988 is presented in Figure 11. Since Elmo is a deep lake and solar radiation cannot penetrate below the euphotic zone, the Chla concentration becomes zero in the deep layers. On 11 April 1988, the lake is more or less well mixed and phosphorus concentration did not vary much throughout the depth but DO concentration gradually declined along with depth due to the con-

Chla and Phosphorus Profiles
Lake Elmo was extensively monitored in 1988 by the Metropolitan Council [65] and had measured phosphorus concentration data at five depths (0, 8, 16, 24, and 32 m) and DO data at 31 depths for open water season. The comparison between simulated and observed concentrations for Chla, phosphorus, and DO on three days in 1988 is presented in Figure 11. Since Elmo is a deep lake and solar radiation cannot penetrate below the euphotic zone, the Chla concentration becomes zero in the deep layers. On 11 April 1988, the lake is more or less well mixed and phosphorus concentration did not vary much throughout the depth but DO concentration gradually declined along with depth due to the contribution of more sink terms (Equation (18)) but for the profile plot, the slope was not steep. The Chla concentration is highest at the surface and did not vary much throughout the depth. On 18 May 1988, the stratification increases and the simulated DO at the bottom is near zero. The Chla concentration is not maximum at the surface, but at 8 m depth from the surface. The phosphorus concentration is higher at the deeper layers because of detrital decay and phytoplankton respiration. On 19 October 1988, the phosphorus near the surface layer is being used by the phytoplankton for growth, the lake became strongly stratified and the bottom layers of the lake become anoxic so that the phosphorus release from sediment contributed to the higher phosphorus concentration at the deeper depths, which increased along with depth from metalimnion and hypolimnion. steep. The Chla concentration is highest at the surface and did not vary much throughout the depth. On 18 May 1988, the stratification increases and the simulated DO at the bottom is near zero. The Chla concentration is not maximum at the surface, but at 8 m depth from the surface. The phosphorus concentration is higher at the deeper layers because of detrital decay and phytoplankton respiration. On 19 October 1988, the phosphorus near the surface layer is being used by the phytoplankton for growth, the lake became strongly stratified and the bottom layers of the lake become anoxic so that the phosphorus release from sediment contributed to the higher phosphorus concentration at the deeper depths, which increased along with depth from metalimnion and hypolimnion. Thrush Lake has measured Chla data at both the epilimnion and hypolimnion. The comparison between simulated and observed concentrations for Chla and DO on three days in 1986 is plotted in Figure 12. Thrush is an oligotrophic lake with lower oxygen demands (low BOD and SOD); therefore, the simulated and observed bottom DO is greater than zero, and there is no phosphorus release from sediment in all three days (Fig-Figure 11. Simulated (line) and observed (dots) chlorophyll-a, phosphorus, and dissolved oxygen profiles for Lake Elmo on (a) 11 April, (b) 18 May, and (c) 19 October 1988. Thrush Lake has measured Chla data at both the epilimnion and hypolimnion. The comparison between simulated and observed concentrations for Chla and DO on three days in 1986 is plotted in Figure 12. Thrush is an oligotrophic lake with lower oxygen demands (low BOD and SOD); therefore, the simulated and observed bottom DO is greater than zero, and there is no phosphorus release from sediment in all three days ( Figure 12). Since no phosphorus data were available, only simulated phosphorus was plotted and has no increase in the hypolimnion. Due to lower light attenuation coefficients, solar radiation can penetrate through the deepest layers of the lake. The calculated euphotic depth (at 1% surface solar radiation) is equal to or greater than the maximum depth (14.6 m) for all of the three selected dates; therefore, the simulated Chla is not zero near the bottom but increases below the mixed layer, especially on 17 June 1986. On 13 May 1986, the Chla and phosphorus concentration did not vary much throughout the depth. The DO concentration gradually decreased below the mixed layer on 13 May and 17 June 1996, due to the contribution of more sink terms. On 14 October 1986, the lake became well mixed due to the fall overturn and had uniform distribution for simulated Chla, phosphorus, and DO. The model underpredicts the mixed layer depth on 13 May and then overpredicts DO at the hypolimnion.   Figures 13 and 14 show examples of simulated time series of P, Chla, and DO with observed data for a deep lake (Elmo) and a shallow lake (Carrie), respectively. Figure 13 shows simulated and observed Chla and phosphorus (P) at three depths (at 1 m, 20 m, 40 m from the surface) of Elmo Lake from 16 April 2007 to 31 December 2009. This is a continuous year-round simulation including two open-water seasons and two ice cover periods, which was not done before using MINLAKE models for Chla and P simulation. The  Figures 13 and 14 show examples of simulated time series of P, Chla, and DO with observed data for a deep lake (Elmo) and a shallow lake (Carrie), respectively. Figure 13 shows simulated and observed Chla and phosphorus (P) at three depths (at 1 m, 20 m, 40 m from the surface) of Elmo Lake from 16 April 2007 to 31 December 2009. This is a continuous year-round simulation including two open-water seasons and two ice cover periods, which was not done before using MINLAKE models for Chla and P simulation.

Chlorophyll-a and Phosphorus Interaction
The first year of simulation is considered as model warm-up period and the results may have more uncertainties due to the assumed initial conditions, for example, phosphorus concentration from the water surface to 20 m was low during the open water season (Figure 13a,b). The simulated ice cover was from 5 December 2007 to 16 April 2008, 7 December 2008 to 10 April 2009 for Elmo Lake, which is marked by blue shaded regions. late spring, the simulated phosphorus concentrations at 1 m increase for a brief period due to spring overturn and then, start to decrease as a result of the phosphorus uptake by phytoplankton. During the early summer (May) of 2008 and 2009, simulated Chla concentration increases gradually from 0.0016 mg/L to a maximum of 0.0033 mg/L (observed on 15 June).
At the deep layers (e.g., 20 and 40 m), DO becomes anoxic during the summer at deep hypolimnion since Elmo Lake is strongly stratified. Anoxic periods at 20 m and 40 m were on average 12 days and 210 days per year from 2007 to 2009, respectively. In the deeper layers, in addition to detrital decay and phytoplankton respiration, sediment release could add up to the available phosphorus. Since there is a long period of anoxic condition at 40 m during the summer, early fall, and some part of the ice cover period, phosphorus release from sediment contributes to a major portion of phosphorus increase. Phosphorus peaks in deepest layers were simulated to occur just after the anoxic condition ends and before the fall or spring mixing/overturns. These overturns sharply reduce phosphorus in very deep layers (e.g., 40 m) and increases phosphorus in other shallower layers. The calculated euphotic depth is 9 m for the simulation period; therefore, there is no photosynthesis below this depth and simulated Chla is zero at 20 m and 40 m. Chla is typically measured near the surface; therefore, there is no measured Chla at the deep depths to compare with simulated values. Near the water surface (1 m), DO concentrations are near saturation as a function of temperature (lowest DO in the middle of summer) and range from 6.69-11.13 mg/L. From late October to late November, before the ice starts to form at the surface of the lake, phosphorus at the surface and other surface layers start to increase due to more mixing and fall overturns. This increase is more evident in 2007 when phosphorus was low in the summer. During the ice cover period, the phosphorus concentration becomes stable at 1 m and 20 m as the organic processes (photosynthesis) become slow due to the near-zero water temperature and/or low light (attenuated by snow cover). After the ice melts out in late spring, the simulated phosphorus concentrations at 1 m increase for a brief period due to spring overturn and then, start to decrease as a result of the phosphorus uptake by phytoplankton. phorus release decreases which results in a gradual decrease in phosphorus concentration. The DO concentration decreases with depth because of no photosynthesis in the deeper layers (below the euphotic zone) plus sedimentary oxygen demands. Simulated DO at 7 m has some fluctuations in the summer of 2008 and 2009 due to short period strong mixing and results in anoxic conditions only in a few days. The simulated high phosphorus concentrations directly correspond to the simulated anoxic DO conditions in the 2008 and 2009 winters (Figure 14b).

Sensitivity Analysis
MINLAKE2020 model is sensitive to several calibration parameters. The model was first calibrated with a regression coefficient (R 2 for measured versus modeled profile data) of 0.8972 for DO simulation; then, only one model parameter was changed at a time and the regression coefficient for each new run of DO simulation was determined. Table 4 lists the calibration parameters with the calibrated and uncalibrated values and the regression coefficients of DO simulations with uncalibrated values for Lake Elmo in 2007-2009. In Figures 15 and 16, sensitivity analysis graphs are plotted for five parameters: maximum photosynthesis rate Gmax (used in Equations (9), (15)- (17) and (19)), sediment phosphorus release rate SP (used in Equation (9)), half-saturation coefficient of phosphorus KP (used in Equations (9), (15)- (17) and (19)), minimum Chla for grazing Chlamin (used in Equation (9)), At the deep layers (e.g., 20 and 40 m), DO becomes anoxic during the summer at deep hypolimnion since Elmo Lake is strongly stratified. Anoxic periods at 20 m and 40 m were on average 12 days and 210 days per year from 2007 to 2009, respectively. In the deeper layers, in addition to detrital decay and phytoplankton respiration, sediment release could add up to the available phosphorus. Since there is a long period of anoxic condition at 40 m during the summer, early fall, and some part of the ice cover period, phosphorus release from sediment contributes to a major portion of phosphorus increase. Phosphorus peaks in deepest layers were simulated to occur just after the anoxic condition ends and before the fall or spring mixing/overturns. These overturns sharply reduce phosphorus in very deep layers (e.g., 40 m) and increases phosphorus in other shallower layers. The calculated euphotic depth is 9 m for the simulation period; therefore, there is no photosynthesis below this depth and simulated Chla is zero at 20 m and 40 m. Chla is typically measured near the surface; therefore, there is no measured Chla at the deep depths to compare with simulated values.

Sensitivity Analysis
MINLAKE2020 model is sensitive to several calibration parameters. The model was first calibrated with a regression coefficient (R 2 for measured versus modeled profile data) of 0.8972 for DO simulation; then, only one model parameter was changed at a time and the regression coefficient for each new run of DO simulation was determined. Table 4 lists the calibration parameters with the calibrated and uncalibrated values and the regression coefficients of DO simulations with uncalibrated values for Lake Elmo in 2007-2009. In Figures 15 and 16, sensitivity analysis graphs are plotted for five parameters: maximum photosynthesis rate G max (used in Equations (9), (15)- (17) and (19)), sediment phosphorus release rate S P (used in Equation (9)), half-saturation coefficient of phosphorus K P (used in Equations (9), (15)- (17) and (19)), minimum Chla for grazing Chla min (used in Equation (9)), and settling velocities of algae (v c ) and BOD (v B ) (used in Equations (9) and (18)) for Elmo and Pearl Lake, respectively. Units of these model parameters are given in Table 4, and the below corresponding equations. In addition to the predicted increase, the Chla concentration decreases in spring of 2011 and 2012, early summer and fall of 2012 as a result of the spring and fall overturn when Gmax is increased from 0.25 to 0.4. The Chla is not very sensitive to the half-saturation coefficient (KP) and the minimum Chla threshold for zooplankton grazing (Chlamin). However, the phosphorus concentration decreased by an average of 60% when KP is decreased from 0.03 to 0.005 mg/L ( Figure 16c). As KP decreases, the Chla increases throughout the time period except for a portion of summer and fall of 2012. This happens due to the increased respiration of algae which also causes a greater reduction is phosphorus. The zooplankton grazing, being negligible for Pearl Lake, does not affect the phosphorus concentration much. Hence, phosphorus is not sensitive of zooplankton grazing in Pearl Lake.   In Figures 15 and 16, the blue line corresponds to simulation results using the calibrated value of the parameter which has produced satisfactory results. The calibrated maximum photosynthesis rate G max is 0.6 for Lake Elmo. When G max is increased to 0.9, the Chla concentration increases throughout the simulation period on average by 30 percent; the phosphorus concentration increases in the ice cover period because of the detrital decay of increased algae. Higher G max corresponds to higher phytoplankton growth which uses more nutrients (phosphorus) except the slow growth period under ice cover. This results in an average 32% decrease in phosphorus concentration in open water season.
When the sediment phosphorus release rate S P was decreased from 0.02 to 0.01, the phosphorus concentration decreased by 55% (on average). The phosphorus starts to decrease during the ice cover period in December 2007 as the sediment phosphorus release start contributing to the phosphorus concentration in anoxic condition. After ice melting, as the algae begin to grow and use up phosphorus, the phosphorus concentration decreases even further in summer. The half saturation coefficient of phosphorus is decreased from 0.07 mg/L to 0.035 mg/L. Since the surface layer is not limited by light, the decrease in the half saturation coefficient results in greater phosphorus limitation which results in increased algal growth. The phosphorus concentration decreases except for the ice cover period which is governed by the sediment phosphorus release in anoxic condition. The increase in the minimum chlorophyll-a threshold for grazing results in less grazing of zooplankton and a subsequent increase in Chla. The settling velocity of algae (both green algae and blue-green algae) is a very sensitive parameter for phytoplankton and BOD simulation. In Figure 15e, as the settling velocity of green algae, blue-green algae, and detritus was decreased to 0.75, 0.5, and 1.0, respectively, the phytoplankton concentration increased due to less loss of phytoplankton from the simulating layer from settling. The phosphorus concentration increases due to the release of phosphorus through respiration of the increased algal population.
For a shallow lake such as Pearl, the effect of maximum photosynthesis rate (Figure 16a) is not as straightforward as for a deep stratified lake such as Elmo Lake (Figure 15a). In addition to the predicted increase, the Chla concentration decreases in spring of 2011 and 2012, early summer and fall of 2012 as a result of the spring and fall overturn when G max is increased from 0.25 to 0.4. The Chla is not very sensitive to the half-saturation coefficient (K P ) and the minimum Chla threshold for zooplankton grazing (Chla min ). However, the phosphorus concentration decreased by an average of 60% when K P is decreased from 0.03 to 0.005 mg/L ( Figure 16c). As K P decreases, the Chla increases throughout the time period except for a portion of summer and fall of 2012. This happens due to the increased respiration of algae which also causes a greater reduction is phosphorus. The zooplankton grazing, being negligible for Pearl Lake, does not affect the phosphorus concentration much. Hence, phosphorus is not sensitive of zooplankton grazing in Pearl Lake.

Long-Term Simulations Using MINLAKE2020
West and Stefan [23] performed a multiple-year simulation (same calibration parameters) using MINLAKE98 for Lake Riley and Elmo. For Lake Riley, a different set of calibration parameters was needed for different years whereas, for Lake Elmo, the model could simulate successfully for 1985-1990 with the regression coefficient for temperature and DO as 0.91 and 0.79, respectively. For a simulation of 1985-1990 using the MIN-LAKE2020 NCDO model, the regression coefficient for temperature and DO are 0.9944 and 0.9715 against 146 profile data points, respectively. Simulated phosphorus and Chla concentrations at different depths are satisfactory as well. MINLAKE2020 performed well for multiple-year simulation allowing the user to simulate 20 consecutive years with the

Long-Term Simulations Using MINLAKE2020
West and Stefan [23] performed a multiple-year simulation (same calibration parameters) using MINLAKE98 for Lake Riley and Elmo. For Lake Riley, a different set of calibration parameters was needed for different years whereas, for Lake Elmo, the model could simulate successfully for 1985-1990 with the regression coefficient for temperature and DO as 0.91 and 0.79, respectively. For a simulation of 1985-1990 using the MIN-LAKE2020 NCDO model, the regression coefficient for temperature and DO are 0.9944 and 0.9715 against 146 profile data points, respectively. Simulated phosphorus and Chla concentrations at different depths are satisfactory as well. MINLAKE2020 performed well for multiple-year simulation allowing the user to simulate 20 consecutive years with the same calibration parameters (Figure 17). For a 20-year simulation (1989-2009) using the MINLAKE2020 NCDO model, the regression coefficients for temperature and DO are 0.9888 and 0.9419 against 864 profile data points, respectively. The simulated Chla and phosphorus match reasonably well with observed values with the same trend. Moreover, the phosphorus and Chla concentration at five simulation depths (1 m, 8 m, 12 m, 20 m, and 30 m) match well with the available observed data [54].

Conclusions
A one-dimensional daily water quality model MINLAKE2020 was developed from the daily temperature and DO MINLAKE2012 model by incorporating phytoplankton, zooplankton, nutrient, and BOD simulation into the model. The inflow-outflow submodel of MINLAKE2020 was disabled for this study to focus on the internal nutrient dynamics inside the lake. The model was applied to six Minnesota lakes with varying characteristics in terms of depth (two shallow lakes, two medium-depth lakes, and two deep lakes) and trophic status (two eutrophic, two mesotrophic, and two oligotrophic lakes). The simulated water temperature, DO, Chla, and phosphorus time series and profiles were compared with available observed data in 15-36 days for two to four years. The model was also applied to long-term simulation over 20 years (1989-2009) for Lake Elmo. Simulation results from the MINLAKE2020 model provide the following conclusions: 1. MINLAKE2020 was calibrated against measured profiles in six Minnesota lakes (Table 4) for the short term (2-4 years) with an average standard error of 1.51 °C for temperature and 2.33 mg/L for DO. The average standard error for DO simulation of From 1989 to 1996, both phosphorus and chlorophyll-a seasonal variations were reasonably stable. Phosphorus started to increase from 1997 and matched with observed data from 2004 to 2007. The average daily phosphorus was 0.0071 mg/L from 1990 to 1996 and 0.0158 mg/L from 1997 to 2009. Comparing these two periods, the average daily phosphorus increased by 0.0087 mg/L. As phosphorus increased, it resulted in some higher peaks in spring algal blooms as shown in Figure 17a. The phosphorus increase trend is caused by the increase in phosphorus release from the lake sediment which is related to the anoxic condition in lakes. Therefore, the phosphorus release for each layer (the last term in Equation (15) times the layer volume) and then each day (sum for all layers) was outputted and added together for the annual phosphorus release amount. From 1990 (excluding 1989 for the initial condition effect) to 1996, the average yearly sediment phosphorus release is 151.8 kg (21.27 kg of standard deviation) but from 1997 to 2009 average yearly sediment phosphorus release is 244.1 kg (53.28 kg of standard deviation). The average annual phosphorus release increases by 60.8%. Figure 17c shows the time series of the simulated DO at 41 m (1 m above the deepest lake bottom) from 1989 to 2009 and clearly shows many anoxic days in the open water seasons and the ice cover periods, which result in phosphorus release from sediment. From 1990 to 1996, the average anoxic days is 228 but from 1997 to 2009 average anoxic days is 253 (34 days of standard deviation). An average of 25 days more of the anoxic condition is one of the causes of the phosphorus increase trend. From DO simulation it was observed that the anoxia started at a lower depth in the water column (16 m) from 1999 to 2009 compared to earlier simulation years. As a result, the anoxia had a greater horizontal-area coverage. Since the sediment release is a function of the bottom area (Equation (15)), the increased bottom area resulted in the increased sediment phosphorus release.

Conclusions
A one-dimensional daily water quality model MINLAKE2020 was developed from the daily temperature and DO MINLAKE2012 model by incorporating phytoplankton, zooplankton, nutrient, and BOD simulation into the model. The inflow-outflow submodel of MINLAKE2020 was disabled for this study to focus on the internal nutrient dynamics inside the lake. The model was applied to six Minnesota lakes with varying characteristics in terms of depth (two shallow lakes, two medium-depth lakes, and two deep lakes) and trophic status (two eutrophic, two mesotrophic, and two oligotrophic lakes). The simulated water temperature, DO, Chla, and phosphorus time series and profiles were compared with available observed data in 15-36 days for two to four years. The model was also applied to long-term simulation over 20 years (1989-2009) for Lake Elmo. Simulation results from the MINLAKE2020 model provide the following conclusions: 1.
MINLAKE2020 was calibrated against measured profiles in six Minnesota lakes ( The addition of phosphorus and Chla simulation in MINLAKE2020 improved model performance in comparison to MINLAKE2012 where Chla was specified input. It greatly affects the DO concentration in some lakes such as Pearl Lake ( Figure 10). Thrush Lake and Carlos Lake also showed significant improvement in DO simulation with MINLAKE2020. The standard error decreased by 2.12 mg/L and 1.76 mg/L for Thrush and Carlos Lake, respectively. 3.
The deep lakes exhibit a certain trend for phosphorus and Chla simulation year by year whereas the shallow lakes might show a significant change in phosphorus and Chla concentration year by year due to two overturn periods (complete mixing) and the complex interactions/connections among phosphorus, Chla, and DO ( Figures 13 and 14), which are evident through governing Equations (9) and (15)- (19) and processes simulated (Figures 1-7). 4.
DO concentration is a primary control of internal loading via anoxic release of phosphorus from the lake sediment. MINLAKE2020 was applied to Lake Elmo for a 20-year (1989  MINLAKE2020 explains the internal link between phosphorus, Chla, and DO. This model can help in choosing/testing effective ways of lake restoration and management. For example, since Lake Elmo experiences internal loading of phosphorus, removing the external nutrient source might not be an effective restoration technique for this lake. Though MINLAKE2020 has simulated six lakes with good correlation to the observed data, more lakes need to be simulated to verify the model in different scenarios before using it professionally for evaluating lake restoration measures. Since the nutrient model is complex and requires a number of calibration parameters, incorporating automatic calibration (by programming software) will make the model more efficient and user-friendly in the future.

K NI
Nitrification rate (d −1 ) K NH(n) Half-saturation constant for preferential uptake of ammonia over nitrate [-] K TN(n) Half-saturation constant of nitrogen for each algal class (mg/L) K P Nocturnal predation rate (d −1 ) ∆T g (I) Time that zooplankton spent in a layer I during the night (d)