Water Ecosystems Tool (WET) 0.1.0 – a new generation of flexible aquatic ecosystem model

. We present the Water Ecosystems Tool (WET) - a new generation of an open source, highly customizable aquatic 10 ecosystem model. WET is a completely modularized aquatic ecosystem model, developed in the syntax of the Framework for Aquatic Biogeochemical Models (FABM), which enables coupling to multiple physical models ranging from zero to three dimensions, and is based on the FABM-PCLake model. The WET model has been extensively modularized, empowering users with flexibility of food web configurations, and also incorporates model features from other state-of-the-art models, with new options for nitrogen fixation and vertical migration. With the new structure, features and flexible customization options, WET 15 is suitable in a wide range of aquatic ecosystem applications. We demonstrate these new features and their impacts on model behavior for a temperate lake for which a model calibration of the FABM-PCLake model was previously published, and discuss the benefits of the new model. of sediment, a fully size-based fish module, and extensive testing and improvements of the model in a 3D application, with the expressed aim of the authors (in their roles as the current main developers of WET), that the model be even more tailorable to a wide array of 475 ecosystems types, across latitudinal, spatial and productivity gradients, simply by turning features on or off, and combining different modules, all configurable at runtime. With this model, which is open source and freely available, we hope to facilitate the consolidation of successful features of many models together in one, with the goal of preventing ‘re-inventions of the wheel’ in the future, and making aquatic ecosystem modelling easier, more flexible and, ultimately, better.


