PCLake+: A process-based ecological model to assess the trophic state of stratified and non-stratified freshwater lakes worldwide

The lake ecosystem model PCLake is a process-based model that was developed to simulate water quality based on ecological interactions in shallow, non-stratifying lakes in the temperate climate zone. Here we present PCLake+, which extends the PCLake model to cover a wide range of freshwater lakes that differ in stratification regime and climate-related processes. To this end, the model was extended with a hypolimnion layer that can be invoked and configured by forcing functions or by simple built-in empirical relationships that impose stratifi- cation. Further adjustments to the original PCLake model have been made with respect to the calculation of 1) light irradiation in the water column, 2) evaporation processes and 3) phenology of macrophytes. The simulation output of PCLake+ for different types of lakes complies well with generally accepted limnological knowledge, thus holding promise for future contributions to ecological theory and application to lakes around the globe.


Introduction
Earth has over 117 million lakes with a surface larger than 2000 m 2 (Verpoorter et al., 2014). These lakes are of major importance for ecosystem services like drinking water, food and human health (Paerl and Otten, 2013). Environmental pressures such as eutrophication, global warming, land use change and overexploitation affect lake ecosystems and endanger the provision of these ecosystem services (Raworth, 2017;Rockström et al., 2009). As a consequence of these pressures, algal blooms are becoming increasingly frequent and severe (Paerl and Huisman, 2009).
Until now, we often fail to predict and prevent algal blooms (Spears et al., 2016(Spears et al., , 2017. Water quality assessment of lake ecosystems is critical in understanding the development of algal blooms and in ultimately preventing them. Global assessment studies are of great importance to achieve this goal, particularly in regions where fresh water is scarce. However, for many lakes there are no data available to assess their water quality in response to environmental change Messager et al., 2016). Moreover, data themselves are not predictive (Cuddington et al., 2013). Therefore, a modelling approach is needed for lake water quality assessment of lake ecosystems at a global scale .
A large diversity of lake ecosystem models to assess water quality exists, each with a level of complexity that matches its own purpose and niche . For instance, there are lake ecosystem models that contain only a few equations that have shaped ecological theory such as 'minimal models' on alternative stable states (Scheffer and Van Nes, 2007). On the other side of the complexity spectrum lay frameworks that allow for a high level of detail and serve the purpose of simulating spatial heterogeneity in lakes (e.g. FABM (Bruggeman and Bolding, 2014), Delft3D-WAQ/ ECO (Los, 2009)). In between, ecological models with an intermediate level of complexity combine relative short runtimes and limited data requirements of simple models with a level of complexity suitable to answer lake water quality related questions. Examples are one-dimensional lake ecosystem models such as AQUASIM (Omlin et al., 2001), MyLake (Saloranta and Andersen, 2007) and PCLake (Janse et al., 2010). These models are in an excellent position to bridge theory with application and are also suitable for global assessment studies by coupling to global assessment models (e.g. global NEWS or IMAGE-GNM, Beusen et al., 2016;Janse et al., 2019; https://doi.org/10.1016/j.ecolmodel.2019.01.006 Received 17 August 2018; Received in revised form 11 January 2019; Accepted 14 January 2019 Kroeze et al., 2016).
Having an intermediate level of complexity, PCLake is potentially suitable for global studies on lake ecosystems. Compared with other models with an intermediate level of complexity, PCLake includes key lake food web components . Ecological feedbacks between these food web components, such as the interaction between macrophytes and algae, are pivotal to critical nutrient loadings at which the ecosystem state shifts from a clear to turbid state or vice versa. As PCLake includes these key feedbacks within the food web, the model is well-suited to sudying impacts of eutrophication and oligorophication on lake ecology and water quality. For instance, the model has been used to assess critical nutrient loading of specific lakes (Hilt et al., 2018;Janse et al., 2008;Janssen et al., 2017;Kong et al., 2016;Kuiper et al., 2017), to simulate the effect of climate warming on these critical nutrient loadings (Mooij et al., 2007;Nielsen et al., 2014) and to assess the stability of ecosystems to disturbances around the critical nutrient loadings Kuiper et al., 2015). The power of PCLake lies in the coherence of the ecological process formulations, minimal overhead and run time and limited required input. Its strong focus on ecology provides an alternative approach to physically oriented lake models that have their basis in hydrodynamic processes and comprehensive energy balances in lakes.
PCLake was originally developed to simulate shallow, non-stratifying lake ecosystems in the temperate climate zone. While applications of PCLake to lakes in different climate zones have shown promising results (e.g. Janssen et al., 2017;Kong et al., 2016;Li et al., 2019), the climate-related processes need improvement. For instance, the phenology of macrophytes in PCLake is too rigid for global applications and the light module does not yet account for latitudinal differences in light intensity. Additionally, the assumption of a fully mixed water column does not hold for many lakes due to thermal stratification (Wetzel, 2001). Applying PCLake to stratifying lakes while assuming a mixed water column will lead to incorrect estimations of relevant lake characteristics such as chlorophyll-a and nutrient values (Janssen, 2017).
Different modes of thermal stratification are defined in limnological theory (Lewis, 1983;Ostë et al., 2010;Wetzel, 2001). We have included Fig. 1 for a graphical representation of these various modes. Shallow lakes are generally classified as polymictic lakes (from Greek polus 'much', and miktos 'mixed') because they are mixed throughout the year (continuous polymictic) or stratify multiple times for only short periods of time (discontinuous polymictic) (Andersen et al., 2017;Lewis, 1996). Deeper lakes are categorized as stratified lakes that mix either one period during winter (warm monomictic lakes), one period during summer (cold monomictic lakes), twice a year during spring and autumn (dimictic lakes) or not at all (oligomictic or amictic lakes) (Wetzel, 2001). During stratification periods the upper and lower layers hardly mix due to temperature dependent density properties of water. Besides a heterogeneous distribution of temperature throughout the water column, also nutrients, suspended solids, oxygen and various biological components distribute differently throughout the water column due to stratification (Jeppesen et al., 2009).
Here we extend PCLake to increase the model's application domain to a wide range of freshwater lakes that differ in thermal stratification modes and climate-related processes (Fig. 1), while staying as close as possible to the original model formulations and scientific niche. In line with the intended focus on ecological processes in PCLake, we deliberately kept the implementation of stratification and climate related processes in PCLake+ as simple as possible, using either forcing functions or empirical relationships.
First, we present the extensions made for thermal stratification and climate-related processes in the model and thereafter elaborate on how to initialize the model. Next, we tested whether the model reproduces expected patterns of oxygen, nutrients and chlorophyll-a based on basic limnological theory. The full set of state variables, parameters and equations can be found in the DATM file at GitHub: https://github. com/pcmodel. DATM files can be automatically translated to various programming languages including Grind for Matlab, R, ACSL, Deft3D, FORTRAN and C++ Van Gerven et al., 2015).

General setup of PCLake+
As in the original PCLake, PCLake+ has been developed to simulate average ecological dynamics of a lake. To this end, the model was extended with a hypolimnion layer that can be invoked and configured by forcing functions or by simple build-in empirical relationships that impose stratification (see for methods below). In case of continuous polymictic lakes, PCLake+ resembles the original PCLake model with a single fully mixed water column (Janse, 2005). In case of a stratifying lake the model was configured according to the schematization in Fig. 2, with the water column being divided into an epilimnion and a hypolimnion separated at mixing depth. Both water layers have their own state variables for nutrients, oxygen, three groups of algae, detritus, inorganic matter and zooplankton. Planktivorous fish, benthivorous fish, predatory fish and macrophytes are physically not restricted by stratification, thus were each captured in a single biomass state variable for the entire water column. River inflow and outflow of water were connected to the epilimnion whereas infiltration and seepage processes were connected to the hypolimnion. Similar to continuous polymictic lakes, stratifying lakes may have a marsh zone with emergent macrophytes, which was connected to the epilimnion. Below we give a thorough description on the implementation of stratification and climate related processes.

Implementation of stratification
The presence of stratification or mixing can be initialized by the user (calcMixDepth = 0) or estimated by the model based on the dimensions of the lake (calcMixDepth = 1). During stratification periods there is limited exchange between the epilimnion and hypolimnion layer. In absence of stratification advective and diffusive mixing processes take place between these layers.

Initialization of stratification by the user (calcMixDepth = 0)
Users may have information on the occurrence of stratification in their lake obtained from measurements or model simulations. Diverse statistical and mechanistic models exist to estimate the occurrence of stratification (Minns and Shuter, 2012). Examples are the statistical model of Ottosson and Abrahamsson (1998) and the mechanistic models Flake (Kirillin et al., 2011) and GLM (Hipsey et al., 2013).
Two important steps need to be made in order to use measured or modelled input data that impose stratification upon the model. First, the mixing depth needs to be set using parameter 'cDepthMix'. The mixing depth is the depth of abrupt change in temperature in the thermal profile of the lake. More complex temporal patterns in thermocline depth can also be imposed as a forcing function.
Second, the timing of stratification needs to be imposed in either of two ways. One option is to set the parameter 'ReadStrat' to 1 and use time series for 'mStrat', with 1 representing stratification periods and 0 representing mixing periods. The other option is to use temperature measurements over depth ('ReadStrat' = 0). The timing of stratification is then based on the time when the absolute difference between the surface and the bottom temperature of the lake exceeds a temperature threshold. This threshold has default value of 1°C according to the minimum temperature difference between the epilimnion and hypolimnion found at which stratification may occur (Woolway et al., 2014).
Similar to the original PCLake model, water temperature for both epilimnion and hypolimnion in PCLake+ can be set based on measured data or simulation output from other models using forcing functions.
Alternatively, they can be calculated according to a cosine function that was already applied in the original version of PCLake (Janse, 2005): where i represents the layer (epilimnion or hypolimnion), aTm i is the daily temperature in the specific layer (ºC), cTmAve i is the average temperature in the specific layer (ºC), cTmVar i is the maximum temperature variation in the specific layer (ºC), sTime is the day of the year, cTimeLag the parameter (days) that determines at which day the minimum temperature is reached and DaysInYear the number of days per year. In case the resulting daily temperature falls below zero, the water temperature will stagnate at zero degrees Celsius, as, realistically, ice formation will commence.

Initialization of stratification by the model (calcMixDepth = 1)
In the case where the mixing depth is unknown, first approximations are made by the model. The approximated mixing depth is based on an empirical relationship which has proven applicable to a great diversity of lakes varying in size, shape and geographical location (Hanna, 1990). This empirical relationship estimates the mixing depth (d mix , m) of a specific lake as function of the fetch (F in m) measured as the maximum length wind can blow over the lake surface in the predominant wind direction (maximum effective length): with a = 0.569 and b = 0.336 (Hanna, 1990). Lakes can stratify if their depth extends below the mixing depth. In all other cases, lakes are classified as continuous polymictic lakes (Fig. 3). Similar to the user initialized stratification, the method to approximate the timing of stratification is based on the difference in temperature between the epilimnion and the hypolimnion, set by default at 1°C (Woolway et al., 2014).

Implementation of climate-related processes
Climate-related processes of PCLake+ differ from the original PCLake in three ways: 1) methods used to calculate light irradiation in the water column have been improved, 2) the calculation of evaporation processes has been adjusted and 3) the representation of the phenology of macrophytes has been adjusted.
First, the method used to calculate light irradiation (W/m 2 ) was changed to broaden the model's applicability to different latitudes. Using the latitudinal position of the lake and the time of year, the daily incoming radiation was calculated according to the equations from Stull (2000). The daily incoming radiation was used to simulate the light profile in the epilimnion and hypolimnion. If the user has information on the incoming light irradiance, this information may be still imposed to the model instead of using these equations.
Second, a rough estimation of the evaporation can now be made using the Thornthwaite equation (Thornthwaite, 1948) using monthly irradiation and temperature. Again, the user can impose time series of evaporation in case these data are available from measurements or model simulations.
The last adjustment made to the model was the timing of energy storage of macrophytes. Like the original model, the state variable of macrophytes in PCLake+ accounts for three functional groups of macrophytes: floating, submerged and emergent macrophytes. In autumn macrophytes store their energy for surviving winter in propagules, seeds or roots. In the original PCLake model, the timing of energy storage was linked to a specific day in the year, yet macrophytes are Fig. 2. Schematization of a stratifying lake in PCLake+ with the water column divided in two layers: epilimnion and hypolimnion. For ease of comparison, this scheme of PCLake+ is designed similar to the scheme of the original PCLake model published in Janse (2005) and later updated by Van Gerven et al. (2015). An important addition to the original PCLake is the inclusion of a hypolimnion layer. Furthermore, all water state variables on the left side of the vertical grey dashed line were duplicated so that they are represented in both the epilimnion and hypolimnion; variables on the right side of the vertical grey dashed line were each captured in a single state variable for the entire water column. For example macrophytes may grow beyond the mixing depth. (For a colour version of this figure, the reader is referred to the web version of this article).
known to respond to physiological and environmental cues to determine timing of energy reallocation (Madsen, 1991;Van Dijk and Janse, 1993). In PCLake+, storage of energy by macrophytes to survive winter is now a function of light availability and the species specific minimum light needed for plant growth as an environmental cue (Hilt et al., 2018).

