© F. Soulignac et al., Published by EDP Sciences 2019

1 Introduction

The efficiency of ecosystem services provided by lakes is strongly influenced by lake water quality. Maintaining the water quality is crucial to guarantee long term availability of ecosystem services (Keeler et al., 2012). Several directives have therefore been issued to protect lake water quality around the world such as the Clean Water Act in the United-States of America and the Water Framework Directive (WFD) in Europe. The main objective of the WFD is to achieve a “good” ecological status for all surface waters (EU, 2000). Ecological status of lakes is classified based on biological, physico-chemical and hydro-morphological indicators. According to the French Ministry of the Environment, Energy and the Sea, the assessment of lake ecological status relies on samples from a unique station and a limited number of measurements throughout one or two year within a six years' management plan. The collected data are assumed to be representative of the lake system.

However, uncertainties in the evaluation of lake ecological status have repeatedly been pointed out (Hering et al., 2010; Reyjol et al., 2014). They may arise from properties of the ecosystem, which are not well addressed by the sampling protocol, such as spatial and temporal heterogeneities, which are observable from multi-sites or more frequent in situ sampling. Satellite observations are also able to enhance the spatial and temporal observations of physical, chemical and biological variables such as phytoplankton. Such spatio-temporal heterogeneities have profound effects on, or result from, the interactions between organisms themselves and their environment. They thus have consequences on a lake ecosystem quality assessment and might lead to misleading results and erroneous assignments to ecological status classes (Bresciani et al., 2011).

The location of the sampling site and the actual definition of quality classes with strict boundaries can also been seen as uncertainties (Kolada et al., 2014). The evaluation of lake ecological status is based on several annual metrics, which have to be combined to provide a single value describing the overall status in which the metric with the lowest score determines the status (the so-called “one out, all out” principle). The way of combining these metrics can also be seen as a further source of uncertainties (Moe et al., 2015). The lake ecological status is generally re-evaluated on basis of a six years' management plan. Therefore, inter-annual variability can also introduce uncertainties (Sondergaard et al., 2016).

Coupled hydrodynamic and ecological modeling did not typically appear in the usual solutions to quantify uncertainties as outlined above. In fact, a few studies applied such an approach in relation to the implementation of the WFD. Trolle et al. (2008) used an ecological model to estimate the level of phosphorus reduction needed to achieve a “good” ecological status, and Girbaciu et al. (2015) used an ecological model to simulate a lake water level and quality. These studies showed that ecological models simulate satisfactorily both physico-chemical and biological variables, the ones which can be used in lake ecological status assessment. Therefore, models are now interesting tools that provide relevant information and supplementary data to implement or test water quality assessment protocols.

We proposed here to use the three-dimensional (3D) modeling approach to quantify the role of lake spatio-temporal heterogeneities in the assessment of WFD ecological status. The aim of this study was to evaluate the impact of spatio-temporal heterogeneities on the evaluation of the water quality of Lake Geneva (France/Switzerland) by using the Delft3D model. The first objective was to estimate the sensitivity of the lake ecological status assessment to the timing of the sampling dates in one year and to the location of the sampling station. The second objective was to describe the spatial heterogeneity of the ecological status within the entire lake in order to evaluate the representativeness of the historical sampling stations, and to locate the more representative sampling area for the assessment of WFD ecological status.

2 Material and methods

2.1 Study site

Lake Geneva (46.45°N, 6.53°E in the World Geodetic System 1984 (WGS84), altitude: 372 m) is a large and deep peri-alpine lake located in the western part of the Alps, on the border between France and Switzerland (Fig. 1). Its mean depth is 153 m, its maximal depth is 309 m and its surface area is 580 km². It is thermally stratified during much of the year, never freezes over and does not undergo complete mixing every year. Its main tributary is the Rhone River, which accounts for 75% of the total inflow (average from 1965 to 2015). The lake is monitored by the French alpine lakes observatory (SOERE-OLA) as part of a long-term in situ monitoring program ordered by the International commission for the protection of Lake Geneva (CIPEL) for water quality and by the Canton of Geneva. Lake Geneva has been eutrophic for several years, and measures to reduce phosphorus loading were first implemented in the 1970s, leading to a decrease in phosphorus concentrations starting in the early 1980s. Nowadays, annual average total phosphorus concentration is 18 μg l⁻¹ (Tran Khac et al., 2018), classifying the lake as mesotrophic (Jacquet et al., 2014).