Introduction
The study and management of aquatic ecosystems have benefitted widely from the ongoing development of various numerical 20 model approaches to a host of ecological questions (Soares and Calijuri, 2021). As the field matures, new and superior approaches and descriptions of individual ecological processes are formulated and improved upon, and management tools must continuously be updated to reflect the current state-of-the-art. However, rather than building new models from scratch, and thus 're-inventing the wheel' over and over (Trolle et al., 2012), another way forward is to consolidate new descriptions of ecological processes into a few proven and well-established biogeochemical models, thereby improving their applicability for 25 management and for the study of ecosystem-wide responses to environmental stressors. There is therefore a call for flexible and configurable models that contain various optional features, allowing them to be tailored to specific uses, without changing the model code or making a new version.
Among the most widely used lake ecosystem models in the world (Mooij et al., 2010;Trolle et al., 2012), the PCLake model was originally developed for the shallow Dutch lake Loosdrecht in the early 90s, under the name PCLoos (Janse and Aldenberg,30 that movement of plankton elements is limited to passive advection and a constant sinking or flotation velocity. These limitations might be acceptable in shallow environments or 0D applications, but quickly become untenable in deeper systems. 65 Here, we present a complete restructuring and -coding of the FABM-PCLake model that adds both flexibility as well as new features to the model. To avoid conflation with the increasing number of PCLake versions available, we have decided to present this new model under its own distinct name, the Water Ecosystem Tool (WET). So far, the name Water Ecosystem Tool (WET) has been associated to the QGIS plugin developed by Nielsen et al. (2017) to setup, configure and run the coupled 70 GOTM-FABM-PCLake model complex. With this paper, we redefine what WET is, and present it as a new generation of an aquatic ecosystem model, originating from the PCLake model, specifically the version by Hu et al. (2016).
In the following sections, we present the suite of new features which have been added to the PCLake framework, together with example dynamics from its first application to a lake ecosystema temperate Danish lake. The new features constitute a complete modularization of the model code, the inclusion of vertical migration algorithms and the addition of a nitrogen 75 fixation option to the phytoplankton module. A plugin (now called QWET) for the GIS software QGIS has also been developed (Nielsen et al., 2021), which provides a graphical user interface to configure and runallows WET in a user-friendly workflowto be configured and run in conjunction with the 1D hydrodynamic model GOTM, and allows (but does not require) linking GOTM-WET to the SWAT and SWAT+ (Soil & Water Assessment Tool) watershed model (Arnold et al., 1998). As for its predecessor, the plugin provides a Graphical User Interface (GUI), allowing the model to be set up and run in a user-friendly 80 workflow. The plugin can be found on the WET website (http://wet.au.dk), along with user guides instructing the user on how to download and set up both WET and QWET.

Model description
Like its predecessor FABM-PCLake, WET can describe interactions between multiple trophic levels and abiotic nutrient 85 dynamics in both the water column and the sediment. The model accounts for the dynamics of dry weight, nitrogen, phosphorous, silica and oxygen, and features bottom-shear-dependent resuspension, as well as two different light-limitation functions for phytoplankton. WET is also implemented within the FABM framework, allowing the model to be coupled to various physical driver models, e.g. GOTM (1D, Burchard et al., 1999) or GETM (3D, e.g. Stips et al., 2004), without changing any of the model code. Within the FABM framework, the physical model takes care of updating and iterating the model state 90 variables forward in time, and the sole responsibility of the WET code is to provide local source and sink terms for its state variables as well as feedback to physical variables such as light or bottom shear stress (Bruggeman and Bolding, 2014).
Here we concentrate on descriptions of the new features and all changes that separate WET from its parent model. We refer to Hu et al. (2016) and Janse (2005;1995) for a detailed description of the basic equations governing the biogeochemical processes and food web dynamics, since these are unchanged, even though the model code has been rewritten and reorganized. 95 A complete list of parameters related to the new features with options and default values can be seen in Table 1.

Modularization of the food web
A major drawback of the FABM-PCLake model is that it retains the rigid food web structure inherited from the original 100 PCLake model and can only run in a fixed food web configuration, with fixed, preordained interactions between food web components. In contrast, WET has been designed to take full advantage of FABM, by being fully modularized. This modularization enables the user to set up an arbitrary number of types of any food web element (e.g. multiple phytoplankton types) within a simulation, or to remove it altogether (e.g. no fish). Thus, it is possible to customize the model to a desired level of complexity with the aim of addressing a specific study system or research question, and without changing any of the 105 model code. Increasing or decreasing simulated food web complexity as the situation requires or data allows is done by simply adding or removing a food web module instances from a single configuration file (fabm.yaml, see Fig. 3). As an example, to add a new zooplankton species to an already existing model, one could first copy the lines for an existing zooplankton species in the configuration file, and then change the name of the copy instance. Secondly one would go through the couplings and parameters sections in the new zooplankton instance, modifying these to fit the desired organism. Finally, one would modify 110 the instances of any predators to include the new zooplankton instance in their diets. Thus, adding or subtrackting instances to a model setup is relatively easy, and testing for the optimal food web configuration in a specific case is possible, if not usually very feasible, by calibrating several different module setups and comparing their performance.
The code base of the WET model consists of eleven FORTRAN files. Six of these are required files, of which two handle model initialization and shared functions, and four constitute a basic chassis of required modules, handling microbial, chemical 115 and physical processes in the water column and upper sediment (see Fig. 2, 'fixed modules'). Besides these core parts, the model consists of five optional modules representing different food web component types (phytoplankton, rooted macrophytes, zooplankton, zoobenthos, and fish, see Fig. 2, 'food web modules'). For all WET modules, fabm.yaml testcase setup files with default parameters can be found in the testcase folder of the source file repository. 120

Primary producers
The WET phytoplankton module is developed with a high degree of flexibility in mind, and contains options to allow it to represent all main phytoplankton groups. These constitute optional dependence on silica for growth (e.g. for modelling diatoms), the option to allow phytoplankton to fix atmospheric nitrogen (e.g. cyanobacteria), and the optional ability to migrate vertically in response to ambient light and nutrient availability (various phytoplankton groups). Both nitrogen fixation and 125 vertical movement algorithms are new features of WET, described in the following sections (Sect. 2.2 and 2.3, respectively).
The option of silica dependence is turned on by setting the a lSi parameter to true in the configuration file (fabm.yaml, see Fig.   3 and Table 1), and its formulation is otherwise identical to the one used in the PCLake model family (Hu et al., 2016;Janse, 2005;Janse and van Liere, 1995). In WET, multiple types of rooted macrophytes can be included. Depending on the requirements of the modeler, all or some of 130 the macrophyte instances can be set up to share a common carrying capacity, by setting changing the nComptsa parameter to >0 forin all competing macrophyte instances, and pointing to the relevant macrophyte instances in the configuration file. Aside from this option and the general modularization, the formulation of the macrophyte module is identical to FABM-PCLake (Hu et al., 2016;Janse, 2005;Janse and van Liere, 1995). 135

Heterotrophic modules
In accordance with the modularization of WET, the original feeding formulations (see Janse, 2005) of heterotrophic modules zooplankton, zoobenthos, and fishhave been adapted to be flexible, allowing predators to feed on multiple prey (including other instances of their own base module). Mixed diets are set up in the configuration file, where each prey is pointed to, and a preference factor (typically between 0 and 1) is specified, following e.g. Fasham et al. (1990). In addition, consumption of 140 particulate organic matter (POM) is likewise an option for the invertebrate modules, i.e. zooplankton and zoobenthos.
For fish, foraging is separated into three foraging modes, planktivory, benthivory, and piscivory. These foraging modes are, assumed separate in time and/or space, such that each take up a fraction (fFishZoo, fFishBen and fFishPisc, which must sum up to 1) of the total foraging effort of the fish population, the foraging mode fraction (fFishZoo, fFishBen and fFishPisc, which must sum up to 1). For each mode, several prey types can be present, each with their own preference factor as for zooplankton 145 and -benthos. Saturation functions are calculated for each foraging mode separately, using a Monod-type formulation, e.g. for piscivory:

150
where aDSatFiPisc is the current saturation level for piscivorous feeding, nPISC is the number of fish prey types, SDPisci is the biomass of fish prey i, PISCprefi is the preference factor for fish prey i, and hDPiscFi is the half-saturation constant for piscivorous feeding. The amount of assimilated biomass from each foraging mode (here tDAssFiPisc, again using piscivory as an example) is then calculated as: where sDFi is the biomass of the predator, kDAssFi is the maximum assimilation rate of fish at 20°C, aFunVegAss is the (optional) dependency of fish predationpiscivory on macrophyte biomass (Janse, 2005;Janse and van Liere, 1995), and uFunTmFish is the temperature correction on fish vital rates (calculated internally in WET). The contributions from the three 160 foraging modes are then summed for the purpose of calculating total assimilation.

Linking fish instances into a pseudo stage structure
WET describes all populations in terms of biomasses, and does not explicitly consider population or age structure of any organism. For fish however, the link between juvenile (zooplanktivorous) and adult (benthivore) fish present in the PCLake 165 model family (Hu et al., 2016;Janse, 2005;Janse and van Liere, 1995), has been generalized to the fish module. Thus, in WET, instances of the fish module can be linked through 'aging' or 'reproduction', where a fixed proportion of biomass is transferred from one instance to another on a fixed date. For both aging and reproduction, this is set up in the configuration file, by setting the qStageOpt parameter, pointing to the recipient fish instance(s), and providing parameters for the date of transfer and the proportion of biomass transferred. In this way, a population of fish can be separated into a stage structure, 170 containing two, three or more stages, with individual parameterizations, diets and predators. Note, however, that this implementation is not truly a stage-structured model, as each instance in the structure can in principle persist indefinitely, regardless of the state of the others.

Phytoplankton nitrogen fixation 175
Depending on the external nutrient inputs, nitrogen fixation by cyanobacteria can be an influential process in freshwater ecosystems (e.g. Paerl et al., 2016). Advancing from the PCLake model family (Hu et al., 2016;Janse, 2005;Janse and van Liere, 1995), WET features the possibility to simulate nitrogen fixation, which can be turned on or off by the user for each individual phytoplankton instance. This is done by setting the lNfix parameter to .true. in the configuration file (see Table 1), and supplying two additional parameters; the maximum fixation rate (cNFixMax, mgN mgDW -1 d -1 ), and the maximum 180 realized fraction of the growth rate at maximum nitrogen fixation rate (fMuNfix, dimensionless). By default, total nutrient limitation for phytoplankton growth is governed by Liebig's Law of the Minimum, and is by default calculated as: where aNutLim is the overall nutrient limitation, and aPLim, aNLim and aSiLim are the Droop functions for phosphorous, nitrogen and silica growth limitation, respectively. In the case of nitrogen fixation being turned on, nitrogen uptake rate is assumed to never be limiting for phytoplankton growth, and consequently, total nutrient limitation is: or simply aPLim, in the absence of silica uptake. This independence of internal nutrient concentration for growth dynamics is balanced by a growth rate reduction due to allocation of energy to nitrogen fixation: Note that this formulation assumes that nitrogen fixation takes place whenever the phytoplankton is less than nitrogen replete (and in spite of other possible limiting nutrients) and that phytoplankton only fixes what nitrogen it cannot uptake through mineral absorption. The final nitrogen fixation rate is calculated as: and has units of mg N mg -1 DW d -1 . This formulation of the nitrogen fixation dynamics in WET is adapted from the CAEDYM model (Hamilton and Schladow, 1997;Hipsey et al., 2005), while this general formulation of nitrogen fixation is common in phytoplankton models (see e.g. Inomura et al., 2020 and references therein). 205

Vertical movement algorithms
Many organisms employ vertical movement (VM) as a means to exploit vertical gradients in e.g. nutrient, oxygen or light availability (e.g. Dini and Carpenter, 1992;Mehner, 2012;Olli, 1999). Existing model variants in the PCLake family were primarily developed for shallow lake applications, and do therefore not consider the ecological ramifications of lake 210 stratification or the movement of organisms in the vertical. While the FABM-PCLake model have been applied to deeper, stratified lake systems in a 1D setup (e.g. Allan, 2018;Chen et al., 2020), its biological descriptions and structure has inherited some limitations in how it deals with spatial heterogeneity from its 0D predecessors. Examples of this includes the fact that FABM-PCLake has no exchange of fish between model depth layers, even when depth resolution is fine (e.g. layers being only a few centimeters thick), and that movement of plankton elements is limited to passive advection and a constant sinking 215 or flotation velocity. These limitations might be acceptable in shallow environments or 0D applications, but quickly become unacceptable in deeper systems.
In WET, all pelagic food web modules (i.e. phytoplankton, zooplankton and fish) now have several options for vertical migration. The purpose of these options is to 1) aAdd flexibility to the types of environment that can be modelled with WET, and 2) as well as tTo increase model accuracy and applicability, by providing more realistic dynamics of all food web elements, 220 in all types of aquatic systems. The various options for vertical migration are presented for each model element below. These options can be configured individually for each instance at runtime by setting the qTrans parameter, and any necessary additional option-specific parameters.

