Model-based data analysis of the effect of winter mixing on primary production in a lake under reoligotrophication

Nutrient loading, in combination with climate change are important drivers of primary productivity in lakes. Understanding and forecasting future changes in primary production (PP) in response to local and global forcing are major challenges for developing sustainable lake management. The objective of this study is to understand and characterize the mechanisms underlying the large differences in observed PP rates and nutrient concentrations between two consecutive years (2012 and 2013) in Lake Geneva, Switzerland. For this purpose, we apply a one-dimensional (1D) physical–biogeochemical model system. The Framework of Aquatic Biogeochemical models (FABM) interface is used to couple the General Ocean Turbulence Model (GOTM) with a biogeochemical model, the Ecological Regional Ocean Model (ERGOM). We calibrated GOTM, by adjusting physical parameters, with the observed temperature profiles. A model calibration method is implemented to minimize model-data misfits and to optimize the biological parameters related to phytoplankton growth dynamics. According to our results, the simulated surface mixed layer depth is deeper and heat loss from the lake and turbulent kinetic energy in the water column are much higher in winter 2012 than that in 2013, pointing to a cooling-driven, deep mixing in the lake in 2012. We found significant differences in internal phosphorus loads in the epilimnion between the two years, with estimates for 2012 being higher than those for 2013. ERGOM predicts weak nutrient limitation on phytoplankton and higher growth rates in 2012. Apparently, the deep mixing event led to high turnover of nutrients (particularly dissolved inorganic phosphate) to the productive surface layers, and a massive algal bloom developed later in the productive season. In contrary, the turnover of nutrients in 2013 was weak and consequently the PP was low. Our findings demonstrate the utility of a coupled physical–biological model framework for the investigation of the meteorological and physical controls of PP dynamics in aquatic systems.


Introduction
Changes in nutrient loads and climate forcing perturb lake ecosystems on different temporal scales. On the one hand, studies have linked long-term responses, e.g. increase in thermal stratification and shifts in phytoplankton and zooplankton community structures, to global warming and different trophic conditions (Magnuson et al., 1990;Watson et al., 1997;Livingstone, 2003;Austin and Colman, 2007;Perroud et al., 2009;Anneville et al., 2019;Lepori, 2019). On the other hand, seasonal and inter-annual variations in dynamics of physical and Vertical mixing in lakes is largely controlled by heat fluxes and wind (Imboden and Wüest, 1995). Typically for temperate lakes, deep mixing occurs in spring and fall, and in summer they remain stratified. The turnover of water masses has important consequences for the development of the trophic state of a lake as it influences the exchange of heat and matter between epilimnion, metalimnion, and hypolimnion (Michalski and Lemmin, 1995). Vertical mixing is one of the important factors controlling PP in a lake as the supply and concentration of nutrients from the deeper parts to the epilimnion depends on mixing strength. Schwefel et al. (2019) investigated the total phosphorus (TP) concentrations in the epilimnion of a meromictic, alpine, lake under present and future climate conditions. Their data show that TP transport from the deep hypolimnion by winter mixing can be the predominant source of TP in the lake epilimnion, much higher than the contributions of nutrient inputs from catchments.
The mixing depth has stochastic control on PP because it changes with meteorological forcing (Lewis, 2011). These changes can occur on yearly scale or can be long-term. Straile et al. (2010) studied the impact of the warm winter of 2006/2007 (W06/07) on nutrient mixing dynamics and phytoplankton biomass in Lake Constance, and compared it with the preceding cold winter of 2005/2005 (W05/06). Their analysis shows that the density gradient between upper (warm) and lower (cold) water layers was strong in W06/07, which prevented deep mixing, and consequently the vertical transport of phosphorus to the upper water layers. In years 2005 and 2006, exceptional deep mixing events stirred up and ventilated the entire water column in Lake Lugano, leading to high PP (Holzner et al., 2009). Goldman et al. (1989) studied yearly fluctuations in PP in the two temperate Lakes Castle and Tahoe. They found that PP patterns in Lake Castle are driven by large-scale climate events, e.g. El Niño or Southern Oscillation, and in Lake Tahoe by strong mixing in winter and spring. As climate warming continues, a reduction in the potential of vertical mixing has been reported for tropical Lake Tanganyika, and thereby limiting the nutrient fluxes to the photic layers (Verburga and Hecky, 2009). Lake Geneva is a peri-alpine lake which renders a variety of ecosystem services (leisure, fishing, transportation etc.) to the region. The lake has experienced severe eutrophication in the last century and slowly, but steadily, recovered from it. Significant efforts (e.g. strict regulation) were undertaken across Switzerland and France, since the 1970s, to improve the trophic status of lakes and water quality (Jacquet et al., 2014). Although phosphorus concentrations in Lake Geneva have significantly decreased (annual averages in the top 30 m in 1980 were 45 μg L −1 , compared to 5 μg L −1 in 2015) due to successful point sources management (Anneville et al., 2019), observed PP and total phytoplankton biomass remain high and strongly variable on interannual and seasonal scales (Tadonléké et al., 2009). These variations are underpinned by dynamic interactions between physical (e.g. winter mixing and winds) and ecological processes (e.g., grazing pressure, nutrient turnover, remineralization etc.).
In this study, we apply a one-dimensional (1D) coupled model system to investigate variations in PP and nutrient dynamics in Lake Geneva between the two consecutive years 2012 and 2013. Interestingly, the maximum value of the observed PP in 2012 was three times higher than that of 2013, whereas dissolved inorganic phosphorus (DIP) concentrations were higher and dissolved inorganic nitrogen (DIN) concentrations were lower in the lake epilimnion in 2012 compared to 2013. We seek to gain insights on the drivers (e.g. winter mixing intensity and riverine inputs) that could generate such changes in PP and nutrient dynamics. To address this, we designed a calibration method that provides the optimized estimates of eco-physiological parameters for 2012, and then we validate the model against the observations from 2013. Together with estimates of physical quantities and biological parameters, we aim to understand how changes in deep mixing dynamics affect PP rates and nutrient concentrations.