Analysis and initialization
Similarly to PCLake, PCLake+ can be used for transient and bifurcation analysis. Transient analysis is used to study temporal dynamics and impacts of the lake's history as captured by the initial conditions. In case of transient analysis, the user needs to carefully determine all initial state variables by using field data or estimates at time zero. A major challenge in transient analysis is that the initial conditions are often partly unknown.
Bifurcation analysis aims to study the equilibrium of the model in response to changes in parameters. This technique is particularly suited for examining the presence of hysteresis due to alternative stable states (e.g. Scheffer and Carpenter, 2003). The procedure of bifurcation analysis in PCLake+ was kept similar to the original version of PCLake. Compared to transient analysis, the exact determination of initial conditions is less important since bifurcation analysis depends on equilibria. Nevertheless, it is important to define the initial conditions carefully, because the initialization must be close enough to the equilibrium of a clear or turbid lake and far enough from the shifting point. Therefore, we advise to first take an initialization step by running a newly parameterized model twice to equilibrium, once with an extremely low (for clear conditions) and once with an extremely high nutrient load (for turbid conditions). Importantly, as the state variables in PCLake+ are dynamic during the year, equilibrium was defined here as recurring identical seasonal patterns for consecutive years. The equilibrium values for all state variables in both clear and turbid conditions by the end of December 31 st can then be used as initial values for the bifurcation analysis of the specific lake. Next, to perform a bifurcation analysis, a range of nutrient loadings needs to be selected. For each nutrient loading value within this range, the model needs to run twice to equilibrium, once initialized for clear and another for turbid conditions. The equilibrium output of interest, commonly chlorophylla, needs then to be averaged over the year or growing season. By plotting the averaged equilibrium chlorohyll-a as function of nutrient load creates the typical final product of the bifurcation analysis: a loadresponse curve. From this load-response curve the critical nutrient load can be deducted as an abrupt change in chlorophyll-a in case of a nonlinear curve, or as a point where a maximum chlorophyll-a level is exceeded in case the curve runs smoothly .
Initial state values for a standard deep lake are supplied in Appendix A. For comparability, the default settings were chosen as close as possible to the default settings as defined by Janse (2005). The default settings were based on a lake with an average depth of 15 m, a fetch of 1000 m, no marsh zone, and with an external nutrient load with an N/P-ratio of 10.