Fig. 1 Lake Geneva contour, isodepths (100, 200 and 300 m), numerical domain (curvilinear grid), location of the five main tributary inflow (Aubonne, Dranse, Rhône, Venoge, Versoix) and outlet and location of the two sampling stations (SHL2 and GE3).

2.2 Data

The variables used in this study were the water temperature, concentrations of chlorophyll-a (Chla), phytoplankton and zooplankton, dissolved oxygen (DO), ammonium (NH₄), nitrate (NO₃), orthophosphates (PO₄), silicate (SiO₂), total phosphorus (TP), particulate organic matter (POC, PON and POP), dissolved organic matter (DOC, DON, DOP), particulate inorganic matter (PIM) and transparency (SDD, Secchi disk depth). The years, 2009, 2010, 2011 and 2012 were chosen for this study because the hydrodynamics model had already been validated by Soulignac et al. (2018). The data were available from the SOERE-OLA database and the Canton of Geneva. Sampling took place at two stations from 2009 to 2012, straight up the deepest point of the lake in the large basin (SHL2: 46.45270°N, 6.58872°E (WGS84), depth: 309 m) and in the small basin (GE3: 46.29721°N, 6.21994°E (WGS84), depth: 72 m) (Fig. 1). At SHL2, sampling was conducted twice a month, except in winter, when it was carried out once a month. At GE3, sampling frequency was once a month between 2009 and 2012. Water samples taken at SHL2 for Chla were collected at ten depths between the surface and 30 m (i.e., at 0, 1, 2, 3.5, 5, 7.5, 10, 15, 20 and 30 m). At GE3, Chla was measured in an integrated sample from the surface to 20 m depth. Chlorophyll-a concentration was analyzed by spectrophotometry (Strickland and Parsons, 1968) after extracting in 90% acetone. Water samples for estimating phytoplankton species were collected by using an integrated bell sampler respectively between the surface and 18 m depth. Zooplankton was sampled by vertical tow from 50 m to the surface using a 200 μm plankton net. Phytoplankton and zooplankton composition were obtained by microscopy. Phytoplankton abundance was expressed as species biovolume per volume unit, and zooplankton abundance in number of individuals for a particular species per surface unit. For DO, nutrients, particulate and dissolved organic and inorganic matters, samples were collected at nine depths (i.e., at 0, 2.5, 5, 7.5, 10, 15, 20, 25 and 30 m), except at GE3 where 25 m was not sampled. Samples were analyzed by standard methods from AFNOR (the French standardization association).

2.3 Modeling

The open source Delft3D software (version 4.03.00) used in this study has been applied to lakes of different sizes all around the world (Chanudet et al., 2012; Kacikoc and Beyhan, 2014; Razmi et al., 2013; Soulignac et al., 2017). Delft3D is composed of a 3D hydrodynamic and a 3D water quality modules. For this study, the hydrodynamic model has been validated, and its configuration was described in Soulignac et al. (2018). The model was calibrated by performing a sensitivity analysis on several physical parameters (background coefficients for horizontal viscosity and diffusion of heat, coefficient for heat transfer at the atmosphere-lake interface), and it was validated based monthly temperature profiles taken at the two monitoring stations, SH2 and GE3. The model was applied on a grid representing Lake Geneva. The surface of the numerical domain was meshed with a curvilinear grid composed of 591 cells of about 1 km² area (Fig. 1). In the vertical direction, 100 Z-layers were used. Their thickness varied from 1 m at the surface to about 7 m at the bottom. The variation in the layer thickness was equal to 1.02.

