The response of small and shallow lakes to climate change: new insights from hindcast modelling

Small and shallow water bodies are a dominant portion of inland freshwaters. However, the effects of climate change on such ecosystems have rarely been quantitatively adressed. We propose a methodology to evaluate the thermal response of a small and shallow lake to long-term changes in the meteorological conditions, through model simulations. To do so, a 3D hydrodynamic model is forced with meteorological data and used to hindcast the evolution of a urban lake in the Paris region between 1960 and 2017. Its thermal response is analyzed through the definition of a series of indices describing its 5 thermal regime in terms of water temperature, thermal stratification and tendency to biomass production. Model results and meteorological forcing are analyzed over time to test the presence of monotonic trends and 3D simulations are exploited to highlight spatial patterns in the dynamics of stratification. The thermal regime of the study site underwent significant changes. Its response was highly correlated with three meteorological variables: air temperature, solar radiation and wind speed. Mean annual water temperature showed a considerable warming trend of 0.6°C/dec, accompanied by longer stratification and by 10 an increase of thermal energy available for biomass production. Water warming was significant during all four seasons, with maxima in Spring and Summer, while stratification and energy for phytoplankton growth increased especially during Spring and Autumn. Stratification only established in the deeper areas of the water body, possibly inducing heterogeneity in the release of nutrient from the sediment and in the development of harmful algal blooms. Numerous similar ecosystems might be experiencing analogous changes, and appropriate management policies are needed to preserve their ecological value. 15


Introduction
Lakes and reservoirs represent 3.7% of the Earth's non-glaciated continental area (Verpoorter et al., 2014), and often act as "sentinels" of climate change . They have experienced considerable warming along the past decades (O'Reilly et al., 2015;Schmid et al., 2014;Schneider and Hook, 2010;Piccolroaz et al., 2020), sometimes even accelerated in respect to the surrounding areas (Schneider et al., 2009). Climate change is expected to further deteriorate the ecological 20 status of a number of lakes worldwide that already suffer from eutrophication. In particular, changes in water temperature and in the patterns of thermal stratification could have a strong influence on the development of harmful algal blooms. Warmer The lake was originated in the 1940s by excavation and represents now a valuable and demanded recreational area. However, it suffers from strong eutrophic conditions and experiences severe harmful algal blooms, especially between late spring and early autumn. In particular, cyanobacteria such as Microcystis and Aphanizomenon, capable to produce toxins, often proliferate and become the dominant species in the lake. This leads regularly to bathing bans and to restrictions in the access to the lake.
For these reasons, the lake is subject to a periodic monitoring. In order to provide a faster and more responsive survey of the 70 main chemico-physical properties of the lake, a high-frequency (every 10 minutes) in situ measuring system has been installed at two different locations (A and B) during the spring 2015. Each measuring site is equipped with sensors at three depths: below the surface at 0.5 m depth, in the middle of the water column at 1.5 m and above the sediment at 2.5 m (Tran Khac et al., 2018). Water temperature measurements are recorded at the surface and bottom layers with a precision of 0.02°C and a resolution of 0.05°C through the thermal sensor SP2T10 (nke INSTRUMENT®), and through the MPx multi-parameter 75 sensor (nke INSTRUMENT®) at the middle of the water column, with a precision and a resolution of 0.05°C (see Fig. 1

-a).
High-frequency water temperature observations are used here for the calibration and validation of the hydrodynamic model. Lake Champs-sur-Marne is polymictic and its thermal behavior is strongly influenced by meteorological conditions. Between March and November periods of thermal stratification alternate with mixing and overturn of the water column. Depending on meteorological conditions, thermal stratification might form during the day and break up at night as well as last up to three or 80 four consecutive weeks.

