ECO: A Generic Eutrophication Model Including Comprehensive Sediment-Water Interaction

The content and calibration of the comprehensive generic 3D eutrophication model ECO for water and sediment quality is presented. Based on a computational grid for water and sediment, ECO is used as a tool for water quality management to simulate concentrations and mass fluxes of nutrients (N, P, Si), phytoplankton species, detrital organic matter, electron acceptors and related substances. ECO combines integral simulation of water and sediment quality with sediment diagenesis and closed mass balances. Its advanced process formulations for substances in the water column and the bed sediment were developed to allow for a much more dynamic calculation of the sediment-water exchange fluxes of nutrients as resulting from steep concentration gradients across the sediment-water interface than is possible with other eutrophication models. ECO is to more accurately calculate the accumulation of organic matter and nutrients in the sediment, and to allow for more accurate prediction of phytoplankton biomass and water quality in response to mitigative measures such as nutrient load reduction. ECO was calibrated for shallow Lake Veluwe (The Netherlands). Due to restoration measures this lake underwent a transition from hypertrophic conditions to moderately eutrophic conditions, leading to the extensive colonization by submerged macrophytes. ECO reproduces observed water quality well for the transition period of ten years. The values of its process coefficients are in line with ranges derived from literature. ECO’s calculation results underline the importance of redox processes and phosphate speciation for the nutrient return fluxes. Among other things, the results suggest that authigenic formation of a stable apatite-like mineral in the sediment can contribute significantly to oligotrophication of a lake after a phosphorus load reduction.


Introduction
Eutrophication models tend to focus on processes in the water column. The benthic-pelagic coupling and processes in the sediment have been neglected or very much simplified in most models [1][2][3][4][5][6][7][8]. However, the longer the residence time in a water system, the more return fluxes of nutrients from the sediment contribute to the availability of nutrients for algae. Because of longer residence times the adequate modeling of return fluxes is generally more important for lakes than for coastal seas. This holds in particular for stratified lakes and reservoirs, in which transient hypoxia affects the return fluxes strongly [9]. Jeppesen et al. [10] argued on the basis of data for 35 lakes that due to internal loading P oligotrophication is delayed 5 to 10 years after external load reduction.
The modeling of benthic-pelagic coupling requires the incorporation of sediment diagenesis processes. Sediment diagenesis models were developed to study processes or to quantify exchange fluxes between sediment and overlying water [11][12][13]. Generic diagenesis models are dynamic and multi-layered. Apart from organic matter and nutrients, these models contain C-Fe-Mn-S-Ca chemistry [14][15][16][17][18][19][20]. Although algal primary production in the water column often formed the broader context, sediment diagenesis models were mostly used with observed water quality as forcing function. Reducing the computational burden, simplified sediment sub-models were developed for incorporation in eutrophication models as process modules, for marine systems [21][22][23][24], [13] as well as for freshwater systems [25][26][27][28]. These models were used to study nutrient budgets [4], [17], eutrophication level [2,8], the effects of hypoxia [29,30] and the effects of mitigative measures [5,6,26,31]. Many of these models are deficient as to specific processes, or cover only one or two of the nutrients (N, P, Si). The simplified substances and processes content of their sediment sub-models often resulted in insufficient temporal variability of the nutrient return fluxes.
The substances usually included in eutrophication models concern the biomasses of phytoplankton species, one or two particulate detrital organic matter fractions and dissolved inorganic nutrients (N, P,Si). Inorganic suspended matter is often included, because it contributes to light limitation of phytoplankton. For the enhanced modeling of phytoplankton, water quality and nutrient return fluxes the addition of some of the following substances is required. Dissolved organic matter strongly affects the under water light regime. Oxygen, nitrate, manganese(IV), iron(III), sulphate and methane participate in redox processes that affect the speciation of phosphate in the sediment. The adsorption of phosphate to iron(III) oxyhydroxides and the precipitation of phosphate minerals determine the concentration of dissolved phosphate. The dissolution of opal silicate rules the concentration of dissolved silicate.
For the formulation of the diagenetic processes in eutrophication models various kinetic concepts are available that differ with regard to reaction mechanism and details taken into account [11,12]. Simplified models mostly use first order kinetics, the process rate being proportional to the concentration of the main reactant. The multi-G-model for the decomposition of organic matter [32] is an example of first-order kinetics. Advanced models may contain more versatile Michaelis-Menten kinetics, that take into account the maximum rate and the inhibition of a microbial enzymatic process [22]. The kinetics for mineral formation take the difference between actual and equilibrium conditions as a driving force, based on a pH-dependent concentration or a solubility product [14]. Theoretically, the more comprehensive a model is and the more realistic its kinetics are, the more generic it may be, because it predicts process rates reliably for a wider range of conditions. However, trade-offs need to be made to minimize the number of poorly quantifiable parameters.
In this article we present the 3D eutrophication model ECO, aiming at the comprehensive, deterministic and accurate simulation of phytoplankton biomass, water quality, sediment diagenesis and sediment-water interaction. We developed its advanced process formulations for substances in the water column and the bed sediment to allow for a much more dynamic calculation of the sediment-water exchange fluxes of nutrients (N, P, Si) as resulting from steep concentration gradients across the sediment-water interface than is possible with other eutrophication models, including its predecessor model DBS [25,26]. ECO is to more deterministically and consequently more accurately calculate the accumulation of organic matter and nutrients in the sediment, and to allow for more accurate prediction of phytoplankton biomass and water quality in response to mitigative measures such as nutrient load reduction. Its generic process formulations should make ECO applicable to both freshwater and marine systems. Due to its modular set-up ECO has minimal restrictions with regard to the simulation of spatial and temporal variability and the modification of its substances and processes content. The calibration of the ECO model for Lake Veluwe in the Netherlands is described. The resulting values of process coefficients are compared with ranges of reported values. The simulation results of the calibrated model, many of which are additional to the results of predecessor model DBS [25,26], are discussed. These results underline the importance of redox processes and phosphate speciation for the nutrient return fluxes, and entail new conclusions concerning the mechanisms that led to the oligotrophication of Lake Veluwe.