Phytoplankton module 225
The WET phytoplankton module contains four modes of vertical movement behavior: Passive transport (no active movement), passive transport plus active chemotaxis (for nutrients), passive transport plus active phototaxis, and passive transport plus combined photo-and chemotaxis (see Table 1 for details on options and parameterization). The VM functions of phytoplankton described below are all based on Ross and Sharples (2007).

230
The first vertical movement option is passive transport (qTrans = 1), and is, identical to what is currently the only available mode in FABM-PClake, where phytoplankton only move vertically as a result of passive advection by the host physical model as well as through a fixed sinking rate (positive, negative or neutral buoyancy). The fixed sinking rate is specified through the cVSet parameter, and is negative in case of sinking.
In addition to passive advection, when chemotaxis is turned on (qTrans = 2 or 4), phytoplankton will swim downwards with 235 constant swimming speed cVSwim (m d -1 ), whenever nutrient limitation growth coefficient aNutLim decreases below a threshold value (fNutLimVMDown, dimensionless). Thus, phytoplankton in WET operate under the assumption that nutrient concentrations are always higher at greater depth.
When phototaxis is turned on (qTrans = 3 or 4), phytoplankton will swim upwards with fixed speed cVSwim, whenever ambient light levels surpass a light threshold value fLVMmin (W m -2 ). 240 The fourth phytoplankton vertical movement option (qTrans = 4) combines chemotaxis and phototaxis, such that chemotaxis takes precedence over phototaxis. Thus, phytoplankton will swim down in nutrient deplete situations (aNutLim < fNutLimVMDown), and up when the cell is nutrient replete (aNutLim > fNutLimVMUp), provided that ambient light levels surpass the threshold value fLVMmin.