The model 2.2.1 Presentation of Delft3D-FLOW
In order to simulate the thermal behavior of the lake and hindcast its evolution in the past decades, the FLOW module from the Delft3D modelling suite was used (Deltares, 2014). Delft3D-FLOW is a well known hydrodynamic model that has been applied 85 in various contexts, from estuaries to rivers, lakes and reservoirs (Piccolroaz et al., 2019;Chanudet et al., 2012;Soulignac et al., 2017). It solves the Reynolds averaged Navier-Stokes equations for an incompressible fluid under the shallow water and the Boussinesq assumptions. The time integration of the partial differential equations is done through an Alternate Direction Implicit method (Deltares, 2014;Leendertse, 1967). For the spatial discretization of the horizontal advection terms the Cyclic scheme was selected (Stelling and Leendertse, 1992). 90 The bathymetry and the two-dimensional mesh of the domain representing the study site are shown in Fig. 1-b. The surface of the lake is divided in 255 20 m × 20 m cells. The Z-model was implemented for the discretization of the vertical axes, with 12 fixed parallel horizontal layers of 30 cm thickness. It is generally accepted that horizontal layers help avoiding artificial mixing, improving model results in terms of thermal stratification (Hodges, 2014). Turbulent eddy viscosity and diffusivity were computed through the k-ε turbulence closure model. Background values were set to zero [m 2 .s −1 ] for vertical viscosity 95 and diffusivity, while they were set to 0.01 m 2 .s −1 , after Soulignac et al. (2017) and according to the grid size, for horizontal viscosity and diffusivity. Bottom roughness was computed through Chézy's formulation with the default value for the Chézy The computation of the heat exchange at the air-water interface is done through Murakami's model (Murakami et al., 1985).  [m]. The heat flux model computes the heat budget at the air-water interface by taking into account the net incident solar radiation (Q s ), the heat losses due to back radiation (long wave, Q b ) and evaporation (latent heat flux, Q e ), and the sensible convective heat flux (Q c ). The total upward heat flux through the air-water interface (Q) is therefore: Finally, evaporative mass flux is here neglected and water volume and depth are therefore considered as constant. This assumption makes it possible to analyze exclusively the impact of changes in the climatic forcing.