Unit of state variables
In the original PCLake, state variables in the water column were expressed as concentrations, in order to easily connect to empirical water quality measurements. To solve the mass balance in PCLake+, however, it proved necessary to express all state variables as masses per m 2 instead of concentrations. This is a direct result of the hydrodynamic components incorporated into the process descriptions of PCLake and PCLake+. The hydrodynamics are calculated in the same calculation step as biogeochemical dynamics. Therefore, the lake water volume required to calculate concentrations is only known after the changes in substances are calculated. Since the concentrations of the next time step depend on the lake volume of the next time step, determining exact concentrations for the next time step is in the effect not possible. To solve the issue, PCLake+ has been rewritten to a mass-based model in which all state variables are expressed in mass per meter squared. In the initialization step, initial values that are expressed in concentrations are converted into mass to accommodate the use of field measurements. As an intermediated step within each time step of the calculation, masses are converted back and forth to concentrations when needed in the equations. Moreover, the user will find the final output as auxiliary values in concentrations to accommodate an easy comparison with field measurements.

Verification
Because the output for non-stratifying lakes (continuous polymictic lakes) remained similar to the values reported by Janse (2005) we verified the model with a focus on stratifying lakes. We qualitatively verified model output in two steps, starting with simulations for deep lakes (15 m) and proceeding with simulations of shallow discontinuous polymictic lakes (1.2 m). We verified the output based on comparison with basic limnological knowledge regarding changes in oxygen, nutrient and chlorophyll-a concentrations corresponding with changes in stratification regimes. Below we describe the model settings used for simulations of both deep and shallow stratifying lakes and specify the expected limnological patterns.

