Onur Kerimoglu
email: kerimoglu.o@gmail.com

Stéphan Jacquet

Brigitte Vinçon-Leite

Bruno J Lemaire

Frédéric Rimet

Frédéric Soulignac

Dominique Trévisan

Orlane Anneville

Modelling the plankton groups of the deep, peri-alpine Lake Bourget *

Keywords: coupled physical-biological model, re-oligotrophication, cyanobacteria, mixotrophy, allometry, priority eect, bistability

Predicting phytoplankton succession and variability in natural systems remains to be a grand challenge in aquatic ecosystems research. In this study, we identied six major plankton groups in Lake Bourget (France), based on cell size, taxonomic properties, food-web interactions and occurrence patterns: cyanobacterium Planktothrix rubescens, small and large phytoplankton, mixotrophs, herbivorous and carnivorous zooplankton. We then developed a deterministic dynamic model that describes the dynamics of these groups in terms of carbon and phosphorus uxes, as well as of particulate organic phosphorus and dissolved inorganic phosphorus. The modular and generic model scheme, implemented as a set of modules under Framework for Aquatic Biogeochemical Models (FABM) enables run-time coupling of the plankton module an arbitrary number of times, each time with a prescribed position across the autotrophy/heterotrophy continuum.

Parameters of the plankton groups were mainly determined conjointly by the taxonomic and allometric relationships, based on the species composition and average cellular volume of each group. The biogeochemical model was coupled to the one-dimensional General Ocean Turbulence Model (GOTM) and forced with local meteorological conditions. The coupled model system shows very high skill in predicting the spatiotemporal distributions of water temperature and dissolved inorganic phosphorus for ve simulated years within the period 2004 to 2010, and intermediate skill in predicting the plankton succession. We performed a scenario analysis to gain insight into the factors driving the sudden disappearance of P. rubescens in 2010. Our results provide evidence for the hypothesis that the abundance of this species before the onset of stratication is critical for its success later in the growing season, pointing thereby to a priority eect.

Introduction

Mechanistic ecosystem models are not only ideal media for synthesizing data and theory and thereby improving our understanding of the functioning of ecosystems, but they are also useful tools for supporting decision-making processes. Implementation of ecosystem models in aquatic systems has been increasingly appearing in the form of coupled hydrodynamic-biogeochemical models, reecting the increasing recognition that biogeochemistry is often strongly driven by hydrodynamics, as well as the advances in computing power [START_REF] Robson | When do aquatic systems models provide useful predictions, what is changing, and what is next[END_REF].

While coupled hydrodynamic-biogeochemical models can often provide reliable estimates of physical parameters such as temperature and salinity, prediction of chemical and biological parameters, and especially occurrence of certain plankton species or functional groups are more dicult to predict [START_REF] Shimoda | Phytoplankton functional type modelling: Running before we can walk? A critical evaluation of the current state of knowledge[END_REF]. Problems start there already at the very preliminary stage of conceptual model building: what constitutes a functional group? From a global perspective, [START_REF] Hood | Pelagic functional group modeling: Progress, challenges and prospects[END_REF] dened a functional group as an entity that plays a particular role in a certain biogeochemical pathway, such as nitrogen xation or silication. However, shifts in the community structure under focus might be driven by competition for the shared resources or trophic interactions, intensity of which may change across seasonal to interdecadal scales, e.g., with thermal stratication dynamics and changes in nutrient loading [START_REF] Sommer | Beyond the Plankton Ecology Group (PEG) Model: Mechanisms Driving Plankton Succession[END_REF][START_REF] Kerimoglu | Seasonal, inter-annual and long term variation in top-down versus bottom-up regulation of primary production[END_REF]. In such cases, traits relating to growth rate, grazer defense, resource utilization, temperature response and motility [START_REF] Litchman | Linking traits to species diversity and community structure in phytoplankton[END_REF] should (also) be considered for identifying the functional groups, for which, ample examples indeed exist [START_REF] Jöhnk | Summer heatwaves promote blooms of harmful cyanobacteria[END_REF][START_REF] Mieleitner | Modelling functional groups of phytoplankton in three lakes of dierent trophic state[END_REF][START_REF] Carraro | Coupling high-resolution measurements to a three-dimensional lake model to assess the spatial and temporal dynamics of the cyanobacterium Planktothrix rubescens in a medium-sized lake[END_REF].

The system under investigation here, Lake Bourget, France, has been recovering from eutrophication since the 1980s [START_REF] Vinçon-Leite | Contribution of mathematical modeling to lake ecosystem understanding : Lake Bourget (Savoy, France)[END_REF]Jacquet et al., 2014a). Starting from 1996, the toxic cyanobacterium, Planktothrix rubescens became a dominant species in the lake [START_REF] Vinçon-Leite | Phytoplankton variability in Lake Bourget: Phytoplankton dynamics and meteorology[END_REF][START_REF] Jacquet | The proliferation of the toxic cyanobacterium Planktothrix rubescens following restoration of the largest natural French lake (Lac du Bourget)[END_REF], see also Fig. 1). P. rubescens is a wide-spread cyanobacterial species especially prevalent in peri-alpine lakes (e.g., [START_REF] Ernst | Abundance and toxicity of Planktothrix rubescens in the pre-alpine Lake Ammersee, Germany[END_REF][START_REF] Salmaso | Eects of nutrient availability and temperature on phytoplankton development: A case study from large lakes south of the Alps[END_REF][START_REF] Dokulil | Deep living Planktothrix rubescens modulated by environmental constraints and climate forcing[END_REF][START_REF] Posch | Harmful lamentous cyanobacteria favoured by reduced water turnover with lake warming[END_REF], but also observed elsewhere (e.g., [START_REF] Konopka | Buoyancy regulation and vertical migration by Oscillatoria Rubescens in Crooked Lake, Indiana[END_REF][START_REF] Halstvedt | Seasonal dynamics and depth distribution of Planktothrix spp. in Lake Steinsfjorden (Norway) related to environmental factors[END_REF][START_REF] Naselli-Flores | Toxic Cyanobacterial Blooms in Reservoirs Under a Semiarid Mediterranean Climate: The Magnication of a Problem[END_REF][START_REF] Padisák | Rarity, ecological memory, rate of oral change in phytoplanktonand the mystery of the Red Cock[END_REF]. Being a potentially microcystin producing species [START_REF] Briand | Variations in the microcystin production of Planktothrix rubescens (cyanobacteria) assessed from a four-year survey of Lac du Bourget (France) and from laboratory experiments[END_REF], its occurrence in lakes and reservoirs has been of major concern for livestock and human health [START_REF] Naselli-Flores | Toxic Cyanobacterial Blooms in Reservoirs Under a Semiarid Mediterranean Climate: The Magnication of a Problem[END_REF][START_REF] Ernst | Abundance and toxicity of Planktothrix rubescens in the pre-alpine Lake Ammersee, Germany[END_REF]. In 2010, P. rubescens suddenly disappeared (Jacquet et al., 2014b) from Lake Bourget, whereas the mixotrophic and small phytoplankton species became relatively more abundant, latter being typical for oligotrophic systems [START_REF] Anneville | Seasonal and inter-annual scales of variability in phytoplankton assemblages: comparison of phytoplankton dynamics in three peri-alpine lakes over a period of 28 years[END_REF][START_REF] Chen | Relationships between phytoplankton growth and cell size in surface oceans: Interactive eects of temperature, nutrients, and grazing[END_REF][START_REF] Mitra | The role of mixotrophic protists in the biological carbon pump[END_REF]. Accordingly, mixotrophy and traits associated with cell size should be taken into account for understanding the mechanisms driving the changes in phytoplankton community composition, and in particular, the disappearance of P. rubescens in Lake Bourget.

Since long, phytoplankton cell size has been recognized to be an important aspect in determining the ecophysiology of phytoplankton (e.g., [START_REF] Finkel | Phytoplankton in a changing world: cell size and elemental stoichiometry[END_REF]Litchman et al., 2010, and references therein), and cell size has been increasingly used as a `master trait' [START_REF] Litchman | Linking traits to species diversity and community structure in phytoplankton[END_REF] in theoretical modelling studies (e.g., [START_REF] Grover | Resource competition in a variable environment: phytoplankton growing according to the variable-internal-stores model[END_REF][START_REF] Armstrong | Grazing limitation and nutrient limitation in marine ecosystems: Steady state solutions of an ecosystem model with multiple food chains[END_REF][START_REF] Litchman | Contrasting size evolution in marine and freshwater diatoms[END_REF]Kerimoglu et al., 2012, submitted), although the integration of size concept in realistic ecosystem models attempting to reproduce mesocosm or eld observations at relevant ecological time scales has been gaining momentum only recently (but see, e.g., [START_REF] Ward | A size-structured food-web model for the global ocean[END_REF][START_REF] Wirtz | Mechanistic origins of variability in phytoplankton dynamics: Part I: niche formation revealed by a size-based model[END_REF][START_REF] Terseleer | Trait-based representation of diatom functional diversity in a plankton functional type model of the eutrophied southern North Sea[END_REF]. Following many decades of dichotomous classication of planktonic organisms as `autotrophs' and `heterotrophs' [START_REF] Flynn | Misuse of the phytoplankton-zooplankton dichotomy: the need to assign organisms as mixotrophs within plankton functional types[END_REF], importance of mixotrophy in ecosystem functioning has been increasingly recognized (Mitra et al., 2014, and references therein). Mixotrophy has been addressed mainly by theoretical work so far (e.g., [START_REF] Thingstad | On the strategy of 'eating your competitor': a mathematical analysis of algal mixotrophy[END_REF][START_REF] Flynn | Building the perfect beast: modelling mixotrophic plankton[END_REF][START_REF] Crane | Coexistence of mixotrophs, autotrophs, and heterotrophs in planktonic microbial communities[END_REF][START_REF] Berge | Modeling succession of key resourceharvesting traits of mixotrophic plankton[END_REF]. The recent work of [START_REF] Ward | Marine mixotrophy increases trophic transfer eciency, mean organism size, and vertical carbon ux[END_REF] constitutes the rst example where mixotrophy is resolved in a global ocean model.

Environmental control of the occurrence of P. rubescens blooms is still under debate. Analysis of the longterm changes in individual lakes and inter-comparison between lakes suggest that phosphorus (and nitrogen, see [START_REF] Posch | Harmful lamentous cyanobacteria favoured by reduced water turnover with lake warming[END_REF] availability is the primary determinant [START_REF] Dokulil | Deep living Planktothrix rubescens modulated by environmental constraints and climate forcing[END_REF]Jacquet et al., 2014b;[START_REF] Anneville | Blue-Green Algae in a Greenhouse Century? New Insights from Field Data on Climate Change Impacts on Cyanobacteria Abundance[END_REF]. P. rubescens is characterized by their slow growth rates [START_REF] Bright | The daily integral of growth by Planktothrix rubescens calculated from growth rate in culture and irradiance in Lake Zürich[END_REF] and tolerance to low light conditions [START_REF] Walsby | Light-dependent growth rate determines changes in the population of Planktothrix rubescens over the annual cycle in Lake Zürich, Switzerland[END_REF]. As a result of the latter, it often develops thin and intense layers within the metalimnion during the growing season (eg., Jacquet et al., 2014b), where they are the rst to harvest the nutrients leaking from the nutrient-rich layers beneath the thermocline. On the other hand, occurrence of P. rubescens displays extreme interannual variability in deep lakes, which is suggested to be driven by meteorological conditions [START_REF] Vinçon-Leite | Phytoplankton variability in Lake Bourget: Phytoplankton dynamics and meteorology[END_REF][START_REF] Salmaso | Long-term phytoplankton community changes in a deep subalpine lake: responses to nutrient availability and climatic uctuations[END_REF]Jacquet et al., 2014b;[START_REF] Anneville | Blue-Green Algae in a Greenhouse Century? New Insights from Field Data on Climate Change Impacts on Cyanobacteria Abundance[END_REF]: after warm winters, P. rubescens have been observed to be at relatively high abundances, which is usually followed by their sustained dominance throughout the growing season [START_REF] Salmaso | Long-term phytoplankton community changes in a deep subalpine lake: responses to nutrient availability and climatic uctuations[END_REF][START_REF] Posch | Harmful lamentous cyanobacteria favoured by reduced water turnover with lake warming[END_REF].

In this study we had the following objectives: 1) developing a plankton model that resolves mixotrophy and relies on allometric relationships for the parameterization of plankton groups; 2) implementing the plankton model in a 1-D coupled hydrodynamical-biogeochemical modeling framework for the simulation of the surface layer dynamics of Lake Bourget; 3) quantifying the performance of the model with respect to both the bulk characteristics of the system and plankton groups 4) gaining insight into the reasons for the relative importance of low phosphorus concentrations and low starting inoculum at the beginning of the growth season for the disappearance of P. rubescens in Lake Bourget.

