Toward an open-access of high-frequency lake modelling and statistics data for scientists and practitioners. The case of Swiss Lakes using Simstrat v2.1

One-dimensional hydrodynamic models are nowadays widely recognized as key tools for lake studies. They offer the 10 possibility to analyse processes at high frequency, here referring to hourly time scale, to investigate scenarios and test hypotheses. Yet, simulation outputs are mainly used by the modellers themselves and often not easily reachable for the outside community. We have developed an open-access web-based platform for visualization and promotion of easy access to lake model output data updated in near real time (simstrat.eawag.ch). This platform was developed for 54 lakes in Switzerland with potential for adaptation to other regions or at global scale using appropriate forcing input data. The benefit of this data platform 15 is practically illustrated with two examples. First, we show that the output data allows for assessing the long term effects of past climate change on the thermal structure of a lake. The study confirms the need to not only evaluate changes in all atmospheric forcing but also changes in the watershed or through-flow heat energy and changes in light penetration to assess the lake thermal structure. Then, we show how the data platform can be used to study and compare the role of episodic strong wind events for different lakes on a regional scale and especially how their thermal structure is temporarily destabilized. With 20 this open-access data platform we demonstrate a new path forward for scientists and practitioners promoting a cross-exchange of expertise through openly sharing of in-situ and model data.