Structure, Schematization, Mass Transport
ECO is based on the open source water quality modeling framework D-Water Quality (DELWAQ), that facilitates the selection of substances and processes from a processes library and has many options for mass conservative numerical integration [24,33]. A DELWAQ model is defined by means of an input file with a system definition (substances, processes and pertinent coefficients), a computational grid, an initial composition, timers, flow fields (a dynamic water balance), dispersion coefficients, loads and meteorological forcing. All input parameters in the model can be specified as constants or temporally and/or spatially varying parameters. The output of a DELWAQ model includes the concentrations and mass balances of all simulated substances, the magnitudes of all imposed and simulated parameters, and the mass fluxes of all simulated processes.
DELWAQ solves the following 3D advection diffusion reaction equation for a finite volumes computational grid: where C: concentration (g.m 23 ) u, v, w: components of the velocity vector (m.s 21 ) Dx, Dy, Dz: components of the dispersion tensor (m 2 .s 21 ) x, y, z: coordinates in three spatial dimensions (m) S: sources and sinks of mass due to loads and boundaries (g.m 23 .s 21 ) P: sources and sinks of mass due to processes (g.m 23 .s 21 ) t: time (s) ECO dynamically simulates a set of substances and processes on a computational grid that is composed of a water grid and a sediment grid having the same horizontal resolution, allowing for a much finer vertical resolution than DBS [25] and a spatial detailedness that is innovative for eutrophication models. Each water compartment has a sediment column, consisting of a number of layers the thicknesses of which may increase with depth. Thin upper sediment layers allow for a realistic representation of the steep concentration gradients across the sedimentwater interface. The number and size of the grid cells is only limited by the maximum acceptable computational burden. The numerical solver applied in the ECO model of Lake Veluwe concerns an explicit upwind scheme in the horizontal direction, combined with an implicit in time scheme with central discretization of advection in the vertical direction. The computational time step is 5 minutes. It was verified that a further reduction of the time step did not affect the simulation results of the Lake Veluwe model significantly. The vertical mass transport in the sediment and across the sediment-water interface resulting from settling, resuspension, dispersion and seepage is generated by a process included in the processes library. This transport leads to the burial of substances below the boundary of the model in the sediment or the addition of substances from below this boundary. Dispersion coefficients and seepage velocity are input parameters. In the Lake Veluwe model, transport in the sediment is based on constant layer thickness and porosity, and on net settling (no resuspension).

Substances
ECO as applied for Lake Veluwe simulates phytoplankton biomass (C/N/P/Si/S) for five species (ALG i ), five fractions of detrital organic carbon, nitrogen, phosphorus and sulphur (particulate POC/N/P/S1, POC/N/P/S2, POC/N/P/S3, POC/N/P/S4 and dissolved DOC/N/P/S), nitrate (NO3), ammonium (NH4), dissolved phosphate (PO4), adsorbed phosphate (AAP), vivianite-P (VIVP), apatite-P (APATP), dissolved silicate (Si), opal silicate (OPAL), dissolved oxygen (DO), sulphate (SO4), dissolved and particulate sulphide (SUD, SUP), methane (CH4), two fractions of inorganic matter (IM1, IM2) and chloride (Cl). Compared to DBS [26], POX3, POX4, DOX, APATP, OPAL, all sulphur components, CH4, IM1 and IM2 have been added. Oxygen, sulphate and methane play a dominant role in the decomposition of organic matter, and together with nitrate determine the vertical gradient of the redox potential in the sediment. This gradient regulates the adsorption or desorption of phosphate to iron(III) oxyhydroxides. VIVP and APATP are phosphate minerals that contain iron or calcium, and in sediment are associated with other iron or calcium minerals. OPAL represents the silicate skeletons resulting from diatom mortality. IM1 and IM2 represent silt (,63 mm) and sand ($63 mm), and are used to fix sediment porosity, to generate burial and to calculate the adsorption capacity for phosphate in the sediment. IM3 is used to force the inorganic suspended sediment concentration, which contributes to the extinction of light and determines the adsorption capacity for phosphate in the water column. Conservative substance Cl was used to check the water balance of the Lake Veluwe model as to accurateness.

Processes Included in ECO
The processes in ECO are selected from DELWAQ's processes library. The processes included in ECO for Lake Veluwe are: N growth, mortality, grazing of phytoplankton N extinction of light N decomposition of detrital organic matter N consumption of electron acceptors, methanogenesis N nitrification, denitrification N adsorption, precipitation of phosphate N dissolution of opal silicate N oxidation, precipitation, speciation of sulphide N oxidation, ebullition, volatilization of methane N reaeration of dissolved oxygen N net settling of particulate components N mass transport in the sediment Figure 1 provides an overview of the substances and processes. Tables 1, 2, 3, 4 specify the mathematical formulations and the symbols. The reaction term P in the above advection diffusion equation is given for all substances in Equations A.1-16. The formulations for the individual processes and parameters are presented in sections B-G of Tables 1, 2. The process formulations are generic for water column and sediment bed, but may turn out differently in individual grid cells depending on local chemical conditions such as the presence or absence of dissolved oxygen.
Compared to DBS [26], the process formulations for organic matter decomposition, electron acceptor consumption, nitrification, phosphate adsorption. phosphate mineral formation and opal dissolution have been improved, and the sulphur and methane processes have been added.

Phytoplankton Processes
Phytoplankton is simulated with BLOOM as was done previously with DBS [26]. This model and the calibration of phytoplankton properties were reported in full detail by Los [34], its concepts and main formulations also by Los & Wijsman [35] and Blauw et al. [24]. Therefore, only the overall term P for algae biomass is included in Table 1 (Eq. A.1). BLOOM allows for the    N the energy type, with relatively high growth rate, low mortality rate and high N/C and P/C ratio; N the nitrogen type, typically with lower internal N/C ratio, lower maximum growth rate, higher mortality rate, higher settling velocity and higher chlorophyll-a content; and N the phosphorus type, typically with lower internal P/C ratio, lower maximum growth rate, higher mortality rate, lower settling velocity and lower chlorophyll-a content.
The biomasses of the phenotypes are modelled as 14 separate variables (ALG i ) with different nutrient and chlorophyll-a contents, extinction coefficient, growth rate, respiration rate, mortality rate and settling velocity (Table 5; [34]). Additional coefficients concern the fractions for autolysis of nutrients and production of detritus (Tables 3, 4; [34]). When ambient conditions change, a phenotype can be instantaneously converted into another phenotype of the same species, representing rapid adaptation of algal cells. In BLOOM the species compete within the constraints for available nutrients (N, P, Si), available light (energy), the maximum growth rate and the maximum mortality rate (both temperature functions). Linear programming is used as an optimization technique to determine the species composition that is best adapted to prevailing environmental conditions. It can be shown mathematically that the principle by which each phytoplankton type maximizes its own benefit effectively means that the total net production of the phytoplankton community is maximized [34]. Process fluxes are calculated on the basis of daily averaged meteorological forcing. Therefore, the simulation time step for the BLOOM phytoplankton processes is 24 hours. Phytoplankton settled in the sediment dies according to BLOOM's mortality rate. Grazing by zooplankton is taken into account by means of a grazing rate that operates on algae biomass (ALG i ) as well as detritus (POC/ N/P). This rate is proportional to imposed zooplankton biomass as described by [24].

Extinction of Light
Light is simulated as photosynthetically active radiation (PAR) [24], [34]. Extinction is modelled as an exponential decrease of light intensity with depth according to the Lambert-Beer formula. The total extinction coefficient is calculated as the sum of the extinction by inorganic suspended matter (IM3 in the Lake Veluwe model), particulate organic matter (POC1-4), dissolved organic matter (DOC), phytoplankton biomass (ALG 1-14 ) and water (Eqs. G1-4 in Table 2). Each substance has a specific extinction coefficient. As opposed to other eutrophication models the extinction due to DOC is modeled explicitly.