Study Site and Data

The study site is Lake Bourget, a peri-alpine ( 45• 44 N, 5

• 52 E, 231 m altitude) lake with a maximum depth of 145 m and a surface area of approximately 45 km 2 . It has a north-south aligned, elongated basin with a length of 18 km and a maximal width of 3 km at the surface. Within the study period (2004)(2005)(2006)(2007)(2008)(2009)(2010), average total phosphorus and nitrogen concentration ranged between approximately 15-45 mg/m 3 and 550-800 mg/m 3 , respectively (Fig. 1), classifying it as a mesotrophic system. Further details about the lake can be found in [START_REF] Vinçon-Leite | Contribution of mathematical modeling to lake ecosystem understanding : Lake Bourget (Savoy, France)[END_REF]Jacquet et al., 2014b).

Meteorological data required to force the hydrodynamical model (see Section 4.3) were taken from the Météo-France station at Vouglans, located at the southern tip of Lake Bourget. In-situ data used in this study were sampled at the deepest location of Lake Bourget. Sampling was performed usually biweekly during the growth season and monthly during winter. Water temperature was measured at high vertical resolution with a conductivity-temperature-depth probe, and interpolated on a regular grid of 1-m intervals.

Nutrient data were collected at several irregular depths (2, 10, 15, 20, 30, 50, 80, 110, 130, 140 m), and interpolated also to a regular 1-m grid. For the phytoplankton species counts, an integrating sampler was used for the 2.5 times Secchi depth in 2005 (typically between 10-20m) and for the top 20 m after 2005.

For the zooplankton species counts, depth-averaged samples were taken with 0-50m hauls. Further details on the regular sampling programme of Lake Bourget can be found in (Jacquet et al., 2014a). Conversion of the phytoplankton species counts to carbon biomass was based on biovolume of each species [START_REF] Druart | Protocoles d'analyse du phytoplancton de l'INRA: prélèvement, dénombrement et biovolumes[END_REF], and assuming 0.2 pgC/µm 3 . Zooplankton counts were rst converted into wet weight, 10% of which was assumed to be dry weight [START_REF] Dumont | The dry weight estimate of biomass in a selection of Cladocera, Copepoda and Rotifera from the plankton, periphyton and benthos of continental waters[END_REF], from which the carbon biomasses were calculated assuming that carbon constitutes 48% of dry weight [START_REF] Andersen | Carbon, nitrogen, and phosphorus content of freshwater zooplankton[END_REF].

3 Identication of functional groups

Algae and Mixotrophs

We dened 3 functional algal groups and a mixotroph group. The rst group consists of only P. rubescens (abbreviated AP R hereafter), reecting the unique eco-physiological properties of this species, and its abundance in Lake Bourget (Table 1). Second is the `small algae' group (AS ), consisting of all species with cell volume up to 10 3 µm 3 , and not being mixotrophs (see below), comprising cyanobacteria (except P.

rubescens), diatoms, chlorophyta and chrysophyta (Table A2). Third is the `large algae' group (AL), those with cell volume larger than 10 3 µm 3 , and not being mixotrophs. Finally, the `mixotroph' group (M ), consists of species with cell volume ≥ 10 2 µm 3 and belonging to Chrysophyta, Dinophyta and Cryptophyta.

In

Zooplankton

We considered two zooplankton groups: herbivores (ZH ),dominated by Daphnia sp., Diaphanosoma and Eudiaptamus species (Table A1) and are considered to feed on phytoplankton and detritus (see the model description, section 4.1); and carnivores (ZC ), dominated by cyclopoid copepodits of stages C5 and C6, and are considered to feed on herbivorous zooplankton, mixotrophs and large algae.

The Model

As a spatially explicit model, temporal and spatial changes in the volumetric concentration of a biological state variable, ci with units mmol m -3 , are described by a system of partial dierential equations of the generic form:

∂ci ∂t + ∂ ∂z ωi(t, z)ci(t, z) -Kz(t, z) ∂ci(t, z) ∂z = s(ci(t, z)) (1)
where, wi the vertical motion, Kz the eddy diusivity and s(ci) the net source term (productiondestruction) resulting from biological interactions.

Structure of the biological model is depicted in Fig. 2. As phosphorus is recognized to be the primary determinant of phytoplankton dynamics in Lake Bourget [START_REF] Jacquet | The proliferation of the toxic cyanobacterium Planktothrix rubescens following restoration of the largest natural French lake (Lac du Bourget)[END_REF](Jacquet et al., , 2014b)) The complete set of equations describing the dynamics of carbon (X C ) and phosphorus (X P ) bound to plankton, for Xi={APR, AL, AS, M, ZH, ZC}, particulate organic phosphorus (P OP ) and dissolved inorganic phosphorus (DI P ) is given in eq. 2a-2d.

s(X

C i ) =min(µ P , µ C )X C i -(ei + Li)X C i - j (1 -ζj)(G j,k=i X C j ) (2a) s(X P i ) = P i ζiUi + (1 -ζi) k G j=i,k Q k X C i -(ei + Li)X P i - j (1 -ζj)G j,k=i Q k X C i (2b) s(P OP ) = i LiX P i - j (1 -ζj)G j,k=P OP X C j -rP OP (2c) s(DIP ) = i (1 -ζi)(1 -P i )X C i k G j=i,k Q k - i P i ζiUiX C i + i eiX P i + rP OP (2d)
Denition and units of state variables, major processes, intermediate quantities and parameters are provided in Tables-2-5. Additional terms are introduced upon their appearance when necessary. In eq.2a-2d and below, index i stands for each of the 6 plankton groups considered here (Fig. 2), and indices j and k stand respectively for the predator and prey items.

µ at Q → ∞ 1.2 1 0.39 1.22 - - α d -1 W -1 m 2 slope
A S A PR A L M Z H M 0.6 0.4 - - - - Z H 0.2 0.5 0.1 0.2 - - Z C - - - 0.2 0.2 0.6
if Qmin < Qmax :