The water quality module solves, for each grid cell and for each state variable, the advection-diffusion-reaction equation which can be simplified as follow:$$M_i^{t+\mathrm{\Delta }t}=M_i^t+\mathrm{\Delta }t\times {\left(\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}\right)}_{Tr}+\mathrm{\Delta }t\times {\left(\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}\right)}_P+\mathrm{\Delta }t\times {\left(\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}\right)}_S$$ $M_i^t$ and $M_i^{t+\mathrm{\Delta }t}$are the mass of a given state variable at the beginning and at the end of a time step Δt. ${\left(\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}\right)}_{Tr}$ represent changes by transport, ${\left(\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}\right)}_P$, changes by physical, chemical or biological processes, and ${\left(\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}\right)}_S$ changes by sources (e.g., river discharges). Ecological processes include the primary production, nitrification of NH₄, denitrification of NO₃, reaeration (transfer of DO from the atmosphere towards the lake surface layer), decomposition of organic matter, oxygen consumption by nitrification and organic matter decomposition, phytoplankton grazing by zooplankton, and particle settling.

The BLOOM module of Delft3D was used to simulate the growth and mortality of phytoplankton depending on nutrient availability, light intensity in the water column, and water temperature. BLOOM uses a routine to calculate the light intensity encountered by the phytoplankton over a 24 h time step (Bigelow et al., 1977; Los, 2009). The four phytoplankton groups the most abundant in Lake Geneva (diatoms, flagellates, green algae and cyanobacteria) were simulated. The calculation of the algal biomass was performed for each group, and the four groups were used to derive the total amount of simulated Chla based on specific C:Chla ratio for each group (Los, 2009). The default parameters of the BLOOM module were used, except for the specific light extinction coefficients of spring algae species (diatoms and flagellates). These coefficients were adjusted in order to remain within the realistic range given by Capblancq (1995) (i.e. about 0.016 m²g Chla⁻¹) and better match with Chla and SDD observations.

The grazing module CONSBL was used to model the grazing of phytoplankton and detritus by zooplankton. Primary consumer's biomass was included in the model as a forcing function, which allowed the control of zooplankton biomass, and the intensity of grazing on phytoplankton. In the case of Lake Geneva, it is feasible because zooplankton densities are available, and it is known that zooplankton has a strong impact on the abundance of phytoplankton in this lake (Anneville et al., 2019).

The simulated Secchi disk depth (SDD) was calculated for each grid cell by the depth at which the simulated light intensity was equal to 15% of the light intensity forced at the surface. The attenuation of light in the water column was calculated as the sum of the background coefficient, plus specific extinction coefficients due to phytoplankton, particulate organic and inorganic matters. The simulated total phosphorus concentration (TP) was calculated by taking into account dissolved phosphorus form (PO₄) and particulate (POP), as well as phosphorus included in phytoplankton.

The years, 2009, 2010, 2011 and 2012, were simulated independently because of computational time. Simulations started on January 1st and finished on December 31st. Annual simulations were initialized by imposing vertical profiles observed at SHL2 and GE3 for the large and small basins. These profiles resulted from a linear interpolation between the profiles measured before and after the 1st of January of the considered year. Initial conditions for NH₄, Chla and PIM were set to zero because observed values are very low. Simulations were realized by using hydrodynamic, meteorological and hydrological forcing. Hydrodynamic forcing came from validated hydrodynamic simulations (Soulignac et al., 2018). They consisted of hourly 3D maps of water velocity, water temperature and vertical dispersion. Meteorological forcing came down to daily incident solar radiation (RadSurf), the same values being applied to all grid cells (Fig. 2). These data came from the Consortium for Small-scale Modeling (COSMO) atmospheric model (http://www.cosmo-model.org). Hydrological forcing consisted of weekly data of the flow rate and ecological variables for the five major tributaries: Rhone, Dranse, Aubonne, Venoge and Versoix (Fig. 2).

Fig. 2 Simulation forcing: Time series of daily incident solar radiation (RadSurf), flow rate (Qin), ammonium concentration (NH4), nitrate concentration (NO3), total phosphorus concentration (TP), silicate concentration (SiO2), particulate inorganic matter (PIM) of the Rhône River, and biomass of zooplankton (Zooplank).

2.4 Model performance indicators

The annual root mean square error (RMSE) and relative root mean square error (RRMSE) were calculated as follows to quantify the difference between field data and modeling results.$$\text{R}\text{M}\text{S}\text{E}={\left(\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\left(S_i-O_i\right)}^2\right)}^{0.5}$$ $$\text{R}\text{R}\text{M}\text{S}\text{E}=\frac{\text{R}\text{M}\text{S}\text{E}}{\max (O)-\text{min}(O)}$$ O and S are synchronous annual vectors corresponding respectively to the depth averaged observed measurements and the depth averaged model results between the surface and 30 m. The length of O (and of S) is n, and fit the number of field campaigns of the year. The analysis was limited to the top 30 m because the euphotic layer was always less than 30 m deep.