Zooplankton & Fish modules
Regular vertical movements between different depths are a common behavior in both fish and zooplankton populations, especially in deeper systems. Such movement behavior is expressed for a variety of reasons, including avoidance of hypoxic regions, predator avoidance, and bioenergetic exploitation of gradients in temperature or food availability (Dini and Carpenter, 1992;Dodson, 1990;Lambert, 1989;Mehner, 2012). 250 An inherent limitation to the FABM framework is that modules are limited in the amount of information they receive about conditions outside the current model cell, e.g. food availability at other depths. Thus, modelled motile organisms are limited to making 'decisions' about movement based on either local conditions or predictable environmental gradients. Due to this limitation, directional vertical movement of zooplankton and fish in WET are restricted to being in response to hypoxia, ambient light levels, or both. The WET zooplankton and fish modules contain four modes of vertical migration behavior: No 255 transport (no transport between layers, even advection), passive transport (no active transport), hypoxia avoidance, and light-Formatted: Normal based diel vertical migration combined with hypoxia avoidance. Of these, the no transport option (qTrans = 0) turns off all exchange between model layers (only relevant in >0D applications), and is mainly a debugging or analytical tool.
As for phytoplankton, when vertical movementqTrans is set to 1 (Ppassive transport), fish and zooplankton are passively advected by the physical model. 260 Hypoxia avoidance is turned on by setting qTrans = 2 (or 3), and restricts the habitat domain to exclude anoxic parts of the water column, a ubiquitous response to hypoxia among zooplankton and fish (Ekau et al., 2010;Vanderploeg et al., 2009).
When hypoxia avoidance is turned on for a WET module, fish or zooplankton swim upwards with swimming speed cVSwim (m d -1 ) whenever the ambient oxygen concentration falls below the critical threshold (cMinO2, g O2 m -3 ) in the current model cell. 265 Zooplankton and fish may employ diel vertical migration for a number of reasons (Dini and Carpenter, 1992;Dodson, 1990;Lambert, 1989;Mehner, 2012;Sainmont et al., 2013), however ambient light levels is most often the proximate trigger for this behavior. Under this setting, iIn addition to passive advection and hypoxia avoidance, when qTrans is set to 3, fish or zooplankton will swim downwards with speed cVSwim, whenever light exceeds a maximum level cMaxLight (W m -2 ), and upwards whenever light decreases below the lower light threshold cMinLight (W m -2 ). Either the downwards or upwards 270 portions of this movement can be turned off, by setting the maximum (minimum) cMaxLight (cMinLight)threshold to a very high (negative) value.