Deep lakes
For deep lakes we focused on lakes in climate regions that could undergo two different stratification regimes depending on the local weather conditions. In particular we studied the model's response in key limnological variables around the following transitions in stratification regime: 1) amictic lake → monomictic cold lake, 2) monomictic cold lake → dimictic lake, 3) dimictic lake → monomictic warm lake and 4) monomictic warm lake → oligomictic lake. These transitions grossly occur -in this order -to decreasing latitude on the globe or along an altitudinal gradient. To this end, we simulated two hypothetical years in four lakes that differ slightly in water temperature such that the lake experiences a transition in stratification regime. We have adjusted the variation in the epilimnetic water temperature (cTmVarEpi) to simulate the four transitions. Since the focus is on the effects of the mixing regime rather than on the effect of temperature, we  (Hanna, 1990). According to this relationship, lakes with a depth in the grey zone lack stratification, whereas lakes in the white zone will stratify. MEL = Maximum Effective Length.
varied the temperature minimally around the transitions (see Table 1 for the temperature settings used). The mixing depths have been determined using Eq. (2).
Then, we confronted the model results with basic limnological knowledge on the change in oxygen, nutrient and chlorophyll-a concentrations when stratification regimes changes within a specific lake (e.g. Lewis, 1983;Wetzel, 2001). According to general theory (e.g. Wetzel, 2001), the transition from an amictic lake to a monomictic cold lake (Lake 1) results in a single mixing event in summer. This mixing event brings anoxic and phosphate-rich hypolimnion water to the surface causing an intense summer phytoplankton bloom (Wetzel, 2001). The transition from a monomictic cold lake to a dimictic lake (Lake 2) allows the lake to mix twice rather than once, with stratification during winter and a short period in summer. The double mixing event could boost both a spring and late summer phytoplankton bloom due to the phosphate-rich hypolimnion water. The transition from a dimictic lake to a monomictic warm lake (Lake 3) allows the lake to mix once instead of twice, with stratification during summer and mixing in winter. The mixing in winter boosts mainly the spring phytoplankton bloom. Finally, the transition from a monomictic warm lake to an oligomictic lake (Lake 4) results in a year-round stratified lake. Since the phosphate-rich hypolimnion water of an oligomictic lake hardly mixes with epilimnion water, phytoplankton concentrations will be lower in the epilimnion compared to when the lake would experience mixing events.