Decomposition of Detrital Organic Matter
The decomposition of detrital organic matter is innovatively formulated as the mineralization and conversion of five fractions, allowing for the simulation of the entire organic matter cycle in surface water and bed sediment. POC1 is the fast decomposing detritus fraction, POC2 the moderately slow decomposing fraction, POC3 the slow decomposing fraction and POC4 the very slow decomposing (refractory) fraction. POC1 and POC2 are the dominant fractions in the water column, whereas POC3 and POC4 are the dominant fractions in the sediment. DOC represents very slow decomposing (refractory) dissolved organic matter being produced in both water column and sediment. A mineralization flux may have proportional conversion fluxes, which represent the gradual turnover of detritus into residual refractory organic matter due to enzymatic stripping and polymerization. In the Lake Veluwe model POC/ N/P/S1 is converted into POC/N/P/S2, POC/N/P/S2 into POC/N/P/S3 and DOC/N/P/S, and POC/N/P3 into POC/ N/P4 and DOC/N/P/S. This decomposition and conversion scheme is an extension of the multi-G model [32] applied in various forms in many models [2,17,20,30,36]. Raick et al. [37], Lancelot et al. [8] and Dale et al. [38] describe different versions of an alternative approach based on the hydrolysis of particulate and dissolved polymeric detritus fractions into monomeric dissolved substances that are taken up and respired by bacteria. This delivers essentially similar results as the approach in ECO, because both concepts have the same rate limiting first decomposition step.
The mineralization rate in ECO conforms to first-order kinetics [25,32,39]. The rate is dependent on the availability of the nutrients (N, P) needed by detritus decomposing bacteria [24], and on the dominant electron acceptor [19,40] (Eqs. B.1, B.3, B.4). The latter dependency was ignored in the Lake Veluwe model. Mineralization in ECO implies the release of ammonium (NH4), phosphate (PO4) and sulphide (SUD), the consumption of electron-acceptors (DO, NO3, SO4) and/or methanogenesis (CH4). Preferential stripping of nutrients from detritus is taken into account by means of an acceleration factor for organic nutrients (Eq. B.2), which is proportional to the difference of actual and objective detritus nutrient contents [24,26]. The objective contents are the low nutrient contents of humic and fulvic matter, the residual products of the decomposition process. The conversion of detritus is proportional to mineralization with ratios specific for each of the detritus fractions (Eq. B.5).

Consumption of Electron Acceptors and Methanogenesis
The dominant electron acceptors consumed for the oxidationmineralization of organic matter are oxygen (DO), nitrate (NO3), sulphate (SO4) and carbon monoxide. In the model, carbon monoxide and hydrogen arising from the fermentation of organic matter at chemically reducing conditions are implicit in the decomposition product methane. The electron acceptors iron(III) and manganese(IV) are assumed to be implicit in the sulphur redox cycle. This can be justified from the poor solubility of their oxy-hydroxides, and the coupling of the metal and sulphur cycles as described by Wang & Van Cappellen [14] and Wijsman et al. [18]. As ruled by redox potential dependent energy yields, different groups of bacteria species consume the electron acceptors in a specific sequence [12,13,41]. When dissolved oxygen is not available, bacteria resort to nitrate for the decomposition of organic matter. If both oxygen and nitrate are depleted, sulphate is used by other bacteria. At the depletion of all three substances D. Sulphur and methane processes E. Exchange of gases with the atmosphere F. Settling and mass transport in the sediment settling G. Extinction of light e t~eb ze poc ze doc ze alg ze im SD~1:7=e t (G.1) where : i = algae types 1-14. j = particulate detrital organic matter fractions 1-4, and dissolved fraction 5. l = sediment layer 1-7, overlying water layer 0 (L 0 = 0.0005 m). OX = particulate (POX 1-4 ) or dissolved (DOX) detrital fraction; X is C, N, P or S. X = any particulate substance; ALGi, POX, AAP, VIVP, APATP, OPAL, IM1, IM2. Y = any dissolved substance. doi:10.1371/journal.pone.0068104.t002 Table 3. Explanation of the symbols in Tables 1-2, with values of the process coefficients in ECO (part 1). methanogenic bacteria convert organic matter into approximately equal amounts of methane and carbon dioxide. The consumption of an electron acceptor is limited by availability and inhibited by the presence of electron acceptors that yield more energy. However, because both sediment and water compartments are not homogenously mixed in the real world, the electron acceptors may be depleted in parts of the compartments, implying that they can be consumed simultaneously in each model compartment. The sequential but spatially overlapping consumption of electron acceptors is reproduced in ECO with Michaelis-Menten kinetics for limitation and inhibition [12,13,[17][18][19]20,22,30], which is a step forward with regard to the formulation of sediment diagenesis in sediment module SWITCH [25] of DBS and in other eutrophication models. Simplified limitation and inhibition functions were described by Aguilera et al. [42] and applied by Brigolin et al. [2]. Dale et al. [38] presented an alternative approach based on double Michaelis-Menten kinetics for electron donors and acceptors and Gibbs free energy yields. In ECO denitrification is inhibited by oxygen, sulphate reduction by nitrate and oxygen, and methanogenesis by sulphate, nitrate and oxygen (Eqs. B.6-11). The Michaelis-Menten factors are scaled to make the sum of the factors equal 1.0. The total consumption rate of electron acceptors is proportional to the mineralization rate of organic matter (Eqs. A3, A12, A23, A16).

Nitrification
Ammonium (NH4) is oxidized to nitrate (NO3) by means of stepwise nitrification, the oxidation of nitrite being the last and rate limiting step [13,18,19]. Because nitrifying bacteria need readily available organic substrates, nitrification proceeds most rapidly in the oxic top sediment layer. Nitrification is often described with first-order kinetics with regard to ammonium [25,26,43] or to both ammonium and oxygen [2,17,30]. In ECO nitrification is formulated according to Michaelis-Menten kinetics with respect to ammonium and dissolved oxygen (Eq. C.1).