Meteorological input data
The meteorological forcing for this study comes from the spatialized SAFRAN (Système d'Analyse Fournissant des Renseignements Atmosphériques à la Neige) meteorological analysis system (Durand et al., 1993;Quintana-Seguí et al., 2008;110 Raimonet et al., 2017). SAFRAN is part of the SAFRAN-ISBA-MODCOU chain of reanalysis that covers the hydrological cycle over France, from meteorology to snow and ice formation to hydrology, respectively (Habets et al., 2008). SAFRAN integrates spatialized data from meteorological models with various sources of observations through data assimilation techniques, in order to create a consistent and spatially detailed record of meteorological data over the french territory. The data are spatialized on a regular square grid (8 km between each cell center) that covers the entire French Territory. The location of 115 Lake Champs-sur-Marne falls midway on the axes connecting the centers of SAFRAN cells number 1457 (North of the lake) and 1566 (South of the lake). The average of these two cells was therefore considered representative of the conditions over the study site and used as input for the hydrodynamic model. the daily variability of the thermal profile and improve the model performance. This is crucial in shallow water bodies, where thermal stratification and mixing can alternate between day and night. Specific humidity (SH) data had to be converted into relative humidity (RH) to match the input data set needed by Delft3D. This was done through the following formula: where w is the mixing ratio of water with dry air [kg.kg −1 ], the subscript s stands for saturation conditions and SH is the 125 specific humidity, numerically very close to the mixing ratio value. The saturation mixing ratio can be calculated as follows: where the atmospheric pressure (p atm ) was considered to be constant and equal to the global average: p atm =1013 hPa. The ratio between the air and vapor ideal gas constants (R a and R v , respectively) is equal to 0.622. The partial vapor pressure at saturation (e s ) is temperature dependent and can be estimated (in hPa) through the Magnus equation: where T is air temperature [°C]. The numerical coefficients in Eq. (4) were issued from (Alduchov and Eskridge, 1997).
Finally, in order to complete the set of meteorological input for Dleft3D, daily wind direction data were downloaded from the closest available MétéoFrance station (ID: 78621001 located in Trappes, roughly 40 km West of the study site), through the INRAE CLIMATIK platform (https://intranet.inrae.fr/climatik/, in French) managed by the AgroClim laboratory of Avignon, 135 France.

Calibration and validation
Delft3D-FLOW stands on a robust mathematical and physical structure and only few parameters have to be calibrated. Here, only those directly involved in the heat-flux model and in the wind module were calibrated: the Secchi depth [m], the mean cloud cover [-] and the wind drag coefficient [-]. The Secchi depth (H S ) is the parameter that defines water transparency. It 140 is correlated with the penetration of solar radiation in water through the light extinction coefficient (γ = 1.7/H S (Poole and Atkins, 1929)) and therefore has a strong influence on the stratification of the water column. In order to get a first estimate for the sky cloudiness parameter, cloud cover data from the MétéoFrance station in Trappes (ID: 78621001) were averaged over the calibration period. The wind drag coefficient was calibrated in order to take into account the presence of tall trees surrounding the contour of the lake, locally reducing wind speed. The calibration was done through a trial and error procedure based on Model results were compared to water temperature data at three depths (surface, middle and bottom of the water column) and two different locations (A and B). Root mean square error (RMSE) and mean bias error (MBE) were calculated to evaluate 150 model performances. For this purpose, high-frequency data were first averaged every hour to match the model output time step and cleaned from the outliers originated by periodic sensor maintenance. The latter were defined as sudden water temperature variations (> 1°C) over the 10 minutes separating two successive measurements, and consequently erased.

Indices for the characterization of the lake thermal regime
The thermal regime of the lake was assessed directly through the analysis of model results in terms of water temperature, as 155 well as through a series of indices that explore the phenology of stratification and highlight the relation between temperature and phytoplankton production. The indices characterizing the thermal regime of the lake are described here-after.

Stratification indices
In order to thoroughly characterize the phenology of stratification in Lake Champs-sur-Marne, two indices for the stability of the water column have been calculated: the Schmidt stability index and an index based on temperature difference between 160 surface and bottom layers. The Schmidt stability index is a well known parameter often used in limnological studies to estimate the resistance of a water body to mixing, and therefore its stability. It has been extensively used in scientific literature to describe the strength of stratification in lakes and, more recently, to analyze its evolution over time in relation to climate change (Vinçon-Leite et al., 2014;Niedrist et al., 2018;Kraemer et al., 2015;Livingstone, 2003) and algal blooms (Wagner and Adrian, 2009). column at one time instant. It has been here calculated following Idso's formulation (Idso, 1973), in which the vertical axes z is considered positive downwards from the surface to the maximum lake depth z M [m]: where: is the depth of the center of volume of the lake, ρ v [kg.m −3 ] is water density at the depth of the center of volume z v , ρ i is the mean uniform density, g [m.s −2 ] is the acceleration of gravity, V [m 3 ] and A 0 [m 2 ] are respectively the volume and the surface area of the lake, and A(z) is the surface of the horizontal section of the lake at depth z. Computed for each time step, the Schmidt stability can also be averaged over each year or season.
Water resistance to mixing as estimated by the Schmidt stability index is closely correlated to temperature stratification.

175
However, universal thresholds for the onset and breakdown of stratification are difficult to define based on this index and cannot be found in the literature, especially for shallow polymictic lakes. For these reasons, a second index based on temperature difference between surface and bottom layers (∆T ) is proposed. In order to assess the succession of stratification events in a polymictic water body, after Kerimoglu and Rinke (2013) and Magee and Wu (2017), the lake was considered to be stably stratified during one day if the minimum of ∆T is greater than 1°C during the whole day. This allows to identify all stably 180 stratified days. The sum over a year of the stably stratified days (SSD) is called the annual number of SSD. The sum can also be evaluated on a seasonal basis (according to the definitions in section 2.4), leading to the seasonal number of SSD.

Growing degree days (GDD) and number of growing days (NGD)
In order to quantify how changes in the thermal regime might impact biomass production and in particular phytoplankton growth, we introduce two indices: the growing degree days (GDD) and the number of growing days (NGD). Growing degree 185 days is a weather based indicator for biological growth, widely used in the field of agronomy. Based on air temperature, it gives an estimate of the rate of development and of the span of its phase for terrestrial plants and insects during the growing season.
where t is the time (here in days) with t 0 the reference day to start the calculation, ∆t is the time step (equal to 1 day), T i is the daily average of the modeled surface water temperature of day i and T base is a physiological threshold below which growth does not occur. GDD can be calculated on an annual or a seasonal basis by adjusting the values of t 0 and t. Annual GDD are 195 calculated starting from the first of January until the 31st of December, while seasonal GDD can be obtained by adjusting t 0 and t to the seasons defined in section 2.4. As in Dupuis and Hann (2009), the value of 4°C was selected for T base . Such value was chosen to be a representative baseline for the growth of the whole phytoplankton community in Lake Champs-sur-Marne, generally composed of cyanobacteria, green algae and diatoms. This results in a succession of algal blooms that spans almost the whole year, usually starting in February when water temperature in the study site surpasses 4°C, until the end of Autumn.

200
The number of growing days (NGD) at day t is defined as the number of days during which It represents a proxy for the duration of the period favorable to growth for the phytoplankton community. Annual NGD and seasonal NGD can be calculated by adjusting the values of t 0 and t as for the GDD. However, in a warm temperate climate water temperature never drops under the threshold of 4°C during Spring and Summer. This impedes any variability in the counts of NGD during these two seasons, that were therefore not addressed.

Long-term analysis
In the present paper we aim at hindcasting the long-term dynamics of a small and shallow urban lake in order to test the influence of climate change on such ecosystems. To do this we focus on 58 years, from 1960 to 2017, for which meteorological data are available from the SAFRAN reanalysis.

Long-term trends 210
The long-term hydrodynamic simulation starts on the first of January 1960. No data were available to set the initial conditions of the model, neither in terms of water temperature, nor in terms of current velocities. However, the model is strongly driven by the meteorological data and the influence of the initial condition vanishes after only a few days (Piccolroaz et al., 2019).
The effect on model results of small perturbations in water temperature initial conditions (± 2°C) was tested and resulted to vanish in 5 to 7 days. The model was therefore initialized with water at rest and with a uniform water temperature of 7 • C. This Model results at measuring site A are analyzed on an annual and seasonal basis for long-term trends. They are exploited directly in terms of water temperature (averaged over the water column) and to calculate the indices defined in section 2.3. Site A is located in one of the deepest part of the lake and therefore considered as representative of the general behavior of the water 220 body. This assumption will be further analyzed through a spatial analysis of the three-dimensional model results (see section 2.4.2). The presence of long term trends is tested (with a threshold for significance α = 0.05) through the Mann-Kendall test (Mann, 1945;Kendall, 1975), a non-parametric test for the individuation of overall monotonic trends performed here through the MATLAB softwere (Burkey, 2020). The Mann-Kendall test is often preferred to simple linear regression in the analysis of meteorological and hydrological time series (Tímea et al., 2017;Wang et al., 2020), as it does not require any assumption on the 225 8 https://doi.org/10.5194/esd-2020-51 Preprint. Discussion started: 24 July 2020 c Author(s) 2020. CC BY 4.0 License.
distribution of the analyzed dataset. Once a trend is detected, its strength is evaluated through the Sen's slope estimator (Sen, 1968), that uses a linear model to evaluate the intensity of the trend. In order to analyze both annual and seasonal trends, the year was divided into four groups of three months each: (i) January, February and March (Winter), (ii) April, May, June (Spring), (iii) July, August, September (Summer), (iv) October, November, December (Autumn). Meteorological forcing is crucial for this work, as it drives the hydrodynamic model and represents the only source of variability in our modelling configuration.

230
The presence of long-term trends in the meteorological dataset was also evaluated by applying the Mann-Kendall test and the Sen's slope estimator to their annual averages.

Spatial analysis of stratification
The phenology of stratification in Lake Champs-sur-Marne was further assessed exploiting the three-dimensional model simulations. The analysis was extended to the whole computational domain, with the objective of investigating the relation between 235 climate change, water depth and thermal structure in a shallow and polymictic water body. The annual number of stably stratified days (SSD) was therefore calculated for every computational cell over the complete simulation period, using temperature values from the surface and local bottom layers. In order to detect the presence of heterogeneity in the overall spatial distribution of stratification, the annual SSD calculated for each computational cell were averaged over the 58 years of the simulation to obtain a space-dependent (but time-averaged) representation of thermal stratification. 3 Results

Model calibration and validation
The model was calibrated on the year 2016 and validated from May 2015 to the end of 2017. Records from field campaigns show that the Secchi depth in Lake Champs-sur-Marne varies between 0.5 and 3 m; using this range of values, it was calibrated 250 and finally set to 1.2 m. Sky cloudiness was calibrated and set to 80%, and a uniform wind drag coefficient was set to 0.005 [-].
During the calibration periods, the model tends to slightly overestimate water temperature during summer heat peaks (with residuals always lower than 2°C) and to slightly underestimate water temperature during winter (again, with residuals always lower than 2°C). However, overall model performances are satisfactory for all three layers, with RMSE values between simulated and observed water temperature of 0.85°C, 0.78°C and 0.81°C at site A respectively for the surface (0.5 m), middle (1.5 RMSE values for sites A and B (ranging between 0.96 and 1.00°C). Fig. 2-a shows simulation results and observations during the whole validation period at 1.5 m depth at site A, where the data series is the longest, together with the related residuals ( Fig. 2-b); results were very similar for the two remaining sites. The model fits the data very well, especially between spring and autumn. During summer (respectively, winter), similarly to what was observed during calibration, it tends to slightly over-260 estimate (respectively, underestimate) water temperature, with residuals always lower in module than 2°C. Over the validation period the model has a low but non negligible bias (MBE) ranging between -0.2 and -0.3°C for sites A and B. Eventually, the number of SSD was calculated as defined in section 2.3.1 during the whole validation period at sites A and B for both the observed and simulated dataset. In case of gaps in the observation series, also simulation output were discarded from the calculation. The number of SSD between 2015 and 2017 calculated using water temperature observations is 122 for site A and 265 128 for site B. Using Delft3D simulations we found very similar results: 125 SSD for site A and 132 for site B. Overall, the model excellently reproduces high-frequency water temperature data, including the diurnal cycle, at both measuring sites.

Meteorological input data
Annual averages were calculated from 1960 to 2017 for the five meteorological variables used as inputs to the hydrodynamic 270 model and tested with the Mann-Kendall test. Strongly significant monotonic trends (p 0.05) were found for the air temperature, solar radiation and wind speed, as shown in Fig. 3. The Sen's slope estimator was used to test the intensity of the significant monotonic trends. Air temperature displays a considerable warming trend of 0.3°C per decade; solar radiation also shows a significant increasing trend, with an intensity of 3.5 W.m −2 per decade. Wind speed decreases quite sharply over time, at an estimated rate of 0.2 m.s −1 per decade. While the increase in air temperature appears extremely linear (see Fig. 3 sharp shift in the behavior of both solar radiation and wind speed appears around the year 1989 ( Fig. 3-b and -c, respectively).
Solar radiation and wind speed appear to be piecewise functions of time that can be roughly divided into two sections with sensibly different mean values, one until 1989 and the second one starting from 1990. Despite this non linear behavior, the presence of monotonic increasing (for solar radiation) or decreasing (for wind speed) trends is confirmed by the very low p-values obtained for these variables through the Mann-Kendall test. Finally, no significant trend was found for relative humidity and 280 wind direction. The two variables appear to be stationary, the former fluctuating around an annual average of roughly 80% and the latter around an annual prevailing wind direction of 200°N (South-West). Three of the five meteorological variables forcing the hydrodynamic model were therefore characterized by strongly significant monotonic trends along the past six decades, corroborating the idea of a changing climate in the region around the study site.

285
Long-term monotonic trends have been researched at site A on an annual and seasonal basis for: water temperature (vertically averaged), number of stably stratified days (SSD), mean Schmidt stability, growing degree days (GDD) and number of growing days (NGD). Figure 4 shows all the significant monotonic trends found from this analysis. On an annual basis, the Mann-Kendall test highlighted the presence of strongly significant increasing trends (p 0.05) for all variables.
Mean annual water temperature shows a very sharp warming tendency of 0.6°C per decade (see Fig. 4-a), even greater than 290 what was found for air temperature (0.3°C). The Pearson correlation coefficient (r) was calculated between water temperature and the five meteorological input variables in terms of annual averages in order to explain this behavior. Air temperature, solar radiation and wind speed were all strongly correlated with modeled water temperature, with correlation coefficients of 0.8 for solar radiation and air temperature and -0.9 for wind speed. Water temperature showed significant increase during all seasons, with higher slopes during Spring and Summer (0.8 and 0.7°C per decade, respectively), and a lower yet considerable intensity 295 during Autumn and Winter (respectively 0.4 and 0.5°C per decade).
The warming trend is accompanied by reinforced stratification. An increase in water column stability was highlighted on an annual basis by both stratification indices: the annual number of SSD increased on average of around three days per decade, while the Schmidt stability index increased of 0.9 J.m −2 per decade ( Fig. 4-b and -c, respectively). Despite a warming trend being present in all seasons, both stratification related indices showed significant increasing trends only during Winter (1 d and 300 0.4 J.m −2 per decade) and Spring (2 d and 2.6 J.m −2 per decade, for the seasonal SSD and the Schmidt index, respectively).
The sharpest increase in stable stratification was therefore found during Spring.
The analysis of the growing degree days and of the number of growing days shows the progressive improvement of conditions for biomass production. The black line in Fig. 4-d shows the evolution over time of annual GDD. Its pattern is highly correlated to that of water temperature and shows an increasing rate of 190°C·d per decade, with a shift around the year 1989. Even The changes found in the meteorological forcing clearly had an impact on the dynamics of the study site. The lake has sensibly warmed and its tendency to thermal stability has increased. Spring showed the sharpest trends in terms of water temperature, water column stability (Schmidt and SSD) and growing degree days, and might ultimately be the season suffering the strongest changes in terms of biomass production and algal blooms.

Spatial analysis of stratification
In order to examine the presence of spatial patterns regarding water column stratification, the three-dimensional model results were exploited for each cell of the domain and as well as for different depth classes. Figure 5-a displays a map of the study site representing for each cell of the domain the average over time (from 1960 to 2017) of the annual number of SSD, similarly to the spatial representation of sediment disturbance by waves found in Bachmann et al. (2000). The map shows a pronounced 320 spatial pattern, with stratification developing only in some regions. The deeper areas were stratified on average for more than 60 days a year, while the shallower northern part never significantly stratified. A strong linear correlation was found between water depth and number of stably stratified days, with a correlation coefficient between these two variables of 0.98.
The evolution over time of the spatial pattern of stratification was also investigated. The computational domain was divided into eight disjointed groups of cells depending on local water depth, which varies between 1.2 m and 3.6 m. For each year in 325 the simulation, the annual number of SSD was then averaged over each of the eight cell-groups, in order to obtain the eight time series of annual SSD shown in Fig. 5-b. Figure 5-b shows that, for this study site, no thermal stratification ever develops in areas with a depth lower than 1.8 m. The depth of the thermocline therefore appears to be always greater than 1.8 m. Figure 5-b also shows an evolution over time of stratification. All cells groups with a depth greater than 1.8 m tend to experience a higher annual number of SSD towards the end of the simulated period, especially starting from 1995. In order to confirm this result, 330 the time series of annual SSD calculated over each same-depth cells group (shown in Fig. 5-b through the color chart) was analyzed with the Mann-Kendall test. All cell-groups deeper than 1.8 m showed statistically significant (p < 0.05) increasing  as part of a large-scale analysis of observations in the northern hemisphere (Vautard et al., 2010), as well as on a global scale 345 through meteorological reanalysis (Torralba et al., 2017). An overall increase in surface solar radiation was recently found for Europe between 1983 and 2015, specifically of 3 W.m −2 .dec −1 for western Europe (Pfeifroth et al., 2018). Solar radiation and wind speed showed here a piecewise non-linear behavior (Fig. 3), with a shift around the late 1980s. SAFRAN is a coherent and spatialized data set, that still partially depends on the surface observation network and might be influenced by changes in the instrumentation (Vidal et al., 2010). However, the existence of a shift in global climate in the 1980s has been highlighted 350 by a number of studies using different data sources (Reid et al., 2016;Mariani et al., 2012;Gallagher et al., 2013).
Meteorological reanalyses usually cover multi-decadal periods and have the great benefit of being spatialized over vast portions of the globe. Even though their use in limnological studies is quite recent, they have already been used to simulate water temperature (Layden et al., 2016;Piccolroaz et al., 2020), stratification dynamics (Frassl et al., 2018) and phytoplankton distribution (Soulignac et al., 2018). As shown in this work, their use as external forcing to hydrodynamic models can yield, 355 provided that data are available for calibration and validation, to accurate simulations of the behavior of water bodies even in the absence of local meteorological observations. This could open to a great range of applications in the field of limnology and paleolimnology (e.g. Jenny et al., 2016). The proposed methodology allows to thoroughly reconstruct the behavior of a water body both in time and space, independently of its proximity to meteorological stations. Furthermore, the use of biogeochemical modules could give additional insights on oxygen, nutrients and phytoplankton dynamics.

360
The use of an extensive data set of high-frequency observations (recorded every 10 min, at three depths and two locations) allowed to test the ability of the model to reproduce not only the general seasonal water temperature pattern, but also daily and sub-daily dynamics. In this respect the model performed very well, with RMSE values always lower than 1°C and comparable to those obtained in other studies with similar applications (Kerimoglu and Rinke, 2013;Magee and Wu, 2017;Lee et al., 2018) through less frequent observations.

365
Significant increasing trends were detected for water temperature both on an annual and seasonal basis. The highest warming was found during Spring and Summer (0.8 and 0.7°C.dec −1 , respectively), while it was weaker during Autumn and Winter (around 0.4°C.dec −1 ). Overall, mean annual water temperature increased at a considerable rate of 0.6°C per decade, greater than that of air temperature (0.3°C.dec −1 ), a behavior also highlighted by other studies (Austin and Colman, 2007;Schneider et al., 2009). A large-scale analysis showed how trends in summer surface lake water temperature are globally highly variable, 370 comprised between -0.7 and 1.3°C per decade (O'Reilly et al., 2015). After this study, lakes subject to similar changes in the meteorological forcing to our study site, namely in terms of air temperature and solar radiation, display an average summer warming trend of 0.53°C per decade, only slightly smaller than the result of our work. Mean annual water temperature showed a non-linear behavior (see Fig. 4-a), probably caused by that of solar radiation and wind speed. In fact, similarly to what was found by Magee and Wu (2017), mean annual water temperature was highly correlated (i.e. r ≥ |0.8|) with air temperature, 375 solar radiation and wind speed. This suggests that meteorological variables might have additive effects that concur to enhance the response of the dependent variable. These effects might be particularly intense for wind over small and shallow lakes, due to their low volume to surface ratio.
Both stratification related indices (Schmidt stability and SSD) showed a mean annual increasing trend, mainly driven by an increase during Spring (respectively 2 d and 2.6 J.m −2 per decade). Similar values were recently found for shallow water 380 bodies in other long-term studies (Magee and Wu, 2017;Moras et al., 2019). The alternation between stratification and mixing in polymictic water bodies strongly influences the distribution and availability of nutrients over the water column (Song et al., 2013). The length of stable stratification and the frequency of mixing events have effects on the dynamics of sedimentation and resuspension of both nutrients and organic matter. This in turn might have strong effects on phytoplankton productivity and on the occurrence of harmful algal blooms, with impacts on the ecological state of the water body. To address this issue, the GDD 385 and the NGD were introduced to link the thermal regime to the biological productivity of the study site. Both the annual GDD and the NGD showed a significant increasing trend. The seasonal GDD were found to significantly increase during all seasons, while the seasonal NGD increased significantly only during Autumn and Winter.
The thermal regime is a key factor in the regulation of the biogeochemical cycle and in the development of algal blooms.
Small and shallow lakes react rapidly to climatic changes, partially due to their low volume and heat capacity. Changes in 390 water temperature can alter the composition of phytoplankton groups and their succession throughout the year. Warmer water temperatures favor species with higher optimum temperatures such as cyanobacteria, capable of producing toxins (Paerl and Huisman, 2008). Stratification induces a separation between the sediment and the surface layers, influencing the distribution of nutrients and biomass over the water column. Due to their ability to migrate upwards towards the water surface, an increase in water column stability might also favor cyanobacteria over other species (Dupuis and Hann, 2009). Furthermore, in polymictic 395 water bodies, a mechanism of accumulation (during stratification) and release (during subsequent mixing) of nutrients from the sediment has been suggested, causing multiple pulses that act as an internal nutrient source (Song et al., 2013;Wilhelm and Adrian, 2008). The combination of increasing trends for water temperature, stable stratification and the widening of the growing season can favor the occurrence of harmful algal blooms, further deteriorating ecosystems that are often already eutrophic (Winder and Sommer, 2012;Jones and Brett, 2014;Noble and Hassall, 2015). A number of small urban lakes 400 similar to our study site might be undergoing similar changes as Lake Champs-sur-Marne. These changes might be especially sharp during Spring, that showed here the greatest increase for water temperature, GDD and SSD. The timing, composition and intensity of spring blooms are extremely important in determining the succession of blooms in the subsequent months (Townsend et al., 1994;Sommer and Lengfellner, 2008;Lewandowska and Sommer, 2010).
Stratification affects many aspects of lake productivity, including nutrient recycling, habitat conditions for microorganisms 405 and their distribution over depth (Hanna, 1998;Wilhelm and Adrian, 2008;Song et al., 2013). The spatial analysis of stratification showed a strong linear correlation between the number of annual SSD and local water depth. The deeper part of the study site does experience considerable stable stratification throughout the year (up to around 60 days on average), even though for non-consecutive periods of time. The increase in the annual number of SSD found for site A is shared with all cell-groups with water depth greater than 1.8 m. However, this is not the case in the shallower northern part of the basin, that never show 410 stable stratification. Depending on the bathymetry, these spatial patterns might be strong in shallow water bodies, inducing heterogeneity in nutrients, phytoplankton and oxygen concentrations. Such heterogeneity might be even stronger in large water bodies, and can only be thoroughly inferred through three dimensional models (Gong et al., 2016;Frassl et al., 2018).

Conclusions
In this work, the long-term hydrodynamics of a shallow urban lake were profitably reconstructed through model simulations 415 from 1960 to 2017. A series of indices were proposed with the objective of thoroughly describing the thermal regime of shallow water bodies, in relation with stratification dynamics and biological productivity. The meteorological data set was derived from the SAFRAN reanalysis and showed a significant increase in air temperature and solar radiation and a significant decrease in wind speed, with a regime shift in the late 1980s. Simulation results show that small urban lakes react rapidly to external meteorological conditions, with only limited resilience to climatic shifts. The increase in water temperature cannot be 420 explained by air warming only. The additive effect of increasing solar radiation and decreasing wind speed acts on different terms of the heat budget at the lake surface, enhancing the changes found in the lake. Water warming (0.6°C/dec) is much quicker than air warming (0.3°C/dec), and especially intense in Spring, as is the lengthening of stratification. This could have favored early phytoplankton blooms, the development of cyanobacteria and ultimately the degradation of the whole aquatic ecosystem. Furthermore, the heterogeneity found in the spatial distribution of thermal stratification might concur to locally 425 create favorable conditions for algal blooms, in terms of nutrient availability or warmer surface water temperature. The use of three-dimensional models is needed to thoroughly infer the dynamics of a water body, including horizontal patterns. Small and shallow lakes are extremely widespread ecosystems, and a number of them might be experiencing analogous changes. Our results suggest the urgent need for appropriate management in order to preserve their ecological value.
Code and data availability.