Shallow lakes
Compared to deep lakes, the period of thermal stratification in shallow discontinuous polymictic lakes is generally frequent but short with a duration of only a few hours to a couple of days (Padisák and Reynolds, 2003;Vilas et al., 2018). Thermal stratification in shallow lakes occurs at sunny and windless days (Padisák and Reynolds, 2003). Colonization of shallow lakes by macrophytes can enhance the period of stratification by dampening turbulence and absorbing heat (Andersen et al., 2017;Vilas et al., 2018).
For the hypothetical discontinuous polymictic shallow lake we simulated an oligotrophic lake colonized by macrophytes. The mixing depth of this lake was manually set to 0.6 m. The epilimnion temperature followed a cosine function (Eq. (1)) with an average of 12°C, a variation of 10°C and small random variations to simulate weather changes. We assumed that during summer months the resulting small temperature changes due to local weather were damped by macrophytes (Andersen et al., 2017;Vilas et al., 2018). Consequently, the hypolimnion temperature was based on the cosine functions without these random temperature variations. When temperature differences between epilimnion and hypolimnion were larger than 1°C stratification occured.
Following limnological theory, stratification in shallow lakes reduces the hypolimnetic oxygen concentrations (Andersen et al., 2017;Vilas et al., 2018). As a result, phosphate is released from sediments and hypolimnetic algal concentrations are expected to become higher than the epilimnetic algal concentrations because of sufficient light and nutrient availability (Andersen et al., 2017;Vilas et al., 2018).