Study site and observations
Lake Geneva, the largest lake in Western Europe, is shared by Switzerland (60%) and France (40%). It is situated at an altitude of 372 m on the north side of the Alps and covers a total area of 580 km 2 . The residence time of the lake water is ∼11 years. Lake Geneva can be divided into three parts: the eastern part close to the Rhône estuary, the largest and deepest central basin (Grand Lac), and the south-western, shallower and narrow part (Petit Lac).
Lake Geneva enjoys the legacy of a long sampling history and has been monitored since early 1950s. Initially by Union Générale des Rhodaniens, and thereafter (1963 onwards) by Commission internationale pour la protection des eaux du Léman (CIPEL). Water samples to study physical and biogeochemical properties of the lake, e.g., temperature, conductivity, transparency, dissolved oxygen and nutrient concentrations, are collected monthly in winter and bimonthly in other seasons from stations GE3 (Petit Lac) and SHL2 (Grand Lac at the deepest location, Perroud et al. 2009). The sampling is done at 19 discrete depth levels: 0, 2.5, 5, 7.5, 10, 15, 20, 25, 30, 35, 50, 100, 150, 250, 275, 290, 300, 305, and 309 m. The database of observations is maintained by the French National Research Institute for Agricultural, Food and Environment (INRAE).
In this study, we apply a 1D coupled physical-biogeochemical model system to analyse and interpret in-situ temperature, PP and nutrient data from SHL2 station for two years (2012 and 2013). PP measurements were done using the 14 C-radiotracer method.

