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

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.

1 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 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 ($\mathrm{9}\times {\mathrm{10}}^{-\mathrm{3}}$ km³) to very large systems such as Lake Geneva (89 km³), 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 Methods

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.

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 online data platform to display these results. The script is controlled by an input file written in JavaScript Object Notation (JSON) format, which specifies the lakes to be modeled 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 Fig. 1.

Figure 1 General workflow diagram. Model input (a) is retrieved and processed by the Python script “Simstrat.py”, which runs the model (Simstrat v2.1) and/or model calibration (using PEST v15.0) (b) and produces output (c). This output is then uploaded to a web interface (https://simstrat.eawag.ch, last access: 29 August 2019) for general use. All scripts and programs are available on https://github.com/Eawag-AppliedSystemAnalysis/Simstrat/releases/tag/v2.1 (last access: 29 August 2019) and https://github.com/Eawag-AppliedSystemAnalysis/Simstrat-WorkflowModellingSwissLakes (last access: 29 August 2019). Simstrat is the one-dimensional hydrodynamic model; CTD is a conductivity–temperature–depth profiler; PEST is the model-independent parameter estimation and uncertainty analysis software; FOEN is the Swiss Federal Office of Environment; MeteoSwiss is the Swiss Federal Office of Meteorology and Climatology; Swisstopo is the Swiss Federal Office of Topography.

Table 1 Input data sources used for the model.

Table 2 Model parameters. The geothermal heat flux is based on existing geothermal data for Switzerland: https://www.geocat.ch/geonetwork/srv/eng/md.viewer\#/full_view/2d8174b2-8c4a-44ea-b470-cb3f216b90d1 (last access: 29 August 2019).

The asterisk (^*) indicates the parameters that were calibrated.

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⁻¹. 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⁻¹) is either obtained from Secchi depth z_Secchi (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: ${\mathit{\epsilon }}_{\mathrm{abs}}=\mathrm{1.7}/z_{\mathrm{Secchi}}$ (Poole and Atkins, 1929; Schwefel et al., 2016). In the second case, ε_abs is set to 0.15 m⁻¹ for oligotrophic lakes, 0.25 m⁻¹ for mesotrophic lakes, and 0.50 m⁻¹ 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⁻¹.

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.

2.3 Calibration

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.

Figure 2 Illustration of the interactive map displayed on the home page of the online platform: https://simstrat.eawag.ch (last access: 29 August 2019). The locations of the lakes discussed in this paper are also indicated with numbers (see Appendix A). Basemap source: Federal Office of Topography ©Swisstopo.

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 m² s⁻¹), the heat capacity of water ($c_{\mathrm{p}}=\mathrm{4.18}\times {\mathrm{10}}^{\mathrm{3}}$ J K⁻¹ kg⁻¹), the volume of the lake V (m³), the area A_z (m²), temperature T_z (^∘C), and density ρ_z (kg m⁻³) at depth z (m), and the mean lake depth $\stackrel{\mathrm{‾}}{z}=\frac{\mathrm{1}}{V}\int z{A}_z\text{d}z$ (m) to calculate time series of derived values:

• mean lake temperature: $\stackrel{\mathrm{‾}}{T}=\frac{\mathrm{1}}{V}\int T_z{A}_z\text{d}z$ (^∘C);

• heat content: $H=c_{\mathrm{p}}\int {\mathit{\rho }}_z{T}_z{A}_z\text{d}z$ (J);

• Schmidt stability: $S_T=\frac{g}{A_{\mathrm{0}}}\int \left(z-\stackrel{\mathrm{‾}}{z}\right){\mathit{\rho }}_z{A}_z\text{d}z$ (J m⁻²);

• 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 $S_T/z_{\text{lake}}\ge \mathrm{10}$ J m⁻³. 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.

Figure 3 Performance of the model for the different lakes, as shown by the root mean square error (RMSE) and the correlation coefficient. Six lakes (with symbol • on the legend) with RMSE > 2 ^∘C are not shown.

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.

3 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, 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).