In ecological modeling studies, Chla has generally the highest error values followed by nutrient concentrations (Arhonditsis et al., 2004). Model performances are considered acceptable if RRMSE is lower than 100% for Chla and lower than 50% for nutrient concentrations.

2.5 Lake ecological status

The French Ministry of the Environment, Energy and the Sea imposes a protocol for the calculation of the ecological status of lakes. Ecological status is calculated annually using data from vertical integrated samples in the euphotic zone. The calculation is based on four samples per year. In situ measurements are performed during two years, and the obtained ecological status states for the entire six years' management plan. This protocol has been designed to allow a comparison between lakes. These samples are taken over the growing season, above the deepest point of the lake and during different temporal windows (TW) (Fig. 3). For Lake Geneva, the first TW corresponds to the end of the winter mixing period, typically from mid-February to end of March. The three other TWs are during the thermal stratification. The second TW is from mid-May to end of June, the third TW extends from July to August and the fourth TW from September to mid-October.

The lake ecological status varies from “very good”, “good”, “medium”, “poor” or “bad”. Five variables used in the lake ecological status evaluation were tested in this study, Chla, NH₄, NO₃, TP and SDD. The annual metric for Chla is calculated as the mean value of three samples from TW2-4. Annual metrics for NH₄ and NO₃ are the maximum values among the four TWs. Annual metrics for TP and SDD are the median values over the four TWs.

The sensitivity of the lake ecological status to the timing of the sampling dates and the location of the sampling station were assessed by using simulations performed for the year, that presented the best fit between model and in situ data (Tab. 1). A set of four sampling dates was randomly selected from the outcomes of the model, one date representing one TW. Based on this selection of four samples, the annual metrics for all the five studied variables, Chla, NH₄, NO₃, TP and SD, were calculated for all grid cells (Fig. 3) following the WFD calculation protocol. This operation was repeated one thousand times. The question about the sensitivity of the assessment result to the timing of the sampling dates was addressed by comparing the one thousand ecological status measured in the cell containing SHL2 and GE3. Analyze of spatial heterogeneity was addressed by comparing the one thousand ecological status obtained in each one on the grid cells.

Finally, the representativeness of SHL2 and GE3 was assessed by using each combination of sampling dates and comparing the ecological status in SHL2 and GE3 to a theoretical ecological status of the whole lake. The theoretical ecological status was set as the median of the ecological status assessed in all the grid cells. A grid cell was estimated as representative of the lake if its status was the same as the theoretical ecological status (Fig. 3). This operation was repeated one thousand times to calculate the probability of representativeness for each grid cell.

Fig. 3 Illustration of the method for calculating the spatial variability of the lake ecological status and the grid cell representativeness from simulation results. Selection of four dates (step 1), calculation of the lake ecological status in each grid cell based on the five variables for one set of randomly selected sampling dates (step 2) and calculation of grid cell representativeness based on each variable (step 3). These steps were repeated for 1000 sets of randomly selected sampling dates.

3 Results

3.1 Performance assessment of Lake Geneva's ecological model