WET testcase -Lake Bryrup 275
To illustrate some of the new features in the WET model, we applied the GOTM-WET model for Lake Bryrup, for which a calibration of the FABM-PCLake model (coupled to the lake version of the 1D hydrodynamic model GOTM, Burchard et al., 1999) was previously published by Chen et al. (2020). Here we have adapted, recalibrated and validated this setup using WET, following the methodology and approaches of Chen et al. (2020), and use the results to illustrate some of the new features of WET, and how these can impact the model behavior. 280

Study site and model configuration
The shallow Lake Bryrup is located in the Central Region of Denmark (56.02• N, 9.53• E), and has a surface area of 37 hectares, a mean depth of 4.6 meters, and a maximum depth of 9 meters. The lake stratifies temporarily during summer, and has a water retention time of 2-3 months. The catchment area is 49.9 square kilometers and heavily farmed, and the lake 285 receives large amounts of nutrients from agricultural and urban drainage (Johansson et al., 2019). Consequently, the lake is eutrophic although management measures have been effective in reducing average total nitrogen and phosphorous concentrations throughout the last decades (see Chen et al., 2020, for a more thorough description of Lake Bryrup).

Formatted: Font: Not Italic
As the chosen physical driver model for the Lake Bryrup test case, we used the Aarhus University fork of lake branch GOTM (version 5.2.2-au, available at https://gitlab.com/WET/gotm), configured with a maximum depth of 9.0 and 18 vertical layers 290 (i.e. vertical grid size of 0.5 m). In order produce high-resolution output for the present figures, the model was also run with a 200 layer resolution. The setup used a lake-specific hypsograph (i.e. the relation between depth and horizontal area) a to capture lake sediment-water column interactions at all depths, by effectively splitting the bottom between model layers, such that each model layer in the 1D setup has an attached bottom layer. Interactions between the water column of a layer, its attached bottom, and the water column layer below is governed by the hypsograph, which specifies the fraction of the bottom area to total layer 295 area. The setup therefore in some aspects functions as a pseudo-2D setup. Each water model layer thus included a sediment layer of 10 cm, similar to the sediment compartment of the PCLake model. See Andersen et al. (2020) and Hu et al. (2016) for detailed descriptions of these aspects of the model, and layer volumes and areas between layers and at the sediment-water interface were derived from a lake-specific hypsograph (i.e. the relation between depth and horizontal area). To simulate lake ice thickness and cover, implementation in GOTM of the ice module from MyLake (Saloranta and Andersen, 2007) was 300 enabled. In WET, the food web comprised three phytoplankton groups (diatoms, cyanobacteria and other algae), rooted, submerged macrophytes, two zooplankton groups (Daphnia and other zooplankton), detritivorous macrozoobenthos, juvenile zooplanktivorous and adult zoobenthivorous fish and piscivorous fish (Fig.ure 1). We applied the European ECMWF ERA5 dataset (Hersbach et al., 2018) at an hourly resolution on air temperature ( o C), air pressure (hPa), dew-point temperature ( o C), cloud cover (%) and wind speed components (m/s) in the north-south and west-east direction as meterological forcing for 305 GOTM. Monthly averages of water inflow (m 3 /s) and nutrient concentrations (NO3, NH4, PO4 and particulate organic matter of nitrogen and phosphorous, mg/L) from the National Monitoring and Assessment Program for the Aquatic and Terrestrial Environment in Denmark (NOVANA) (Kronvang et al., 1993;Lauridsen et al., 2007) were used as boundary conditions (see Chen et al., 2020, for details), applied in the topmost model layer (both in-and outflow). The model was executed with an hourly time step and daily (midday) output of model results, except for the high resolution runs, which ran with a ten minute 310 time step and output.
Model results were compared against monthly to semi-monthly data for water temperature ( o C), dissolved oxygen (mg/L), inorganic nutrient concentrations (NO3, NH4 and PO4, mg/L), total nitrogen and phosphorous concentrations (TN and TP, mg/L) and chlorophyll a concentrations (chl. a, µg/L) obtained from the NOVANA program available at www.miljoportal.dk. Spatial resolution of in-lake dataset for calibration and validation varied with 1-12 measurements per sample date across the 320 water column for water temperature and DO (median of 7 measurements per date) and 1-3 samples per date for water nutrient concentrations (median of 2). All chl a. concentrations were sampled in the surface (between 0-3 m depth) once per sample date. We evaluated model performance by computing root-mean-square-error (RMSE) and the coefficient of determination (R 2 ) for the daily output of each state variable. Model results have been processed and visualized in Python with the packages xarray version 0.15.1 (Hoyer and Hamman, 2017), matplotlib version 3.1.2 (Hunter, 2007), and seaborn version 0.11.1 325 (Waskom, 2021), as well as with the open-source Python program PyNcView (available via The Python Package Index, pip install pyncview).

Results of the testcase with example dynamics of new features
The new WET version of Lake Bryrup model differs from the FABM-PCLake version in a few areas, first and foremost by 330 taking advantage of the new options for vertical mobility for motile cyanobacteria and zooplankton, and by allowing dispersion of fish between layers (see Sect. 2.34). In addition, the new modularization have allowed us try out an alternate food web configurations, namely the differentiation of zooplankton into two boxes (mesozooplankton and daphnia), although in this instance we are limited by the availability of data for calibration.
The new model yields comparable or slightly improved overall performance metrics to the ones obtained by Chen et al. (2020), 335 with slightly lower or equal performance statistics for oxygen and nitrogen state variables, but improved performance for phosphorous, chlorophyll a, and temperature (see Table 21).

Phytoplankton modularity
As in Chen et al. (2020), the WET recalibrated Lake Bryrup model contains three phytoplankton groups, diatoms, green algae and cyanobacteria (Fig. 3). The general seasonal phytoplankton succession in Lake Bryrup involves a spring diatom bloom, an early summer bloom of green algae, and a cyanobacteria bloom late summer, or when the water column stratifies for a 345 prolonged period of time. The central panels of Fig. 3 illustrates the separation of multiple phytoplankton groups into separate versions of the same general module, making it easy to add, switch out or remove phytoplankton groups from the model, and to parameterize these individually. The chosen set of phytoplankton categories reflects the available data, and matches the previous FABM-PCLake version of the model.

Phytoplankton nitrogen fixation
To illustrate the interplay/dynamics between nutrient limitation, nitrogen fixation and growth of cyanobacteria with and without nitrogen fixation ability, we included an additional cyanobacteria group with identical parameters to the calibrated cyanobacteria group besides turned on N fixation with default values (lNfix = true and fMuNFix = 0.9). Although phytoplankton groups in Lake Bryrup are P limited in the period where cyanobacteria blooms frequently occur (approx. mid-June to mid-355 August with simulated P limitation between 0.8 to 0.5), N-fixing cyanobacteria dominated the community in contrast to the non N-fixing cyanobacteria ( Fig. 4A and B). Both cyanobacteria groups were N limited in Spring, which allowed N-fixing cyanobacteria to get an advantage by decreasing N limitation via N-fixation (Fig. 4C) and increase biomass concentrations, thereby having a head start for the period with low mixing and warmer surface waters. As expected, N-fixation rates increased with increased N limitation and N-fix. Cyanobacteria biomass concentration ( Figure 4C, secondary y-axis). The relatively low 360 N limitation experienced by the N-fixing cyanobacteria group (N limitation factor between 0.95 to 1.0) resulted a low growth rate penalty during periods with N-fixation between 0 to 1.5%. In scenarios with altered external nutrient loads to switch from a P-limited to a N-limited system in late summer to fall, N-fixing cyanobacteria still dominate the phytoplankton community in late summer with N-fixation rates now significantly increased (for instance 20-fold in one scenario) as the cyanobacteria are now more N limited. 365 in the water column, and panels C-F illustrating the dynamics associated with each of the four vertical movement settings (here 370 only for cyanobacteria). Of these, panel C corresponds to the old FABM-PCLake condition (Hu et al., 2016). Much of the seasonal succession of phytoplankton, and especially the shift from denser-than-water and non-motile phytoplankton (pre-July dynamics of Fig. 5C-F, see also Fig. 3A and B) to a mid-to-late-summer cyanobacterial bloom are driven mainly by the presence/absence and location of lake stratification (as indicated by the temperature gradients in Fig. 5A). (A) and associated presence/absence of lake stratification (B), which, together with sinking, drive the spatial distributions of denser-than-water 375 and non-motile phytoplankton (pre-July dynamics of Fig. 4, panel C and D, see also Fig. 3, panels A and B). The early part of the growing season are dominated by relatively fast-sinking non-motile diatoms and green algae, which have qTrans set to 1 (passive advection and fixed sinking rate, see Sect. 2.34.1.1), in the configuration file (Fig. 3).

Phytoplankton vertical mobility
With the formation of a late summer stable stratification, relatively buoyant or motile cyanobacteria (Fig. 3C, Fig. 5C-FqTrans = 3, phototaxis, see Sect. 4.1.3) become dominant, in part due to their ability to stay in the epilimnion (Fig. 4, panel C). Even 380 under the passive transport setting (Fig. 5C), cyanobacteria are able to sustain a bloom under stratified conditions, due to their slow sinking rate (cVSet = -0.022 m d -1 ). However, when phototaxis is turned on (Fig. 5E), higher biomasses are reached by the cyanobacteria. When nutrient taxis without phototaxis is turned on (Fig. 5D), cyanobacteria aggregate at the bottom of the mixed layer when nutrients at the surface are very scarce (Fig. 5B)Panel D of Fig. 4 illustrate the same dynamics as panel C, with the exception of qTrans now being set to 4 (combined chemo-and phototaxis, see Sect. 4.1.4). Under the combined 385 nutrient and phototaxisthis setting (Fig. 5F, cyanobacteria aggregate around the bottom of the thermocline, forming a deep chlorophyll maximum (DCM), where nutrients are more available, and cyanobacteria are less nutrient limited (i.e. higher simulated nutrient limitation factor) , but are also able to exploit higher light intensities at the surface, resulting in higher biomasses, compared with nutrient taxis alone. (Fig. 4, panel E).

Zooplankton vertical mobility
In WET, heterotrophic pelagic modules such as fish and zooplankton can exhibit vertical movement, in the forms of hypoxia avoidance behavior and light-triggered diel vertical migration (see Sect. 4.22.3.2). Figure 56 illustrates both of these three options for zooplankton and fish vertical mobility behavior in WETs (panel A),. Ttriggered by the up-and downwards migration of the formation of a hypoxic deep zonelayer (panel Fig. 6D), and the daily light cycle (panel Fig. 6FB), the Daphnia 395 in Fig. 6A conducts diel vertical migrations when qTrans is set to 3. The active movement of zooplankton is modulated by advection (panel CFig. 6E), which diffuses part of the migrating zooplankton during days of high turbulent mixing (e.g. on the 18 th 3 rd and 20 th 7 th of July in Fig. 65). In panel B of Fig. 6 (qTrans = 2), upwards swimming is triggered when hypoxic conditions extends up into the water column, which does not happen when vertical movement of Daphnia is limited to passive advection (panel C of Fig. 6, qTrans = 1). 400

Observations on the performance of the new vertical movement and nitrogen fixation features.
Model configuration and complexity is often constrained by data availability. Symptomatically, although Lake Bryrup is 405 included in the Danish long-term ecological monitoring program (NOVANA), there is no data against which to validate, for instance, spatial distributions of higher organisms. Here, we have been mostly concerned with demonstrating the features of WET, however the ability to model diel vertical migration (DVM) will be essential in many applications worldwide, especially large and deep lakes across the globe, as well as in many marine or estuarine environments., and iI In such cases the modeler will often have to rely on studies that specifically target zooplankton (or fish) DVM in similar environments. For fish however, 410 a large step (if somewhat of a low-hanging fruit) towards more realistic model representation has been taken in WET, with the removal of the unrealistic absence of movement between model layers in its predecessor FABM-PCLake model.
While DVM in zooplankton and fish is an elusive dynamic to observe and understand, the importance of vertical movement processes for the composition and seasonal succession of phytoplankton communities is more easily recognized., with mMotile Formatted: Superscript Formatted: Normal phytoplankton having have a distinct advantage in highly stratified conditions, where the ability to stay in the euphotic 415 epilimnion (Wentzky et al., 2020), or to balance opposing gradients of light and nutrients by migrating to the hypolimnion (Leach et al., 2018) can be important. We demonstrated that the choice of DVM method has profound an impacts on the model behavior (Fig. 5), e.g. whether mobile phytoplankton will concentration at the surface layer or near the bottom of the thermocline (i.e. forming a DCM).
The option of nitrogen fixation in phytoplankton is another feature which might improve model performance in stratified or 420 oligotrophic lakes, where nitrogen limitation can be important (Reinl et al., 2021), although the ability of N-fixation to fully compensate for nitrogen limitation has recently been called into question (Shatwell and Köhler, 2019). In the present study, the impact of nitrogen fixation was expected to be low, as Lake Bryrup is predominantly limited by phosphorous, particularly in the main part of the growing season. Nevetheless, by adding a nitrogen-fixing cyanobacterial competitor, we shoved how shifting nutrient conditions throughout the growth season shaped the relative competitive landscape, as well as the nitrogen 425 fixation rate of the N-fixing cyanobacteria throughout the model period.

Advantages of model modularization
As demonstrated here, the WET model can reproduce the behavior of FABM-PCLake. But this food web configuration might not be optimal for every, nor even necessarily in this, use case. While beyond the scope of this paper, the WET model 430 application to Lake Bryrup could potentially be calibrated in several food web configurations to find the optimum conceptual representation. This procedure could in principle be applied in all model applications as an extra layer in the calibration process, but will not be feasible in many cases, in the face of the already daunting task of calibrating a model with hundreds of parameters. Unfortunately, we believe that this is a case where there is no replacement for experience with both models and the study system, and the scientist will have to rely on their expertise However, this also means that there are more choices to 435 make when setting up a model for a new system, and so the modeler is rewarded for detailed to configure the most accurate and realistic food web for any given research questiona priori knowledge of the modelled system and its dynamics. However, on a smaller scale, it will often be possible to test different food web compositions, by e.g. adding or subtracting an extra phytoplankton, zooplankton or fish module to an already calibrated model, and doing simpler recalibrations.
Differences in ecosystem structure and functions between lakes in different climatic regions would also most likely warrant 440 changes to food web configurations. For example, fish populations in (sub)tropical shallow lakes have in general smaller body size, shorter life span, faster growth, multiple reproduction events and stronger preference for the littoral zone compared to temperate lakes (Meerhoff et al., 2012, and references herein). In combination with increased proportion of herbivorous and omnivorous fish species Iglesias et al., 2017), this difference would most likely weaken trophic cascades and hereby diminish the impact of several lake restoration strategies (Jeppesen et al., 2010). So to reproduce 445 for instance warm water shallow lake dynamics and responses to potential restoration efforts, configuration of the food web to the specific lake is likely needed.

Concluding remarks regarding the new features
Overall, the new changes that separate WET from its predecessors provide the model with a high degree of flexibility and 450 adaptability, with the distinct advantage of allowing one model code base to handle many different application cases, instead of requiring many distinct models for different purposes. By taking advantage of the modularization, distinct food web configurations can be set up for different systems. Meanwhile, the new vertical movement and nitrogen fixation algorithms allow the model to be applicable in a much wider array of physical settings, across gradients in latitude, depth or salinity. This flexibility of the new model may contribute to limiting cases of parallel development and 'reinventing the wheel', while 455 promoting comparability between different model implementations. As stated in the introduction, we believe that the consolidation of various model features and approaches into a few flexible and customizable models is a crucial process for the overall progress of the field of ecological modelling. However, for such models to be truly and successfully flexible, such customization must be possible with relative ease, or risk going unused. The modularized structure of WET under the FABM framework supports the user by making everything configurable in a single setup file, by making it easy to switch modules on 460 or off, and by having all options relevant to an organism type available from the base module. By switching out module sections or an entire configuration file, a new model setup can be run with the same executable, without having to change or recompile any source code. For an even simpler workflow, running the model through the QWET plugin for the QGIS software provides the user with a simple step-by-step GUI based process, for which a tutorial is available on the WET homepage.

Future work
WET is under active and continuous development. Currently, effort is being applied to improving improve the applicability of the model to subtropical and tropical regions, which will include improvements to the macrophyte module and support for herbivory in fish. Another work in progress is a complete overhaul of the heterotrophic modules, replacing old feeding 470 formulations with more realistic descriptions, and introducing options for dynamic diets, feeding strategies and foraging efforts of fish, based on optimal foraging theory.
Apart from these new features, other areas for future work include improved handling of resuspension of sediment, a fully size-based fish module, and extensive testing and improvements of the model in a 3D application, with the expressed aim of the authors (in their roles as the current main developers of WET), that the model be even more tailorable to a wide array of 475 ecosystems types, across latitudinal, spatial and productivity gradients, simply by turning features on or off, and combining different modules, all configurable at runtime. With this model, which is open source and freely available, we hope to facilitate the consolidation of successful features of many models together in one, with the goal of preventing 're-inventions of the wheel' in the future, and making aquatic ecosystem modelling easier, more flexible and, ultimately, better.

Author Contributions
The model was coded and developed mainly by NASM, with significant contributions and aid from DT, TKA, FRSH and KB, based on a previous model code by FRSH. The testcase was adapted and carried out by TKA, and NASM prepared the manuscript with contributions from all co-authors. 495

Competing Interests
The authors declare that they have no conflict of interest.