Figure 4 Evolution of several indicators for Lake Brienz over the period 1981–2018; all linear regression have p values ≪ 0.001: (a) yearly mean lake surface temperature (0.69 ^∘C decade⁻¹), yearly mean air temperatures (0.49 ^∘C decade⁻¹), yearly mean tributary temperatures (0.26 ^∘C decade⁻¹), yearly mean lake temperatures (0.22 ^∘C decade⁻¹), and yearly mean bottom temperatures (0.16 ^∘C decade⁻¹), with linear regression, (b) contour plot of the linear temperature trend through depth and month, (c) yearly start (+3.7 d decade⁻¹) and end (−7.5 d decade⁻¹) days of summer stratification, with linear regression, (d) yearly mean (line), min, and max (shaded area) Schmidt stability, with linear regression, (e) yearly maximum Brunt–Väisälä frequency ($\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻² decade⁻¹), with linear regression, (f) yearly mean (line), min, and max (shaded area) heat content.

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⁻¹ (Fig. 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 (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⁻² decade⁻¹, which 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 neighboring 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, 2016; 1955–2013, Livingstone, 2003), and +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 ∼ 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).

Figure 5 Comparison of timing of stratification and ice cover for the considered lakes. The colored areas represent the mean periods of summer stratification (red) and ice cover (blue); the vertical lines represent the last year (here 2017). The transparency for the ice cover indicates the freezing frequency: full transparency means that ice was never modeled, while no transparency means that ice was modeled every winter. Lakes are ordered from left (low elevation) to right (high elevation). The time period of data used is indicated in Appendix A.

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 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 throughflow of $\sim -\mathrm{40}$ W m⁻² 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⁻¹ (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⁻¹ (p<0.001) and a later breakdown of the stratification by +3.7 d decade⁻¹ (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 $\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻² decade⁻¹. 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⁻¹ in June), resulting in a cooling layer between 20 and 60 m (−0.2 ^∘C decade⁻¹) 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.

Figure 6 (a) Mean wind field on 28 June 2018 (data source: MeteoSwiss, COSMO-1 model, coordinate system CH1903+) and delay in Schmidt stability increase for the modeled lakes: from no delay (white) to a delay of more than 5 d (red). (b) Schmidt stability (daily average) in Lake Neuchâtel during the period of the storm.

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 $\sim \mathrm{1.4}\times {\mathrm{10}}^{\mathrm{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 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.

4 Conclusion

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.

Appendix A: Properties of the modeled lakes

Table A1 This table summarizes the main properties of the 54 lakes we model in this work. The full dataset is available as a JSON file. The superscript “a” after the lake name indicates that this lake was not calibrated due to the lack of observational data. MeteoSwiss is the (Swiss) Federal Office of Meteorology and Climatology. FOEN is the (Swiss) Federal Office for the Environment. The superscript “b” indicates lakes where Secchi disk depths are available. For lakes with clearly defined multiple basins such as Lake Lucerne, Lake Zurich, Lake Constance and Lake Lugano, each basin is considered as a separated lake connected to the other basins by inflows/outflows.

Table A2 Meteorological stations.

Appendix B: Ice module

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 T_w (^∘C) drops below the freezing temperature T_f (^∘C). T_f 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 E_f is calculated as

$$\begin{array}{}\text{(B1)}& E_{\text{f}}={\mathit{\rho }}_{\text{w}}{c}_{\text{pw}}{z}_{\mathrm{1}}(T_{\text{f}}-T_{\text{w}}).\end{array}$$

Here, ρ_w (1000 kg m⁻³) is the density of freshwater, c_pw the heat capacity of water (4182 J kg^−1 ∘C⁻¹), and z₁ the height of the topmost grid cell. E_f and the latent heat of freezing l_h (3.34×10⁵ J kg⁻¹) as well as the density of black ice ρ_ib (916.2 kg m⁻³) are used for calculating the initial height of black ice h_ib (m) in Eq. (B2); thereafter, T_w is set equal to T_f:

$$\begin{array}{}\text{(B2)}& h_{\text{ib}}=E_{\text{f}}/\left(l_{\text{h}}{\mathit{\rho }}_{\text{ib}}\right).\end{array}$$

If ice cover is present and if the atmospheric temperature T_a (^∘C) is smaller than or equal to T_f, the growth of black ice dh_ib∕dt continues as described in Saloranta and Andersen (2007).

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\text{d}{h}_{\text{ib}}}{\text{d}t}}=\sqrt{\mathrm{2}{k}_{\text{i}}/\left({\mathit{\rho }}_{\text{ib}}\times l_{\text{h}}\right)\times \left(T_{\text{f}}-T_{\text{i}}\right)}\end{array}$$

Here, k_i (2.22 W K⁻¹ m⁻¹) is the thermal conductivity of ice at 0 ^∘C and T_i (^∘C) the ice temperature calculated as

$$\begin{array}{}\text{(B4)}& {\displaystyle }& {\displaystyle }{T}_{\text{i}}={\displaystyle \frac{P{T}_{\text{f}}+T_{\text{a}}}{\mathrm{1}+P}}\text{(B5)}& {\displaystyle }& {\displaystyle }P=max\left({\displaystyle \frac{k_{\text{i}}{h}_{\text{s}}}{k_{\text{s}}{h}_{\text{ib}}}},{\displaystyle \frac{\mathrm{1}}{\mathrm{10}{h}_{\text{ib}}}}\right).\end{array}$$

There, k_s (0.2 W K⁻¹ m⁻¹) is the thermal conductivity of snow and h_s (m) the height of the snow layer. When T_a is smaller than the snow temperature (default set to 2 ^∘C), water-equivalent precipitation p_r (m h⁻¹) is turned into fresh snow h_{s_new} (m) as

$$\begin{array}{}\text{(B6)}& h_{\text{s\_new}}=p_{\text{r}}{\displaystyle \frac{{\mathit{\rho }}_{\text{w}}}{{\mathit{\rho }}_{\text{s0}}}},\end{array}$$

where ρ_s0 (250 kg m⁻³) is the initial snow density. The existing snow cover h_s (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).

$$\begin{array}{}\text{(B7)}& {\displaystyle }& {\displaystyle \frac{\text{d}{h}_{\text{s}}}{\text{d}t}}=-h_{\text{s}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\text{s}}}{\left[{\mathit{\rho }}_{\text{s}}+\text{d}{\mathit{\rho }}_{\text{s}}\right]}}\right)+h_{\text{s\_new}}\text{(B8)}& {\displaystyle }& {\displaystyle \frac{\text{d}{\mathit{\rho }}_{\text{s}}}{\text{d}t}}={\mathit{\rho }}_{\text{s}}{C}_{\mathrm{1}}{w}_{\text{s}}{e}^{-C_{\mathrm{2}}{\mathit{\rho }}_{\text{s}}}-\left({\mathit{\rho }}_{\text{s}}-{\displaystyle \frac{{\mathit{\rho }}_{\text{s}}{h}_{\text{s}}+{\mathit{\rho }}_{\text{s0}}{h}_{\text{s\_new}}}{h_{\text{s}}+h_{\text{s\_new}}}}\right)\end{array}$$