Simulation results at the two monitoring stations, SHL2 and GE3, showed that the ecological model reproduced well the seasonal variability of the five studied variables, Chla, NH₄, NO₃, TP and SDD, over the four simulated years (Fig. 4 and Tab. 1). The observed phytoplankton spring bloom, was mainly composed of diatoms in both in situ observations and model. The following decrease in biomass was correctly reproduced by the model (RMSE values for Chla ranged from 1.59 to 5.50 μg l⁻¹). The RMSE values for TP ranged from 2 to 7 μP l⁻¹ and the model captured quite well the increase and decrease in TP during the mixing (January to March) and production (March to November) periods, respectively. The model correctly predicted lower TP concentrations at GE3 compared to SHL2. The seasonal variability in NH₄ and NO₃ were well reproduced by the model, and the RMSE values ranged from 6 to 15 μg NH₄ l⁻¹ and 0.19 to 0.51 μg NO₃ l⁻¹. The spring decrease in transparency was well reproduced, and RMSE values for SDD ranged from 1.7 to 3.1 m. Simulation results presented an early onset of algal development in 2011, the model did not reproduce the peak of the Mougeotia gracillima in late summer 2009 and 2011, and the bloom of Aphanizomenon was also not predicted in 2012.

Year 2010 presented the lowest RMSE for Chla, TP and SDD (Tab. 1). Simulation results showed an over-consumption of NO₃ early in the year compared to the in situ measurements, especially at GE3. But in 2010, the ecological model reproduced well the vertical distribution of the four studied variables, Chla, NH₄, NO₃ and TP (Fig. 5). The model predicted correctly the spring subsurface algal development at SHL2. The model reproduced well NH₄ at SHL2 and GE3. For TP, simulations predicted well the observed increase in TP, which was uniformly distributed between the surface and 30 m depth, and the later decrease. Given the good match between in situ observations and simulated data, the analysis of the sensitivity of the lake ecological status were performed for the year 2010.

Fig. 4 Comparison between depth averaged observed and simulated data. Values are depth averaged between the surface and 30 m except for Chla at GE3 which was calculated between the surface and 20 m. Simulation results include the nine pixels around the cell corresponding to the given monitoring station.

Fig. 5 Comparison between the vertical distributions of observed and simulated data between the surface and 30 m depth in 2010. Simulation results are plotted by considering the grid cell corresponding to the given monitoring station.

3.2 Effect of sampling dates on the ecological status assessment at Lake Geneva's monitoring stations

Nitrate and ammonium concentrations as well as SDD indicated a “good” or a “very good” status respectively, for all the combinations of dates at SHL2 and GE3 (Fig. 6). However, an uncertainty remains for NO₃ because the concentrations are not well reproduced by the model (Sect. 3.1).

The Chla-based status predicted by the model was “medium” for 54% of the combinations of dates and “poor” for 46% at SHL2 (Fig. 6). Regarding GE3, the status was “medium” for 86% and “poor” for 14% of the sampling date combinations. Based on TP, the calculated status was “very good” for 28% of the combinations and “good” for 72% at SHL2. It was “very good” for 77% and “good” for 23% at GE3. These results showed that the date of sampling strongly influence the assessment of Lake Geneva ecological status when considering Chla and TP.

Fig. 6 Distribution of the lake ecological status among five classes, “very good”, “good”, “medium”, “poor” and “bad”, and calculated with 1000 sets of randomly selected sampling dates in 2010. The line thought the boxes indicates the median, the bottom and the top of the box correspond to the first and third quartile respectively.

3.3 Spatial variability of Lake Geneva's ecological status

The assignment of an ecological status to Lake Geneva in 2010 was sensitive to the location of the sampling station. This assignment is particularly determined by Chla, TP and SDD (Fig. 7). Simulation results for Chla, TP and SDD showed a high probability to reach a better status in the western part of the lake, which includes GE3, compared to the eastern part, which includes SHL2.

