Toward an open access to high-frequency lake modeling and statistics data for scientists and practitioners – the case of Swiss lakes using Simstrat v2.1
- Surface Waters Research and Management, Eawag, Swiss Federal Institute of Aquatic Sciences and Technology, Kastanienbaum, Switzerland
- deceased, 2019
Correspondence: Damien Bouffard (email@example.com)
One-dimensional hydrodynamic models are nowadays widely recognized as key tools for lake studies. They offer the possibility to analyze processes at high frequency, here referring to hourly timescales, 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 (http://simstrat.eawag.ch, last access: 29 August 2019). 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 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 throughflow 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 this open-access data platform, we demonstrate a new path forward for scientists and practitioners promoting a cross exchange of expertise through openly sharing in situ and model data.
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 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., 2018). Today's most widely referenced one-dimensional (1-D) models include (in alphabetical order) DYRESM (Antenucci and Imerito, 2000), FLake (Mironov, 2005), General Lake Model (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.
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 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 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 ( km3) to very large systems such as Lake Geneva (89 km3), over an altitudinal gradient (Lake Maggiore at from 193 m a.s.l. to Daubensee at 2207 m a.s.l.) and over all trophic states (14 eutrophic 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 updates the online platform weekly (https://simstrat.eawag.ch, last access: 29 August 2019) 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.
2.1 Model and workflow
We use the 1-D lake model Simstrat v2.1 to model 54 Swiss lakes or reservoirs (see Appendix A for details of modeled 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 https://github.com/Eawag-AppliedSystemAnalysis/Simstrat/releases/tag/v2.1 (last access: 29 August 2019). 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 (ii) an ice and snow module. The ice and snow module employed in the model is based on the work of Leppäranta (2014, 2010) and Saloranta and Andersen (2007), and is further described in Appendix B.
2.2 Input data
Table 1 summarizes the type and sources of the data fed to Simstrat. For meteorological forcing, homogenized hourly air temperature, wind speed and direction, solar radiation, and relative humidity from the Federal Office of Meteorology and Climatology (MeteoSwiss, Switzerland) 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 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 watershed data are a source of error in small alpine lakes, yet such error can be compensated during the calibration process. The light absorption coefficient εabs (m−1) is either obtained from Secchi depth zSecchi (m) measurements (for Inkwilersee, Lake Biel, Lake Brienz, Lake Geneva, Lake Neuchâtel, lower Lake Zurich, Oeschinen Lake, upper Lake Constance, and Sihlsee) or set to a constant value based on the lake trophic status. In the first case, the following equation is applied: (Poole and Atkins, 1929; Schwefel et al., 2016). In the second case, εabs is set to 0.15 m−1 for oligotrophic lakes, 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 (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 1 d, the dataset is linearly interpolated. In total, < 1 % of the dataset is corrected using this approach. Longer data gaps of up to 20 d 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.
Model parameters are set to default values, and four of them are calibrated (see Table 2). The parameters p_radin and f_wind scale the incoming longwave 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 lake's morphology and its exposure 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, last access: 29 August 2019), a model-independent parameter estimation software (Doherty, 2016). As a reference for calibration, temperature observations from CTD profiles are used. Calibration is performed on a yearly basis, unless significant changes are made either to 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.
2.4 Output/available data on the online platform
The online platform (accessible at https://simstrat.eawag.ch, last access: 29 August 2019) is automatically fed every week with model results, metadata, and plots for all the 54 modeled lakes (see Fig. 2). It allows for efficient display and open sharing of the model results for 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 datasets 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 (g=9.81 m2 s−1), the heat capacity of water ( J K−1 kg−1), the volume of the lake V (m3), the area Az (m2), temperature Tz (∘C), and density ρz (kg m−3) at depth z (m), and the mean lake depth (m) to calculate time series of derived values:
mean lake temperature: (∘C);
heat content: (J);
Schmidt stability: (J m−2);
timing of summer stratification: we use a threshold based on the Schmidt stability to determine the beginning and end of summer stratification. The lake is assumed to be stratified for 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 https://plot.ly/python, last access: 29 August 2019). 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.
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, analyze, 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 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 eight lakes and between 2 and 3 ∘C for six lakes (Fig. 3). There were too few in situ observations on eight lakes to perform a proper calibration and all parameters were thereby set to default values. Overall, the performance is comparable to the RMSE range of ∼ 0.7–2.1 ∘C reported in a recent global 32-lake modeling 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 is often missing.
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 (Sect. 3.1), and secondly by investigating the differential response of lakes across Switzerland to episodic forcing (short-term extremes; Sect. 3.2).
3.1 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 (p<0.001) trend of +0.69 ∘C decade−1 (Fig. 4a). For the same period, monthly in situ observations indicate a similar trend of 0.72 ∘C decade−1 (p∼0.07), while the trend of air temperature at the meteorological station in Interlaken is lower (+0.50 ∘C decade−1, p<0.01). Based on physical principles, lake surface temperature is expected to increase less than air temperature (Schmid et al., 2014); 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 of 1981–2015, the ascending trend in solar radiation is 5 W m−2 decade−1, which corresponds to an equilibrium temperature increase of about 0.2 ∘C decade−1. The warming rate at the surface of Lake Brienz is larger than observed trends in neighboring lakes with reported increases of +0.46 ∘C decade−1 for upper Lake Constance (1984–2011; Fink et al., 2014), +0.41 ∘C decade−1 for lower Lake Zürich (1981–2013, Schmid and Köster, 2016; 1955–2013, Livingstone, 2003), and +0.55 ∘C decade−1 for lower Lake Lugano (1972–2013; Lepori and Roberts, 2015). This can be explained by the lower light penetration in Lake Brienz (ranging from ∼ 1 to ∼ 10 m) compared to other lakes, the increase in solar radiation being distributed into a shallower layer and thereby warming the lake surface slightly more. This low light penetration results from upstream hydropower operation on the glacier-fed river (Finger et al. 2006).
The temperature increase was significantly smaller in the hypolimnion, with a minimum trend at the lake bottom of 0.16 ∘C decade−1 (p<0.001), leading to a depth-averaged rate of temperature increase of 0.22 ∘C decade−1 (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 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−1. The riverine temperature remains colder than the lake surface temperature, leading to a yearly average loss of energy by throughflow of 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 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−1 (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.
The vertically heterogeneous warming modeled in Lake Brienz is consistent with previous observations showing that the difference in warming between the surface and the bottom increases the strength and duration of the stratified period (Zhong et al., 2016; Wahl and Peeters, 2014). We simulate an earlier onset of the stratification in spring of −7.5 d decade−1 (p<0.001) and a later breakdown of the stratification by +3.7 d decade−1 (p<0.001) (Fig. 4c). Both the warming trend and the increase in length of the stratified period increase the Schmidt stability (Fig. 4d) and heat content (Fig. 4f). Finally, the yearly maximum stratification strength (Brunt–Väisälä frequency; Fig. 4e) gradual increases over the investigated period with a rate of s−2 decade−1. The simulated increase in overall stability (Fig. 4d–f) reduces vertical mixing and affects the vertical storage of heat with less heat transferred immediately below the thermocline causing a slight decrease in temperature observed in autumn at ∼ 30 m depth (Fig. 4b). This effect is even more clearly seen in other lakes like Lake Geneva (https://simstrat.eawag.ch/LakeGeneva, last access: 29 August 2019) with the surface waters warming strongly (+1 ∘C decade−1 in June), resulting in a cooling layer between 20 and 60 m (−0.2 ∘C decade−1) in late summer. Such a reduction of vertical exchange is self-strengthening and enhances the differential vertical warming.
Such analyses can be extended to all modeled lakes. An intercomparison of the temporal extent of summer stratification and winter ice cover period is illustrated in Fig. 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. (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 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.
3.2 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. This resolution only allows to detect long-term trends when measurements are conducted over an extended period typically longer than 30 years. However, traditional monitoring programs cannot resolve the impact of short-term events and their consequences for the ecosystem. This is a strength of high-frequency (hourly timescale) lake modeling, 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 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 28 June 2018, when Switzerland experienced a strong but by no means exceptional storm with northeasterly winds mainly affecting the northwestern part of the country – the mean wind speed during that day is shown spatially in Fig. 6a. The evolution of the stratification strength, illustrated here by the Schmidt stability, is given in Fig. 6b for one of the most affected lakes, Lake Neuchâtel (https://simstrat.eawag.ch/LakeNeuchatel, last access: 29 August 2019; Fig. 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 ∼ 5 d (Fig. 6b). This also resulted in a total increase of the lake heat content by 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 modeled lakes. In Fig. 6a, lakes are colored 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 northwestern parts of Switzerland were more affected than large lakes or 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 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.
The workflow presented in this paper allows openly sharing 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. First, we showed that the high-frequency modeled temperature data allow 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 throughflow heat energy, and in light penetration to accurately assess the evolution of the lake thermal structure. Then, we showed that the high-frequency modeled 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 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 and (ii) simple statistical analyses such as monthly temperature profiles and long-term temperature trends.
For researchers, this work can facilitate (i) scenario modeling of any of the lakes, as the basic model setup is ready to 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 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 open-access data platform is a new path forward for scientists and practitioners.
The workflow was developed for Swiss lakes but can be easily extended to other geographical area or at global scale by using other meteorological input data. Simstrat and the Python workflow are available on https://github.com/Eawag-AppliedSystemAnalysis/Simstrat/releases/tag/v2.1 (last access: 29 August 2019, https://doi.org/10.5281/zenodo.2600709, Bärenbold et al., 2019) and https://github.com/Eawag-AppliedSystemAnalysis/Simstrat-WorkflowModellingSwissLakes (last access: 29 August 2019, https://doi.org/10.5281/zenodo.2607153, Gaudard, 2019). Meteorological data are available from MeteoSwiss (https://gate.meteoswiss.ch/idaweb/, last access: 29 August 2019), hydrological data are available from FOEN (https://www.hydrodaten.admin.ch, last access: 29 August 2019), and CTD data were provided by various sources listed here: https://simstrat.eawag.ch/impressum (last access: 29 August 2019). The calibration software PEST is available on http://www.pesthomepage.org/ (last access: 29 August 2019, Doherty, 2010).
The ice and snow module employed is based on the work of Leppäranta (2014, 2010) and Saloranta and Andersen (2007), and includes the following physical processes:
air-temperature-dependent formation and growth of black ice, including the insulating effect of a snow cover;
snow layer build-up, including the compression effect due to the weight of fresh snow;
buoyancy-driven formation of white ice;
shortwave irradiance reflection and penetration into the underlying water column; and
melting of snow, white and black ice due to both the direct heat flux through the atmospheric interface and the absorption of shortwave irradiance.
Three layers are used to represent black ice, white ice, and snow. An instant supply of water through cracks in the black ice is assumed to occur in order to form white ice. The water stored in ice and snow is neither withdrawn during ice formation nor added during melting to the water balance. Furthermore, the effect of liquid water pools on top of or between the layers is neglected.
B1 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 freshwater, 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×105 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 ice cover is present and if the atmospheric temperature Ta (∘C) is smaller than or equal to Tf, the growth of black ice dhib∕dt 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 h−1) is turned into fresh snow hs_new (m) as
where ρs0 (250 kg m−3) is the initial snow density. The existing snow cover hs (m) undergoes compression (first terms of Eqs. B7 and B8) by the new layer as described in Yen (1981); thereafter, the new and existing layers are combined in both height and density (second terms of Eqs. B7 and B8).
Here, ρs (kg m−3) is the snow layer density kept within with the maximum snow density set to 450 kg m−3, C1 (5.8 m−1 h−1) and C2 (0.021 m3 kg−1) are snow compression constants, and ws (m) is the total weight above the layer under compression expressed in water-equivalent height.
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 shortwave irradiance penetration.
B2 Above the freezing point (melting)
If ice cover is present and if Ta>Tf, melting starts. Each layer melts from above through the atmospheric interface and by penetrating shortwave radiation:
where Hx_y (W m−2) is the layer-dependent heat flux (in the following, subscript x represents the subscript letters “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 kJ kg−1. For the uppermost layer (y=top; Eq. B12), the heat flux includes layer-dependent uptake of shortwave radiation Hs_x, longwave absorption Ha, or layer-dependent emission Hw_x, as well as sensible Hk and latent Hv heat. If the layer is not in direct contact with the atmosphere, only Hs is used for melting from above (y=under; Eq. B13).
Here, we follow Leppäranta (2014, 2010) for determining the heat flux terms in Eq. (B12). The transmittance of shortwave irradiance through each layer depends on each layer's thickness hx as well as on the layer-specific bulk attenuation coefficient λx(m−1; default λs=24, λiw=3, and λib=2; Leppäranta, 2014).
There, Hs_w is the radiation penetrating through the ice cover to the water below and Is (W m−2) the incoming shortwave irradiance. We introduce the albedo parameter Ap, which tunes shortwave irradiance in order to match observed water temperatures, thus adjusting the melting and indirectly the duration of the ice cover. Furthermore, depending on which layer is in contact with the atmosphere, we use a layer-dependent constant albedo Ax (default As=0.7, Aiw=0.4, and Aib=0.3; Leppäranta, 2014).
Calculating Ha requires the longwave emission parameters ka=0.68, kb=0.036 (mbar−1), and kc=0.18 (Leppäranta, 2010), atmospheric water vapor pressure ea (mbar), cloud cover C, and the Stefan–Boltzmann constant W m−2 K−4). For Eqs. (B19) and (B20), the temperature Tx is given in Kelvin. Hw_x is layer dependent for the emissivity Ex with and Es(ρs) from 0.8 at ρs=250 kg m−3 to Es=0.9 for ρs=450 kg m−3. Calculating Hk and Hv requires the atmospheric density ρa=1.2 kg m−3, the heat capacity of air cpa=1005 J kg−1 K−1, the wind speed at 10 m height w10, the convective (bc) and latent (bl) bulk exchange coefficients both set to 0.0015 (Leppäranta, 2010; Gill, 1982), as well as the specific humidity of both measured qa (mbar) and at saturation q0. There, , where pa is the air pressure and at Ta=0 ∘C (Leppäranta, 2014).
As Hs_w warms the water under the ice, melting takes place from underneath with the energy Hbottom (W m−2):
After obtaining Hbottom, the temperature of the first cell is set to Tf and the decrease of ice cover from below becomes
Equation (B24) is only applied to hib and hiw. In principle, hs melts completely from above using Eq. (B11) before hib and hiw reach zero; however, if no ice is present, hs is set to zero. By combining Eqs. (B11) and (B24), the total melting of each ice layer is calculated as
When hx<0 due to melting, the surplus energy is used for melting neighboring layers according to the following procedure: if the melting is initiated from above, the surplus energy is used to melt the layer directly underneath; if the melting is caused by the water below, the layer directly above receives the surplus melting energy; if , the water in the topmost grid cell is heated with the remaining energy.
B3 Ice model performance
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. Modeled and monthly measured total ice cover from 2012 to 2018 is shown in Fig. B1 (RMSE of 0.078 m). The modeled 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 cm). During these years, the quality of the forcing dataset becomes crucial. In the case of Sihlsee, the timing and duration of snowfall prolong the duration of the ice-covered period. We use the meteorological station at Samedan (SAM) located 4 km from the lake in a region with rapid topographical change. This, in combination with monthly ice thickness measurements, results in the divergence during 2016 and 2018.
The algorithm below is based on the equations from the lake HeatFluxAnalyzer (see http://heatfluxanalyzer.gleon.org/, last access: 29 August 2019), following the methods of Meyers and Dale (1983):
declination of the Sun (rad): , where DOYs is the day of year after the winter solstice (21 December);
cosine of the solar zenith angle (–): , where φ is the latitude in radians and H is the hour of the day, assuming the solar noon is at 12:30 LT;
air mass thickness coefficient (–): m=35cos Z ;
dew point temperature (∘C): , where pw (mbar) is the water vapor pressure;
precipitable water vapor (cm): , where G is an empirical constant dependent on latitude and day of year (see tables from Smith, 1966);
attenuation coefficient for water vapor (–): ;
attenuation coefficient for aerosols (–): λa=0.935m;
attenuation coefficient for Rayleigh scattering and permanent gases (–): (m (0.000949pa+0.051))0.5, where pa (mbar) is the air pressure;
effective solar constant (W m−2): , where DOY is the day of year; and
clear-sky solar radiation (W m−2): .
The new version of Simstrat was developed by FB, AG, and LRV. The workflow was developed by AG. The ice model was 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.
The authors declare that they have no conflict of interest.
We thank Davide Vanzo for helping with the docker and the scripts, and Michael Pantic for helping restructuring version 2.1 of Simstrat. We finally thank Marie-Elodie Perga for her comments on a preliminary version of the paper. The full list of acknowledgements regarding in situ observations can be found here: https://simstrat.eawag.ch/impressum (last access: 29 August 2019).
This paper was edited by Min-Hui Lo and reviewed by three anonymous referees.
Antenucci, J. and Imerito, A.: The CWR dynamic reservoir simulation model DYRESM, Science Manual, The University of Western Australia, Perth, Australia, 2000.
Bärenbold, F., Gaudard, A., and Raman Vinna, L.: Simstrat v2.1 (Version v2.1), Zenodo, https://doi.org/10.5281/zenodo.2600709, 2019.
Bruce, L. C., Frassl, M. A., Arhonditsis, G. B., Gal, G., Hamilton, D. P., Hanson, P. C., Hetherington, A. L., Melack, J. M., Read, J. S., Rinke, K., Rigosi, A., Trolle, D., Winslow, L., Adrian, R., Ayala, A. I., Bocaniov, S. A., Boehrer, B., Boon, C., Brookes, J. D., Bueche, T., Busch, B. D., Copetti, D., Cortés, A., de Eyto, E., Elliott, J. A., Gallina, N., Gilboa, Y., Guyennon, N., Huang, L., Kerimoglu, O., Lenters, J. D., MacIntyre, S., Makler-Pick, V., McBride, C. G., Moreira, S., Özkundakci, D., Pilotti, M., Rueda, F. J., Rusak, J. A., Samal, N. R., Schmid, M., Shatwell, T., Snorthheim, C., Soulignac, F., Valerio, G., van der Linden, L., Vetter, M., Vinçon-Leite, B., Wang, J., Weber, M., Wickramaratne, C., Woolway, R. I., Yao, H., and Hipsey, M. R.: A multi-lake comparative analysis of the General Lake Model (GLM): Stress-testing across a global observatory network, Environ. Model. Softw., 102, 274–291, https://doi.org/10.1016/j.envsoft.2017.11.016, 2018.
Burchard, H., Bolding, K., and Villarreal, M. R.: GOTM, a general ocean turbulence model: theory, implementation and test cases, Space Applications Institute, 1999.
Carlson, R. E.: A trophic state index for lakes, Limnol. Oceanogr., 22, 361–369, 1977.
Doherty, J.: PEST, Model-independent parameter estimation – User manual (5th edn., with slight additions): Brisbane, Australia, Watermark Numerical Computing, available at: http://www.pesthomepage.org/ (last access: 29 August 2019), 2010.
Doherty, J.: PEST: Model-Independent Parameter Estimation, 6th edn., Watermark Numerical Computing, Australia, 2016.
Fink, G., Schmid, M., Wahl, B., Wolf, T., and Wüest, A.: Heat flux modifications related to climate-induced warming of large European lakes, Water Resour. Res., 50, 2072–2085, 2014.
Gaudard, A.: Simstrat-WorkflowModellingSwissLakes (Version v1.0), Zenodo, https://doi.org/10.5281/zenodo.2607153, 2019.
Gaudard, A., Schwefel, R., Vinnå, L. R., Schmid, M., Wüest, A., and Bouffard, D.: Optimizing the parameterization of deep mixing and internal seiches in one-dimensional hydrodynamic models: a case study with Simstrat v1.3, Geosci. Model Dev., 10, 3411–3423, https://doi.org/10.5194/gmd-10-3411-2017, 2017.
Gelman, A., Carlin, J., Stern, H., and Rubin, D.: Bayesian Data Analysis, 3rd edn., Chapman and Hall/CRC, New York, 2013.
Gill, A. E.: Atmosphere-Ocean Dynamics, Academic Press, San Diego, California, USA, ISBN 0-12-283520-4, 1982.
Goudsmit, G.-H., Burchard, H., Peeters, F., and Wüest, A.: Application of k−ϵ turbulence models to enclosed basins: The role of internal seiches, J. Geophys. Res.-Oceans, 107, 3230, https://doi.org/10.1029/2001JC000954, 2002.
Gray, D. K., Hampton, S. E., O'Reilly, C. M., Sharma, S., and Cohen, R. S.: How do data collection and processing methods impact the accuracy of long-term trend estimation in lake surface-water temperatures?, Limnol. Oceanogr.-Meth., 16, 504–515, https://doi.org/10.1002/lom3.10262, 2018.
Hamilton, D. P., Carey, C. C., Arvola, L., Arzberger, P., Brewer, C., Cole, J. J., Gaiser, E., Hanson, P. C., Ibelings, B. W., Jennings, E., Kratz, T. K., Lin, F.-P., McBride, C. G., Marques, M. D. de, Muraoka, K., Nishri, A., Qin, B., Read, J. S., Rose, K. C., Ryder, E., Weathers, K. C., Zhu, G., Trolle, D., and Brookes, J. D.: A Global Lake Ecological Observatory Network (GLEON) for synthesising high-frequency sensor data for validation of deterministic ecological models, Inland Waters, 5, 49–56, https://doi.org/10.5268/IW-5.1.566, 2015.
Hipsey, M. R., Bruce, L. C., and Hamilton, D. P.: GLM – General Lake Model. Model overview and user information, Technical Manual, The University of Western Australia, Perth, Australia, available at: http://swan.science.uwa.edu.au/downloads/AED_GLM_v2_0b0_20141025.pdf (last access: 29 August 2019), 2014.
Jennings, E., Eyto, E., Laas, A., Pierson, D., Mircheva, G., Naumoski, A., Clarke, A., Healy, M., Šumberová, K., and Langenhaun, D.: The NETLAKE Metadatabae – A Tool to Support Automatic Monitoring on Lakes in Europe and Beyond, Limnol. Oceanogr., 26, 95–100, https://doi.org/10.1002/lob.10210, 2017.
Kiefer, I., Odermatt, D., Anneville, O., Wüest, A., and Bouffard, D.: Application of remote sensing for the optimization of in-situ sampling for monitoring of phytoplankton abundance in a large lake, Sci. Total Environ., 527–528, 493–506, https://doi.org/10.1016/j.scitotenv.2015.05.011, 2015.
Lepori, F., and Roberts, J. J.: Past and future warming of a deep European lake (Lake Lugano): What are the climatic drivers?, J. Great Lakes Res., 41, 973–981, 2015.
Leppäranta, M.: Modelling the Formation and Decay of Lake Ice, in: The Impact of Climate Change on European Lakes, edited by: George, G., Springer Netherlands, Dordrecht, 63–83, 2010.
Leppäranta, M.: Freezing of lakes and the evolution of their ice cover, Springer, New York, 2014.
Livingstone, D. M.: Impact of secular climate change on the thermal structure of a large temperate central European lake, Clim. Change, 57, 205–225, 2003.
Livingstone, D. M., Lotter, A. F., and Kettle, H.: Altitude-dependent differences in the primary physical response of mountain lakes to climatic forcing, Limnol. Oceanogr., 50, 1313–1325, 2005.
Meyers, T. P. and Dale, R. F.: Predicting Daily Insolation with Hourly Cloud Height and Coverage, J. Climate Appl. Meteor., 22, 537–545, https://doi.org/10.1175/1520-0450(1983)022< 0537:PDIWHC> 2.0.CO;2, 1983.
Mironov, D. V.: Parameterization of lakes in numerical weather prediction – Part 1: Description of a lake model, German Weather Service, Offenbach am Main, Germany, 2005.
Perga, M.-E., Bruel, R., Rodriguez, L., Guénand, Y., and Bouffard, D.: Storm impacts on alpine lakes: Antecedent weather conditions matter more than the event intensity, Glob. Chang. Biol., 24, 5004–5016, https://doi.org/10.1111/gcb.14384, 2018.
Perroud, M., Goyette, S., Martynov, A., Beniston, M., and Anneville, O.: Simulation of multiannual thermal profiles in deep Lake Geneva: A comparison of one-dimensional lake models, Limnol. Oceanogr., 54, 1574–1594, 2009.
Poole, H. H. and Atkins, W. R. G.: Photo-electric measurements of submarine illumination throughout the year, J. Mar. Biol. Assoc. UK, 16, 297–324, 1929.
Råman Vinnå, L., Wüest, A., Zappa, M., Fink, G., and Bouffard, D.: Tributaries affect the thermal response of lakes to climate change, Hydrol. Earth Syst. Sci., 22, 31–51, https://doi.org/10.5194/hess-22-31-2018, 2018.
Riley, M. J. and Stefan, H. G.: Minlake: A dynamic lake water quality simulation model, Ecol. Model., 43, 155–182, https://doi.org/10.1016/0304-3800(88)90002-6, 1988.
Sadro, S., Sickman, J. O., Melack, J. M., and Skeen, K.: Effects of Climate Variability on Snowmelt and Implications for Organic Matter in a High-Elevation Lake, Water Resour. Res., 54, 4563–4578, https://doi.org/10.1029/2017WR022163, 2018.
Saloranta, T. M.: Modeling the evolution of snow, snow ice and ice in the Baltic Sea, Tellus A, 52, 93–108, https://doi.org/10.1034/j.1600-0870.2000.520107.x, 2000.
Saloranta, T. M. and Andersen, T.: MyLake – A multi-year lake simulation model code suitable for uncertainty and sensitivity analysis simulations, Ecol. Model., 207, 45–60, 2007.
Schmid, M., Hunziker, S., and Wüest, A.: Lake surface temperatures in a changing climate: a global perspective, Clim. Change, 124, 301–305, 2014.
Schmid, M. and Köster, O.: Excess warming of a Central European lake driven by solar brightening, Water Resour. Res., 52, 8103–8116, https://doi.org/10.1002/2016WR018651, 2016.
Schwefel, R., Gaudard, A., Wüest, A., and Bouffard, D.: Effects of climate change on deep-water oxygen and winter mixing in a deep lake (Lake Geneva)–Comparing observational findings and modeling, Water Resour. Res., 52, 8811–8826, https://doi.org/10.1002/2016WR019194, 2016.
Smith, W. L.: Note on the Relationship Between Total Precipitable Water and Surface Dew Point, J. Appl. Meteor., 5, 726–727, https://doi.org/10.1175/1520-0450(1966)005< 0726:NOTRBT> 2.0.CO;2, 1966.
Stepanenko, V., 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, 1977–2006, https://doi.org/10.5194/gmd-9-1977-2016, 2016.
Thiery, W., Stepanenko, V. M., Fang, X., Jöhnk, K. D., Li, Z., Martynov, A., Perroud, M., Subin, Z. M., Darchambeau, F., Mironov, D. V., and Van Lipzig, N. P. M.: LakeMIP Kivu: evaluating the representation of a large, deep tropical lake by a set of one-dimensional lake models, Tellus A, 66, 21390, https://doi.org/10.3402/tellusa.v66.21390, 2014.
Wahl, B. and Peeters, F.: Effect of climatic changes on stratification and deep-water renewal in Lake Constance assessed by sensitivity studies with a 3D hydrodynamic model, Limnol. Oceanogr., 59, 1035–1052, https://doi.org/10.4319/lo.2014.59.3.1035, 2014.
Woolway, R. I., Simpson, J. H., Spiby, D., Feuchtmayr, H., Powell, B., and Maberly, S. C.: Physical and chemical impacts of a major storm on a temperate lake: a taste of things to come?, Clim. Change, 151, 333–347, https://doi.org/10.1007/s10584-018-2302-3, 2018.
Yen, Y. C.: Review of thermal properties of snow, ice and sea ice, Cold Regions Research and Engineering Laboratory, Hanover, New Hampshire, 1981.
Zhong, Y., Notaro, M., Vavrus, S. J., and Foster, M. J.: Recent accelerated warming of the Laurentian Great Lakes: Physical drivers, Limnol. Oceanogr., 61, 1762–1786, https://doi.org/10.1002/lno.10331, 2016.