Here, ρ_s (kg m⁻³) is the snow layer density kept within ${\mathit{\rho }}_{\text{s0}}<{\mathit{\rho }}_{\text{s}}<{\mathit{\rho }}_{\text{sm}}$ with the maximum snow density set to 450 kg m⁻³, C₁ (5.8 m⁻¹ h⁻¹) and C₂ (0.021 m³ kg⁻¹) are snow compression constants, and w_s (m) is the total weight above the layer under compression expressed in water-equivalent height.

If the snow mass m_s (kg m⁻²) becomes heavier than the upward acting buoyancy force B_i(kg m⁻²), white ice with height h_iw (m) and density ρ_iw (875 kg m⁻³; Saloranta, 2000) is formed between the snow and the black ice layers to achieve equilibrium between B_i and m_s.

$$\begin{array}{}\text{(B9)}& {\displaystyle }& {\displaystyle }{B}_{\text{i}}=h_{\text{ib}}\left({\mathit{\rho }}_{\text{w}}-{\mathit{\rho }}_{\text{ib}}\right)+h_{\text{iw}}\left({\mathit{\rho }}_{\text{w}}-{\mathit{\rho }}_{\text{iw}}\right)\text{(B10)}& {\displaystyle }& {\displaystyle \frac{\text{d}{h}_{\text{iw}}}{\text{d}t}}={\displaystyle \frac{m_{\text{s}}-B_{\text{i}}}{{\mathit{\rho }}_{\text{s}}}}\end{array}$$

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 T_a>T_f, melting starts. Each layer melts from above through the atmospheric interface and by penetrating shortwave radiation:

$$\begin{array}{}\text{(B11)}& {\displaystyle \frac{\text{d}{h}_{x\mathit{\_}\text{upper}}}{\text{d}t}}=-{\displaystyle \frac{H_{x\mathit{\_}y}}{\left(l_{\text{h}}+l_{\text{e}}\right){\mathit{\rho }}_x}},\end{array}$$

where H_{x_y} (W m⁻²) 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 l_e (J kg⁻¹). Non-sublimation melting is default with l_e set to zero; for sublimation melting, the user can set l_e to 2265 kJ kg⁻¹. For the uppermost layer (y=top; Eq. B12), the heat flux includes layer-dependent uptake of shortwave radiation H_{s_x}, longwave absorption H_a, or layer-dependent emission H_{w_x}, as well as sensible H_k and latent H_v heat. If the layer is not in direct contact with the atmosphere, only H_s is used for melting from above (y=under; Eq. B13).

$$\begin{array}{}\text{(B12)}& {\displaystyle }& {\displaystyle }{H}_{x\mathit{\_}\text{top}}=H_{\text{s}\mathit{\_}x}+H_{\text{a}}+H_{\text{w}\mathit{\_}x}+H_{\text{k}}+H_{\text{v}}\text{(B13)}& {\displaystyle }& {\displaystyle }{H}_{x\mathit{\_}\text{under}}=H_{\text{s}\mathit{\_}x}\end{array}$$

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 h_x as well as on the layer-specific bulk attenuation coefficient λ_x(m⁻¹; default λ_s=24, λ_iw=3, and λ_ib=2; Leppäranta, 2014).

$$\begin{array}{}\text{(B14)}& {\displaystyle }{H}_{\text{s\_s}}=I_{\text{s}}{A}_{\text{p}}\left(\mathrm{1}-A_x\right)\left(\mathrm{1}-e^{\left(-{\mathit{\lambda }}_{\text{s}}{h}_{\text{s}}\right)}\right)\text{(B15)}& {\displaystyle }{H}_{\text{s\_iw}}=I_{\text{s}}{A}_{\text{p}}\left(\mathrm{1}-A_x\right)\left(e^{\left(-{\mathit{\lambda }}_{\text{s}}{h}_{\text{s}}\right)}-e^{\left(-{\mathit{\lambda }}_{\text{s}}{h}_{\text{s}}-{\mathit{\lambda }}_{\text{iw}}{h}_{\text{iw}}\right)}\right)\text{(B16)}& {\displaystyle }\begin{array}{rl}{H}_{\text{s\_ib}}=& I_{\text{s}}{A}_{\text{p}}\left(\mathrm{1}-A_x\right)\left(e^{\left(-{\mathit{\lambda }}_{\text{s}}{h}_{\text{s}}-{\mathit{\lambda }}_{\text{iw}}{h}_{\text{iw}}\right)}\right.\\ & \left.-e^{\left(-{\mathit{\lambda }}_{\text{s}}{h}_{\text{s}}-{\mathit{\lambda }}_{\text{iw}}{h}_{\text{iw}}-{\mathit{\lambda }}_{\text{ib}}{h}_{\text{ib}}\right)}\right)\end{array}\text{(B17)}& {\displaystyle }{H}_{\text{s\_w}}=I_{\text{s}}{A}_{\text{p}}\left(\mathrm{1}-A_x\right)\left(\mathrm{1}-e^{\left(-{\mathit{\lambda }}_{\text{s}}{h}_{\text{s}}-{\mathit{\lambda }}_{\text{iw}}{h}_{\text{iw}}-{\mathit{\lambda }}_{\text{ib}}{h}_{\text{ib}}\right)}\right).\end{array}$$