For Chla, the calculated status was mostly “medium” in the western part of the lake, “poor” in the largest part of the eastern part of the lake, and “bad” in the far eastern part of the lake in front of the Rhône River inflow (Fig. 7). The monitoring station SHL2 was located at the frontier between “medium” and “poor” status areas. The status based on TP presented also an east-west gradient of probability, ranging from “very good” to “poor” with a higher probability to present a worse quality status in the eastern part of the lake as compared to the western part. A few grid cells also presented a “bad” status at the site of the Rhône River inflow. The monitoring station SHL2 was located at the frontier between “very good” and “good” status. The status based on SDD was “very good” for the largest part of the lake, which includes SHL2 and GE3, but several grid cells presented a high probability for a “good” status in the eastern part of the lake. In opposite, for NH₄ and NO₃, results indicated a “very good” and a “good” status respectively for all the grid cells. For NO₃, an uncertainty remains because the concentrations are not perfectly reproduced by the model (Sect. 3.1).

Fig. 7 Spatial variability of the lake ecological status calculated with 1000 sets of randomly selected sampling dates in 2010. The colors indicate the probability of each grid cell of being from “bad” to “very good” ecological status.

3.4 Representativeness of Lake Geneva's monitoring stations

The western part of the lake, which includes SHL2 and GE3, presented high representativeness values for Chla (Fig. 8). The representativeness value was 0.84 for SHL2 and 0.77 for GE3. All grid cells were fully representative based on NH₄ and NO₃. This was because ecological status based on these parameters presented few spatial variabilities (Fig. 7, Sect. 3.3). Based on TP, SHL2 was more representative than the eastern and western parts (Fig. 8). The representativeness values were respectively 0.85 and 0.56 at SHL2 and GE3. Based on SDD, SHL2 and GE3 were representative. Results showed that Lake Geneva monitoring station SHL2 was more representative than GE3.

Simulation results showed that the central region of the lake was the more representative for the consideration of the minimal representativeness of the five variables (Fig. 9). The most representative grid cell, given by maximal value of the minimal representativeness, was located 4 km west of the actual monitoring station SHL2.

Fig. 8 Lake representativeness in 2010. The colors indicate the representativeness of each grid cell vs. the whole lake.

Fig. 9 Lake minimal representativeness in 2010. The colors indicate the representativeness of each grid cell vs. the whole lake.

4 Discussion

4.1 Model performance

Model performances are considered acceptable if the relative root mean square error (RRMSE, Sect. 2.4) is lower than 100% for Chla, and lower than 50% for nutrient concentrations (Arhonditsis et al., 2004). In the study of Lake Eğirdir, Turkey, by Kacikoc and Beyhan (2014) using Delft3D, RRMSE for Chla was between 30% and 100%, and lower than 50% for nutrient concentrations. In our study, Chla RRMSE was 30% at SHL2 and GE3 for year 2010 which was used for lake ecological status sensitivity analysis (Tab. 1). The relative root mean square error for NH₄, NO₃ and TP were lower than 50%, except for NO₃ at GE3. An explanation could be that the phytoplankton module of our ecological model calculate the light intensity encountered by the algae over a 24 h time step without considering the horizontal dimension (Bigelow et al., 1977; Los, 2009). As a consequence, phytoplankton located in shallow littoral and transported horizontally to a deeper part of the lake would experience much more light than in reality. This might explain why phytoplankton grew too intense in the littoral area near GE3, and consumed too much NO₃. Our model, thus, presented good performances and could be safely used but with caution regarding NO₃ and Chla in shallow littoral zone.

Simulation results did not predict adequately M. gracillima and Aphanizomenon development. These limits in our results may results in the fact that the model was computed using four phytoplankton groups (diatoms, flagellates, green algae and cyanobacteria) and may fail in reproducing the dynamic of specific taxa, which may have characteristic properties and life-traits.

4.2 Variability of the lake ecological status