Results
We developed PCLake+ by expanding the properties of the existing ecological lake model PCLake (Table 2). These changing properties have expanded the application domain of the model by including additional processes and a hypolimnion layer to the model. We simulated four deep lakes (15 m) and one shallow discontinuous polymictic lake (1.2 m) to verify the model output of PCLake+ against expected general limnological patterns.

Deep lakes
The comparison of the simulation output for the four deep lakes with different stratification regimes is shown in Fig. 4. A lake in a climate region which can make the lake either amictic or monomictic cold (Fig. 4, Lake 1) differs in vertical patterns in oxygen, phosphate and chlorophyll-a depending on the regime it experiences. Due to the Table 1 Temperature settings of 4 hypothetical lakes that change in stratification regime due to a small change in the yearly epilimnetic temperature variation (cTmVarEpi). Note that water is warmer in the hypolimnion than in the epilimnion for the entire year in cold amictic lakes and for part of the year in monomictic cold lakes. This is because water with a temperature of 4°C is denser than water with a temperature of 1°C and therefor sinks to deeper layers.

Table 2
Comparison of properties for PCLake and PCLake+ with respect to the application domain, main processes included and technical settings.

Processes
Trophic Interactions X X Suspension and resuspension X X Burial X X Marsh zone related processes X X Flexible stoichiometry X X Grazing X X Lake measures (mowing, flushing, fishing, etc) X X Phenology of macrophytes related to available light -X Dependence of light irradiation on latitude -X Climate-dependent evaporation -X Thermal stratification ( mixing event in summer when the lake was monomictic, hypolimnetic phosphate rich water was mixed with the epilimnion. This boosted a late summer chlorophyll-a peak in autumn that was absent when the lake was amictic. Additionally, the spring peak of algae was also higher when the lake is monomictic. When lakes reside in a climate region that can make the lake either monomictic cold or dimictic (Fig. 4, Lake 2) the difference between these regimes was made by the presence of a short stratification event during summer. In the dimictic situation, the additional stratification period during summer allowed for the hypolimnion to become anoxic for a short period, leading to higher hypolimnetic phosphate concentrations. These increased concentrations led to a higher chlorophylla peak during summer compared to the monomictic situation.
A lake located in a climate region which can make the lake either dimictic or monomictic warm (Fig. 4, Lake 3) results in respectively the presence or absence of a short stratification event during winter. The difference in phosphate and chlorophyll-a concentrations between both stratification regimes in such a lake was small as stratification occurs in winter when the biological activity is reduced. Chlorophyll-a was slightly higher when the lake is monomictic warm, but as this was observed throughout the year rather than a response directly after a mixing event, this is likely a result of the small increase in annual temperature rather than stratification dynamics.
Lastly, we consider a lake in a climate region which can undergo either a short period of mixing (monomictic warm) or no mixing event at all (oligomictic) (Fig. 4, Lake 4). As the anoxic and phosphate rich hypolimnion water does not mix with the epilimnion, the chlorophyll-a concentration was lower in the oligomictic lake compared to the monomictic warm lake.

Shallow lakes
The discontinuous polymictic shallow lake is stratified between the end of March when macrophytes start to grow until September (Fig. 5). Especially during spring, the stratification periods were characterized by a slightly lower oxygen concentration in the hypolimnion. The drop in oxygen levels was not very obvious since most stratification periods last for a few days maximum, after which the hypolimnion was mixed again with the oxygen-richer epilimnion. The hypolimnetic chlorophylla levels were higher than in the epilimnion, even though the increase in phosphate levels was small. High chlorophyll-a values were possible due to the higher nutrient release from the hypolimnion while light availability was high. Later in the year, macrophytes were so abundant that the light availability was not sufficient to sustain high hypolimetic algae concentrations.

Discussion
With the further development of PCLake into PCLake+ we aim to increase the application domain of the model by including additional processes and a hypolimnion layer to the model. Our qualitative verification showed that PCLake+ is able to simulate dynamics of stratifying lakes known from limnological theory. In the original PCLake model, these dynamics are missing, leading to unrealistic patterns in oxygen, nutrients and chlorophyll-a concentrations (Fig. 4 compared to Appendix B) when applying PCLake to stratifying lakes. Our qualitative verification against general limnological patterns is a first confirmation of the validity of the model; a quantitative validation using a dataset of multiple lakes around the globe will be a next step in confirming the applicability of PCLake+ to thermally stratified lake systems.
PCLake+ offers a number of advantages over the original PCLake model. First, anoxia caused by stratification and the associated increase of epilimnion phosphorus after mixing events were absent in the original PCLake (Appendix B, Figure B.1) but is now an emergent property of PCLake+. This facilitates a better simulation of the subsequent algae blooms after mixing (Appendix B, Figure B.2). Additionally, the incoming nutrients via rivers in PCLake+ enter the system in the epilimnion layer instead of instantaneously mixing throughout the entire water column, leading to more realistic estimates of nutrient availability in the epilimnion. Lastly, in the original model, all state variables, including chlorophyll-a were averaged over the entire water column whereas measurements are typically taken in the epilimnion. Especially in deep lakes, large part of the water column is photosynthetically unproductive because of light limitation. As a result, taking an average of the entire column may lead to unrealistically low values of chlorophyll-a compared to measurements if the lake is deep (Appendix B, Figure B.3-B.5). In PCLake+, this is solved by splitting the water layer into two distinct compartments, the epilimnion and the hypolimnion in order to compare emperical measurements from the epilimnion with epilimnion model output.
While the development of PCLake+ has led to an increased application domain, there are still open challenges. First, PCLake+ accounts Results for PCLake+: simulations for a shallow discontinuous polymictic lake with an average depth of 1.2 m and a mixing depth of 0.6 m. The fetch of the hypothetical lake was set to 1000 m, there was no marsh zone and we used an external nutrient load with a N/P-ratio of 10 (P-load was 0.002 g P.m −2 .d -1 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
for stratification periods that last for at least one full day, prescribed by the model's temporal resolution of one day. Simulation of stratification periods of a few hours that may occur in small lakes are thus beyond the scope of the model. Additionally, we used empirical relationships to determine stratification depth and timing in case the user holds no data on stratification. For global assessment studies, this simplification is in most cases acceptable and even desirable.
We can envision future developments of the model to include brackish and salt lakes, for which the model is currently unsuitable. By adding chemical stratification and a good process description of the biochemical processes that strongly depend on sulphur reduction, PCLake+ would be suitable to also model brackish and salt lakes. For future applications of PCLake+ in global carbon studies, the carbon cycles have to be extended, including the addition of pH as state variable.
Another obvious limitation of PCLake+ is that feedbacks from ecological processes to hydrological and climate-related processes are missing. For instance, hydrodynamic resistance of water flow by macrophytes or the influence of macrophytes on stratification is not covered by the model (Andersen et al., 2017;Van Nes et al., 2002;Vilas et al., 2018). In applications where such feedbacks are deemed to be important, the biogeochemical components of PCLake+ should be coupled with more sophisticated hydrodynamical frameworks such as FABM and DELWAQ (Hu et al., 2016;Janssen et al., 2017). Coupling of PCLake+ to other models has additional advantage as simulations can then be run in various spatial settings. To facilitate such coupling, PCLake+ is available in DATM (Database Approach To Modelling) Van Gerven et al., 2015). With DATM, the model code can be translated to various programming languages including Grind for Matlab, R, ACSL, Delft3D, FORTRAN and C++.
Summarizing, in the development of PCLake+ we embraced the philosophy of the PCModel family that ecology should be the core of the model while keeping the formulation of e.g. hydrodynamics and energy balances as simple as possible. We foresee three application domains for PCLake+: theoretical studies, global assessment studies and case studies for specific lakes. Experiences with PCLake have shown that the model is well suited for studying fundamental ecological mechanisms Kong et al., 2016;Kuiper et al., 2017). Within the application domain of global assessment studies, where often large data gaps exist, the generic and simplified assumptions in the default model settings allow broad estimates of a wide range of lakes. In its application to specific lakes, the model provides much flexibility to the user to impose temporal input data on e.g. temperature, evaporation and stratification process originating from measurements or other models (cf. Janssen et al., 2017;Kong et al., 2016). As with any model, the proof of the pudding is in the eating. In this paper, we assessed the model by 1) comparing the output with expert knowledge and 2) examining the model's sensitivity to different temperatures and lake depths. But only in its application will PCLake+ reveal its real strengths and weaknesses. Experiences with PCLake showed that the domain of application of the model has changed over time. The intended application domain of the original PCLake model was eutrophication in temperate shallow lakes (Janse et al., 2008(Janse et al., , 2010. Thereafter, however, the model was successfully appliedamong others -in climate studies (Mooij et al., 2007;Nielsen et al., 2014), food web studies , coupled with a spatially heterogeneous model (Hu et al., 2016;Janssen et al., 2017) and showed its value when applied for temporal dependent input in the subtropical lakes Kong et al., 2016;Li et al., 2019). Our next step with PCLake+ is a calibration and quantitative validation of the model's output. In order to test for the generality of the model and to prevent overfitting, we advocate a similar method as used to verify PCLake; by testing with a dataset containing multiple lakes around the globe rather than to verify for a single lake. A subsequent step would be to apply PCLake+ within an ensemble study in which multiple models are run based on a unified dataset. Analyzing the ensemble output would not only put PCLake+ into perspective with other models but will also be useful for global assessments of lake ecosystems and related ecosystem services

Conclusion
With the development of PCLake+ we have a model that is applicable to a wide range of both deep and shallow lakes around the world that exhibit various stratification regimes and experience different climatic conditions. PCLake+ maintains the philosophy of the PCmodel family with its focus on ecological dynamics. The extended model produces the same output as the original PCLake in case of continuous polymictic lakes which makes the model backward compatible. The first output of PCLake+ complies well with basic limnological knowledge for deep and discontinuous polymictic shallow lakes. These results holds a promise for future model development and application to lakes around the globe.

Conflicts of interest
There is no competing interest to declare.