There, H_{s_w} is the radiation penetrating through the ice cover to the water below and I_s (W m⁻²) the incoming shortwave irradiance. We introduce the albedo parameter A_p, 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 A_x (default A_s=0.7, A_iw=0.4, and A_ib=0.3; Leppäranta, 2014).

$$\begin{array}{}\text{(B18)}& {\displaystyle }{A}_x\left\{\begin{array}{l}{A}_{\text{s}},h_{\text{s}}>\mathrm{0}\\ A_{\text{iw}},h_{\text{s}}=\mathrm{0}\mathit{\&}{h}_{\text{iw}}>\mathrm{0}\\ A_{\text{ib}},h_{\text{s}}+h_{\text{iw}}=\mathrm{0}\end{array}\right.\end{array}$$

Calculating H_a requires the longwave emission parameters k_a=0.68, k_b=0.036 (mbar⁻¹), and k_c=0.18 (Leppäranta, 2010), atmospheric water vapor pressure e_a (mbar), cloud cover C, and the Stefan–Boltzmann constant $\mathit{\sigma }(\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴). For Eqs. (B19) and (B20), the temperature T_x is given in Kelvin. H_{w_x} is layer dependent for the emissivity E_x with $E_{\text{iw}}=E_{\text{ib}}=\mathrm{0.97}$ and E_s(ρ_s) from 0.8 at ρ_s=250 kg m⁻³ to E_s=0.9 for ρ_s=450 kg m⁻³. Calculating H_k and H_v requires the atmospheric density ρ_a=1.2 kg m⁻³, the heat capacity of air c_pa=1005 J kg⁻¹ K⁻¹, the wind speed at 10 m height w₁₀, the convective (b_c) and latent (b_l) bulk exchange coefficients both set to 0.0015 (Leppäranta, 2010; Gill, 1982), as well as the specific humidity of both measured q_a (mbar) and at saturation q₀. There, $q_{\text{a}}=\mathrm{0.622}{e}_{\text{a}}/p_{\text{a}}$, where p_a is the air pressure and $q_{\mathrm{0}}=\mathrm{0.622}\times \mathrm{6.11}/p_{\text{a}}$ at T_a=0 ^∘C (Leppäranta, 2014).

$$\begin{array}{}\text{(B19)}& {\displaystyle }& {\displaystyle }{H}_{\text{a}}=\left(k_{\text{a}}+k_{\text{b}}\sqrt{e_{\text{a}}}\left[\mathrm{1}+k_{\text{c}}{C}^{\mathrm{2}}\right]\right)\mathit{\sigma }{T}_{\text{a}}^{\mathrm{4}}\text{(B20)}& {\displaystyle }& {\displaystyle }{H}_{\text{w}\mathit{\_}x}=E_x\mathit{\sigma }{T}_{\text{f}}^{\mathrm{4}}\text{(B21)}& {\displaystyle }& {\displaystyle }{H}_{\text{k}}={\mathit{\rho }}_{\text{a}}{c}_{\text{pa}}{b}_{\text{c}}\left(T_{\text{a}}-T_{\text{f}}\right)w_{\mathrm{10}}\text{(B22)}& {\displaystyle }& {\displaystyle }{H}_{\text{v}}={\mathit{\rho }}_{\text{a}}{l}_{\text{h}}{b}_{\text{l}}\left(q_{\text{a}}-q_{\mathrm{0}}\right)w_{\mathrm{10}}\end{array}$$

As H_{s_w} warms the water under the ice, melting takes place from underneath with the energy H_bottom (W m⁻²):

$$\begin{array}{}\text{(B23)}& H_{\text{bottom}}=\left(T_{\text{w}}-T_{\text{f}}\right)c_{\text{pw}}{\mathit{\rho }}_{\text{w}}{z}_{\mathrm{1}}.\end{array}$$

After obtaining H_bottom, the temperature of the first cell is set to T_f and the decrease of ice cover from below becomes

$$\begin{array}{}\text{(B24)}& {\displaystyle \frac{\text{d}{h}_{x\mathit{\_}\text{lower}}}{\text{d}t}}=-{\displaystyle \frac{H_{\text{bottom}}}{l_{\text{m}}{\mathit{\rho }}_x}}.\end{array}$$

Equation (B24) is only applied to h_ib and h_iw. In principle, h_s melts completely from above using Eq. (B11) before h_ib and h_iw reach zero; however, if no ice is present, h_s is set to zero. By combining Eqs. (B11) and (B24), the total melting of each ice layer is calculated as

$$\begin{array}{}\text{(B25)}& {\displaystyle \frac{\text{d}{h}_x}{\text{d}t}}={\displaystyle \frac{\text{d}{h}_{x\mathit{\_}\text{lower}}}{\text{d}t}}+{\displaystyle \frac{\text{d}{h}_{x\mathit{\_}\text{upper}}}{\text{d}t}}.\end{array}$$

When h_x<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 $h_{\text{ib}}<=h_{\text{iw}}<=\mathrm{0}$, 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.

Figure B1 Ice model performance in Sihlsee (2012 to 2018) showing modeled white ice (orange), black ice (green), and total ice cover (white and black ice combined, in blue) against measurements (black).

Appendix C: Estimation of clear-sky solar radiation

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): $\mathit{\delta }={\sin }^{-\mathrm{1}}$ $\left(-\mathrm{0.39779}\cos \frac{\mathrm{2}\mathit{\pi }{\mathrm{DOY}}_{\mathrm{s}}}{\mathrm{365.24}}\right)$, where DOY_s is the day of year after the winter solstice (21 December);