µ P = µ∞fQ = µ∞ 1 - Qmin Q (3)
where fQ describes quota-dependent nutrient limitation of biomass growth. For pure heterotrophs (ζ = 0, Table 3), which are assumed to be homeostatic (Qmax = Qmin, Table 3), nutrient limitation is assumed to be absent, i.e. only the carbon limited growth is taken into account (µ

P = ∞).
Carbon limited growth rate is calculated as the sum of specic carbon xation rate through photosynthesis and carbon assimilation rate through heterotrophy:

µ C = ζiPi + (1 -ζi) C i k G j=i,k (4)
Photosynthetic carbon xation rate is calculated as a function employed by [START_REF] Schwaderer | Eco-evolutionary dierences in light utilization traits and distributions of freshwater phytoplankton[END_REF] that accounts for inhibition at high irradiances:

P = µmaxfI = µmaxIP AR I 2 µmax αI 2 opt + I(1 -2 µmax αI opt ) + µmax α (5)
where, µmax = µP (Qmax) (eq.3), fI describe light limitation and IP AR stands for photosynthetically active radiation.

Heterotrophic carbon assimilation rate in eq. 4 is given by the product of total ingestion rate and carbon assimilation eciency, C (Table 3 andsection 4.2.2). Ingestion rate of the each prey item is described by:

G j=i,k = Gi,max p j=i,k X C k KG,i + k p j=i,k X C k (6)
where, p j,k is the preference of zooplankton j for the food item k (Table 5), and X C k=P OP =POC was calculated as P OC = P OP/QP OM , assuming Redeld ratio for QP OM (Table 4). Copepods are known to switch between ltering and raptorial feeding behaviors based on the food availability [START_REF] Kiørboe | How zooplankton feed: mechanisms, traits and trade-os[END_REF], which is here accounted by adjusting the feeding preferences of ZC based on relative abundance, as in [START_REF] Fasham | A nitrogen-based model of plankton dynamics in the oceanic mixed layer[END_REF]:

p j,k = p j,k X k k p j,k X k (7)
where p j,k replaces p j,k in eq.6.

Analogous to eq.2a, the rst term in eq.2b describes the cumulative phosphorus binding rate through uptake of dissolved nutrients and nutrient assimilation rate of the ingested food. Uptake rate of dissolved nutrients is represented according to the Michaelis-Menten kinetics:

U = UmaxfN = Umax DIP DIP + KU (8)
where, fN describes nutrient limitation of nutrient uptake rate. Regulation of nutrient assimilation diers between homeostatic and non-homeostatic organisms in our model: In the non-homeostatic case, (i.e., Qmin < Qmax) as in phytoplankton and mixotrophs, nutrient uptake rate is adjusted by P , as a linear function of nutrient quota as described in eq.9:

if Qmin < Qmax :

P = Qmax -Q Qmax -Qmin (9)
which stems from the Droop model describing down regulation of nutrient uptake [START_REF] Morel | Kinetics of nutrient uptake and growth in phytoplankton[END_REF] the lower boundary of Umax being 0), and generalized by [START_REF] Crane | Coexistence of mixotrophs, autotrophs, and heterotrophs in planktonic microbial communities[END_REF] to apply also for assimilation of ingested nutrients (eq.2b). The nal, quota-adjusted nutrient uptake rate appears as a sink term in eq. 2d, summed across all plankton (multiplied by their autotrophic fractions).

In the case of homeostatic regulation of carbon and phosphorus assimilation (i.e., Qmin = Qmax), the guiding assumption is that stoichiometry of elemental gains should be equal to the stoichiometry of the organism (eq.10). There are then two possibilities: If the carbon gain would exceed phosphorus gain with the default assimilation eciencies, carbon assimilation eciency is re-adjusted (eq. 11). This implies decreasing specic growth rate with decreasing Q k (of prey) when Q k < Qj for a xed prey concentration, which is in line with observations [START_REF] Hessen | Ecological stoichiometry: An elementary approach using basic principles[END_REF]. On the other hand, if phosphorus gain exceeds carbon gain, P-assimilation eciency is re-adjusted (eq. 12), similar to [START_REF] Grover | Stoichiometry, Herbivory and Competition for Nutrients: Simple Models based on Planktonic Ecosystems[END_REF].

if Qmin = Qmax = Q : P i k G j=i,k Q k C i k G j=i,k = Qi (10) if Qi C i k G j=i,k > P i k G j=i,k Q k : C i = P i k G j=i,k Q k Qi k G j=i,k (11)
else :

P i = Qi C i k G j=i,k k G j=i,k Q k (12)
In both homeostatic/non-homeostatic cases, the un-assimilated portion of the ingested or taken-up phosphorus ((1

-P i ) k G j=i,k Q k ) is added back to the DIP pool (eq.2d).
Last two terms in eq. 2a-2b correspond to excretion, mortality and losses caused by predation by other plankton groups. Excretion is assumed to be associated with heterotrophy (0 for purely autotroph plankton), occurs at a xed specic rate (Table 3) and is recycled back as DIP (eq.2d). Planktonic losses due to mortality is recycled as P OP , and is the sum of a constant rate and a specic rate that linearly scales with X C , accounting for losses to parasites and unresolved higher predators [START_REF] Steele | The role of predation in plankton models[END_REF]:

Li = li + l q i X C i ( 13
)
Gas vesicles of P. rubescens are known to collapse at high pressures (70-90 m in the strains found in Lake Zurich [START_REF] Walsby | The critical pressures of gas vesicles in Planktothrix rubescens in relation to the depth of winter mixing in Lake Zurich, Switzerland[END_REF]. In reality, this only aects their buoyancy, and consecutively leads to an increased sedimentation of P. rubescens, but as we do not explicitly account for buoyancy in this study, we mimic this eect by parameterizing the mortality rate of AP R as a sigmoidal function of depth, that has the iniction point at z=40 m, converging to l i=P R (=0.02 d -1 , as all other plankton, see Table3) at the surface, and 5*l i=P R (=0.1 d -1 ) towards 80 m:

li=P R = l i=P R 1 + 4 1 1 + e 0.1(40-z) (14)
Ingestion appears as a sink term for P OP (eq.2c), which is assumed to be a food source for mixotrophs (Table 5). Last terms in 2c-2d describe remineralization of organic phosphorus O P into dissolved inorganic form DIP at a xed rate r.

Finally, temperature dependence of all reaction rates was included in the model through the Q10 rule:

rate(T ) = rate(T = T ref )fT (15) fT (Q10) = Q10 (T -T ref )/T ref (16)
where, T is the ambient water temperature in °C, which is provided by the physical model. For remineralization (r), eq.15 directly applies, where fT in eq16 is computed with the Q10B for bacteria (Table 4), whereas for plankton, fT in eq.15 is replaced by a response function f T obtained by weighing the au- totrophic and heterotrophic response functions with the corresponding autotrophy (ζ ) and heterotrophy (1-ζ ) fractions:

f T = ζfT (Q10A) + (1 -ζ)fT (Q10H ) (17)
where, Q10A and Q10B are the Q10 values for autotrophs and heterotrophs (Table 4).

Parameterization 4.2.1 Algae and mixotrophs

Parameterization of processes for the algal and mixotrophic groups was based on recent trait-based studies, which consider size as a fundamental trait, but recognize also the taxa-specic dierences.

Umax, Qmin

These parameters were estimated allometrically using the biovolume-weighted average cell volume calculated for each group (Table 6). For AL, which is mainly composed of diatoms (Table A2), we used the scaling coecients specically for freshwater diatoms found by [START_REF] Litchman | Contrasting size evolution in marine and freshwater diatoms[END_REF]. For the other groups, we use the coecients given by [START_REF] Edwards | Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton[END_REF]. Cell-specic values were converted to C-specic values by allometrically scaling the carbon contents using the coecients provided by [START_REF] Menden-Deuer | Carbon to volume relationships for dinoagellates , diatoms , and other protist plankton[END_REF]. [START_REF] Crane | Coexistence of mixotrophs, autotrophs, and heterotrophs in planktonic microbial communities[END_REF] suggested Qmin of mixotrophs should be 20% higher due to costs of maintaining the organismal apparatus required for herbivory. Here we generalize this by assuming

Q real min = Qmin + (Qmax -Qmin)(1 -ζ) (after estimating Qmax from Qmin, see below), such that at ζ = 0, Q real min = Qmax, i.e.
, pure herbivores are homeostatic. Qmax

Compared to other parameters, signicant allometric relationships for phosphorus-Qmax are scarce. Data set collected by [START_REF] Litchman | Contrasting size evolution in marine and freshwater diatoms[END_REF] for marine diatoms suggests scaling coecients for phosphorus-Qmax to be almost identical to that of Qmin, i.e., Qmax being proportional to Qmin with a proportionality constant of Qmax/Qmin = 10 -9.32+10.6 = 19.05. Considering typical values used in literature (e.g., [START_REF] Gal | Implementation of ecological modeling as an eective management and investigation tool: Lake Kinneret as a case study[END_REF], we assumed a more modest storage capacity of Qmax/Qmin = 10 . µ∞

Given the average volumes and the taxa of the species involved in these groups (3), rst the µmax values were visually determined from [START_REF] Edwards | Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton[END_REF], Fig3C, to be 10 0.1 (=1.1), 10 -0.1 (=0.9) and 10 -0.25 (=0.55) respectively for As, A l , and M (for the reference temperature: 20°C). Maximum growth rates for Apr was determined to be 0.35, considering relatively low growth rates [START_REF] Bright | The daily integral of growth by Planktothrix rubescens calculated from growth rate in culture and irradiance in Lake Zürich[END_REF]. Then, µ∞ values were calculated using these µmax, allometrically calculated Qmin and Qmax, and the property µmax = µ∞(1 -Qmin/Qmax).

KU

We could not nd any signicant allometric relationship for KU (regarding phosphorus uptake), therefore we adjusted these parameters taking into consideration the taxonomical averages reported by [START_REF] Edwards | Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton[END_REF], Fig. 4E.

α, Iopt

There does not seem to exist a signicant allometric relationship for the parameters regarding light utilization, therefore corresponding parameters were determined based on taxonomic statistics provided by [START_REF] Schwaderer | Eco-evolutionary dierences in light utilization traits and distributions of freshwater phytoplankton[END_REF], Fig. 1.

Zooplankton

Gmax

Average biovolume-weighted cell volume of the mixotrophs, M, is 1.3 10 3 µm 3 (Table 1). Considering the species composition of zooplankton groups (Table A1) and body volume for individual species (e.g., [START_REF] Hansen | Zooplankton grazing and growth: Scaling within the 2-2,000-µm body size range[END_REF], we estimate the average volume of ZH and ZC respectively to be 10 7 and 10 8 µm 3 .

We then allometrically scale Gmax of mixotrophs and zooplankton, respectively using the coecients for agellates and for all herbivore groups provided by [START_REF] Hansen | Zooplankton grazing and growth: Scaling within the 2-2,000-µm body size range[END_REF], as listed with converted units in Table 6. (2000), where : protists<3000µm 3 ; : diatoms >3000µm 3 ; ª: dinoagellates; G max from [START_REF] Hansen | Zooplankton grazing and growth: Scaling within the 2-2,000-µm body size range[END_REF] for the original units of /h, where ð: all agellates, : all groups. C was set for each heterotroph from the relationship [START_REF] Straile | Gross growth eciencies of protozoan and metazoan zooplankton and their dependence on food concentration, predator-prey weight ratio, and taxonomic group[END_REF], where GGE refers to gross growth eciency and NGE refers to net growth eciency. GGE's were taken from [START_REF] Straile | Gross growth eciencies of protozoan and metazoan zooplankton and their dependence on food concentration, predator-prey weight ratio, and taxonomic group[END_REF] to be 0.40, 0.3 and 0.3 respectively for M, ZH and ZC; and NGE (growth/(growth+metabolic losses))

U max Q min Q C G max X i
C =GGE/NGE
was estimated to be N GE = Gmax/(Gmax + m + e). p j,k Preferences p j,k were set (Table 5) based on knowledge about the feeding interactions in general and in Lake Bourget when available.

Qmin, Qmax

Unlike the autotrophs, zooplankton are homeostatic [START_REF] Hessen | Ecological stoichiometry: An elementary approach using basic principles[END_REF]. The plankton model described above acts homeostatically when Qmin= Qmax, values of which, for herbivores and carnivores are respectively based on D. longispina and as an approximate average of Heterocope and Acanthodiaptomus [START_REF] Andersen | Carbon, nitrogen, and phosphorus content of freshwater zooplankton[END_REF]. -1 for non-algal particles, i.e., P OC here (which was estimated from P OP assuming the Redeld ratio). The same study provides a range from 0.012 to 0.045 m 2 mmolC -1 for various phytoplankton species, but for the sake of simplicity, we assume here a value of 0.03 m 2 mmolC -1 for all phytoplankton species.

Q10A, Q10H , Q10B

In order to account for the temperature dependence of reaction rates, Q10-rule was employed in this study, for which, the coecients were taken from [START_REF] Eppley | Temperature and phytoplankton growth in the sea[END_REF] for autotrophs and [START_REF] Hansen | Zooplankton grazing and growth: Scaling within the 2-2,000-µm body size range[END_REF] for heterotrophs (table 4). For bacteria, Q10 coecient (used only for remineralization of organic matter r) was assumed to be identical to that of autotrophs. Especially phytoplankton is known to be inhibited at high temperatures [START_REF] Butterwick | Diversity in the inuence of temperature on the growth rates of freshwater algae, and its ecological relevance[END_REF], but we had to ignore temperature inhibition in this study, as we are lacking the data required to parameterize such functions (e.g., [START_REF] Jöhnk | Summer heatwaves promote blooms of harmful cyanobacteria[END_REF].

Remaining parameters were either set to values commonly used in similar modelling studies (r,e,m) or manually adjusted based on the qualitative knowledge (w,p jk ,m q ).

Physical model

GOTM is a 1-D hydrodynamic model of a water column based on the Reynolds-averaged Navier-Stokes equations and a repository of turbulence closure models. Details about GOTM can be found in [START_REF] Burchard | Description of a exible and extendable physical-biogeochemical model system for the water column[END_REF]. For the application of GOTM to Lake Bourget, eddy diusivity was calculated using kclosure using default model parameters. Air-water uxes were calculated according to the bulk formula of [START_REF] Kondo | Air-sea bulk transfer coecients in diabatic conditions[END_REF] and back-radiation was calculated according to [START_REF] Bignami | Longwave radiation budget in the Mediterranean Sea[END_REF]. All the required atmospheric forcing data were provided in hourly resolution. Short wave radiation was calculated by GOTM as a function of cloud cover, solar altitude and albedo correction [START_REF] Payne | Albedo of the Sea Surface[END_REF] as described by [START_REF] Rosati | A general circulation model for upper ocean simulation[END_REF].

In GOTM, the only source term for heat arises from the attenuation of the shortwave radiation I across the water column, described as the sum of photosynthetically non-active and active radiation:

I(z) = I0ae -z η 1 + I0(1 -a)e -z η 2 -0 z i k c,i c i (z )dz ( 18
)
where I0 is the short wave radiation at the surface, a is the weighting parameter for the dierential attenuation of the red and blue-green wavelength components of the light spectrum, η1 and η2 are the corresponding absorption length scales, and ci (as in eq. 1) and kc,i are, respectively the concentration and the specic extinction coecient of item i, constituting the feedback of the coupled biological model to the physical model. The second term in eq.18 corresponds to IP AR (in eq.5). To account for the background turbidity due to suspended material not explicitly accounted for by the model, Jerlov-IA optical class was assumed, corresponding to a = 0.62, η1 = 0.6 and η2 = 20 [START_REF] Paulson | Irradience measurements in the upper ocean[END_REF].

Model coupling and operation

The biological model was developed within the Framework of Aquatic Biogeochemical Models (FABM [START_REF] Bruggeman | A general framework for aquatic biogeochemical models[END_REF]. FABM acts as an interface between the biogeochemical model and physical model, GOTM which acts as the host of the biological model, by performing the numerical integration of the advection-diusion-reaction equations (eq.1). According to the on-line coupling scheme employed; state variables of the biogeochemical model are advected based on the the settling rates (ω in eq.1) specied by the biogeochemical model and vertical diusion rates as a function of eddy diusivity (Kz in eq.1) and vertical gradients (eq.1); ambient water temperature (T ) and (IP AR) estimated by GOTM is provided to the biological model; and in return, absorption of heat (eq.18) is inuenced by the concentration of the biological state variables. Through the GOTM-FABM coupler interface (a Fortran namelist le), discretization scheme for advection was chosen as the ULTIMATE QUICKEST algorithm and the source-sink dynamics was chosen to be discretized by the 4 th order Runge-Kutta scheme. Simulations were ran with an integration time step of 300 seconds and a vertical resolution of 1 m. Changing the vertical resolution to 0.5 m and time step to 100 seconds was observed to make no dierence.

A number of processes that are required for an accurate simulation of the winter mixing are not resolved by the coupled physical-biological model system. These include the variation of the mass and energy content of each vertical layer with depth (i.e., lake hypsography), benthic-pelagic exchange at the lake bottom, various overwintering strategies of plankton, and the pressure sensitivity of the gas vesicles of P. rubescens, which determine their ability to regulate their buoyancy later in the season (which itself is also not resolved anyway).Therefore, simulations were started each year shortly before the onset of stratication, precisely on a date when eld data that can be used as initial conditions exist. For the initial conditions for temperature, vertically resolved proles were used. For the initial phytoplankton and zooplankton concentrations, i.e., X C i,0 , for which only integrated data respectively from top 20 m and 50 m are available, we assumed that the measured concentrations were homogeneous throughout the water column. For the initial concentration of dissolved nutrients, i.e., DIP0, no dierence was observed between using vertical proles and homogeneous distribution (obtained by averaging the vertical prole), so we assumed homogeneous distributions, which facilitated testing the model sensitivity to phosphorus availability in the system. The concentration of initial particulate organic phosphorus, P OP0 was calculated as P OP0 = T P -DIP0i X P i,0 , where TP stands for measured total phosphorus and X P i,0 = X C i,0 * Qmax,i implying that the plankton were at their maximum quota, as a result of being exposed to high nutrient concentrations throughout the winter. For all the biological variables, no-ux boundary condition was assumed both for the surface and bottom of the water column.

Skill assessment metrics

For evaluating the success of the physical model, we used Taylor & Target diagrams [START_REF] Jolli | Summary diagrams for coupled hydrodynamic-ecosystem model skill assessment[END_REF].

In Target diagrams, standard deviation normalized model bias (B*), and standard deviation normalized unbiased (calculated from the anomalies around the means) root mean square deviation (RMSD), multiplied by the sign of the dierence between the standard deviation of model and observation are mapped on Cartesian coordinates. In Taylor diagrams, correlation coecients between the observations and simulations, and the standard deviation of the model, normalized to that of the observations (equal at 1) are mapped as, respectively, the angle and radius on polar coordinates.

Results

Temporal occurrence and taxonomic composition of each phytoplankton group for the simulated years are shown in Fig. 3. Occurrence of the small algae, AS do not seem to follow a systematic pattern: while in 2006 and 2009, this group, represented mainly by small diatoms, was the major contributor of the spring blooms, in 2005, mainly consisting of cyanobacteria, they were responsible for an intense autumn bloom and nally in 2004 and 2010 this group made up a relatively low and stable background algal biomass throughout the year. AL were in most years responsible for the spring blooms with a major contribution by diatoms, but in 2004, they were abundant also throughout the summer and autumn. Majority of mixotrophs, M , appeared in summer, without any obvious systematic pattern with regard to their taxonomic constituents.

Phytoplankton growth and limitation functions with the parameters specied for each group listed in Table . 3 are shown in Fig. 4, facilitating a comparison between the groups and understanding the succession patterns shown in Fig. 3. AS is characterized by highest growth rate and nutrient uptake rates, partially explaining its abundance in spring and autumn. AL also has relatively high growth rates, but it has signicantly slower nutrient uptake rates, explaining its absence in summer. AP R is characterized by the lowest growth rates, explaining its delayed growth in the season in most years, and lowest light requirements and greatest high-irradiance intolerance explaining its growth in deeper layers (Jacquet et al., 2014b). M grow almost as slow as the AP R, and is poor in both light and nutrient acquisition. 34. f T , f Q , f I and f N stand for temperature limitation (eq. 16), nutrient limitation of biomass growth (eq. 3), light limitation of photosynthesis (eq. 5) and nutrient limitation of uptake (eq. 8), respectively. The match between the measured and simulated water temperatures is outstanding. For DIP in the top 20 meters, the spring draw-down, and late-spring replenishment in some years was reproduced, although the latter was underestimated in 2010. The gradual withdrawal of TP throughout the season is also well reproduced. Slightly higher model estimations at the beginning of the simulations reect uncertainties regarding the phosphorus content of the cells, which we assumed to be at Qmax, and the value of Qmax especially for P. rubescens, which is usually the most abundant species in the system by the end of winter.

Skill metrics for both DIP and TP calculated for the upper 20 m are good (Fig. 8), characterized by near-zero model bias and correlation coecients at around 0.8.

Model estimates of the lumped plankton groups, namely the phytoplankton, mixotrophs and zooplankton are slightly better than those of individual groups, indicated by lower normalized bias, and modeled variability closer to the observations (Fig. 8). Annual average abundance of all groups are reasonably predicted by the model, indicated by low model bias (Fig. 8). The model also provides a decent estimate of seasonal succession, although some events are reproduced with delays, resulting in phase errors, hence relatively low correlation coecients and large unbiased RMSD. This is especially true for AS. An ubiquitous dierence between the measured and simulated plankton abundances is the frequency of temporal uctuations: the simulated trajectories are substantially smoother than the observations, such as those for P. rubescens.

Inter-annual variability is present in model estimations for some variables, but for some others, like the mixotrophs and the ZH , the model seems to be repeating an identical seasonal cycle each year. Seasonal amplitudes of ZH , M and AP R were underestimated, whereas that of AS was overestimated, leading to lower standard deviations in comparison to the observations. Vertical distribution of particularly water temperature, but also the DIP and to a lesser extent, TP are realistically represented, as indicated by the relatively high skill scores attained by the point-wise comparisons (Fig. 9). Depth of the mixed layer, and its seasonal changes throughout the year are accurately represented by the model (Fig. 10). Although the very high nutrient concentrations within the last 10 meters of the water column are not captured (not shown), and the nutrient concentrations at the surface are slightly underestimated, the general features, like the depth of the nutrient-depleted layer and the dierence between the DIP concentrations at deeper layers in 2004 and 2010 are reproduced reasonably well (Fig. 10). Relatively low skill score attained by TP seem to be caused by an unrealistic accumulation of organic material near the thermocline (Fig. 10), which is mainly driven by AP R, AL and M (Fig. 5).

In Jacquet et al. (2014b), it was hypothesized that, among others, two important factors for the disappearance of P. rubescens in Lake Bourget in 2010 was decreasing winter-phosphorus concentrations and too small winter inoculum. Despite the fact that no P. rubescens was detected in the phytoplankton assemblages throughout the year 2010, the model was initiated with a winter inoculum of 10 -2 mmolC/m 3 , which is close Figure 9: Target-Taylor diagrams for temperature, dissolved phosphorus and total phosphorus within the upper 80 m, depth-resolution retained, unlike in Fig. 8 to the minimum observed concentration in 2004. Starting from this small inoculum, P. rubescens biomass did not grow substantially (Fig. 6)). In order to gain some insight into the reasons that eliminated P. rubescens from the system, a scenario analysis was performed, where the initial, winter concentrations of P. rubescens and DIP were replaced with the values observed in 2004 (Fig. 11). Higher initial DIP (0.91 instead of 0.46 mmol/m

3 ) resulted in slightly higher P. rubescens concentrations, although the dierence was not appreciable. On the other hand, higher winter inoculum (0.62 instead of 0.01 mmol/m 3 ) resulted in larger than 10 mmol/m 3 of P. rubescens during autumn. Setting both the DIP and AP R concentrations to 2004 values resulted in dynamics very similar to those obtained in 2004, suggesting that other dierences between these two years, e.g., with respect to the meteorological conditions and initial concentration of other plankton groups, were not critical in determining P. rubescens dynamics.

Discussion

In this study, we described a novel biogeochemical model coupled to a 1-D hydrodynamical model, and its application to Lake Bourget. We demonstrate that the model is able to provide a reasonable representation of the seasonal phosphorus cycle and to some extent, the succession of the plankton groups in this ecosystem.

Our scenario analysis suggests that the winter concentration of P. rubescens has a determining role for its dynamics in the following growth season.

The model is implemented as three FABM [START_REF] Bruggeman | A general framework for aquatic biogeochemical models[END_REF] modules, which can be coupled at run-time (i.e., without the need for changing and re-compiling the model code): a nutrient module, a detritus module, and a plankton module that can be placed anywhere along the autotrophyheterotrophy continuum [START_REF] Flynn | Misuse of the phytoplankton-zooplankton dichotomy: the need to assign organisms as mixotrophs within plankton functional types[END_REF] with the 'autotrophy fraction' parameter ζ, enabled by a generalized formal description of autotrophic and heterotrophic processes (Section 4.1). In the present implementation, the plankton module is coupled six times to describe interactions of three purely autotrophic, two purely heterotrophic and a mixotrophic plankton groups. Realistic representation of spatio-temporal distribution of inorganic and particulate phosphorus, as well as the plankton groups constitute a proof of the concept. These generic modules thus enable description of various food web congurations in run-time and therefore can facilitate site-specic congurations in future studies.

Incorporation of mixotrophy is another novel aspect of the current study, given that mixotrophy has been,

with exceptions (e.g., [START_REF] Ward | Marine mixotrophy increases trophic transfer eciency, mean organism size, and vertical carbon ux[END_REF], neglected in ecosystem scale model applications. Here, the description of mixotrophy is quite simplistic, relying on the assumption that there exists a linear trade-o between autotrophic and heterotrophic abilities [START_REF] Crane | Coexistence of mixotrophs, autotrophs, and heterotrophs in planktonic microbial communities[END_REF], and the assumption that the parameters of heterotrophic and autotrophic traits reect the conditions that enforced pure heterotrophy or autotrophy in the meta-analysis data that formed the basis for parameterization in this study. It should be further noted that, in reality, a major carbon and nutrient source for mixotrophs is often bacteria [START_REF] Mitra | The role of mixotrophic protists in the biological carbon pump[END_REF], which is not explicitly modelled in this study, but is implicitly represented by the POM pool, P content of which is known, but C content is assumed to be proportional to P by the Redeld ratio. Further resolution of the microbial loop and the organic matter pools in Lake Bourget would be out of the scope of the current study, but constitutes a potential future goal, as with ongoing oligotrophication, relative importance of mixotrophs in the system can be anticipated to increase [START_REF] Kamjunke | Phosphorus gain by bacterivory promotes the mixotrophic agellate Dinobryon spp. during re-oligotrophication[END_REF][START_REF] Mitra | The role of mixotrophic protists in the biological carbon pump[END_REF].

In this study, mixotrophs and herbivorous and carnivorous zooplankton were assigned as individual functional groups based on their position in the food web. P. rubescens was assigned as a separate algal group because of its distinctive eco-physiology characterized by low-light tolerance [START_REF] Walsby | Light-dependent growth rate determines changes in the population of Planktothrix rubescens over the annual cycle in Lake Zürich, Switzerland[END_REF], slow growth [START_REF] Bright | The daily integral of growth by Planktothrix rubescens calculated from growth rate in culture and irradiance in Lake Zürich[END_REF], and grazing defense [START_REF] Kurmayer | Strategies for the co-existence of zooplankton with the toxic cyanobacterium Planktothrix rubescens in Lake Zürich[END_REF] as well as the fact that it represented about half of the phytoplankton assemblage in some years (Fig. 1,3). Previous modelling studies focusing on similarly P. rubescens-infested lakes unequivocally assigned P. rubescens an own group [START_REF] Omlin | Biogeochemical model of Lake Zürich: model equations and results[END_REF][START_REF] Copetti | A biogeochemical model of Lake Pusiano (North Italy) and its use in the predictability of phytoplankton blooms: rst preliminary results[END_REF][START_REF] Mieleitner | Modelling functional groups of phytoplankton in three lakes of dierent trophic state[END_REF][START_REF] Carraro | Coupling high-resolution measurements to a three-dimensional lake model to assess the spatial and temporal dynamics of the cyanobacterium Planktothrix rubescens in a medium-sized lake[END_REF]. Remaining phytoplankton were separated by size, forming 'small' and 'large' algal groups, given the importance of size in determining light and light utilization traits, among others [START_REF] Finkel | Phytoplankton in a changing world: cell size and elemental stoichiometry[END_REF]. There exist modelling studies where taxonomic classes were parameterized directly as functional groups (e.g., [START_REF] Jöhnk | Summer heatwaves promote blooms of harmful cyanobacteria[END_REF][START_REF] Carraro | Coupling high-resolution measurements to a three-dimensional lake model to assess the spatial and temporal dynamics of the cyanobacterium Planktothrix rubescens in a medium-sized lake[END_REF], but categorization based on size, thereby combining dierent taxonomic classes within a given functional group like in this study is not uncommon (e.g., [START_REF] Mieleitner | Modelling functional groups of phytoplankton in three lakes of dierent trophic state[END_REF][START_REF] Gal | Implementation of ecological modeling as an eective management and investigation tool: Lake Kinneret as a case study[END_REF][START_REF] Rinke | A simulation study of the feedback of phytoplankton on thermal structure via light extinction[END_REF].

Parameterization of the plankton groups was largely based on allometric relationships and taxonomic properties, inspired by recent meta-analyses [START_REF] Litchman | Contrasting size evolution in marine and freshwater diatoms[END_REF][START_REF] Schwaderer | Eco-evolutionary dierences in light utilization traits and distributions of freshwater phytoplankton[END_REF][START_REF] Edwards | Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton[END_REF]. This approach was to some extent necessary, due to the lack of information on physiological parameters of the specic plankton species or groups observed in Lake Bourget. On the other hand, we believe that this approach should produce a more consistent and transferable (i.e., applicable to other systems) parameter set than that could be obtained by a pure calibration exercise, which is problematic due to poor parameter identiability of plankton functional group models [START_REF] Anderson | Plankton functional type modelling: running before we can walk[END_REF][START_REF] Mieleitner | Modelling functional groups of phytoplankton in three lakes of dierent trophic state[END_REF]. In this study, parameter tuning was restricted to parameters like zooplankton preferences, sinking rates (of AL and detritus) and remineralization rates. Setting zooplankton feeding preferences is a common problem in ecosystem models, as model results are sensitive to dierent parameterizations [START_REF] Gentleman | Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF], and as the experimental data so far are based on specic prey-predator pairs that are hardly generalizable. For the parameterization of sinking and remineralization rates, we acknowledge potential biases in the form of an underestimation of sinking rates and an overestimation of remineralization rates to compensate the lack of river nutrient uxes.

Based on the classical variable internal stores model schemes, the storage capacity, Qmax/Qmin has been shown to have a critical importance in the presence of temporal [START_REF] Grover | Resource competition in a variable environment: phytoplankton growing according to the variable-internal-stores model[END_REF] and/or spatial [START_REF] Kerimoglu | Role of phytoplankton cell size on the competition for nutrients and light in incompletely mixed systems[END_REF] heterogeneity. Dierent allometric scalings observed in Qmax and Qmin for nitrogen, contrasting with similar coecients in those for phosphorus have been suggested to explain marine diatoms being larger than freshwater diatoms [START_REF] Litchman | Contrasting size evolution in marine and freshwater diatoms[END_REF]. On the other hand, it is questionable whether Qmax should be considered as a xed physiological parameter, as the observed maximum nutrient contents in phytoplankton cells might as well be reecting the acclimative down-regulation of the uptake rate, which depends on the specic settings, like limitation by other nutrients, supply rates and light limitation, as shown by optimality based model schemes [START_REF] Pahlow | Optimality-based model of phytoplankton growth and diazotrophy[END_REF][START_REF] Wirtz | Autotrophic Stoichiometry Emerging from Optimality and Variable Colimitation[END_REF]. This might be the reason for lack of robust patterns in distributions of Qmax across cell size or taxa, hence, its exclusion from analyses aiming to identify allometric and taxonomic variations in physiological traits (e.g., [START_REF] Litchman | The role of functional traits and trade-os in structuring phytoplankton communities: scaling from cellular to ecosystem level[END_REF][START_REF] Edwards | Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton[END_REF]. In this study, acknowledging the poor understanding of the issue, we assumed a xed proportionality between Qmax and Qmin, in order to avoid introducing a potentially strong but unjustied trait that can eect the competitive outcomes.

Overall, the model performance is good, given the limitations of the idealized 1-D approach. Model estimates for water temperature are almost excellent (Fig. 7-10). Average surface concentrations of DIP and TP can considered to be very good (r>0.85, Fig. 8). The measured DIP concentrations during summer tend to be somewhat higher than the simulations, which might be, to some extent, due to the river uxes that not represented by the model, however, considering the long water residence time of the lake (about 10 years) there is probably some other mechanism involved. Vertical proles of DIP, and to a lesser extent, TP within the upper 80 m are quite realistically reproduced by the model. Annual average concentration of plankton groups are well reproduced, indicated by low bias, but the seasonal dynamics is poorly represented, resulting in phase errors, and low correlation coecients (Fig. 6). Interannual variability of the data is also more pronounced than that of the simulations. Especially P. rubescens follows contrasting patterns in dierent years: in the form of smooth patterns in 2004 and 2009, and very rapid changes in the years 2005 and 2006, which are likely to be related with some vertical or lateral transport mechanism like those induced by internal waves or upwelling events that occur frequently in Lake Bourget [START_REF] Cuypers | Impact of internal waves on the spatial distribution of Planktothrix rubescens (cyanobacteria) in an alpine lake[END_REF], given their low growth rates [START_REF] Bright | The daily integral of growth by Planktothrix rubescens calculated from growth rate in culture and irradiance in Lake Zürich[END_REF] and limited losses imposed by zooplankton [START_REF] Kurmayer | Strategies for the co-existence of zooplankton with the toxic cyanobacterium Planktothrix rubescens in Lake Zürich[END_REF]. On the other hand, their growth at 10-20 m depth (Fig. 5) notably coincides well with the observations in Lake Bourget [START_REF] Jacquet | The proliferation of the toxic cyanobacterium Planktothrix rubescens following restoration of the largest natural French lake (Lac du Bourget)[END_REF](Jacquet et al., , 2014b;;[START_REF] Vu | Highfrequency monitoring of phytoplankton dynamics within the European water framework directive: application to metalimnetic cyanobacteria[END_REF]. In reality, this phenomena is usually attributed to buoyancy regulation [START_REF] Walsby | Comparison of the depth where Planktothrix rubescens straties and the depth where the daily insolation supports its neutral buoyancy[END_REF], whereas in this study, no active mobility mechanism has been considered and their sub-surface growth largely reects the parameterized low-light tolerance and high-light inhibition of the species [START_REF] Walsby | Light-dependent growth rate determines changes in the population of Planktothrix rubescens over the annual cycle in Lake Zürich, Switzerland[END_REF]. The depth-dependent mortality we assumed for P. rubescens (eq.14) does not play a signicant role in the occurrence of such thin layers, whereas it certainly pulls the depth of maximum concentration a few meters upwards. Our preliminary work (not shown) suggest that the thin-layers can be represented by specifying an optimum growth temperature or irradiance, however, using a process-based buoyancy regulation model (e.g., [START_REF] Kromkamp | A computer model of buoyancy and vertical migration in cyanobacteria[END_REF] would be more useful for gaining a better understanding of the relevance of the formation of thin layers.

Finally, high frequency uctuations in other variables are not captured, which might be related with extreme river discharge events [START_REF] Vinçon-Leite | Contribution of mathematical modeling to lake ecosystem understanding : Lake Bourget (Savoy, France)[END_REF] in the system. A 3-D modelling study revealed that in a relatively smaller (surface area: 5 km 2 ) Lake Pusiano, a major ood event caused a signicant variation in the horizontal distribution of P. rubescens [START_REF] Carraro | Coupling high-resolution measurements to a three-dimensional lake model to assess the spatial and temporal dynamics of the cyanobacterium Planktothrix rubescens in a medium-sized lake[END_REF]. Given that the two major tributaries of Lake Bourget are in the South [START_REF] Bryhn | Predicting future eects from nutrient abatement and climate change on phosphorus concentrations in Lake Bourget, France[END_REF], and the outows are in the North, such can be expected in Lake Bourget as well to some extent, but the eects are presumably smaller due to the much smaller quantity of uxes in relation to the water volume in Lake Bourget in comparison to Lake Pusiano, which has a water residence time of about 1 year [START_REF] Carraro | Coupling high-resolution measurements to a three-dimensional lake model to assess the spatial and temporal dynamics of the cyanobacterium Planktothrix rubescens in a medium-sized lake[END_REF], i.e., one tenth that of Lake Bourget.

Occurrence of P. rubescens is usually attributed to the nutrient availability, when long-term dynamics [START_REF] Dokulil | Deep living Planktothrix rubescens modulated by environmental constraints and climate forcing[END_REF][START_REF] Posch | Harmful lamentous cyanobacteria favoured by reduced water turnover with lake warming[END_REF]) and lakes at dierent trophic states are considered [START_REF] Salmaso | Eects of nutrient availability and temperature on phytoplankton development: A case study from large lakes south of the Alps[END_REF][START_REF] Anneville | Blue-Green Algae in a Greenhouse Century? New Insights from Field Data on Climate Change Impacts on Cyanobacteria Abundance[END_REF]. Within our study period, winter DIP fell from 0.9 to 0.46 mmolP/m 3 (corresponding to 27.9 to 14.3 mgP/m 3 ) between 2004 and 2010 (Fig. 1), which does not encompass a large range relative to the above examples. Indeed, our scenario analysis where we used the winter DIP concentration of 2004 for simulating 2010 did not result in a signicant dierence from the reference simulation of 2010 (Fig. 11). In contrary, when the winter concentration of P. rubescens in 2004 was used for simulating 2010, a large P. rubescens bloom was formed, unlike the reference simulation of 2010 (Fig. 11). These results suggest that the initial abundances by the end of the winter are highly important for the dynamics of P. rubescens within the growth season, which is in line with the ndings of Jacquet et al. (2014b), who showed a signicant relationship between the winter-spring, spring-summer and summerautumn biomasses. Qualitatively the same results were obtained when the scenario analysis was repeated with an altered model parameterization (not shown), where, P. rubescens was described as a species with slightly stronger resource competition traits (50% lower KU and 25% higher µ∞, given the uncertainties of the allometric scaling of these parameters), and slighty weaker defense against predation (mixotrophs feeding on them with p M,k =0.5,0.4 and 0.1 for k=POM, AS and AP R, given the lacking evidence for the feeding interactions between P. rubescens and mixotrophic species in Lake Bourget). The sensitivity of growing season P. rubescens abundance to the initial population size after winter does not simply result from a delayed geometric growth: when the scenario analysis is repeated with a structurally dierent model where P. rubescens is the only autotroph, i.e., not having any competitors, even the scenario with low initial inoculum and phosphorus concentrations results in a substantial bloom with concentrations reaching up to 50 mmolC m -3 (Appendix B), which is close to the concentrations observed in the years with intense P. rubescens occurrences. These ndings therefore suggest that the disappearance of P. rubescens from Lake Bourget was actually due to competitive exclusion, facilitated by the phenomena called the priority eect [START_REF] Hodge | Eects of Temporal Priority on Interspecic Interactions and Community Development[END_REF], in which, the species with higher initial concentrations outcompete their competitors with lower initial concentrations by making the abiotic environment inhabitable by the others. The priority eect has been previously reported from eld measurements [START_REF] Tapolczai | Occurrence and mass development of Mougeotia spp. (Zygnemataceae) in large, deep lakes[END_REF] and mesocosm experiments [START_REF] Sommer | Climate change and the phytoplankton spring bloom: Warming and overwintering zooplankton have similar eects on phytoplankton[END_REF]. Vertical distributions obtained under dierent scenarios (not shown) suggest that the priority eect is manifested as a bistability of both phytoplankton composition and vertical distribution of phytoplankton biomass (as in [START_REF] Ryabov | Vertical distribution and composition of phytoplankton under the inuence of an upper mixed layer[END_REF]. When P. rubescens start the season from relatively higher concentrations, they end up being the dominant species starting from summer by forming a dense layer at the metalimnion and thereby prevent the growth of mainly the small species through intensifying the nutrient limitation in the upper layers. To the contrary, when they start from low concentrations, small phytoplankton persist in the surface layers throughout the summer, making it thereby dicult for P. rubescens to grow through intensifying light limitation. The self-stabilizing nature of these two alternative states may partially explain the inter-annual variability of the P. rubescens abundances commonly observed in various deep lakes [START_REF] Vinçon-Leite | Phytoplankton variability in Lake Bourget: Phytoplankton dynamics and meteorology[END_REF][START_REF] Salmaso | Long-term phytoplankton community changes in a deep subalpine lake: responses to nutrient availability and climatic uctuations[END_REF]Jacquet et al., 2014b;[START_REF] Anneville | Blue-Green Algae in a Greenhouse Century? New Insights from Field Data on Climate Change Impacts on Cyanobacteria Abundance[END_REF].

It should be noted that, although the additional P. rubescens growth caused by higher initial phosphorus was relatively small, the higher P. rubescens concentration by the beginning of winter might be of substantial importance for the dynamics that follow later in the next season, as scenarios with higher initial P. rubescens concentrations suggest (Fig. 11). Moreover, the range of phosphorus concentrations in this scenario analysis is relatively small relative to the historical observations in Lake Bourget itself and in other lakes, so it does not imply that phosphorus availability is unimportant for the occurrence of this P. rubescens. In contrary, we believe that the long-term decrease of phosphorus in Lake Bourget was a prerequisite for the nal disappearance of P. rubescens as suggested by (Jacquet et al., 2014b). A model-aided test of this hypothesis can be achieved with a model validated for contrasting conditions with respect to the trophic state and the availability of P. rubescens that are observed either in a single lake across a long-term transition, or multiple lakes of dierent trophic status.

In Lake Bourget [START_REF] Vinçon-Leite | Phytoplankton variability in Lake Bourget: Phytoplankton dynamics and meteorology[END_REF]Jacquet et al., 2014b;[START_REF] Anneville | Blue-Green Algae in a Greenhouse Century? New Insights from Field Data on Climate Change Impacts on Cyanobacteria Abundance[END_REF], as well as in other deep lakes like Lake Zurich [START_REF] Walsby | The critical pressures of gas vesicles in Planktothrix rubescens in relation to the depth of winter mixing in Lake Zurich, Switzerland[END_REF][START_REF] Posch | Harmful lamentous cyanobacteria favoured by reduced water turnover with lake warming[END_REF], Lake Garda [START_REF] Salmaso | Long-term phytoplankton community changes in a deep subalpine lake: responses to nutrient availability and climatic uctuations[END_REF] and Mondsee [START_REF] Dokulil | Deep living Planktothrix rubescens modulated by environmental constraints and climate forcing[END_REF], the abundance of P. rubescens by the end of winter is largely determined by meteorological conditions during winter. In such deep, warm monomictic lakes, cooling in mild winters may end up being insucient to cause a full overturn, leading to a reduced dilution of the P. rubescens concentrations, as well more P. rubescens cells survive the winter [START_REF] Walsby | The critical pressures of gas vesicles in Planktothrix rubescens in relation to the depth of winter mixing in Lake Zurich, Switzerland[END_REF]. Beyond the fact that phosphorus reached low concentrations in Lake Bourget, it seems like the coincidence of such a cold winter with a weak autumn population, might be the explanation for the complete elimination of the P. rubescens in 2010. Indeed, in 2010, both the autumn and winter air temperatures were among the coldest in the region [START_REF] Anneville | Blue-Green Algae in a Greenhouse Century? New Insights from Field Data on Climate Change Impacts on Cyanobacteria Abundance[END_REF], and P. rubescens had a relatively low autumn concentration (Jacquet et al., 2014b, and this study, Fig 6). In this context, it is worth noting that the warming trend observed in Lake Bourget [START_REF] Vinçon-Leite | Long-term temperature evolution in a deep sub-alpine lake, Lake Bourget, France: how a one-dimensional model improves its trend assessment[END_REF] can be expected to further continue given the anticipated warming of especially the winters in the alpine region [START_REF] Beniston | Mountain weather and climate: A general overview and a focus on climatic change in the Alps[END_REF], possibly promoting the re-establishment of P. rubescens in Lake Bourget, given the evidence that warming promoted P. rubescens in Lake Zurich [START_REF] Anneville | Seasonal and inter-annual scales of variability in phytoplankton assemblages: comparison of phytoplankton dynamics in three peri-alpine lakes over a period of 28 years[END_REF] and Lake Geneva [START_REF] Gallina | Estimating future cyanobacterial occurrence and importance in lakes: a case study with Planktothrix rubescens in Lake Geneva[END_REF].

In Lake Bourget, mixing does not reach the lake bottom during warm winters, such as in 2007 and 2008, which leads to signicantly lower phosphorus concentrations and high P. rubescens abundances later in the season season (Jacquet et al., 2014b). Therefore an extended model based analysis of plankton succession encompassing the winter dynamics would be highly relevant. Although the model presented in this study can serve as a starting point, several other processes need to be taken into account, including the lake hypsography, benthic-pelagic exchange, overwintering strategies of plankton, and the buoyancy regulation of P. rubescens, as well its inhibition through collapsing of gas vesicles under high pressure. Other rare behavior of P. rubescens worth considering in future modelling studies include colony formation behavior and its eect on predation by zooplankton. In this study, we assumed a non-zero but low preference of P. rubescens by the herbivorous zooplankton, but the evidence suggests that the grazing rate depends on the colony size [START_REF] Oberhaus | Grazing of two toxic Planktothrix species by Daphnia pulicaria: potential for bloom control and transfer of microcystins[END_REF], and there is a signicant relationship between the nutrient status and colony size (Jacquet et al., 2014b), suggesting that this preference should not be constant in time and space.

B

Sensitivity of scenario analysis to the model structure If the disappearance of P. rubescens in 2010 is really associated with the `priority eect' as claimed in the main text, in the absence of its competitors (other phytoplankton groups), P. Rubescens should reach to substantial abundances in 2010, even if starting from a small initial abundance after winter. In order to test this prediction, we performed the scenario analysis with an alternative model structure, in which P. Rubescens is the only autotroph available in the system, herbivorous zooplankton (ZH) feeds on P.

Rubescens and carnivorous zooplankton (ZC) feeds on ZH. Results of the scenario analysis obtained with this simple model structure is qualitatively dierent than those obtained with the reference model described in the text: initial inoculum of P. Rubescens indeed only determines the timing of the bloom, and not the presence/absence (Fig. B1). This analysis suggests that the disappearance of P. Rubescens from Lake Bourget is not a simple linear response to the inoculum size, but a case of competitive exclusion facilitated by the priority eect.

Figure 1 :

 1 Figure1: Observed concentrations of total phosphorus and total nitrogen (0-140 m average), and total phytoplankton and contribution by P. rubescens in Lake Bourget (0-20 m average) for the period 2004-2010.

KG

  [START_REF] Hansen | Zooplankton grazing and growth: Scaling within the 2-2,000-µm body size range[END_REF] did not nd systematic relationships between KG and body size or taxanomic group and report an average of 240 mgC m -3 = 20 mmolC m -3 which we take for all heterotrophs (M , ZH , ZC ).

  the water column increases the absorption of light. Dierent phytoplankton species and detritus have been observed to have dierent specic-light extinction coecients. Oubelkheir et al. (2005) suggests 0.02 m 2 mmolC

Figure 3 :

 3 Figure 3: Measured contribution of major taxonomic classes to each phytoplankton group throughout the simulated years.

Figure 4 :

 4 Figure 4: Limitation functions and their eects on autotrophic functions for each plankton group with the parameter values listed in Tables3-4. f T , f Q , f I and f N stand for temperature limitation (eq. 16), nutrient limitation of biomass growth (eq. 3), light limitation of photosynthesis (eq. 5) and nutrient limitation of uptake (eq. 8), respectively.

  Seasonal and spatial distributions of the state variables estimated by the model are in general reasonable, as exemplied for the year 2004 in Fig.5. As the surface waters warm and thermal stratication develops at around mid-March, the mixed layer becomes conned to the surface layers. This relieves phytoplankton from light limitation, leading to blooms of rapid growing AS and AL, and consumption of DIP within the mixed layer. Phytoplankton growth fuels mixotrophs and herbivorous zooplankton, leading to POP production and consequently further growth of mixotrophs. Increasing abundances of herbivorous zooplankton and mixotrophs result in a summer-bloom of carnivorous zooplankton, suppressing, in turn, herbivory and hence, preventing the collapse of AS. Autumn is dominated by P. rubescens, which has been growing since midsummer and at increasingly deeper layers, reaching to about 15-20 meters towards the end of autumn. The season ends with stratication gradually weakening and mixed-layer extending back to deeper layers.

Figure 5 :

 5 Figure 5: Spatio-temporal distribution of simulated major variables in 2004.

Figure 6 :

 6 Figure 6: 0-20m averaged concentration of simulated (lines) and measured (circles) concentration of plankton groups.

Figure 7 :

 7 Figure 7: 0-20m averaged concentration of simulated (lines) and measured (circles) concentration of biological variables.

Figure 8 :

 8 Figure 8: Target-Taylor diagrams for 0-20 m average concentration of (top) individual plankton groups and (bottom) dissolved phosphorus and total phosphorus, total phytoplankton, total mixotrophs and total zooplankton carbon biomass.

Figure 10 :

 10 Figure 10: Simulated (lines) and measured (circles) Temperature, DIP and TP proles from 2004 and 2010 within the upper 80 m. Only the rst measurements in each month are shown. Proles from identical months in dierent years are aligned.

Figure 11 :

 11 Figure 11: Surface (0-20 m) average water column average state variables with initial DIP concentration of DIP 0 = 0.46 mmolP m -3 as in 2010 (solid lines), DIP 0 = 0.91 mmolP m -3 as in 2004 (dashed lines), in combination with initial P. rubescens concentration of P R 0 = 0.01 mmolC m -3 as (was assumed) in 2010 marker) and P R 0 = 0.62 mmolC m -3 as in 2004 (triangles). Initial conditions for all other variables and meteorological forcing was as in 2010.

  DIP 0 =0.46, PR 0 =0.01 DIP 0 =0.91, PR 0 =0.01 DIP 0 =0.46, PR 0 =0.62 DIP 0 =0.91, PR 0 =0.62

Figure B1 :

 B1 FigureB1: Surface (0-20 m) average water column average state variables with initial DIP concentration of DIP 0 = 0.46 mmolP m 3 as in 2010 (solid lines), DIP 0 = 0.91 mmolP m 3 as in 2004 (dashed lines), in combination with initial P. rubescens concentration of P R 0 = 0.01 mmolC m 3 as (was assumed) in 2010 (no marker) and P R 0 = 0.62 mmolC m 3 as in 2004 (triangles). Initial conditions for all other variables and meteorological forcing was as in 2010.

Table 1

 1

, 2004-2010 average biovolume fraction and biovolume-weighted cell volume of each algal and mixotroph group are provided.

Table 1

 1

	: Properties of algal/mixotroph groups, averaged for the period 2004-2010. % Contr.: percentage
	contribution to the total biovolume of phytoplankton and mixotrophs; V: biovolume-weighted average volume
	of the group
	Class	% Contr. V [µm 3 ]
	Small	22.9	236.9
	Large	15.5	13717.1
	P. rubescens	40.1	84.8
	Mixotrophs	21.3	12978.8

Table 2 :

 2 Denition and units of state variables, processes and intermediate quantities

	Symbol	Unit	Denition
	X C	mmolC m -3	C bound to plankton
	X P	mmolP m -3	P bound to plankton
	P OP	mmolP m -3	Particulate organic phosphorus
	DIP	mmolP m -3	Dissolved inorganic phosphorus
	Q	molP molC -1	Phosphorus quota
	µ C(P )	d -1	C (P) limited growth rate
	P	d -1	Photosynthesis rate
	G jk	d -1	Ingestion rate of prey k by predator j
	U	molP molC -1 d -1	Nutrient uptake rate
	L	d -1	Mortality rate
	the `mixotroph species' in Crane and Grover (2010). A central parameter in this unied representation is
	the fraction of photosynthethic autotrophy, ζ, set to ζ = 1.0 for 3 autotroph groups (Xi={APR, AL, AS}),
	ζ = 0.5 for the mixotrophs (Xi={M}) and ζ = 0.0 for the zooplankton (Xi={ZH, ZC}) (Table 3, see section
	3 for the identication of functional groups).

Table 3 :

 3 Model parameters specic to plankton groups. Horizontal lines separate the parameters relevant for all plankton, only autotrophs and only heterotrophs.

	Symbol	Unit		Denition	AS	AL	APR	M	ZH	ZC
	ζ	-				fraction of autotrophy	1	1	1	0.5	0	0
	w	m d	-1		sinking velocity	0	-0.1	0	0	0	0
	l q	d	-1 m 3 mmolC -3	quadratic mort. rate	0.001	0.001	0.001	0.01	0.02	0.04
	l	d	-1			linear mortality rate	0.02	0.02	0.02	0.02	0.02	0.02
	e	d	-1			excretion rate	0	0	0	0.01	0.02	0.02
	Qmax	molP molC	-1	upper bound of P quota	0.0145	0.04	0.0145	0.0082	0.0117	0.005
	Qmin	molP molC	-1	subsistence quota	0.00145	0.004	0.00145	0.0045	0.0117	0.005
	µ∞	d	-1

Table 4 :

 4 Model parameters that are not specic to plankton groups.

	Symbol Unit	Denition	value
	r	d -1	remineralization rate	0.15
	k c,phy	m 2 mmolC -1	specic light ext. coef. 0.03
	k c,det	m 2 mmolC -1	specic light ext. coef. 0.02
	w P OP	m d -1	sinking rate	-1.0
	Q P OM	molP molC -1 P:C ratio of POM	1/116
	T ref	°C	reference temperature	20
	Q10 B	-	Q10 for bacteria	1.5
	Q10 A	-	Q10 for autotrophs	1.5
	Q10 H	-	Q10 for heterotrophs	2

Table 5 :

 5 Model parameters: grazing preferences p j,k of predator j (rows) for prey k (columns). -stands for p j,k = 0. for Z C indicates that preferences are dynamically adjusted (see eq. 7). POM

Table 6 :

 6 Allometric scaling coecients for U max , Q min , G max (of form 10 a V b ) and Q C (of form aV b ). U max and Q min from :Edwards et al. (2012) þ:Litchman et al. (2009); Q C from Menden-Deuer and Lessard

Acknowledgements

This study is a contribution to the PhyCl project and to the observatory of alpine lakes (OLA). Data were made available by SOERE OLA-IS, INRA Thonon les Bains, CISALB and ecoinformatics ORE-INRA team.

Authors are grateful to all the technicians, engineers and researchers who have contributed to the monitoring of Lake Bourget. We thank two anonymous referees for their comments. Funding for OK was provided by the University of Savoie Mont Blanc (AAP 2012). OK was additionally supported by the German Research Foundation (DFG) through the priority program 1704 `DynaTrait'.

Appendices

A

Composition of functional groups