Phosphate Processes
The adsorption of phosphate (PO4) to iron(III) oxyhydroxides and other inorganic components of sediment is a pH dependent equilibrium process [44]. Desorption occurs when pH increases. The adsorption capacity of sediment depends on the presence of oxygen. At anoxic conditions the capacity is far smaller, mainly due to the reduction of iron(III) oxyhydroxides, but some adsorption capacity is maintained in sediment [9,23,28,45,46]. The adsorption of phosphate can be described in various ways with linear or Langmuir equilibrium sorption [25,26,28,47]. In ECO, phosphate adsorption is formulated according to a pH dependent Langmuir model according to Eqs. C.2-5 [44]. The pH is a forcing function. Innovatively, the adsorption capacity is deduced from the reactive iron content of the sediment corrected for the oxidized iron(III) fraction. The correction factor is equal to 1.0 when the oxygen concentration is above a critical value, but smaller at lower oxygen concentrations. The sorption flux is calculated from the difference of the actual and equilibrium adsorbed phosphate concentrations multiplied with a sorption rate. The equilibrium adsorbed phosphate concentration follows from Equation C.3 (Table 1).
In the sediment vivianite (iron(II) phosphate) or apatite-(calcium phosphate)-like minerals may precipitate at the supersaturation of pore water [48]. Boers & de Bles [49] studied the possible precipitation of vivianite in Lake Loosdrecht (peat lake), but did not find conclusive evidence. However, Dittrich et al. [19] included vivianite in their sediment diagenesis model for Lake Zug. The formation of Ca-P minerals in marine sediment has been stipulated by De Jonge [50], Slomp et al. [23], Slomp & Van Cappellen [51] and Anschutz et al. [45]. However, the coprecipitation of phosphate with calcite is not always found for marine sediment [28]. Unlike vivianite that is only stable under anoxic conditions, the stability of calcium phosphates is not affected by the redox potential. Therefore, the precipitation of apatite-like minerals may lead to the more or less permanent storage of phosphate in sediment [23]. The mechanisms and controlling factors of P-mineral formation are not well understood [46]. Generally, the precipitation and dissolution rates are modeled proportional to the difference of an ion activity product and the   solubility product, whereas the dissolution rate is also proportional to the concentration of the mineral [12][13][14]18,19]. In ECO the precipitation of phosphate is formulated in a simplified way with first-order kinetics [23,25], the difference between the actual and equilibrium dissolved phosphate concentrations being the driving force (Eq. C.6). This is justified by the assumption that the pH in the sediment and the pore water concentration profiles of iron(II) and calcium are rather constant, although iron(II) varies in the oxidizing-reducing transition zone. Unlike other eutrophication models ECO distinguishes between vivianite-P and apatite-P. In ECO, the precipitation of vivianite takes place only when the average dissolved oxygen concentration is below a critically low value (0.25 mg O 2 /L). The dissolution rate of the apatite-like mineral depends on the extent of undersaturation as well as the concentration of the mineral (Eq. C.7). The dissolution of the vivianite-like mineral is proportional to the concentrations of the mineral and dissolved oxygen.

Dissolution of Silicate
In many models the dissolution of opal silicate (OPAL) is dealt with as first-order degradation [8,24,37]. However, this is a chemical-physical process proceeding in an undersaturated solution [11,52,53]. In ECO, OPAL is slowly converted into dissolved silicate (Si) proportional to the concentration of OPAL and the difference between the saturation and actual dissolved silicate concentrations (Eq. C.8).