• cosine of the solar zenith angle (–): $\cos Z=max\left(\sin \mathit{\phi }\sin \mathit{\delta }+\cos \mathit{\phi }\cos \mathit{\delta }\cos \frac{\mathit{\pi }\left(\mathrm{H}h-\mathrm{12.5}\right)}{\mathrm{12}},\mathrm{0}\right)$, 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 $(\mathrm{1244}{\cos }^{\mathrm{2}}Z+\mathrm{1}{)}^{-\mathrm{0.5}}$;

• dew point temperature (^∘C): $T_{\text{d}}=\mathrm{243.5}\log \frac{p_{\text{w}}}{\mathrm{6.112}}/$ $\left(\mathrm{17.67}-\log \frac{p_{\text{w}}}{\mathrm{6.112}}\right)+\mathrm{33.8}$, where p_w (mbar) is the water vapor pressure;

• precipitable water vapor (cm): $w_{\text{p}}=e^{\mathrm{0.1133}-\log \left(G+\mathrm{1}\right)+\mathrm{0.0393}(\mathrm{1.8}{T}_{\text{d}}+\mathrm{32})}$, where G is an empirical constant dependent on latitude and day of year (see tables from Smith, 1966);

• attenuation coefficient for water vapor (–): ${\mathit{\lambda }}_{\text{w}}=\mathrm{1}-\mathrm{0.077}(w_{\text{p}}m{)}^{\mathrm{0.3}}$;

• attenuation coefficient for aerosols (–): λ_a=0.935^m;

• attenuation coefficient for Rayleigh scattering and permanent gases (–): ${\mathit{\lambda }}_{\text{Rg}}=\mathrm{1.021}-\mathrm{0.084}$ (m (0.000949p_a+0.051))^0.5, where p_a (mbar) is the air pressure;

• effective solar constant (W m⁻²): $I_{\text{eff}}=\mathrm{1353}(\mathrm{1}+$ $\left.\mathrm{0.034}\cos \frac{\mathrm{2}\mathit{\pi }\mathrm{DOY}}{\mathrm{365.24}}\right)$, where DOY is the day of year; and

• clear-sky solar radiation (W m⁻²): $H_{\mathrm{cs}}=I_{\text{eff}}\times \cos Z\times {\mathit{\lambda }}_{\text{w}}{\mathit{\lambda }}_{\text{a}}{\mathit{\lambda }}_{\text{Rg}}$.