Yearly phosphorus budget
In order to investigate the potential impact a deep mixing event may have on PP rates, a full budget of bioavailable phosphorus was conducted in the production zone for the two years considered in our study. The net supports the phytoplankton biomass production in the trophogenic layer of a lake. It corresponds to the sum of all watershed contributions entering and leaving the system  −  , with the net difference of P stock in the epilimnion before and after the production period. This stock difference corresponds to the total mass of TP in the epilimnion after the winter mixing minus the remaining mass in autumn at the end of the PP period.
The productive zone was defined with = −30 m based on the analysis of long-term phytoplankton biomass concentration and PP profiles (not shown here, Tadonléké et al. 2009). The area as a function of depth, ( ), was taken from hypsometric curves (Federal Office of Topography, 2017).
Public data of in-stream P and discharge measurements in lake tributaries were accessible for the Rhône River only, which delivers ≈ 85% of the total river inflow (https://www.cipel.org). Total dissolved P load originating from the Rhône watershed was obtained from flowweighted concentrations extracted from the weekly sampling at station ''Porte-du-Scex'' (upstream at the river-lake confluence, https://www. bafu.admin.ch) combined with the daily discharge data measured at the same site (https://www.hydrodaten.admin.ch). Fluxes of total P discharged from waste water treatment plants (WWTPs) were collected annually by CIPEL and publicly available in the annual reports (https: //www.cipel.org). The loss of P at the outlet of the lake  was calculated from surface TP concentration measured near Geneva (station GE3) averaged over the top 10 m, multiplied by daily flows measured in Geneva at the gauge station ''Les Halles'' (same data source as above).

Meteorological data
COSMO (Consortium for Small-scale Modelling) models are developed to forecast weather in the complex alpine region which includes Switzerland (Baldauf et al., 2011). We obtained hourly meteorological COSMO model data (e.g., short wave radiation, wind, air pressure, cloud cover and relative humidity), for both years, from MeteoSwiss. The spatial resolution of the model is 1.1 km.

Modelling framework
In this study, we apply a coupled physical-biogeochemical model system to investigate changes in dynamics of lake processes, such as winter mixing, evolution of temperature profiles, nutrient, PP. We use the FABM interface (Bruggeman and Bolding, 2014) to couple GOTM (Burchard et al., 2006) with ERGOM (Neumann et al., 2002).

Physical model
GOTM is a 1D hydrostatic model for studying thermo-hydrodynamic processes of marine and lake ecosystems. It contains state-of-the-art parameterizations for describing the vertical turbulent fluxes of heat and salinity. The model is based on Reynolds-averaged, Navier-Stokes, equations of momentum, temperature and salinity (Burchard et al., 2006). The temperature equation is given as: where the coordinate is taken to point upwards with the origin z = 0 at the lake surface. T is the temperature in units • C . The specific heat capacity ( ) and the molecular diffusivity of heat ( ′ ) are given in units J kg −1 • C −1 and m 2 s −1 .
is the kinematic radiative flux at a given depth and time and has units • C m s −1 . It is to note that ( , ) in GOTM is estimated from the surface irradiance ( 0 ), attenuation lengths at the surface ( 1 ) and over 10-40 m depth ( 2 ), and depth (Paulson and Simpson, 1977), where R is a weighting parameter, and the depth at which becomes 10% of 0 is assumed to be the Secchi depth ( ) in the model. Hence, Eq. (3) becomes We assume 1 = 0.7 m, as this value is representative for Lake Geneva and corresponds to the mean value of 1 for Jerlov water types IA and IB (Jerlov, 2014). We substitute the value for 1 and measurements of , obtained from SHL2 station, in Eq. (4) to get the estimates of 2 . The momentum and salinity equations are formulated equivalently, and to avoid repetitions they are presented in Appendices A.1 and A.2.
In this study we have applied GOTM-5.3 version, available at https: //github.com/gotm-model/code. A second order turbulent model has been considered for computing vertical turbulence fluxes (Umlauf and Burchard, 2005), and longwave back radiation is estimated by the method described in Clark et al. (1974). In addition, the model grid resolution at the surface is higher (grid zooming parameter, ddu = 3) than the bottom (ddl=0). For calibration, we adapt values for GOTM parameters from their typical ranges listed in Wüest and Lorke (2003) for lakes. In the next step, we adjust scaling parameters for winds and short wave radiations to obtain good agreement with observed temperature profiles. The detailed descriptions of GOTM are given in Burchard et al. (2006),  and .

Characterizing lake thermal stratification
Using the COSMO data, in-situ data from Lake Geneva, and the GOTM model output, we estimated the net heat flux for the years under study: ( , ) = lat ( ) + sen ( ) + + lw ( ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ Here, lat is the latent heat flux, sen is the sensible heat flux, − lw and + lw are the upward and downward long-wave radiation heat fluxes, while sw is the short-wave radiation heat flux, which penetrates across the water column. In particular, we denote the surface net heat flux as ( , = 0) = 0 ( ). GOTM provides the evolution of the thermal structure in Lake Geneva. Based on the temperature field, , and the salinity reference observed in Lake Geneva, 0 ≈ 0.2 g L −1 , we determined the density field using the equation of state ( , ) in Millero and Poisson (1981). The density distribution allowed estimation of the water stability, The position of the pycnocline, pyc , corresponds to the depth of the maximum stability, 2 max . For simplicity, we defined the surface mixedlayer (SML) as pyc . To avoid capturing the diurnal surface mixedlayer (Imberger, 1985), we look for 2 max in 5 m ≤ depth ≤ 309 m. The turbulent kinetic energy, , is a physical quantity modelled by GOTM. The root-mean-square of the depth-averaged in the SML, √ ⟨ ⟩, provides a level of velocity fluctuations in the SML and a quantity to assess the energetic state of the water column. To quantify the kinetic energy generated by cooling-driven convection in the SML, we estimated the convective velocity scale, , known as well as the Deardorff scale (Deardorff, 1970), , the thermal expansion coefficient, the gravitational acceleration, and 0 the water reference density. Having √ ⟨ ⟩ ∼ implies that convection is the dominant mechanism generating kinetic energy in the SML.

Biogeochemical model
We apply the ecological model to simulate nutrient and phytoplankton growth dynamics in Lake Geneva. Although the model was originally developed for studying biogeochemistry in marine ecosystems, its complexity in terms of trophic structure and the description of physio-ecological processes make ERGOM applicable for lakes as well. The model resolves growth dynamics of three functional groups of phytoplankton: diatoms ( 1), flagellates ( 2) and cyanobacteria ( 3), that account for nutrient and light limitations. The nitrogen limitation ( ) for the uptake of is described by a squared form of the Michaelis-Menten function, as it stabilizes the model solutions through reducing phytoplankton uptake at low nutrient concentrations (Neumann et al., 2002).
where and are concentrations of nitrate and ammonium in units of mmol N m −3 and is the half-saturation constant for nitrogen uptake in mmol N m −3 . Likewise, the phosphorus limitation function ( ) is given as: where is given in units of mmol P m −3 and is the P:N ratio. The model has adapted the light limitation ( ) parameterization from Steele (1962)  which applies to all three phytoplankton groups. Here, is the photosynthetically active radiation in units W m −2 and it is estimated from the surface irradiance ( 0 ), , and the attenuation coefficients of water and planktonic particles ( and ). The light optimum ( ) is the maximum value between 50% of 0 and the minimum irradiance ( ). In contrast to phytoplankton, ERGOM includes one bulk formulation for zooplankton that represents all species of heterotrophs. All dead particulate organic matter are accounted in bulk detritus. The model also resolves consumption and production of oxygen in the water column and sediment dynamics explicitly. The detailed first description of the model is given in Neumann (2000) and Neumann et al. (2002). Below, we briefly summarize the source parameterizations for phytoplankton and zooplankton.
Phytoplankton. The total phytoplankton biomass is given by the sum of the concentrations of diatoms, flagellates and cyanobacteria. The growth rate of diatoms ( 1 ) depends on , and and is given as: where 1 is the maximum potential growth rate of diatoms in units d −1 . The model considers smooth increase in the growth rate of flagellates ( 2 ) with temperature, in addition to dependencies on , and : where 2 is the maximum potential growth rate of flagellates in units d −1 .
and are the half-saturation temperature for flagellates and water temperature in • C. For cyanobacteria, ERGOM assumes an exponential increase in the growth rate ( 3 ) with temperature and no nitrogen limitation: where 3 is the maximum potential growth rate of cyanobacteria in d −1 . and are the half-saturation temperature for cyanobacteria and a temperature control parameter in units • C and [ • C] −1 , respectively.
Zooplankton. The grazing by bulk zooplankton, , on phytoplankton is described by the Ivlev formulation (Ivlev, 1966) and assumed to be dependent on temperature with lower preference for cyanobacteria: where has units d −1 . and are the optimal temperature for zooplankton and the total biomass of phytoplankton in units • C and mmol C m −3 , and is the Ivlev constant. The maximum grazing rates for phytoplankton are denoted as , in units d −1 , with the assumption that the rates on diatoms and flagellates ( 1 , 2 ) are higher than that on cyanobacteria: 1 , 2 > 3 .

Detritus.
The dead phytoplankton and zooplankton are accounted for in bulk detritus. Besides this, the resuspended matter from sediments also becomes a part of detritus. Eq. (B.1) in Appendix B shows the source terms for detritus.
Nutrients. Dynamics of nutrients ( , ) are represented by losses of nitrogen and phosphorus by physiological processes (e.g. respiration and excretion) and biogeochemical processes (e.g. remineralization of detritus and sediment). Eqs. (B.2) and (B.3) in Appendix show the rates of ammonification and nitrification. Whereas Eq. (B.5) denotes the source-minus-sink equation of DIP. The ordinary differential equations (ODEs) for nutrients are given in Neumann et al. (2002).
Sediment . The rate of sedimentation ( ) in ERGOM depends on the sedimentation rate of detritus ( ) and the critical stress ( ) at the bottom of the lake, see Eq. (B.6) in Appendix B.
Oxygen. The sources for oxygen ( ) in the model are PP by phytoplankton and the gas exchange through lake surface ( 2 ). Eq. (B.7) in Appendix B shows the ordinary differential equation for dissolved oxygen. Fig. C.8 in Appendix C shows the schematics of ERGOM, and its source code is available at https://github.com/fabm-model/fabm/ blob/master/src/models/gotm/ergom.F90.

Model calibration
From the results of our preceding test simulations (not shown here), we have identified that to realistically simulate the nutrient dynamics in the Lake Geneva, it is critical that the model accounts for temporal variations in the riverine input loads. This is particularly relevant for 2012 and 2013 because the observed monthly nutrient fluxes (especially DIN) from the rivers and tributaries are significantly different between the two years ( Fig. 1).
Our data assimilation (DA) approach is to calibrate the model for 2012, and validate the optimized solution with the measurements of PP, DIN, and DIP from 2013, while considering the observed monthly variations in the areal riverine influxes of DIN and DIP (in units mmol N m −2 d −1 and mmol P m −2 d −1 ) in the simulations. The river fluxes are one of the source terms in the simulation of the state variables DIP and DIN (see Eq. (B.5)). To quantify the discrepancy between measurements and model outputs, we designed a misfit cost function.

Cost function
We implement a probabilistic-based, multi-variate, misfit function ( ) which is similar to that of Schartau and Oschlies (2003). On a given date, the mean squared error (MSE) is given as the square of the difference between the depth-average value of an observed profile and the corresponding model output for a measurement type. We consider three Chi-squared cost functions (Ward et al., 2013), defined as sums of weighted MSEs, to quantify the total temporal data-model misfits in the lake epilimnion ( ), metalimnion ( ) and hypolimnion ( ℎ ). It is important to account for misfits in these layers of the water column individually because dynamics of physical and biogeochemical processes may vary between them, particularly during the stratified period. For any observation type, the misfit costs in the three layers are given as: where ( ) , ( ) and ( ℎ ) are the mean observations of epilimnion (0 to 10 m), metalimnion (10 to 50 m) and hypolimnion (> 50 m) on a date , and ( ) , ( ) and ( ℎ ) are the equivalent model outputs. Likewise, , and ℎ are the observed uncertainties (standard deviations) in epilimnion, metalimnion and hypolimnion on a given date. Note, , and ℎ describe the conditional probability of explaining the observed data, ℎ , given model and its parameters, and can be derived from their Gaussian probability density functions using Bayes' theorem (Sivia and Skilling, 2006). Three ( = 1..., 3) observation types ( , , ) are assimilated in the cost function and are assumed to be independent. The overall misfit cost ( ) is given as the sum of , and ℎ normalized to the number of observed profiles ( ), over the entire period, for all three variables ( ):

Parameter selection
The selection of parameters for optimization is based on their relevance to PP and the degree of sensitivity of the results to the parameter values. An ideal approach would be to assign the maximum growth rates parameters ( 1 , 2 , 3 ) the values identified in Neumann (2000) for Baltic Sea, and run the simulations for Lake Geneva. We unsuccessfully tried this approach in the preceding test (hand-tuned) simulations. Therefore, we have decided to optimize 1 , 2 , 3 for the respective year, and on purpose grant an extra degree of freedom to the model. Through a sensitivity analysis, we found , (mortality rate of zooplankton), (remineralization rate of detritus), and (sinking speed of detritus) parameters to be critical in resolving nutrient dynamics. Hence, we selected the above six physio-ecological parameters, also listed in Table 1. Other model parameters, that are not subjected to optimization, are assigned values that are identified in previous studies, Neumann (2000), Neumann et al. (2002), Fennel and Neumann (2014), and Kerimoglu et al. (2017). These parameters are listed in Table 2. Typically in lakes P:N ratios of seston are highly variable, and are lower than the Redfield ratio for marine waters as phosphorus limitation of algal growth is stronger than nitrogen limitation (Healey and Hendzel, 1980). Lakes with residence time greater than six months, generally, have P:N < 0.03 (Wetzel, 2001). Keeping all this in mind and after assessing performance of our hand-tuned solutions, we decided to fix P:N of seston ( parameter) in the model to a value of 0.02.

Optimization procedure
In this study, we adopted the optimization procedure of Krishna et al. (2019). We use functions from the R package FME (Soetaert et al., 2010) to implement, a two-step, optimization method that converges to its global minimum, and yields us the corresponding parameter estimates ( ). The first step is to hand-tune the control parameters to obtain a model solution (reference solution) that gives a reasonably good agreement with the observed profiles of , and . The reference solution is provided as initial parameter values for the start of optimizations. In the next step, we apply a simulated annealing algorithm (SANN), to perform stochastic search (not using gradient information) in the parameter space to find a solution in the vicinity of (Bélisle, 1992). This is done to avoid the algorithm getting trapped in the local minima of . Finally, we apply the gradient-based algorithm, by Broyden-Fletcher-Goldfarb-Shanno (Broyden, 1970;Fletcher, 1970;Goldfarb, 1970;Shanno, 1970), to refine the model solution.

Measured temperatures and nutrient concentrations
The measured epilimnion temperatures at SHL2, in summer and autumn, were higher in 2013 than 2012 (Fig. 2). The maximum surface water temperature was 28 • C in 2013 (in August), whereas for 2012 it was 23 • C. Also, the summer stratification was stronger in 2013 than 2012 . These two years are interesting in terms of PP dynamics, as in 2012 the maximum observed PP (140 mmol C m −3 d −1 ) was three times higher than the one for 2013 (Fig. 2). In general, DIN concentrations in the epilimnion and metalimnion were higher in 2013 than 2012. The average DIN concentration over the winter and early spring period was ∼40 mmol N m −3 in 2013 and ∼30 mmol N m −3 in 2012 (Fig. 2). In contrast to DIN, DIP concentrations in the winter and spring of 2012 were higher than those in 2013 (Fig. 2).

Phosphorus budgets for 2012/2013
The budgets of the bioavailable P for the two productive seasons of 2012 and 2013 are summarized in Fig. 3. As explained in the Method section above, we estimate the net P-uptake as the net-ecosystemproduction (gross production minus respiration), shown as the right bars in Fig. 3, by balancing the P input, outflow, and the stock difference between spring and autumn in the productive surface layer. The bioavailable P input, relevant for the primary production in Lake Geneva, consists of: (i) DIP from the Rhône River, (ii) bioavailable P from all WWTPs of the entire remaining catchment, except those from the Rhône River catchment that are already included in term (i) and finally, (iii) all other bioavailable P from various sources.
The Rhône River DIP input is continuously monitored 6 km upstream of the lake in Porte-du-Scex by the Swiss NADUF program (https: //www.hydrodaten.admin.ch/en/2009.html). The DIP annual inputs were 40 and 36 t P yr −1 for the years 2012 and 2013, respectively, and were 29 and 26 t P during the seven months of the two productive seasons (mid-March to mid-October) of 2012 and 2013, respectively. Term (ii), the WWTP discharges, contributed 34 and 35 t P during the productive seasons in 2012 and 2013, respectively. In summary, the first two monitored contributions of bioavailable P input results to 63 and 61 t P per productive season of 2012 and 2013, respectively.
The estimation of term (iii) via monitoring is almost impossible, as this non-point sources term consists of various contributions from agricultural land, rain depositions, and mineralization of flushed-in allochthonous organic matter. For this reason, we use the lake overall P balance to estimate the annual P input and scale it, using the water inflow, to the seven months production period. The annual P input is given by the sum of the P outflow (151 t P yr −1 in 2012; 121 t P yr −1 in 2013), the net-sedimentation (150 t P yr −1 ) equivalent to 11 g C m −2 yr −1 (Span et al., 1990;Loizeau et al., 2012), and the decline of the P content (-55 t P yr −1 ) calculated as the 5-year average . This balance yields annual P inputs of 246 and 217 t P yr −1 for 2012 and 2013, respectively. We feel comfortable with these estimates, as they agree well with the linear P budget model (Müller et al., , 2014 which would yield 228 and 215 t P yr −1 , respectively, and we consider the small discrepancies as well within the range of errors. Scaling to the seven productive months would lead to bioavailable P inputs of 168 and 146 t P for the two summers of 2012 and 2013, respectively. The P net uptake during the productive season can now be calculated from the input (168 t P in 2012 and 146 t P in 2013), the spring-to-autumn stock difference in the productive layer (222 t P in 2012 and 118 in 2013) and the outflow (-88 t P in 2012 and -71 t P in 2013). This leads to a P-uptake of 302 and 193 t P in 2012 and 2013, which lets us expect a difference in net-ecosystem production of 36%.

Model-data comparison
From our DA method, we obtain an optimized set of parameters ( Table 1) that simulates the observed PP and nutrient dynamics in 2012. Although the value of 3 is higher than those listed in Neumann (2000) and Neumann et al. (2002), it is well within the range of the maximum potential growth rates of cyanobacteria (Shimoda and Arhonditsis, 2016;Lee et al., 2018). Furthermore, we validate the model solution with the observations from 2013.

Temperature profiles
By nuanced variations (< 10%) in the scaling factors for wind and short wave radiation (GOTM parameters), we obtain the simulated temperature profiles for both years that compare well to the observed ones (Fig. 4, panels A and B). Also, the summer stratification is well resolved by GOTM. The maximum epilimnion temperatures in 2012 and 2013 are similar at 20 • C in late summer. L is the riverine dissolved P input (including discharge from sewage treatment plants) during the production season (March to October), and L accounts for the lake outflow. TP is TP available after the winter mixing (in March) for PP, and TP is the remaining total phosphorus (in October) after PP. Finally, Uptake is the phosphorus load taken up by phytoplankton and converted to particulate form during the productive period. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Nutrients and PP
The observed and simulated DIN, averaged over the top 20 m, show good agreement (Fig. 4, panels C and D). Although the model is calibrated only for 2012, the seasonal patterns of the observed DIN dynamics, e.g. spring drawdown, summer depletion and winter replenishment, are well reproduced for both years. However, the model appears to be less sensitive in resolving the monthly variations. Both the model and observations show higher concentrations of DIN in 2013 than 2012. In general, the model overestimates the observed DIP in winter of both years (Fig. 4, panels E and F). However, the time of DIP depletion matches with the observations. The model is also able to capture higher DIP concentrations in 2012 than 2013.  Fig. 5), the lower the external concentration of the respective nutrient and the stronger its limitation on primary producers. Strikingly, the period of P-limitation is longer in 2013 than 2012. In other words, our model suggests that DIP remained depleted in the epilimnion from mid June to end of September in 2013, thereby limiting phytoplankton growth over a longer part of the productive period. Fig. 6, panels A and B, shows the simulated and observed PP averaged over the top 20 m of the water column. There is a clear difference in the rates of PP between the two years. The maximum rate of PP in 2012 is nearly three times higher than that of 2013. The estimate of integrated PP (over the productive period) in 2012 is 567 g C m −2 , where as for 2013 it is 430 g C m −2 . In 2012, a large contribution to the total PP comes from the summer bloom, constituting mainly cyanobacteria as suggested by the model (Fig. 6, panel C). Conversely, in 2013 PP occurred mostly in spring, dominated by diatoms.

Deep mixing event of 2012
February 2012 was characterized by a strong cooling event that lasted almost two weeks (see the green shaded area in the estimated net surface heat flux shown in Fig. 7a). During the first week of February, the air temperature dropped significantly in a couple of days. This abrupt atmospheric cooling led to a massive heat release from the lake's water towards the atmosphere, reaching about 0 ≈ −650 W m −2 , and drove a sharp increase of the deepening rate of the estimated SML, ℎ SML , due to vigorous convection and entrainment of deep waters into the SML. This process is shown in the time series of the SML thickness (see blue triangles and green shaded area in Fig. 7b). Indeed, during the first two weeks in February 2012, the SML deepened to almost 200 m. The meteorological event was attributed to a cold wave associated with a large-scale atmospheric disturbance connected to the Artic Oscillation that hit from East Asia to Western, Southern Europe (WMO-RCC, 2013). Such an extreme event contrasts with the mild meteorological conditions registered in winter 2013 (compare blue line with the red line in Fig. 7a). Schwefel et al. (2019) proposed an empirical model to estimate the relative P flux into the epilimnion as a function of the surface mixed layer depth for Lakes Zug, Zürich, Geneva, and Tanganyika. Our results show that in the winter of 2012, SML deepened to 180 m, which is about 60% of the maximum depth of Lake Geneva (309 m). If we compare this to the model of Schwefel et al. (2019), it would correspond to 90% relative to its maximum P upwelling. According to our P budgets (based on observations), ∼ 350 tons of TP have been transported to the epilimnion (Fig. 3), which is ∼ 80% of the maximal P upwelling from the deeper layers of the lake in winter 2012.

Simulated temporal changes in phytoplankton
The average phytoplankton biomass (for epilimnion) over the productive period in 2012 (450 mg C m −3 ) is three times higher than the one for 2013 (Fig. 6, panels C and D). This could be because simulated nutrient limitation on growth of phytoplankton is weak in 2012 (Fig. 5).   Furthermore, there are significant differences in simulated dynamics of diatoms, dinoflagellates, and cyanobacteria between the years. Seasonal succession pattern of phytoplankton growth is clearly evident in 2012. Our results show that dinoflagellate bloom occurs earlier in 2012 than 2013. To understand this, we have to consider, at first, the timing of the onset of thermal stratification. Although the winter of 2012 is characterized by deep mixing events, seasonal stratification arrives much earlier (by April) which indicates an increase in water temperature (Fig. 7). As warm water stimulates the growth rate of dinoflagellates and also enables them to switch trophic strategy, they may have dominated in summer 2012. Wilken et al. (2018) observed a faster increase in abundance of heterotrophic flagellates, while decline in autotrophic biomass, under the climate change scenario during spring. According to the model, a massive cyanobacteria bloom may have developed in late summer 2012. Warm surface waters, stratified conditions, and abundant DIP in the epilimnion, could have elevated the growth rates of cyanobacteria. Modelling experiments of Gray et al. (2019) show that the increase in water temperature, shallower mixed layer depth, and meso-eutrophic conditions result in higher biomass and dominance of cyanobacteria.
Conversely in 2013, the simulated average concentrations of cyanobacteria in the epilimnion (below 5 mg C m −3 ) are much lower than 2012. This is because DIP gets depleted by early summer, and consequently limiting the growth of phytoplankton. However, there is no significant difference in concentrations of diatoms between the years.

Model biases
In general, ERGOM underestimates carbon fixation by phytoplankton in the deeper epilimnion layers (Fig. C.9 in Appendix C). This could be because of the dynamics of deep living phytoplankton, e.g. Planktothrix rubescens, are not resolved by ERGOM. It is to note that P. rubescens are known to dwell close or below the thermocline and are tolerant to low light conditions (Kerimoglu et al., 2017). Furthermore, the abundance of filamentous species like P. rubescens has increased, 2005 onwards, in Lake Geneva (Anneville et al., 2019). When we perform optimization, the missing dynamics of deep phytoplankton is, apparently, compensated by introducing uncertainty in the estimate of 3 .
Our results also reveal a few biases that may introduce uncertainties in the model solutions. Fig. C.10, in Appendix C, shows the temporal evolution of the simulated and observed temperature values averaged over the layer between 20 to 50 m. For both years, the model underestimates the observed temperature by 1-2 • C during late summer and autumn periods. River intrusions have significant effect on physical and biological dynamics in Lake Geneva (Giovanoli, 1990;Loizeau and Dominik, 2000). As river water plunges into the lake, it moves along the delta-channel as buoyancy-driven flows and perturbs the metalimnion temperature (Lambert and Giovanoli, 1988;Wüest et al., 1988). Bouffard and Perga (2016) performed a spatial sampling of Lake Geneva to study the effects of flood-driven river intrusions on the temperature, turbidity and oxygen gradients in the lake hypolimnion. They found that the turbidity peak between 30-50 m, associated to river inflow, coincided with a peak in the temperature profile at SHL2 station. This indicates that the turbid layer in the metalimnion changes, locally, the density of the water and increases the temperature. It is to note that our physical model do not include river intrusion. This could be a possible explanation for the discrepancy in the observed and simulated temperature profiles in the metalimnion.
Internal seiches in Lake Geneva is a known physical process (Imboden and Wüest, 1995;Lemmin et al., 2005). These are wind-caused, long-standing waves strongly constrained by lake bathymetry and summer stratification. They have modes of oscillation, the positive phase is associated with shallowing of the thermocline and the negative phase with deepening of the stratified layer (Lemmin et al., 2005). Perroud et al. (2009) applied different physical models to study multi-annual thermal profiles in Lake Geneva and found that the models yield rootmean-square error value of 2-3 • C in the layers with the greatest thermocline variability. However, the models that included a parameterization for internal seiche performed better than the ones without. Since the dynamics of internal seiches are not explicitly resolved by GOTM, this may contribute to uncertainty in the simulated deep epilimnion temperature profiles, particularly after strong wind events. In addition, the temperature of the surface mixed-layer can be affected by differential heating or cooling between the nearshore shallow regions and offshore interior waters. This process could be particularly relevant in fall and winter periods during which the shallow regions cool faster than interior waters, flushing cool waters towards the interior basin, which is not accounted in a 1D model (Fer et al., 2002;Ulloa et al., 2019).

Conclusion
We have applied a coupled model system to investigate the differences in PP and nutrient dynamics in Lake Geneva for two consecutive years (2012 and 2012). According to the physical model, the winter of 2012 (particularly in February) was exceptionally cold, and it had triggered a deep-convection in the lake that resulted in almost complete churning of the water column. This has redistributed nutrients to the productive surface layers. Consequently, throughout the production period, DIP replete conditions prevailed in 2012 and led to high algal biomass build-up. Furthermore, as water temperature increased in summer and autumn of 2012, abundances of dinoflagellates and cyanobacteria were stimulated. Conversely, in 2013 there was a strong P limitation of phytoplankton growth, and hence low biomass yield. ERGOM suggests that as DIP gets exhausted early in 2013, it may have limited the uptake of DIN by phytoplankton in the lake epilimnion. In a nutshell, an increase in TP by ∼52% in 2012 relative to 2013 due to deep-mixing event has contributed to ∼ 32% higher PP in 2012 (567 g C m −2 ) than 2013 (430 g C m −2 ).
To evaluate the model skills, we calibrated it with the observations from 2012 and validated with those of 2013. We found that the model is quite robust in simulating PP, phytoplankton growth, and nutrient dynamics as long as we account for temporal changes in riverine DIP and DIN fluxes. However, for more realistic modelling of PP in oligotrophic lakes, we propose to resolve growth dynamics of deep dwelling phytoplankton, e.g. Planktothrix rubescens, in biological models. In addition, changes in stoichiometry of seston with depth should also be considered in the model. Our study strengthens the utility of models to analyse and interpret field observations, and highlights the significance of deep mixing events in driving PP in lakes under reoligotrophication as climate warming continues.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. In this framework, and are the mean velocities in eastward and northward directions in units m s −1 , whereas ′ and ′ represent the fluctuation components resulting from the Reynolds' averaging, respectively. The quantities ′ , are are vertical eddy velocity in

A.2. Salinity equation
The vertical distribution of salinity in the water column is determined by: where , ′′ and −1 are salinity in units PSU, molecular diffusivity of salt in m s −2 and relaxation time given in second. The term −1 ( − ) describes relaxation of simulated to the observed salinity ( ).

Appendix B. Ecological model
Detritus. The rate of detritus production ( ) is given by: where , and are losses of phytoplankton, zooplankton and sediment to detritus in units d −1 .
Nutrients. The rate of ammonification ( ) is given by the sum of rates of excretion by zooplankton ( ), respiration by phytoplankton ( ) and remineralization of detritus ( ) and sediment ( ). Where is the areal phosphate flux in mmol P m −2 d −1 , and , 1 and 2 are constants. The fluxes of phosphate are taken in account as an additional source of in the model surface layer, and temperature and oxygen dependencies of phosphate release from sediment are explicitly resolved for the bottom layer. These are indicated by kronecker deltas ( , ) in Eq. (B.5). Sediment. The rate of sedimentation ( ) in ERGOM depends on the sinking speed of detritus ( ) and the critical stress ( ) at the lake bottom. where is the rate of change of oxygen production in units mmol O 2 m −3 d −1 and ( 1 + 2 + 3) is the total phytoplankton biomass in mmol C m −3 . The parameters and are O 2 -to-NH 3 and O 2 -to-NO 3 ratios. The term, ( , accounts for variations in O 2 -to-DIN stoichiometry with changes in and . As 2 is the source for only at the surface model layer with index , it is indicated by the kronecker delta ( , ) in Eq. (B.7). And is the thickness of the surface model layer in units m. The piston velocity ( ) is given in units m d −1 , and and are the saturation and water concentrations of oxygen in mmol O 2 m −3 . Note, depends on water temperature (Livingstone, 1993).