Introduction
Aquatic research is particularly oriented towards providing relevant tools and expertise for practitioners. Understanding and monitoring inland waters is often based on in situ observations. Today, the physical and biogeochemical properties of many lakes are monitored using monthly to bi-monthly vertical discrete profiles. Yet, part of the dynamics is not captured at this temporal scale (Kiefer et al., 2015). An emerging alternative approach consists in deploying long-term moorings with sensors 5 and loggers at different depths of the water column. However, this approach is seldom used for country-level monitoring, although it is promoted by research initiatives such as GLEON (Hamilton et al., 2015) or NETLAKE (Jennings et al., 2017).
It is common to parameterize aquatic physical processes with mechanistic models, and ultimately use them to understand aquatic systems through scenario investigation or projection of trends in for example a climate setting. In the last decades, many lake models have been developed. They often successfully reproduce the thermal structure of natural lakes (Bruce et al., 10 2018). Today's most widely referenced one-dimensional (1D) models include (alphabetic order) DYRESM (Antenucci and Imerito, 2000), FLake (Mironov, 2005), GLM (Hipsey et al., 2014), GOTM (Burchard et al., 1999), LAKE (Stepanenko et al., 2016), Minlake (Riley and Stefan, 1988), MyLake (Saloranta and Andersen, 2007), and Simstrat (Goudsmit et al., 2002). The results from these models are mainly used by the modellers themselves, and often not easily accessible for the outside community. 15 The performance of lake models is determined by the physical representativeness of the algorithms and by the quality of the input data. The latter include (i) lake morphology, (ii) atmospheric forcing, (iii) hydrological cycle (e.g. inflow, outflow and/or water level fluctuations), and (iv) light absorption. In situ observations, such as temperature profiles, are required for calibration of model parameters. To support this approach, it is important to promote and facilitate the sharing of existing datasets of observations among scientists and practitioners. Conversely, scientists and practitioners should benefit from the 20 model output, which is often ready-to-use, high-frequency and up-to-date. Yet, model output data should not only be seen as a tool for temporal interpolation of measurements. Models also provide data of hard to measure quantities which are helpful for specific analyses (e.g., the heat content change to assess impact of climate change, or the vertical diffusivity to estimate vertical turbulent transport). Models finally support the interpretation of biogeochemical processes which often depend on the thermal stratification, mixing and temperature. In a global context of open science, collaboration between the different actors 25 and reuse of field and model output data should be fostered. Such win-win collaboration serves the interests of lake modellers, researchers, field scientists, lake managers, lake users, and the public in general.
In this work, we present a new automated web-based platform to visualize and distribute the near real time (weakly) output of the one-dimensional hydrodynamic lake model Simstrat through an user-friendly web interface. The current version includes 54 Swiss lakes covering a wide range of characteristics from very small volume such as Inkwilersee (9 x 10 -3 km 3 ) to very 30 large systems such as Lake Geneva (89 km 3 ), over an altitudinal gradient (Lago Maggiore at from 193 m. a.s.l. to Daubensee at 2207 m. a.s.l.) and over all trophic states (14 euthrophic lakes, 10 mesotrophic lakes and 21 oligotrophic lakes, Appendix A). We focus here on describing the fully automated workflow, which simulates the thermal structure of the lakes and weekly updates the online platform (https://simstrat.eawag.ch) with metadata, plots and downloadable results. This state-of-the-art framework is not restricted to the currently selected lakes and can be applied to other systems or at global scale.

Model and workflow 5
We use the 1D lake model Simstrat v2.1 to model 54 Swiss lakes or reservoirs (see Appendix A for details of modelled lakes) in an automated way. Simstrat was first introduced by Goudsmit et al. (2002) and has been successfully applied to a number of lakes (Gaudard et al., 2017;Perroud et al., 2009;Råman Vinnå et al., 2018;Schwefel et al., 2016;Thiery et al., 2014).
Recently, large parts of the code were refactored using the object-oriented Fortran 2003 standard. This version of Simstrat provides a clear, modular code structure. The source code of Simstrat v2.1 is available via GitHub at: 10 https://github.com/Eawag-AppliedSystemAnalysis/Simstrat/releases/tag/v2.1. A simpler build procedure was implemented using a docker container. This portable build environment contains all necessary software dependencies for the build process of Simstrat. It can therefore be used on both Windows and Linux systems. A step-by-step guide is provided on GitHub.
In addition to the improvements already described by Schmid and Köster (2016), Simstrat v2.1 includes (i) the possibility to use gravity-driven inflow and a wind drag coefficient varying with wind speed -both described by Gaudard et al. (2017), and 15 (ii) an ice and snow module. The ice and snow module employed in the model is based on the work of Leppäranta (2014Leppäranta ( , 2010 and Saloranta and Andersen (2007), and is further described in Appendix B.
A Python script was developed to (i) retrieve the newest forcing data directly from data providers and integrate them into the existing datasets, (ii) process the input data and prepare the full model and calibration setups, (iii) run the calibration of the model for the chosen model parameters, (iv) provide output results, and (v) update the simstrat.eawag.ch online data platform 20 to display these results. The script is controlled by an input file written in JSON format, which specifies the lakes to be modelled together with their physical properties (depth, volume, bathymetry, etc.) and identifies the meteorological and hydrological stations to be used for model forcing. The overall workflow is illustrated in Figure 1. Table 1 summarizes the type and sources of the data fed to Simstrat. For meteorological forcing, homogenized hourly air 25 temperature, wind speed and direction, solar radiation and relative humidity from the Federal Office of Meteorology and Climatology (MeteoSwiss, CH) weather stations are used. For each lake the closest weather stations are used. Air temperature is corrected for the small altitude difference (see Appendix A) between the lake and the meteorological station, assuming an adiabatic lapse rate of -0.0065 °C m -1 . This correction is a source of error in high altitude lakes like Daubensee for which dedicated meteorological station would be needed. The cloud cover needed for downwelling longwave radiations are estimated 30 by comparing observed and theoretical solar radiation (Appendix C). For hydrological forcing, homogenized hourly data from the stations operated by the Federal Office for the Environment (FOEN) are used. For each lake, the data from the available stations at the inflows are aggregated to feed the model with a single inflow. The aggregated discharge is the sum of the discharge of all inflows, and the aggregated temperature is the weighted average of the inflows for which temperature is measured. Inflow data are often missing for small or high altitude lakes (Appendix A). Missing inflows and more generally 5 watershed data is a source of error in small alpine lakes, yet, such error can be compensated during the calibration process.

Input data
The light absorption coefficient abs [m -1 ] is either obtained from Secchi depth Secchi [m] measurements (for Inkwilersee, Lake Biel, Lake Brienz, Lake Geneva, Lake Neuchatel, Lower Lake Zurich, Oeschinensee, Upper Lake Constance, and Sihlsee), or is set to a constant value based on the lake trophic status. In the first case, the following equation is applied: abs = 1.7/ Secchi (Poole andAtkins, 1929, Schwefel et al. 2016). In the second case, abs is set to 0.15 m -1 for oligotrophic lakes, 10 0.25 m -1 for mesotrophic lakes, and 0.50 m -1 for eutrophic lakes. The values correspond to observations of Secchi depths in Swiss lakes (Schwefel et al. 2016) and fall into the decreasing range of transparency from an oligotrophic to eutrophic system (Carlson 1977). For glacier-fed lakes (typical above 2000 m) rich in sedimentary material, abs is set to 1.00 m -1 .
The timeframe of the model is determined by the availability of the meteorological data (air temperature, solar radiation, humidity, wind, precipitation). Initial conditions for temperature and salinity are set using conductivity-temperature-depth 15 (CTD) profiles or using the temperature information from the closest lake. We apply different data patching methods to remove data gaps from the forcing depending on the length of the data gap. For small data gaps with duration not exceeding one day, the dataset is linearly interpolated. In total < 1 % of the dataset is corrected using this approach. Longer data gaps of up to 20 days are replaced by the long-term average values for the corresponding day of the year. Only ~ 1.5 % of the dataset is corrected using this approach. 20

Calibration
Model parameters are set to default values, and four of them are calibrated (see Table 2). The parameters p_radin and and f_wind scale the incoming long-wave radiation and the wind speed, respectively, and can be used to compensate for systematic differences between the meteorological conditions on the lake and at the closest meteorological station. The parameter a_seiche determines the fraction of wind energy that feeds the internal seiches. This parameter is lake-specific, as it depends on the 25 lake's morphology and it's exposition to different wind directions. Finally, the parameter p_albedo scales the albedo of ice and snow applied to incoming shortwave radiation, which depends on the ice/snow cover properties. The calibration parameters were selected according to their importance for the model (e.g. based on previous sensitivity analysis), and their number was deliberately kept small in order to keep the calibration process simple and focused. Calibration is performed using PEST v15.0 (see http://pesthomepage.org), a model-independent parameter estimation software (Doherty, 2016). As a reference for 30 calibration, temperature observations from CTD profiles are used. Calibration is performed on a yearly basis, unless significant changes are made to either the model, the forcing data, or the observational data (e.g. release of a new version of Simstrat or delivery of a large amount of new observational data). For the eight lakes without observational data, parameters are set to their default value (see Table 2) with no calibration performed, and the lack of calibration is indicated on the online platform.

Output / Available data on the online platform
The online platform (accessible at https://simstrat.eawag.ch) is automatically fed every week with model results, metadata and plots for all the 54 modelled lakes (see Figure 2). It allows for efficient display and open sharing of the model results for 5 interested users. While the framework is here restricted to Swiss lakes, the code could be easily adapted to other lakes outside Switzerland and used at the global scale. From the model results, we directly obtain time series of several model output variables. Those dataset include temperature, salinity, Brunt-Väisälä frequency, vertical diffusivity, and ice thickness. In addition, we use the following known physical and lake-related properties: the acceleration of gravity = 9.81 m 2 s -1 , the heat capacity of water p = 4.18 • 10 3 J K -1 kg -1 , the volume of the lake [ • Schmidt stability: • Timing of summer stratification: we use a threshold based on the Schmidt stability to determine beginning and end 15 of summer stratification. The lake is assumed to be stratified for / ≥ 10 J m -3 . Using a different criterion (e.g., temperature difference between surface and bottom water) results in variations in the calculated stratification period; however, the general pattern among lakes remains similar.
• Timing of ice cover: we use the existence of ice to determine beginning and end of ice covered period.
From these results, we create static and interactive plots. The latter are created using the Plotly Python Library (see 20 https://plot.ly/python). The plots can be categorized as follows: • History (e.g., contour plot of the whole temperature time series, line plot of the whole time series of Schmidt stability); • Current situation (e.g., latest temperature profile); • Statistics (e.g., average monthly temperature profiles, long-term trends).
All Output and processed data are directly available from the online platform. 25

Results and discussion
Analysis of model output allows to compare the response of the different systems to specific events or to long term changes.
The Simstrat model web interface provides regional long-term high-frequency data updated in near real-time as output. This represents a novel way to monitor, analyse and visualize processes in aquatic systems, and, most importantly, grant the entire community direct access to the findings. The coupling between Simstrat and PEST provides an effective way to calibrate 30 model parameters. The uncertainty quantification finally allows an appropriate informed use of the output data. Yet, more advanced methods for both parameter estimation and uncertainty quantification such as Bayesian inference (Gelman et al., 2013) should be applied to Simstrat.
Out of the 46 calibrated lakes, the post-calibration root mean square error (RMSE) is < 1 °C for 17 lakes, between 1 and 1.5 °C for 15 lakes, between 1.5 and 2 °C for 8 lakes and between 2 °C and 3°C for 6 lakes ( Figure 3). There were too few in-situ observations on 8 lakes to perform a proper calibration and all parameters where thereby set to default values. Overall, the 5 performance is comparable to the RMSE range of ~0.7-2.1 °C reported in a recent global 32-lake modelling study using GLM (Bruce et al., 2018) also including Lake Geneva, Lake Constance and Lake Zurich. The correlation coefficient remains always higher than 0.93 suggesting also that the model successfully reproduce the thermal structure of the investigated lakes Overall, the quality of the results is better for lowland lakes than for high altitude lakes where local meteorological and watershed information are often missing. 10 We illustrate the potential of high-frequency lake model data with two examples: first by briefly showing the long-term changes caused by climate change in Lake Brienz (section 3.1), and secondly by investigating the differential response of lakes across Switzerland to episodic forcing (short-term extremes, section 3.2).

Long-term evolution of the thermal structure of lakes in response to climate trends
Over the period 1981-2015, yearly averaged simulated surface temperatures in Lake Brienz increased with a significant 15 (p<0.001) trend of +0.69 °C/decade ( Figure 4a). For the same period, monthly in situ observations indicate a similar trend of 0.72 °C/decade (p~0.07), while the trend of air temperature at the meteorological station in Interlaken is lower (+0.50 °C/decade, p<0.01). Based on physical principles, lake surface temperature is expected to increase less than air temperature , however Schmid and Köster (2016) also observed a higher trend in lake surface temperature than in air temperature for Lower Lake Zurich and assigned the excess warming to a positive trend in solar radiation. For the period 1981-20 2015, the ascending trend in solar radiation is 5 W/m 2 /decade that corresponds to an equilibrium temperature increase of about 0.2°C/decade. The warming rate at the surface of Lake Brienz is larger than observed trends in neighbouring lakes with reported increases of +0.46 °C/decade for Upper Lake Constance (1984( -2011( , Fink et al. 2014), +0.41°C/decade for Lower Lake Zürich (1981-2013, Schmid and Köster, 20161955-2013, Livingstone, 2003, +0.55°C/decade for Lower Lake Lugano (1972-2013, Lepori and Roberts, 2015. This can be explained by the lower light penetration in Lake Brienz (ranging from 25 ~1 m to ~10 m) compared to other lakes; the increase in solar radiation being distributed into a shallower layer and thereby warming slightly more the lake surface.
The temperature increase was significantly smaller in the hypolimnion, with a minimum trend at the lake bottom of 0.16 °C/decade (p<0.001), leading to a depth-averaged rate of temperature increase of 0.22 °C/decade (p<0.001). The temperature difference between the inflow and the outflow also contributes to the heat budget. While no significant change in the yearly 30 total discharge was observed at the gauging stations of FOEN for the inflows of the Aare and of the Lütschine rivers for the period 1981 -2015, the weighted inflow temperature increased by 0.26 °C/decade. The riverine temperature remains colder than the lake surface temperature leading to a yearly average loss of energy by through-flow of ~ -40 W/m 2 for 2015. This result is consistent with the recent observations of Råman Vinnå (2018) suggesting that tributaries significantly affect the thermal response of lakes with residence time up to 2.7 years (as Lake Brienz). The contribution of the river to the heat budget of Lake Brienz is also ~ 4 times larger than that previously estimated for Upper Lake Constance (Fink et al. 1994), a lake with 5 a longer residence time. The increasing difference over time between the inflow temperature and the outflow temperature (taken as the lake surface temperature) leads to a non-negligible cooling contribution from the river of ~ 0.14 °C/decade (p<0.05). The temporal change in the discharge and its temperature resulting from climate change should therefore be taken into account in studies attempting to predict the change in lake thermal structure. Such analyses can be extended to all modelled lakes. An inter-comparison of the temporal extent of summer stratification and winter ice cover period is illustrated in Figure 5. An altitude-dependent decrease of the duration of summer stratification is observed, along with a stronger corresponding increase in the duration of the inverse winter stratification from 1200 m. a.s.l. This is possibly linked to an altitude dependency of climate-driven warming in Swiss lakes, first reported by Livingstone et al. 25 (2005), which may be caused by a delay in meltwater runoff (Sadro et al., 2018). Here this process is not directly resolved but incorporated through the calibration procedure spanning all seasons.
In conclusion, the online platform provides all the data to estimate the past warming rate of lakes and evaluate how the different external processes contribute to their heat budgets. The change in the thermal structure depends mostly on the change in atmospheric forcing, yet, other factors such as the changes in discharge and temperature from the tributaries and the light 30 absorption into the lake should also be taken into account. We specifically show that the warming rate of the lake surface temperature significantly differs from that of depth-averaged temperature, thereby highlighting the benefit of using either in-situ observations resolving the thermal structure over the water column or hydrodynamic model output for assessing climate change impacts on lake thermal structure.

Event based evolution of the lake thermal structure.
A major drawback of traditional lake monitoring programs in Switzerland is the coarse temporal resolution, with measurements often performed on a monthly basis. Thought sufficient for direct long-term trend studies as shown in section 3.1 when 5 conducted over an extended period typically longer than 30 years (Gray et al., 2018). However, traditional monitoring programs cannot resolve the impact of short-term events and their consequences for the ecosystem. This is a strength of highfrequency (hourly time scale) lake modelling, which allows for simulation and comparison of the effects associated with rapid and often severe events such as storms. Based on high-frequency observations, Woolway et al. (2018) showed the effects of a major storm on Lake Windermere. They observed a decrease in the strength of the stratification, a deepening of the thermocline 10 and the onset of internal waves oscillations ultimately upwelling oxygen depleted cold water into the downstream river.
Furthermore, Perga et al. (2018) illustrated how storms could be just as important as gradual long-term trends for changes in light penetration and thermal structure in an Alpine lake.
Here we demonstrate how high-frequency model output can be used to study the influence of specific events on the thermal dynamics of lakes. As an example, we focus on the 28 th of June 2018 when Switzerland experienced a strong but by no means 15 exceptional storm with Northeasterly winds mainly affecting the North-Western part of the country -the mean wind speed during that day is shown spatially in Figure 6a. The evolution of the stratification strength, illustrated here by the Schmidt stability, is given in Figure 6b for one of the most affected lakes, Lake Neuchâtel (https://simstrat.eawag.ch/LakeNeuchatel, Figure 2). This lake, with the main axis well-aligned to synoptical winds, experienced a ~8 % decrease in the Schmidt stability over this half-day event. Yet, the effects were not long-lasting and the Schmidt stability reverted to its pre-storm value within 20 ~5 days (Figure 6b). This also resulted in a total increase of the lake heat content by ~1.4·10 16 J from the start of the storm to the time of recovery. We used the Schmidt stability recovery duration as a way to assess the short-term effect of the storm on the different modelled lakes. In Figure 6a, lakes are coloured based on the delay in Schmidt stability increase (in days) caused by the storm. The impact of the storm was not limited to Lake Neuchâtel but rather showed a regionally-varying pattern.
Particularly small-to medium-sized lakes in the North-Western parts of Switzerland were more affected than large lakes or 25 lakes located in the Southern part of Switzerland. However, the thermal structure of these lakes quickly reverted to the seasonal early summer warming trend.
So far, climate-driven warming has been recognized to cause an overall increase in lake stratification strength and duration, and a gradual warming of the different layers (Schwefel et al., 2016;Zhong et al., 2016;Wahl and Peeters, 2014). Air temperature trend was the most studied forcing parameter. Yet, the dynamics of extreme events (such as heat waves, drought 30 spells, storms), including their changes in strength and distribution, has been comparatively overlooked. Scenario exploration, climate change studies, or historical forcing reanalysis should be integrated in such web-based hydrodynamic platforms to assess their roles in modifying the lake thermal structures and heat storage.

Conclusion
The workflow presented in this paper allows openly sharing of high-frequency, up-to-date and permanently available lake model results for multiple users and purposes. We demonstrated the benefit of the platform through two simple case studies. 5 First, we showed that the high frequency modelled temperature data allows a complete assessment of the effect of climate change on the thermal structure of a lake. We specifically show the need to evaluate changes in all atmospheric forcing, in the watershed or through-flow heat energy and in light penetration to accurately assess the evolution of the lake thermal structure.
Then we showed that the high frequency modelled data can be used to investigate special events such as wind storms, there in-situ measurements under current temporal resolution are failing. More generally, these results are well suited for the 10 following applications and target groups: • For the public, the platform serves as an informative website enabling easy access to broad quantities of regional scientific results, with the intention of raising interest about lake ecosystem dynamics.
• For lake managers, the platform makes relevant information available, such as (i) near real time temperature and stratification conditions of the lakes, (ii) simple statistical analyses such as monthly temperature profiles and long-15 term temperature trends.
• For researchers, this work can facilitate (i) scenario modelling of any of the lakes, as the basic model setup is readyto-use, (ii) improvement of the lake model with addition of previously unresolved processes (e.g., resuspension with changed light properties), (iii) access to variables that were previously not or irregularly available (e.g., vertical diffusivity, heat content, stratification and heat fluxes), and (iv) specific comparative analyses, whereby a given 20 question can be investigated simultaneously over many lakes (e.g., the impact of climate change or a regional storm). By promoting a cross-exchange of expertise through openly sharing of in-situ and model data at high frequency, this openaccess data platform is a new path forward for scientists and practitioners.

Author contribution
The new version of Simstrat was developed by FB, AG and LRV. The workflow was developed by AG. The ice model was 5 developed by LVR. The concept of the workflow was defined by DB. All authors contributed to the validation of the model and interpretation of the results. AG and DB wrote the manuscript with contributions from FB, LVR and MS.

Competing interests
The authors declare that they have no conflict of interest

Acknowledgment 10
We thank Davide Vanzo for helping with the Docker and the scripts and Michael Pantic for helping restructuring the version 2.1 of Simstrat. We finally thank Marie-Elodie Perga for her comments on a preliminary version of the manuscript. The full list of acknowledgements regarding in-situ observations can be found here: https://simstrat.eawag.ch/impressum .  , 527-528, 493-506, doi:10.1016/j.scitotenv.2015.05.011, 2015: Effects of climate change on deep-water oxygen and winter mixing in a deep lake (Lake Geneva)-Comparing observational findings and modeling, Water Resources Research, 52(11), 8811-8826, doi:10.1002/2016WR019194, 2016 Smith, W. L.: Note on the Relationship Between Total Precipitable Water and Surface Dew Point, J. Appl. Meteor., 5(5), 726-727, doi:10.1175/1520-0450(1966) 005<0726:NOTRBT>2.0.CO;2, 1966. 5 Stepanenko, V. M., Mammarella, I., Ojala, A., Miettinen, H., Lykosov, V. and Vesala, T.: LAKE 2.0: a model for temperature, methane, carbon dioxide and oxygen dynamics in lakes, Geosci. Model Dev., 9(5), 1977-2006, doi:10.5194/gmd-9-1977-2016, 2016.       -1981-2018 Below the freezing point (ice formation) The ice module is activated as the water temperature in the topmost grid cell Tw (°C) drops below the freezing temperature Tf (°C). Tf can be set to zero for a vertical grid size ≤ 0.5 m, the user can adapt (raise) this value to fit coarser grids. If temperature is below the freezing point, the energy incorporated into the change of state Ef is calculated as here ρw (1000 kg m -3 ) is the density of fresh water, cpw the heat capacity of water (4182 J kg -1 °C -1 ) and z1 the height of the topmost grid cell. Ef and the latent heat of freezing lh (3.34·10 5 J kg -1 ) as well as the density of black ice ρib (916.2 kg m -3 ) are used for calculating the initial height of black ice hib (m) in Eq. B2, thereafter Tw is set equal to Tf If an ice cover is present and if the atmospheric temperature Ta (°C) is smaller or equal to Tf, the growth of black ice dhib/dt 10 continues as described in (Saloranta and Andersen, 2007).
Here ki (2.22 W K -1 m -1 ) is the thermal conductivity of ice at 0 °C and Ti (°C) the ice temperature calculated as There ks (0.2 W K -1 m -1 ) is the thermal conductivity of snow and hs (m) the height of the snow layer. When Ta is smaller than the snow temperature (default set to 2 °C) water equivalent precipitation pr (m hour -1 ) is turned into fresh snow hs_new ( here ρs (kg m -3 ) is the snow layer density kept within ρs0 < ρs < ρsm with the maximum snow density set to 450 kg m -3 , C1 (5.8 m -1 hour -1 ) and C2 (0.021 m 3 kg -1 ) are snow compression constants, and ws (m) is the total weight above the layer under compression expressed in water equivalent height. 5 If the snow mass ms (kg m -2 ) becomes heavier than the upward acting buoyancy force Bi (kg m -2 ), white ice with height hiw (m) and density ρiw (875 kg m -3 Saloranta, 2000) is formed between the snow and the black ice layers to achieve equilibrium between Bi and ms.
In this model, we assume continuous supply of water through cracks in the black ice to form white ice. The formation of white ice takes place instantaneously each time step and we do not consider the influence of pools under the snow for melting or short wave irradiance penetration.

Above the freezing point (melting)
If an ice cover is present and if Ta > Tf melting starts. Each layer melts from above through the atmospheric interface and by 15 penetrating short wave radiation where Hx_y (W m -2 ) is the layer-dependent heat flux (in the following, x represents s, iw or ib). The model supports melting through both sublimation (solid to gas) and non-sublimation (solid to liquid) with the inclusion/exclusion of the latent heat of evaporation le (J kg -1 ). Non-sublimation melting is default with le set to zero, for sublimation melting the user can set le to 2265 20 kJ kg -1 . For the uppermost layer (y = top, Eq. B12) the heat flux includes layer dependent uptake of short wave radiation Hs_x, To test the ice module, Simstrat was calibrated in Sihlsee with PEST using monthly resolved vertical temperature profiles (2006 to 2008, RMSE 1.2 °C) for four parameters including the new p_albedo parameter for scaling snow/ice albedo. Modelled and monthly measured total ice cover from 2012 to 2018 is shown in Fig. B1 (RMSE 0.078 m). The modelled thickness agrees well with measurements during years with an extensive ice covered period (2013, 2014 and 2017, max height > 5 cm). The model performance is not ideal for years with short temporal ice duration and thin ice thickness (2016 and 2018, max height 5 < 5 cm). During these years, the quality of the forcing dataset becomes crucial. In the case of Sihlsee, the timing and duration of snowfall prolongs the duration of the ice-covered period. We use the meteorological station at Samedan (SAM) located four kilometres from the lake in a region with rapid topographical change. This in combination with monthly ice thickness measurements result in the divergence during 2016 and 2018.
10 Figure B1. Ice model performance in Sihlsee (2012Sihlsee ( to 2018 showing modelled white ice (orange), black ice (green) and total ice cover (white-and black ice combined, in blue) against measurements (black).

C. Estimation of clear-sky solar radiation
The algorithm below is based on the equations from the Lake HeatFluxAnalyzer (see http://heatfluxanalyzer.gleon.org/), following the methods of Meyers and Dale (1983). Clear-sky solar radiation [W m -2 ]: cs = eff • cos •