Sulphide Processes
The oxidation of sulphide is established by fast chemical reaction and slow microbial conversion [14,18]. The oxidation of dissolved sulphide is modeled as a first-order process with regard to both total dissolved sulphide (SUD) and dissolved oxygen (DO) [17,19,30]. At reducing conditions sulphide may precipitate as amorphous iron(II) sulphide, which is thermodynamically unstable and dissolves at oxidizing conditions. The precipitation and dissolution rates are usually formulated on the basis of the difference of the ion activity and solubility products [13,14,19,44], but also with first order kinetics with respect to dissolved Fe 2+ and H 2 S [15,30]. In ECO however, the precipitation of sulphide (SUD) and the dissolution of particulate sulphide (SUP) is described with first-order kinetics (Eqs. D. [2][3][4]. The difference between the actual and equilibrium free sulphide concentrations is the driving force. It is assumed that sulphide forming metals like iron(II) are abundantly present in the sediment. The free dissolved sulphide concentration (S 22 ) is calculated from the total dissolved sulphide concentration and the pH. The chemical equilibrium reactions and pertinent equilibrium constants were derived from Stumm & Morgan [44].

Methane Oxidation
Bacteria oxidize methane with dissolved oxygen (DO) or sulphate (SO4). Often first-order kinetics with regard to both methane and oxygen are applied [14,17,18,30]. In ECO, methane oxidation is modeled according to Michaelis-Menten kinetics with regard to both methane and the oxidizing agent (Eqs. D.5-6), whereas the oxidation with dissolved oxygen excludes the oxidation with sulphate [13].

Exchange of Gases with the Atmosphere
The formulation of the diffusion flux of oxygen (DO: reaeration) and methane (CH4: volatilization) across the water surface is based on the double film concept. The flux is the temperature dependent mass transfer coefficient in water multiplied with the difference of the saturation and actual concentrations, the concept of which was described by Wanninkhof [54]. In the Lake Veluwe model, the wind speed dependent formulation for the transfer coefficient according to Banks & Herrera [55] is used (Eqs. E.1-3). The saturation concentration of oxygen is calculated according to the temperature and chlorinity dependent formulation described by Weiss [56]. Additionally, methane may escape to the atmosphere via gas bubbles, which are formed in the sediment at the supersaturation of dissolved methane. Similar to the approach of Canavan et al. [20], it is assumed that all methane formed at supersaturation is removed instantly (Eq. E.4). The saturation concentration depends on pressure (water depth) and follows from the formulation described by Di Toro [13].

Settling and Mass Transport in the Sediment
The settling flux of a particulate component is proportional to its concentration and settling velocity (Eq. F.1). In the Lake Veluwe model we consider time average net settling (no resuspension), and the settling velocities for algae are equal to 0.0.
Mass transport in the sediment results from advection and dispersion [11,13]. The dispersion and advection of dissolved substances across the sediment-water interface imply return fluxes of nutrients to the water column and a sediment oxygen demand. Advection arises from upward or downward seepage of solutes, and from settling and resuspension of particulates. Settling leads to burial, resuspension to upward transport. The latter requires that sediment and water quality at the lower boundary of the sediment are defined. In the Lake Veluwe model burial results from settling at constant porosity and constant sediment layer thickness. Burial and seepage fluxes follow from Equations F.2-3. Dispersion is caused by seasonally varying activity of benthic organisms resulting in bioturbation of particulates and bio-irrigation of solutes. Solutes are also transported due to flow induced micro turbulence and molecular diffusion [57]. In ECO, dispersion is calculated as diffusive transport according to Fick's second law [11] (Eqs. F. [4][5]. The dispersion coefficient of solutes is the sum of the bio-irrigation coefficient, the flow induced dispersion coefficient and the molecular diffusion coefficient. Except for the latter which is corrected for tortuosity, all coefficients decrease exponentially over depth in the Lake Veluwe model.  The temperature constants for the decomposition of organic matter and other microbial processes are 1.047 and 1.07, respectively. Apart from the dissolution of opal silicate (kt = 1.07), the temperature constants of precipitation and dissolution processes are equal to 1.0. Stochiometric constants for the various redox processes are specified in Equations A.3, A.12-14 and A. 16. The stochiometric constants s x,i for algae biomass in Equations A.2-6 and A10-11 are given in Table 5.

The Calibration of ECO
The Lake Veluwe case Lake Veluwe (The Netherlands) is a shallow lake which was formed when the Flevoland polder was created in 1956 ( Figure 2; [58]). The lake has a total surface area of 32.3 km 2 . 60% of the lake has a sandy sediment with approximately 1.5% calcite at an average depth of 0.75 m. The other 40% has a silty sediment with approximately 8% calcite at an average depth of 2.5 m [59]. Bordering higher land, the shallow part has seepage. The deeper part has infiltration due to the low Flevoland polder in the northwest [60]. In the second half of the 1960s the water quality of the lake deteriorated, because of increased nutrient loading. As of 1970, this resulted in high chlorophyll-a concentrations and an almost permanent bloom of cyanobacteria, predominantly Oscillatoria agardhii. From February 1979 onwards, external phosphorus loading was reduced from 3 to 1-1.5 g P.m 22 .yr 21 by means of phosphorus removal at the sewage treatment plant discharging its effluent to the lake. Additionally, the lake was flushed in winter with water from the Flevoland polder (poor in algae and phosphorus, but rich in calcium and nitrate), which reduced the average winter half year residence time of dissolved compounds from 0.35 to 0.15 year. In the second half of the 1980s, the lake was flushed in summer as well, reducing the average summer half year residence time of dissolved compounds from 0.5 to 0.25 year. Consequently, total phosphorus and chlorophyll-a concentrations decreased as of 1980, reaching strongly reduced levels in the second half of the 1980s. Whereas diatoms and green algae became dominant over cyanobacteria, submerged macrophytes began to re-appear as of 1985 [60,61].

Model Input, Schematization and Forcing
The DBS model of Lake Veluwe had one water compartment and 4 sediment layers [25,26]. For the ECO model, Lake Veluwe was divided into a ''shallow'' compartment and a ''deep'' compartment, that have the above mentioned depths and areas, and different net settling velocities and sediment compositions. Each compartment has 10 sediment layers, the thicknesses of which increase with sediment depth: 460.1 cm (1-4), 0.2 cm (5), 0.4 cm (6), 1.0 cm (7), 2.0 cm (8), 6.0 cm (9), 10.0 cm (10). Layer 10 serves as a buffer layer. Sediment porosity and the initial composition of the sediment with regard to organic matter and inorganic P were based on measurement data for 1979-1981 [59,60]. The porosities for the sandy sediment in the shallow compartment and the silty sediment in the deep compartment are 0.4 and 0.7, respectively. Starting from measured organic matter and total phosphorus contents, the initial sediment composition was adjusted until the model showed a stable composition. The initial composition of the water column followed from the fortnightly routine monitoring data from Rijkswaterstaat, that were also used for the calibration of the model. Inflows and outflows of water were derived from ten days water balances for 1976-1992 based on constant water level (unpublished data from Rijkswaterstaat). The water balance includes inflows from adjacent Lake Wolderwijd, the Flevoland polder, several brooks, the Harderwijk sewage treatment plant, seepage and precipitation. Outflows concern flow to adjacent Lake Drontermeer, infiltration and evaporation. The two compartments exchange water and substances by means of a high dispersion coefficient that warrants very small concentration differences between the two water compartments. Loads of chloride, phosphorus, nitrogen, dissolved silicate and organic matter for each of the inflows are based on observed concentrations (unpublished data from the district water boards Vallei en Veluwe and Zuiderzeeland, made available by Rijkswaterstaat). The sulphate loads on the model were determined proportional to those of chloride. The model is forced with weekly data for solar radiation (PAR; W.m 22 ) and wind speed (m.s 21 ) (data from Royal Netherlands Meteorological Institute), and fortnightly data for water temperature ( o C), suspended inorganic sediment (IM3; g.m 23 ) and pH (monitoring data from Rijkswaterstaat).
All model input data is available on the OpenEarth model repository. OpenEarth is a free and open source initiative to deal with Data, Models and Tools in water related science & engineering projects [62]. To access the data one needs to sign up on the OpenEarth internet site: http://publicwiki.deltares.nl/ display/OET/OpenEarth/. The internet address to the model data is: https://svn.oss.deltares.nl/repos/openearthmodels/ trunk/deltares/lake_veluwe_eco/.

Calibration Procedure
Ethics statement: N/A as to human research, clinical trials, animal research, cell line research, and observational and field studies, because these have not been carried out for this study. No permission for field experiments was required because for this study no field experiments have been carried out. The water quality and water balance data used were collected within the framework of the regular monitoring program of Rijkswaterstaat, who is the management authority for national water systems in the Netherlands such as Lake Veluwe. Many of ECO's process coefficients were calibrated for the Lake Veluwe case, whereas stochiometric ratios were fixed at values derived from literature. Initial values were taken from calibrated previous versions of ECO, GEM [24] and DELWAQ-BLOOM-SWITCH [25,26], or from other reported models. The calibration was done by adjusting the values of process coefficients and by adopting or rejecting new values depending on the improvement or deterioration of model performance, which was assessed by graphical comparison of predicted and observed water and sediment quality parameters. Additional criteria were used for individual coefficients as described below for the resulting coefficients. The calibration was carried out in a step-wise and iterative manner, process by process, going from processes affecting independent state variables to processes affecting dependent state variables, starting from the coefficients for which the model appeared most sensitive. Given the many interactions of substances and processes in the model, this procedure was repeated until no further improvement was observed.

Resulting Process Coefficients (Tables 3, 4)
In general the individual process coefficients of a calibrated water quality model should have values that are in line with values reported for other similar models. This holds in particular for generic multi coefficient models such as ECO. In view of the generic nature of ECO, the process coefficient values resulting from the calibration are compared with ranges of values reported for other models, both freshwater and saline water models because of the relative scarcity of data for freshwater models.
We calibrated the rates of nitrification (Eq. C.1) to reproduce observed nitrate and ammonium concentrations and to obtain an appropriate denitrification flux in the sediment. ECO's nitrification rate for the sediment is within the ranges found in the literature, whereas its value for the water column is somewhat higher. These ranges are for marine water k ni = 0.014-0.07 g N.m 23 .d 21 , for marine sediment 2.56-20.0 g N.m 23 .d 21 [1,2,4,8,17,22,24,29,30,36,37,63], for fresh water k ni = 0.01 g N.m 23 .d 21 , for fresh sediment 9.79-22.4 g N.m 23 .d 21 [6,19,20]. To allow comparison with first order, second order and single Michaelis-Menten rates we re-calculated these rates on the basis of temperature (see above; kt = 1.07) and representative Lake Veluwe concentrations of ammonium (1.0 mg N. L 21 in sediment and water) and oxygen (4 mg O 2 .L 21 in sediment). The ranges of half saturation concentrations for nitrification found in literature are: for marine water Ks an = 0.07-0.28 mg N.L 21 and Ks do = 0.16-0.032 mg O 2 .L 21 [1,8,22,29,36,63] for fresh water Ks an = 0.7 mg N.L 21 and Ks do = 0.64 mg O 2 .L 21 [19]. ECO's values for these coefficients are within the overall range (Ks an ) or slightly higher (Ks do ).
We found few data in the literature for the coefficients of phosphate processes (Eqs. C.2-7). ECO's adsorption coefficients were calibrated to reproduce summer dissolved P concentrations in the water column. The adsorption capacity of the sediment must be large to store enough adsorbed phosphate in the top layer to support the high return fluxes in the summers of 1976-1980. The fraction reactive iron fe i in the suspended sediment and the bed sediment was derived from measurement data [59]. The bulk of the sediment has a reactive iron fraction of 0.025 g Fe.g IM1 21 , the top 4 mm a fraction of 0.075 g Fe.g IM1 21 assuming accumulation due to diagenetic processes: Iron is reduced below the top layer and precipitates in the top layer. The fraction of oxidized reactive iron at reducing conditions f ox is 0.2 for the top 4 mm and 0.1 in the deeper sediment. Anschutz et al. [45] found for lagoon sediment, that due to iron reduction the reactive iron(III) concentration decreases to approximately 50% at a depth of 5 cm and to approximately 25% at a depth of 10 cm. ECO's lower values of f ox for Lake Veluwe sediment are justified by more intensive redox processes due to a much higher organic matter load.
The precipitation of vivianite and apatite-like phosphate minerals is dependent on local chemical conditions. Dittrich et al. [19] determined that the vivianite precipitation rate in a fresh sediment was given by 0.0096 (IAP/SOL-1) g.m 23 .d 21 . This can Figure 6. Simulated sediment quality parameters in Lake Veluwe. The parameters represent the bulk concentrations of particulate organic carbon and total phosphorus, and the pore water concentration of dissolved inorganic phosphorus in the sediment of Lake Veluwe for the period 1976-1985. The silty sediment is in the ''deep'' compartment, the sandy sediment is in the shallow compartment. doi:10.1371/journal.pone.0068104.g006 Table 6. Simulated and observed total phosphorus and organic carbon contents (g.kg DM 21 ) in sediment. North Atlantic platform sediment, and an equilibrium concentration PO4 ae = 0.11 mg P.L 21 . Anschutz et al. [45] reported dissolved P concentrations in marine pore water between 0.05 and 0.25 mg P.L 21 at the precipitation of apatite. From the data presented by Mort el al. [46] it can be inferred that PO4 ae for authigenic Ca-P minerals in the pore water of Baltic Sea sediment could be as high as 4.5 mg P.L 21 . We selected the equilibrium concentrations PO4 ve and PO4 ae applied in ECO in line with concentrations observed for pore water in Lake Veluwe sediment [25]. These agree well with the reported ranges. We adjusted the dissolution rates to keep predicted pore water concentrations in line with the observed concentrations. The values for the opal dissolution rate found in literature all concerned first order rates for marine systems. For comparison with ECO's calibrated rate (Eq. C.8) we re-calculated them for reference temperature 20 o C (see above) and a dissolved silicate concentration 1.0 mg Si.L 21 : for water k od = 0.000625-0.005 (g Si.m 23 ) 21 .d 21 , for sediment k od = 0.0005-0.00048 (g Si.m 23 ) 21 .d 21 [8,24,37]. ECO's rate calibrated for fresh sediment is significantly lower. The equilibrium dissolved concentration Si e in ECO followed from maximum observed pore water concentrations in the sediment of Dutch shallow lakes [64].
Several publications provide data on the oxidation rate of sulphide (Eq. D.1): for marine water k so = 51.3 (mg O 2 .L 21 ) 21 .d 21 , for marine sediment k so = 2.06-23.0 (mg O 2 .L 21 ) 21 .d 21 [1,2,17,22,29,30,63], for fresh sediment k so = 85.0 (mg O 2 .L 21 ) 21 .d 21 [20]. Because of insufficient data for fresh sediment we selected ECO's oxidation rate of sulphide as an approximate average of the range for marine sediment. We gave the precipitation and dissolution rates of particulate sulphide k sp and k sd in ECO (Eqs. D.2-4) high values to ensure that predicted dissolved sulphide stays close to the equilibrium concentration Cs e . We derived the equilibrium constants for the sulphide speciation in solution from Stumm & Morgan [44].
Low values for the net settling velocities of particulate matter in several shallow lakes in the Netherlands that are exposed to rather continuous windy conditions resulted from unpublished model calibrations. However, the net settling velocities v x (Eq. F.1) recalibrated by us for very shallow Lake Veluwe are slightly larger than reported previously [26], respectively 0.23 m/day and 0.115 m/day for the deep and the shallow compartments. The difference is justified by less resuspension occurring in the deep compartment.
We derived the seepage and infiltration velocities (Eq. F.1) from an unpublished hydrological study for Lake Veluwe [26,60]. v spg is 20.015 and 0.0057 m.d 21 in the relatively deep and shallow compartments, respectively.
The dispersion coefficients for solutes and particulates in the sediment D s and D p (Eqs. F. [4][5] have maximal values at the sediment-water interface and decrease exponentially with sediment depth. The present values in ECO for Lake Veluwe are The exponential decrease implies a decrease to 10% at a sediment depth of 5 cm, below which the coefficient is almost equal to the corrected molecular diffusion coefficient. These values agree with values found in the literature that imply a maximal bioirrigation enhancement factor 5 [25,29,36]. The bioturbation dispersion coefficient D p near the sediment-water interface in ECO varies linearly between 2.8610 26 m 2 .d 21 (June-September) and 2.8610 27 m 2 .d 21 (January-February). Exponential decrease implies that D p decreases to 10% at a sediment depth of 4 cm and is virtually equal to zero at a depth of 10 cm. This agrees with values reported in literature for other models: D p = 0-2.0 10 26 m 2 .d 21 (seasonal variation fresh water system; [20,25]), D p = 2.36-8.0 10 26 m 2 .d 21 (marine water; [17,29,38]).
The extinction coefficients for particulate substances in ECO

Water Quality and Phytoplankton
The simulated and observed concentrations indicated in Figure 3 (phytoplankton, organic matter, secchi depth, oxygen, chloride and sulphate) and Figure 4 (nutrients) represent the average water quality of the deep and shallow compartments. Substantially more parameters are shown than for the DBS Lake Veluwe model [26], which is possible due to ECO's comprehensiveness. The observed concentrations represent the data collected for two sampling locations, one in each compartment.
The results for chloride show that the water balances imposed on the model for 1978-1985 are quite accurate. The rather large differences of simulated and observed chloride suggest that the water balances for 1976 and 1977 are not very accurate, which has a bearing on the accuracy of the nutrient loads. No data are available for sulphate in Lake Veluwe, but the predicted ratio of sulphate and chloride is in line with data for other similar lakes in the Netherlands.
Simulated and observed chlorophyll-a match quite well for most of the simulated period. The steep transition caused by the reduction of the nutrient loads and flushing is reproduced by the model. However, chlorophyll-a is overpredicted for the spring peaks in 1976-78, whereas the summer peak in 1979 and the spring peak in 1981 are underpredicted. Both underpredictions are essentially caused by P-limitation as appears from the lack of dissolved phosphate in the water column during these periods. This points at temporary underprediction of the internal P-load (Preturn flux from the sediment) and/or the temporary underestimation of the external P-load. Structural overprediction is found for 1985, most probably due to the absence in the model of macrophytes and microphytobenthos that began to re-appear in Lake Veluwe in 1985 due to its much enhanced transparency [60,61]. Nutrient uptake by macrophytes and microphytobenthos decreases the total nitrogen and total phosphorus concentrations in the lake, so that less nutrients are available to phytoplankton and phytoplankton biomass (chlorophyll-a) is reduced.
The simulated phytoplankton species composition (Fig. 5) is very similar to the composition simulated with the model described by Van der Molen et al. [26], which also contained BLOOM. The dominance of cyanobacteria(mainly Oscillatoria agardhii) in 1976-1979 changes into dominance of diatoms in spring and of green algae in most of the summer as of 1980. Cyanobacteria are almost entirely absent as of 1985. Van der Molen et al. [26] stipulated that this agreed with observed species composition.
Total organic carbon (TOC) is the sum of the carbon stored in algae biomass, the particulate detritus fractions (POC1-4) and dissolved organic carbon (DOC). The simulation results and the observed concentrations for TOC, DOC and dissolved oxygen agree quite well, which supports the claim that primary production by algae is simulated rather accurately by the model. Over-and underpredictions are consistent with over-and underpredictions of chlorophyll-a. The model reproduces observed low secchi depth in 1976-1979 with deviations that are also consistent with this. Simulated secchi depth agrees quite well with the observed increase of the secchi depth in 1980-1985. The light limitation of phytoplankton that occurred in the deep compartment of Lake Veluwe in 1976-1978 was replaced by P-limitation except for periods in winter.
Whereas simulated total phosphorus agrees well with measured total phosphorus as of 1980, summer peaks are underpredicted for 1976-1979. This may be due to underpredicted P-return fluxes from the sediment, as dissolved phosphate is similarly underpredicted, but the phosphorus loads may have been underestimated too. Possibly not enough adsorbed phosphate is stored in the top sediment in the model, or adsorbed phosphate is not released fast enough during the summers of 1977-1979. We could not tune the adsorption related coefficients of the model such that underprediction in 1977-1979 was removed and that overprediction as of 1980 was prevented at the same time. Enhanced precipitation of phosphate minerals in the sediment, the apatitelike mineral in particular, was needed to keep dissolved phosphate as low as observed and to bring about sufficient P-limitation of the algae. Trial simulations showed overprediction of dissolved phosphate as of 1980 when we decreased the ratio of apatite over vivianite precipitation r av or the precipitation rate k vp . Although no specific evidence exists for the presence of an apatite-like mineral in the sediment of Lake Veluwe, the assumption of the formation of an authigenic apatite-like mineral is justifiable by the presence of calcite in this sediment and the possible coprecipitation of apatite with calcite.
Some of the summer dips and the winter peaks of total nitrogen are overpredicted by the model. The overprediction of winter peaks shows in nitrate as well. Nitrogen return fluxes might be overestimated due to insufficient denitrification, but nitrogen loads may have been overestimated too. Trial simulations pointed out that enhanced nitrification and denitrification caused under-prediction of summer chlorophyll-a in 1977-1979, because phytoplankton is N-limited in summer and most of the autumn.
The observed dips and peaks of dissolved silicate as resulting from the spring diatom blooms are reproduced by the model, which supports the validity of the predicted algae species succession.

Sediment and Pore Water Quality
Simulated concentrations of POC, total P, and dissolved inorganic P in the sediment of Lake Veluwe show accumulation in the hypertrophic period 1976-1979 (Fig. 6). Whereas ECO provides much more vertical detail than the SWITCH module of DBS [25], a selection of the results for three or four out of ten layers is presented. After the large decrease of the P-load on the lake in 1979, this is followed by a decline of concentrations in the sediment during 1980-1985, which for the top sediment of the deep compartment is mainly due to the bioturbation of particulate phosphate into deeper sediment and the loss of dissolved phosphate into the groundwater below the ''deep'' compartment (see also [26,60]). The decline of total phosphorus in the top sediment of the shallow compartment is much weaker due to the import of phosphorus from upwelling groundwater. As can been seen from Table 6, the simulated total phosphorus and organic carbon contents in the sediment agree quite well with observed contents, that were derived from Brinkman & Van Raaphorst [59] (see [25]). We derived simulated contents from bulk concentrations by division by (1-porosity) and sediment density (2.6 kg.L 21 ). The concentrations predicted by the model for the individual phosphorus components point out that at the end of the 1976-1985 simulation approximately 40% of total P in the sediment is stored in an apatite-like mineral in both compartments. Observed dissolved phosphorus in pore water of the upper 10 cm sediment for 1983 shows rather irregular patterns, but a distinct maximum occurred between 1-2 cm sediment depth and minimum values were found below 3 cm depth [59]. Concentrations are higher in summer than in spring. These trends are reproduced by the model. The mean observed concentration is 0.198 g P.m 23 (st. dev. = 0.192 g P.m 23 , n = 94). The simulated average dissolved inorganic phosphorus concentrations (PO4) in the top 10 cm of sediment for 1983 vary around 0.17 g P.m 23 in the deep compartment, and around 0.14 g P.m 23 in the shallow compartment. Simulated dissolved organic phosphorus varies around 0.035 g P.m 23 .
The simulated pore water concentrations of the electronacceptors dissolved oxygen, nitrate, and sulphate in the sediment of Lake Veluwe clearly show seasonal variation ( Fig. 7; sandy sediment). The oxic layer is much thinner in the summer half year than in the winter half year, when the decomposition of organic matter is slow. During the summers of 1976-78 the oxic top layer becomes approximately 1 mm thick in the model, which implies the collapse of the adsorption capacity for phosphate and the strong increase of dissolved phosphorus. A structural change occurs in 1979, when the penetration of dissolved oxygen, nitrate and sulphate into the sediment in summer increases and methane decreases because of highly reduced detritus settling on the sediment, which leads to the enhanced adsorption of phosphate and strongly decreased dissolved phosphorus.
Typical summer situations for the simulated fractional contributions of oxygen consumption, denitrification, sulphate reduction and methanogenesis to the mineralization of organic matter in the sediment of Lake Veluwe are shown in Figure 8 for sandy sediment. The contributions are partially overlapping, with a near 100% contribution of oxygen in the upper sediment layer (1 mm) and a near 100% contribution of methanogenesis in the lower sediment layer (4-10 cm). The contributions of oxygen consumption and denitrification increase significantly after 1979, whereas the contribution of sulphate reduction decreases strongly.

Mass Balances and Return Fluxes
The mass balances for phosphorus and nitrogen of Lake Veluwe for 1978 and 1983 that result from the model are displayed in Figure 9. The total external load is the sum of the external load (from surface water and point sources), the seepage load and atmospheric deposition. The internal load is the diffusive return flux of ammonium and nitrate or of phosphate from the sediment to the overlying water. Outflow is the net transport of nutrients to adjacent Lake Dronten north east of Lake Veluwe. Infiltration in the ''deep'' compartment and seepage in the ''shallow'' compartment represent the transport to and from groundwater. Burial is the removal of substances into the sediment below the simulated layers resulting from the net settling of sediment. Storage is the difference between initial and final composition, which in a ''dynamic steady state'' should be equal to zero. Positive storage and burial imply accumulation in the sediment.
The mass balance data produced by ECO are more detailed compared to the DBS data [26], and allow for a more detailed assessment of the various nutrient fluxes and removal processes. As appears from Figure 9, from 1978 to 1983 the total external load of phosphorus dropped 68%, whereas the total external load of nitrogen increased 38% due to flushing with nitrogen rich water. The internal load of phosphorus decreased only 15%, and becomes much larger than the external load in 1983, despite the enhanced adsorption of phosphorus in the top sediment due to the strongly reduced input of detritus. This shows that internal eutrophication can continue years after an external load reduction. The average simulated summer half year return flux of dissolved phosphate for 1978 is 11.6 g P.m 2 .d 21 for the silty sediment and 7.5 g P.m 2 .d 21 for the sandy sediment. For 1983 these fluxes are 6.2 g P.m 2 .d 21 and 5.9 g P.m 2 .d 21 , respectively. The mean phosphate return flux observed in-situ with a benthic chamber in the summer half year of 1983 is 4 mg P.m 2 .d 21 (s.d. = 5, n = 2) for the silty sediment [59,25]. The predicted phosphate return fluxes for 1983 have similar magnitudes as the return fluxes observed for the silty sediment.
The internal load of nitrogen is 87% of the total external load in 1978. It becomes much smaller and drops to only 18% of the total external load in 1983 because coupled nitrification-denitrification is enhanced strongly in a more oxidizing top sediment layer. The denitrification flux becomes larger but remains approximately 50% of the total external load at a strongly reduced residence time.
In 1978 the net outflow of phosphorus is only 3% of the total external phosphorus load, the accumulation in the sediment 76%, and loss by infiltration into groundwater 14%. Due to flushing the net outflow increases to 36% of the total external load in 1983. The accumulation is down to 22%, whereas infiltration increases to 37%. Similarly the net outflow of nitrogen increases from 4% to 35% of the total external loads. The loss of nitrogen by infiltration goes down from 38% to 19%, whereas accumulation is unimportant. The large removal fluxes of nutrients by infiltration indicate its importance for the oligotrophication of Lake Veluwe. With respect to phosphorus, flushing was equally important.

Discussion and Conclusions
The present article describes the content and the calibration of the eutrophication model ECO as based on DELWAQ that combines a number of significant modeling innovations. Com-pared to other eutrophication models, DBS [26] included, ECO contains the following innovations: N Comprehensive water and sediment quality are simulated in a generic way on the basis of a computational grid covering both water column and sediment bed. Process formulations are deterministic and identical for water column and sediment bed. Chemical conditions determine how formulations turn out locally.  N The mass balances for organic carbon, nutrients (N, P, Si), sulphur and oxygen are fully closed for water and sediment. All relevant organic and inorganic components are included.
N Covering the entire organic matter cycle in water and sediment, detrital organic matter is simulated with four particulate fractions and one dissolved fraction, having different decomposition rates, the one fraction being produced from the other. N Light extinction is calculated from contributions of all major components, including dissolved organic matter.
N The 1D, 2D or 3D computational grid for water and sediment is defined by input. Number and thickness of layers are only limited by the maximum acceptable computational burden.
The process coefficients that resulted from the calibration of the Lake Veluwe model are in line with the ranges of values that were reported for other models. The simulation results show ECO's capacity to simulate water and sediment quality and sedimentwater exchange fluxes dynamically and quite accurately. The results underline the importance of redox processes and phosphate speciation for the nutrient return fluxes.
Differences between simulated and observed water quality mainly concern the timing and the magnitude of concentration peaks and dips. Parts of these differences are due to inevitable inaccuracies in model forcing, in water and nutrient loads and in suspended inorganic sediment in particular. The tendency of the model to underestimate the variability of concentrations is partially caused by the application of net settling in stead of alternating settling and resuspension. However, the stochastic patchiness of lake water quality cannot be reproduced by a model with two compartments and forcing that is averaged over time and space.
The underprediction of the nutrient return fluxes in Lake Veluwe in the summers of 1976-1979 may be due to the set-up of the present model. The simulation shows that the oxic top layer of the sediment can become very thin in the model, but can never disappear when the overlying water still contains substantial concentrations of dissolved oxygen. Near bottom thermal stratification may have occurred in the ''deep'' compartment during summer periods with low wind speed, leading to near zero oxygen concentrations just above the sediment, and to stronger chemical reduction, further loss of adsorption capacity and reduced nitrification-denitrification in the upper sediment layer. The present model with fully mixed water column cannot generate these low oxygen concentrations, and therefore tends to underpredict the phosphate and ammonium return fluxes in summer. The near bottom anoxic conditions due to stratification might be reproduced with a 3D model based on stratified water flow.
The overprediction of nutrients and phytoplankton for 1985 is probably due to the absence of submerged macrophytes in the model. These macrophytes that re-appeared in substantial quantities in Lake Veluwe as of 1985 store significant amounts of nutrients that are therefore not available for phytoplankton. The growth of microphytobenthos also not included in the model may have contributed too to observed lower nutrient concentrations in the water column. The overprediction of phytoplankton in Lake Veluwe for 1985 and years after (not shown here) shows that to more accurately simulate phytoplankton and nutrients in shallow lakes with nutrient limitation and a large population of macrophytes these macrophytes and possibly also microphytobenthos should be included in a eutrophication model.
The simulation of the water quality of Lake Veluwe with comprehensive ECO provides additional insight in how lakes may respond to nutrient load reduction and flushing. The following conclusions are drawn additional to the conclusions resulting from previous modeling of Lake Veluwe [25,26]. In the model a 68% phosphorus load reduction in 1979 resulted in a 32% reduction of the summer half year average return flux of phosphate from the sediment in Lake Veluwe four years after the load reduction. This shows that internal eutrophication can continue years after an external load reduction. Although the oligotrophication of Lake Veluwe resulted primarily from the load reduction, the removal of nutrients by infiltration and flushing contributed substantially. Moreover, the observed decrease of total phosphorus in the lake could only be reproduced by the model by including the precipitation of phosphate in an apatite-like mineral in the sediment, that binds phosphate more or less permanently. This suggests that authigenic formation of a stable apatite like mineral can contribute significantly to oligotrophication of a lake after a nutrient load reduction. The model also shows that the removal of nitrogen in the sediment by coupled nitrification-denitrification may be enhanced when the thickness of the oxidizing top sediment layer increases due to reduced detritus settling.