The majority of the studies linked to the WFD deals with the creation of metrics used to classify lakes by ecological status (Padisak et al., 2006; Pasztaleniec, 2016), the redundancy issue between several metrics (Kelly et al., 2016) and the intercalibration between different methods of lake water quality assessment (Poikane et al., 2014). The objective of our study is original because it aimed at exploring the representativeness of the sampling protocol in a context of spatio-temporal heterogeneity of the lake basin. Our results showed that the evaluation of Lake Geneva ecological status was sensitive to the timing of the four sampling dates (Fig. 6) and to the location of the sampling station (Fig. 7) when using metrics based on Chla, TP and SDD. Sensitivity to the timing of the sampling is due to strong short-term fluctuations in phytoplankton abundance. This variability has been underline by high frequency monitoring or satellite data (Kiefer et al., 2015), and for the first time we provide evidence that it can strongly impact results of the water quality assessment. Spatial heterogeneity in phytoplankton abundance are well known, and here results indicated that the location of the sampling station could bias the estimation of lake ecological status. Rimet et al. (2016) identified two types of areas with different status in Lake Geneva. A more oligotrophic status was measured near the shoreline with low population densities and more eutrophic areas situated in more densely inhabited areas.

Our results showed that the western part of Lake Geneva presented a better quality compared to the eastern part. Such a heterogeneity in water quality at the scale of the lake basin is in agreement with previous descriptions that depicted high Chla concentrations in the eastern part of the Lake (Kiefer et al., 2015). According to CIPEL, Lake Geneva is considered as oligotrophic in the eastern part and mesotrophic in the western part. Such a spatial gradient of water quality across the lake could be explained by two reasons. Firstly, the inflow of the main tributary, the Rhône River, in the eastern part of the lake brought phosphorus and particulate inorganic matter, which decreased the transparency. Secondly, the algal development was promoted in the eastern part of Lake Geneva close to the nutrients inputs by the Rhône, which was more sheltered, and where the thermal stratification persisted (Soulignac et al., 2018).

This study showed that an annual metric can be close to a boundary between two classes, thus the assessment of the ecological status becomes very sensitive to data uncertainties and to protocols for the data sampling (Kolada et al., 2014) with major implications in terms of management, especially if the evaluation was astride the “good” and the “medium” status following the WFD. In fact, “Member States should aim to achieve the objective of at least good water status by defining and implementing the necessary measures within integrated programs of measures” (EU, 2000). Consequently, an uncertain evaluation of the lake status as a “good” status will limit effort of ensuring the good water quality or, in the contrary, an uncertain evaluation of the lake status as a “medium” status will induce the implementation of unnecessary programs of measures, which could be expensive.

An ongoing question in the management of aquatic systems is whether or not the number of samples and/or sampling stations should be increased in order to reduce uncertainties. The monitoring plan is clearly a compromise between the cost of supplementary measures and the related benefits expected in terms of the quality of lake status estimation (Porst et al., 2016). In this study, we suggest that a well-designed 3D ecological model is a time and cost efficient way to complement existing monitoring plans. We highlighted the contribution of 3D ecological modeling not only to propose a better time and frequency sampling protocol, but also to interpolate the daily variability of the physico-chemical lake characteristics, and objectively estimate the annual and pluri-annual mean lake characteristic through a six-year management plan, in order to appreciate the whole lake physico-chemical characteristics and review the process of the status evaluation based on a single station, particularly in a large lake.

However, 3D coupled hydrodynamic and ecological modeling is applicable in the situations in which lakes are monitored, and for which in situ observations are available. In this study, in situ measurements provided were used to validate the modeling approach. Without real observations, this type of study could not have been performed.

4.3 Data representativeness for ecological status monitoring

Our simulation results showed that the central part of Lake Geneva was the most representative of the spatially averaged ecological status (Fig. 8). These results confirmed that the historical monitoring station of Lake Geneva (SHL2) was well chosen. Conversely, the easternmost part of the lake appeared to be quasi-uncorrelated with the rest of the lake. This can be linked to the occurrence of strong spatio-temporal heterogeneities in Chla as observed from satellite imagery (Kiefer et al., 2015). A complementary monitoring station located in the eastern part of the lake would therefore be an important asset to monitor transient events. This consideration based on the Chla parameters also holds for SDD and TP. Nevertheless, it is worth remembering that NH₄ and NO₃ parameters were well described by a single sampling station with a representativeness close to one whatever the station location considered. Therefore, the need of supplementary spatial information is critical for Chla, SDD and TP but not, or to a lesser extent, for NH₄ and NO₃.

Recent developments in optical remote sensing enabled to provide spatial information on Chla and SDD from satellite observations at a decametric scale (Lee et al., 2016; Pahlevan et al., 2017). Following the approach of our study, the exploitation of those remote sensing data along with the modeling tools should provide significant additional information to exploit better the sparse in situ time series developed in the framework of water quality monitoring.

Acknowledgments

This work was supported the French Agency for Biodiversity (AFB, previously the French National Agency for Water and Aquatic Environments (ONEMA)), ModeL contract N° 15000239, and by the European Space Agency (ESA) Scientific Exploitation of Operational Missions element (SEOM), CORESIM contract N° A0/1-8216/15/I-SBo. We wish to thank the French Alpine Lakes Observatory (SOERE-OLA) and the International Commission for the Protection of Lake Geneva (CIPEL). Data were from © SOERE OLA-IS, AnaEE-France, INRA Thonon-les-Bains, CIPEL (2016), developed by Eco-Informatics ORE INRA Team. We also wish to thank the Department of environment, transport and agriculture (DETA) of the Geneva water ecology service for providing additional in situ data as well as Isabel Kiefer for satellite data. Finally, we would like to thank Tineke Troost and Hans Los (Deltares) for productive discussions and Pierre Keraudren for English editing.

References

All Tables

All Figures

	Fig. 1 Lake Geneva contour, isodepths (100, 200 and 300 m), numerical domain (curvilinear grid), location of the five main tributary inflow (Aubonne, Dranse, Rhône, Venoge, Versoix) and outlet and location of the two sampling stations (SHL2 and GE3).
In the text

	Fig. 2 Simulation forcing: Time series of daily incident solar radiation (RadSurf), flow rate (Qin), ammonium concentration (NH4), nitrate concentration (NO3), total phosphorus concentration (TP), silicate concentration (SiO2), particulate inorganic matter (PIM) of the Rhône River, and biomass of zooplankton (Zooplank).
In the text

Fig. 3 Illustration of the method for calculating the spatial variability of the lake ecological status and the grid cell representativeness from simulation results. Selection of four dates (step 1), calculation of the lake ecological status in each grid cell based on the five variables for one set of randomly selected sampling dates (step 2) and calculation of grid cell representativeness based on each variable (step 3). These steps were repeated for 1000 sets of randomly selected sampling dates.

In the text

	Fig. 4 Comparison between depth averaged observed and simulated data. Values are depth averaged between the surface and 30 m except for Chla at GE3 which was calculated between the surface and 20 m. Simulation results include the nine pixels around the cell corresponding to the given monitoring station.
In the text

	Fig. 5 Comparison between the vertical distributions of observed and simulated data between the surface and 30 m depth in 2010. Simulation results are plotted by considering the grid cell corresponding to the given monitoring station.
In the text

	Fig. 6 Distribution of the lake ecological status among five classes, “very good”, “good”, “medium”, “poor” and “bad”, and calculated with 1000 sets of randomly selected sampling dates in 2010. The line thought the boxes indicates the median, the bottom and the top of the box correspond to the first and third quartile respectively.
In the text

	Fig. 7 Spatial variability of the lake ecological status calculated with 1000 sets of randomly selected sampling dates in 2010. The colors indicate the probability of each grid cell of being from “bad” to “very good” ecological status.
In the text

	Fig. 8 Lake representativeness in 2010. The colors indicate the representativeness of each grid cell vs. the whole lake.
In the text

	Fig. 9 Lake minimal representativeness in 2010. The colors indicate the representativeness of each grid cell vs. the whole lake.
In the text