Articles | Volume 15, issue 22
Model description paper
18 Nov 2022
Model description paper |  | 18 Nov 2022

SIMO v1.0: simplified model of the vertical temperature profile in a small, warm, monomictic lake

Kristina Šarović, Melita Burić, and Zvjezdana B. Klaić

A simple 1-D energy budget model (SIMO) for the prediction of the vertical temperature profiles in small, monomictic lakes forced by a reduced number of input meteorological variables is proposed. The model estimates the net heat flux and thermal diffusion using only routinely measured hourly mean meteorological variables (namely, the air temperature, relative humidity, atmospheric pressure, wind speed, and precipitation), hourly mean ultraviolet B radiation (UVB), and climatological yearly mean temperature data. Except for the initial vertical temperature profile, the model does not use any lake-specific variables. The model performance was evaluated against lake temperatures measured continuously during an observational campaign in two lakes belonging to the Plitvice Lakes, Croatia (Lake 1 and Lake 12). Temperatures were measured at 15 and 16 depths ranging from 0.2 to 27 m in Lake 1 (maximum depth of 37.4 m) and 0.2 to 43 m in Lake 12 (maximum depth of 46 m). The model performance was evaluated for simulation lengths from 1 to 30 d. The model performed reasonably well, and it was able to satisfactorily reproduce the vertical temperature profile at the hourly scale, the deepening of the thermocline with time, and the annual variation in the vertical temperature profile, which shows its applicability for short-term prognostic simulations. A yearlong simulation initiated with an approximately constant vertical profile of the lake temperature ( 4 C) was able to reproduce the onset of stratification and convective overturn. The epilimnion temperature was somewhat overestimated, especially with the onset of the convective overturn. The upper limit of the metalimnion was well captured, while its thickness was overestimated. Nevertheless, the values of the model performance measures obtained for a yearlong simulation were comparable with those reported for other, more complex models. Thus, the presented model can also be used for long-term simulations and the assessment of the onset and duration of lake stratification periods when water temperature data are unavailable, which can be useful for various lake studies performed in other scientific fields, such as biology, geochemistry, and sedimentology.

1 Introduction

Water temperature is a critical factor that directly influences a whole range of lake properties. It controls the solubility of gases and minerals, the rate of chemical reactions, and biological activity and diversity (e.g., Benson and Krause, 1980; Rasconi et al., 2017; Krumgalz, 2018). In addition, the vertical temperature profile in a lake (and consequent lake stratification and water column stability) and the length of the stratification period play a vital role in the transport pathways of gases and nutrients and, consequently, their distribution within a lake (e.g., Vachon et al., 2019; Ladwig et al., 2021). Furthermore, there is a two-way interaction between lakes and the atmosphere. While the thermodynamic behavior of lakes is mainly driven by meteorological conditions, the distinct physical features of lakes (such as surface roughness, albedo, heat capacity and/or temperature, and evaporation rate) introduce surface heterogeneity in the domain of interest. Thus, their presence modifies surface–atmosphere fluxes and local and regional weather and climate (e.g., MacKay, 2012; Klaić and Kvakić, 2014; Bryan et al., 2015; Kristovich et al., 2017; Wu et al., 2019). Thus, over the last couple of decades, increasing scientific interest has been focused on both modeling the thermal regime of lakes (e.g., Stepanenko et al., 2013, 2016; Thiery et al., 2014; MacKay et al., 2017) and its sensitivity to climate change (e.g., Råman Vinnå et al., 2021).

Due to their relative simplicity and computational efficiency, there is a widespread use of one-dimensional (1-D) water temperature prediction models. There are different types of 1-D models of varying complexity, although they can generally be divided into three groups: (1) mixed layer models based on the energy budget approach, (2) differential models based on solving the 1-D heat transfer equation (thermal diffusivity models), and (3) second-order turbulence closure models. Energy budget-based models assume series of well-mixed layers (often just two, namely, the epilimnion and hypolimnion), and they use the kinetic energy produced by wind shear on the surface to account for the mixing dynamics within these layers and/or to estimate their depths (e.g., Bell et al., 2006; Mironov et al., 2010; Hipsey et al., 2019). Thermal diffusivity-based models usually consist of many well-mixed layers for which the heat transfer equation is solved (e.g., Hostetler and Bartlein, 1990; Liston and Hall, 1995; Stefan et al., 1998; Sun et al., 2007). The second-order turbulence closure models are also known as k-ε, where k is the turbulent kinetic energy per unit mass and ε is the turbulent kinetic energy dissipation rate (e.g., Goudsmit et al., 2002). They solve the turbulent kinetic energy transport equation and are, computationally, considerably more expensive than the previous two types (e.g., Goudsmit et al., 2002; Stepanenko et al., 2011, 2014).

Except in the basic underlying approach, lake models differ in the processes they include, such as wind sheltering, sediment heat flux, attenuation of light, phase change, convective mixing, and others. Direct implementation of a particular process in a model or the simplification or even omission of the process is usually justified by the model purpose. Lake models are developed for various purposes, including improvement of numerical weather prediction and climate models (e.g., Mironov et al, 2010; MacKay, 2012), evaluation of the effects of climate change (e.g., Stefan et al., 1998; Wu et al., 2020; Råman Vinnå et al., 2021), or facilitation of specific limnological studies. Some of these specific studies address gas (e.g., methane and/or CO2) emissions (e.g., Stepanenko et al., 2011), oxygen and nutrient levels (e.g., Bell et al., 2006), heat and mass exchange between the atmosphere and a water body (Sun et al., 2007), and evaporation and lake level fluctuation (Hostetler and Bartlein, 1990).

To run lake models, input data, which are generally not available from routine meteorological measurements, are needed. Specifically, these data include both shortwave and longwave radiation component data. The goal of this study is to formulate a simplified model for predicting the vertical temperature profile in a small, warm, monomictic lake, which, except for the ultraviolet B radiation (UVB), is forced solely by routinely available observed surface meteorological data (namely, the air temperature, relative humidity, atmospheric pressure, wind speed, and precipitation). Conversely, other lake temperature models that are forced with observational data (e.g., Bell et al., 2006; Sun et al., 2007; Martynov et al., 2010; MacKay, 2012, 2017) require both shortwave and longwave radiation component data and do not provide further details on determining them. The proposed model employs carefully chosen parameterizations of longwave and shortwave radiation. Although these parametrizations are well known, in the present study, they are built into a lake temperature model for the first time. Furthermore, in comparison with the model of Sun et al. (2007), the present model does not neglect the turbulent diffusion for small lakes and uses a different approach for calculating the light attenuation with depth. In addition, we examined the sensitivity of the proposed model performance to the length of the simulated time interval. To the best of our knowledge, such a detailed evaluation has not been reported in previous lake temperature modeling studies. Since vertical temperature profiles in lakes are not routinely measured, we also addressed the ability of the proposed model to simulate the onset and termination of lake stratification by a yearlong simulation initiated with a uniform temperature over a completely mixed water column. A similar study was performed by Martynov et al. (2010) for two small dimictic lakes in the USA using an eddy diffusivity model and a two-layer model; Goudsmit et al. (2002) analyzed the performance of a k-ε model in a two-year length simulation, while Bruce et al. (2018) analyzed a set of 32 lakes all over the globe using the General Lake Model (GLM). All of these models are more complex and/or require more extensive input data than the one proposed in this study.

The model proposed here is evaluated using lake temperature experimental data measured at two lakes of the Plitvice Lakes, Croatia. Details about the study area and data collection are presented in Sect. 2. The model's governing equations and parametrizations used are described in Sect. 3. Measures of the model performance and evaluation approach are described in Sect. 4. The results are presented and discussed in Sect. 5, and a comparison with other models is presented in Sect. 6. Finally, a short summary and conclusions are given in Sect. 7.

2 Study area and measurements

2.1 Study area

Plitvice Lakes is a karstic lake system situated in the mountainous region of Croatia (Fig. 1). The system consists of 16 named and several smaller unnamed lakes. The lakes are interconnected with cascades and waterfalls, making a chain approximately 9 km long and extending in a roughly south–north direction. With its unique geomorphology and exceptional biodiversity, the area has been a subject of scientific research dating as early as 1850 (NPPL, 2021). An extensive, multidisciplinary overview of abiotic studies focusing on the Plitvice Lakes area is provided by Klaić et al. (2018).

The numerical model proposed in this paper was applied to the two largest lakes of the system, Prošće and Kozjak Lakes (Fig. 1c and d). Prošće Lake (hereafter Lake 1) is the southernmost and the first lake in the system, while Kozjak Lake (hereafter Lake 12) is the 12th lake in the chain and the largest and deepest lake in the system. The characteristics of each lake are given in Table 1. Based on their surface areas, both lakes can be considered small (e.g., Forcat et al., 2011).

Figure 1Location of Plitvice Lakes (red bubble; source: © Google Maps) (a); closer look at the entire lake system (b); Lake 1 (c), and Lake 12 (d). Locations of the lake temperature measuring points P1 (φ=44.8676 N, λ=15.5981 E, height of the lake surface 636 m a.s.l.) and K1 (φ=44.8902 N, λ=15.6038 E, 535 m a.s.l.), and meteorological measuring site M (φ=44.8811 N, λ=15.6197 E, 579 m a.s.l.) are shown with yellow circles. Panels (b)(d) show composite pictures of the lake bathymetries and the digital orthophoto images (DOF:, last access: 3 December 2020, Print Rights – Under the Microsoft® Bing™ Maps Platform APIs' Terms of Use).

Table 1Characteristics of the studied lakes.

* a.s.l. – above sea level.

Download Print Version | Download XLSX

2.2 Observational data

2.2.1 Lake temperatures

This study uses lake temperatures measured at two different points (Fig. 1), one in Lake 1 (point P1, φ=44.8676 N, λ=15.5981 E, 636 m a.s.l.) and the other in Lake 12 (point K1, φ=44.8902  N, λ=15.6038 E, 535 m a.s.l.). Each point was positioned in the deepest part of the corresponding lake. Lake temperatures were measured and logged with HOBO TidBiT MX Temp 400, as previously described for Lake 1 in Klaić et al. (2020b) and for Lake 12 in Klaić et al. (2020a). The accuracy of the sensors is ±0.20C for temperatures between 0 and 70 C and ±0.25C for temperatures between −20 and 0 C. The initial sampling frequency of lake temperatures was 1 Hz, while 2 min means were stored. However, since meteorological data were available at a resolution of one hour, we used hourly mean lake temperatures in the present study.

At site P1 (Lake 1), 15 factory-calibrated sensors were positioned at fixed depths of 0.2, 0.5, 1, 1.5, 3, 5, 7, 9, 11, 13, 15, 17, 20, 23, and 27 m. As Lake 12 is deeper than Lake 1, an additional sensor was placed at a depth of 43 m at site K1 together with 15 sensors at the same depths as at site P1.

The temperature recording started on 7 July 2018 at K1 and 6 July 2019 at P1 (Table 2). Temperatures were recorded continuously, except during several short periods ( 1–2 d, once in approximately four months) when the sensors were pulled out of the lakes for the purpose of data acquisition. These periods without measurements are shown as thin, vertical white lines in Fig. 2a, c, and e. Due to the malfunction of some sensors during the first year of the measurement campaign, data for some observational depths at K1 are missing.

Table 2Availability of measured data. The positions of the measuring points are shown in Fig. 1.

Download Print Version | Download XLSX

Figure 2Measured water temperatures in Lake 1 (a and b) and Lake 12 (c and d) and water temperature in Lake 12 after interpolation of the measured data (e and f).


Missing data are shown as white areas from July 2018 to July 2019 in Fig. 2c or as intermitted lines in Fig. 2d. The inoperative sensors were later replaced. Missing data at specific depths were subsequently replaced by data calculated by spatial linear interpolation from the two adjacent depths using existing data (Fig. 2e and f). However, temporal interpolation was not performed, since it would fail to reproduce the temporal variability in lake temperature at particular depths during periods of data acquisition. Interpolated lake temperatures were used solely to illustrate the evolution of Lake 12 stratification (Fig. 2e); they were omitted in the calculations of the model performance measures (Sect. 4).

2.2.2 Meteorological data

Meteorological data were measured at the automatic meteorological station Plitvička Jezera (point M in Fig. 1, φ=44.8811 N, λ=15.6197 E, altitude 579 m a.s.l.). The station belongs to the network of the Croatian Meteorological and Hydrological Service (CMHS). The CMHS also provided quality control of these data. In the present study, we used hourly mean values of the air temperature, atmospheric pressure, UVB radiation, atmospheric relative humidity, hourly precipitation amount measured at 2 m above ground level, and wind speed measured at 10 m above ground level (Fig. 3). Wind direction data were also available but were not used in the study.

Figure 3Available meteorological data from the automatic meteorological station Plitvička Jezera (φ=44.8811 N, λ=15.6197 E, 579 m a.s.l.): (a) precipitation amount (P) and relative humidity (RH), (b) air temperature (Ta) and UVB radiation (UVB), and (c) wind speed (V) and atmospheric pressure (p).


The station is approximately 2 km northeastward of the P1 site and 1.6 km southeastward of the K1 site. Despite the comparable distance from both the P1 and K1 sites, the meteorological conditions observed at point M are expected to be more representative for Lake 12 than for Lake 1, because this point is located at the slope adjacent to Lake 12 at approximately 200 m away from its shoreline. In addition, topographic obstacles are found between points P1 and M (Fig. 1b), and the altitude difference between P1 and M is higher than the difference between K1 and M (Table 2).

3 Model description and governing equations

The model is based on the one-dimensional energy balance equation used in similar liquid water models (e.g., Hostetler and Bartlein, 1990; Liston and Hall, 1995; Sun et al., 2007). Because ice was not observed on the two lakes during the measurement campaign (Fig. 2a and c), ice formation was not addressed in the present study. Thus, a simplified approach using water temperature instead of enthalpy is used. Considering that, more often than not, the lake bathymetry is not available, in addition to our goal to keep the model as simple as possible and to limit the input data, it is assumed that the water body has a constant horizontal cross-sectional area (which can be of any shape). Thus, we come to the following equation:

(1) c p ρ T t = z k m + k t T z - Φ z + M conv ,

where cp is the water specific heat capacity (J kg−1 K−1), ρ is the water density (kg m−3), T is the water temperature (C), t is time (s), z is depth (m), km and kt are the molecular and turbulent thermal conductivity (W m−1 K−1), Φ is the heat flux (W m−2), and Mconv is the convective mixing term (W m−3).

The water density is calculated from the Chen and Millero (1986) formula, assuming zero salinity:

(2) ρ = 999.8395 + 6.7914 × 10 - 2 T - 9.0894 × 10 - 3 T 2 + 1.0171 × 10 - 4 T 3 - 1.2846 × 10 - 6 T 4 + 1.1592 × 10 - 8 T 5 - 5.0125 × 10 - 11 T 6 .

The molecular thermal conductivity of water is 0.6 W m−1 K−1 (e.g., Sun et al., 2007). The turbulent thermal conductivity is a function of time and depth, because it depends on meteorological forcing. Here, we also follow the method of Henderson-Sellers (1985), where the turbulent thermal conductivity is calculated as follows:

(3) k t z = c p ρ k u * z / Pr 0 exp - k * z 1 + 37 Ri 2 - 1 ,

where k=0.4 is the von Karman constant, u* is the friction velocity at the surface (m s−1), k* is the latitude-dependent parameter of the Ekman profile, Pr0=1 is the neutral value of the turbulent Prandtl number, and Ri is the Richardson number. The Ekman profile parameter and the Richardson number are calculated as in Sun et al. (2007):

(4) k * = 6.6 sin φ 1 / 2 U 2 - 1.84 ,

where φ is the latitude and U2 is the wind speed at 2 m above the water surface (m s−1), and

(5) Ri = - 1 + 1 + 40 N 2 k 2 z 2 / u * 2 exp - 2 k * z 1 / 2 20 ,

where N is the Brunt–Väisälä frequency (s−1):

(6) N = - g / ρ ρ / z 1 / 2 .

The wind speed U2 is determined from the logarithmic formula:

(7) U 2 = u * log 2 / z 0 / k ,

where z0 is the roughness length (m). The air shear velocity u* and the roughness length z0 are calculated as in Verburg and Antenucci (2010).

Although Sun et al. (2007) suggest that, for shallow lakes (less than 50 m deep), the turbulent thermal conductivity is negligible, this is not in accordance with findings of numerous other studies, which suggest that the turbulent thermal conductivity can be much larger than the molecular thermal conductivity, even for shallow lakes (e.g., Jassby and Powell, 1975; Quay et al., 1980; Vachon et al., 2019). It should be kept in mind that these studies often determine the turbulent diffusion coefficient based on measured change rate of lake water temperature vertical distribution. This means that the contributions of all present mixing processes are included (i.e., shear-induced turbulence, breaking internal waves, boundary layer turbulence). However, the mixing processes and their contributions to turbulent mixing may differ from lake to lake. In the present study, turbulent thermal diffusion was taken into account using Eq. (3).

3.1 Energy budget and boundary conditions

In addition to turbulent thermal diffusion, the only other term in Eq. (1) accounting for meteorological forcing is the heat source term. The surface net heat flux consists of the net shortwave radiation (Sn), net longwave radiation (Ln), sensible heat flux (Hs), latent heat flux (Hl), and heat flux brought by precipitation (Hp). The surface boundary condition can be written as follows:

(8) Φ ( 0 ) = S n + L n + H s + H l + H p .

At the bottom, it is assumed that there is no heat flux and that the temperature gradient equals zero, meaning there is no heat diffusion either. Thus, the bottom boundary conditions can be written as follows:


All heat flux terms in Eq. (8) are defined to be positive when downward. Shortwave and longwave radiation measurements are not very common, and sensible and latent heat fluxes cannot be measured directly (Brunel, 1989; Bahr et al., 2012). Thus, obtaining the heat flux terms in Eq. (8) is expensive and complicated. Therefore, methods for calculating each term using commonly available meteorological data only are proposed in sections 3.1.1. to 3.1.4.

3.1.1 Shortwave radiation

As previously indicated by other authors (e.g., Bell et al., 2006; Martynov et al., 2010; MacKay, 2012), sufficient radiation data (both shortwave and longwave) are not generally available from routine meteorological measurements, and this is also the case for meteorological station M, where only UVB radiation was measured. A number of studies provide correlations among UVA, UVB, total UV, or global solar radiation (G) (Kudish and Evseev, 2000; Kudish et al., 2005; Podstawczynska, 2009; Pokhrel and Bhattarai, 2012; Pashiardis et al., 2017) and show that significant variability occurs in the UV/G ratio between sites, which is mainly due to local atmospheric conditions. Podstawczynska (2009) indicated that air turbidity and cloudiness are the two main factors that determine the variability of daily solar energy transmission through the atmosphere. Pashiardis et al. (2017) found that the UV/G ratio increases with solar elevation and that the presence of clouds reduces the UV component less than the global solar radiation due to the strong absorption of water in the near-infrared spectrum.

Winslow et al. (2001) proposed a model for estimating the total daily solar irradiance from daily precipitation and minimum and maximum temperatures, along with latitude, elevation, and mean annual temperature. This model showed significant improvement over the widely used empirical Bristow and Campbell (1984) model and was applicable for a wide range of climates. Therefore, it is also used in this study.

According to Winslow et al. (2001), the daily solar irradiance at the Earth's surface is equal to

(11) S surf = τ cf D 1 - β s rh T max S top ,

where Ssurf is the total daily solar irradiance at the surface (J m−2 d−1), τcf is the cloud-free atmospheric transmittance, βs is an additional parameter required to introduce variation between sites, rhTmax is the relative humidity at the moment when daily maximum air temperature (Tmax) is reached, and Stop is the total daily solar irradiance at the top of the atmosphere (J m−2 d−1). The total daily solar irradiance is calculated following Wald (2019):

(12) S top = S 0 1 + ε ecc 3600 × 24 π ( cos φ cos δ sin ω s + ω s sin φ sin δ ) ,

where S0=1362 W m−2 is the solar constant, εecc is the eccentricity of Earth's orbit, δ is the solar declination, φ is the location latitude, and ωs is the half-day length (time between sunrise and noon or noon and sunset) in radians. εecc and ωs are functions of the day in the year only, while δ also depends on the location longitude, since its noon value is used, which yields more precise results. Details on calculating these parameters are included in Wald (2019).

The cloud-free atmospheric transmittance in Eq. (11) accounts for the transmittance of dry clean air (τ0) and the transmittance due to absorption by aerosols (τa) and water vapor (τv), and it also incorporates a correction for elevation (celev):

(13) τ cf = τ 0 τ a τ v c elev .

To calculate τ0, τa, τv, celev, D, and βs, we follow Winslow et al. (2001). The transmittance of dry clean air is dependent only on the latitude (φ) and is calculated as follows:

(14) τ 0 = 0.947 - 1.033 × 10 - 5 φ 2.22 for φ 80 τ 0 = 0.774 for φ > 80 .

The absorption by aerosols is extremely variable. Similar to Winslow et al. (2001), we set τa=1 (i.e., no absorption).

The absorption by water vapor is calculated from the following:

(15) τ v = 0.9636 - 9.092 × 10 - 5 T mean + 30 1.8232 ,

where Tmean is the mean annual air temperature (C). On wet days, when the daily precipitation is above 1 mm, τv is reduced by 0.13. The site elevation correction factor (celev) is calculated as follows:

(16) c elev = 1 - 2.2569 × 10 - 5 z a . s . l . 5.2553 ,

where za.s.l. is the site elevation (m).

From Eq. (11), τcfStop is the maximum cloud-free value of Ssurf. The effect of cloudiness is indirectly taken into account by introducing the factor D(1−βrhTmax). This is based on the finding that the solar irradiation from sunrise, when minimum humidity is expected (rhTmin1), until the maximum daily air temperature (and minimum humidity rhTmax) is reached, is proportional to the decline of the relative humidity, Ssurf_Tmax (1−βrhTmax). The factor D=Ssurf/Ssurf_Tmax is introduced to account for the surface solar irradiation from the moment when the air temperature reaches its daily maximum until sunset. D is calculated assuming that the air temperature reaches its maximum around 3 pm:

(17) D = 1 - ω s - π / 4 2 / 2 ω s 2 - 1 .

The factor βs in Eq. (11) is mainly constant, except for regions with very large daily temperature ranges:

(18) β s = max 1.041 , 23.753 Δ T m / T mean + 273.16 ,

where ΔTm is the mean annual temperature range between the daily air temperature maximum and minimum.

Hourly shortwave radiation data were generated from the calculated daily solar irradiance by using the measured UVB radiation data as a weight function:

(19) S ( h ) = UVB ( h ) S surf UVB day = UVB ( h ) S surf 3600 h = 1 24 UVB ( h ) ,

where Ssurf and UVBday are the daily values (J m−2 d−1) and S(h) and UVB(h) are the mean values (W m−2) for the hth hour of the total and UVB solar radiation, respectively. When UVB radiation data are unavailable, the standard daily radiation profile can be used.

Unlike the other terms in Eq. (8), shortwave radiation is not completely absorbed in the lake surface layer but partially passes through the water. The net shortwave radiation reaching a particular depth is calculated using the arctangent model, which was chosen for its simplicity for implementation, as suggested by Henderson-Sellers (1986), but also for its better representation of the light attenuation in the shallow layers, which are usually a lot thinner than the deeper ones:

(20) S n ( z ) = ( 1 - α ) S exp - K 1 z 1 - K 2 tan - 1 K 3 z ,

where Sn(z) is the net shortwave radiation at water depth z (W m−2), α=0.06 is the water surface albedo, and K1, K2, and K3 are empirical constants. K1 corresponds to the light extinction coefficient λe=0.1 (value of 0.1 is appropriate for clear oligotrophic lakes). K2 is calculated as

(21) K 2 = 2 1 - 1 - β exp λ e z A / π ,

where β=0.4 accounts for the absorption in the surface layer, and zA=0.6 m is the depth of the surface absorption layer, where the exponential decay starts. The third parameter, K3=4, is not a direct function of λe and β but is rather a measure of the rapidity of falloff with depth in the upper layers.

3.1.2 Longwave radiation

The net longwave radiation is the difference between the incoming downward atmospheric longwave radiation (La) and the outgoing upwards radiation from the lake surface (Ls). As direct measurement data of longwave radiation by pyrgeometers are not routinely available, longwave radiation may be calculated using the following formula:

(22) L n = ( 1 - r ) L a - L s = ( 1 - r ) ε a σ T a + 273.15 4 - ε σ T s + 273.15 4 ,

where r is the water reflectivity for longwave radiation, ε and εa are the emissivities of the lake surface and the atmosphere, respectively; Ts is the water surface temperature (C), Ta is the air temperature at 2 m height (C), and σ=5.67× 10−8 W m−2 K−4 is the Stefan–Boltzmann constant. The reflectivity and emissivity of water are assumed to be 0.04 and 0.96, respectively (e.g., Sun et al., 2007). The emissivity of the atmosphere depends on the water vapor and atmospheric temperature profile. Assuming a standard atmosphere, Brutsaert (1975) derived a formula for calculating the atmospheric emissivity under clear-sky conditions:

(23) ε ac = 1.24 e a / T a + 273.15 1 / 7 ,

where ea is the water vapor pressure (hPa), which is related to the relative humidity (rh) and saturation vapor pressure (es):

(24) e a = e s ( T a ) rh .

To calculate the saturation water pressure, we use the formula from Winslow et al. (2001):

(25) e s ( T a ) = 6.11 exp m T a / n + T a for T a > 0 C m = 17.269 n = 237.7 for T a < 0 C m = 21.875 n = 265.3 .

Although other empirical formulas for calculating atmospheric emissivity are available, Brutsaert's (1975) expression (Eq. 23) was reported as the best in many studies of different climates (Wang and Dickinson, 2013). Because Eq. (23) refers to clear-sky conditions, it is necessary to additionally account for cloud effects. Assuming that the emissivity of the water droplets in the clouds is approximately equal to 1, Crawford and Duchon (1999) calculate the total atmospheric emissivity as follows:

(26) ε a = 1 - f ε ac + f ,

where f is the cloud fraction term defined using the ratio of the previously estimated surface shortwave radiation and surface clear-sky shortwave radiation:

(27) f = 1 - S surf / τ cf S top .

For clear-sky conditions, the cloud fraction term equals 0. However, since the ratio of the surface solar irradiance to the clear-sky irradiance never reaches 0, the cloud fraction term never reaches the theoretical maximum of 1, even in total cloud cover conditions. Note that, even though the model will be run with a time resolution of one hour, the daily mean atmospheric emissivity will be used.

Equation (26) is considered the best formula in many studies (Wang and Dickinson, 2013). By substituting Eqs. (23) and (26) in Eq. (22), we obtain the expression for calculating the net longwave radiation:

(28) L n = ε 1 - f ε ac + f σ T a + 273.15 4 - ε σ T a + 273.15 4 .

3.1.3 Latent and sensible heat flux

To calculate the latent and sensible heat flux, we use a slightly modified algorithm provided by Verburg and Antenucci (2010). Their code, which is publicly available at the National Institute of Water and Atmospheric Research (NIWA) website (NIWA, 2021), uses the bulk aerodynamic method based on the Monin–Obukhov similarity theory (Monin and Obukhov, 1954). According to this method, the sensible and latent heat fluxes can be calculated as follows:


where CH and CE are the transfer coefficients for sensible and latent heat flux, respectively; ca=1005 J kg−1 K−1 is the specific heat of air, Lv≈2500 kJ kg−1 is the latent heat of evaporation, ρa is the air density (kg m−3), and qs and qa are the specific humidities (kg kg−1) at the water surface and measuring levels, respectively. Air density and specific humidity were determined from the ideal gas law equation and from the observed relative humidity, respectively.

The transfer coefficients were calculated in an iterative procedure, initially assuming neutral atmospheric conditions:


where CD is the drag coefficient, h is the height above ground (m), z0 and zE are the roughness lengths (m), and ψM and ψE are the stability functions for momentum and vapor, respectively. The stability functions are defined through the stability parameter ζ=h/L, where L is the Monin–Obukhov length:

(34) L = - ρ a u 3 T V k g H s c a + 0.61 T a + 273.15 H l L v ,

where Tv is the virtual temperature. Obviously, L depends on Hs and Hl, while Hs and Hl depend on the stability of the atmosphere. Therefore, to calculate Hs and Hl, an iterative procedure has to be used. The procedure is initiated by assuming neutral conditions (ψM=ψE=1). Further details on the calculation of roughness lengths, stability functions, and the iterative process itself can be found in Verburg and Antenucci (2010).

3.1.4 Heat brought by precipitation

Assuming the first lake layer in the numerical model gets completely mixed with the precipitation that falls during a time period Δt (s), then the temperature of that layer would equal

(35) T 1 + p = Δ z 1 T 1 + P / 1000 × 3600 Δ t T prec Δ z 1 + P / 1000 × 3600 Δ t ,

where T1 and T1+p represent the water temperature of the first layer before and after the precipitation has been introduced (C), Tprec is the precipitation temperature (C), Δz1 is the thickness of the first layer (m), and P is the hourly precipitation (mm h−1). The heat flux brought in by precipitation Hp (W m−2) can then be calculated as

(36) H p = 1 Δ t Δ z 1 + P 1000 × 3600 Δ t ρ c p T 1 + p - T 1 = ρ c p P 1000 × 3600 T prec - T 1 .

Since Tprec was not available, the air temperature was used instead.

3.2 Convective mixing

During the night, the net heat flux at a lake surface is generally negative. Consequently, unstable lake stratification is established. However, this unstable stratification is short lived, because the higher density water forming on top of the lake quickly sinks and mixes with the lower density water below it, thus restoring equilibrium (i.e., minimum potential energy).

As Sun et al. (2007) pointed to the importance of introducing a convective mixing mechanism in a water temperature model, we also incorporated this mechanism in the present model. Namely, after each time step of integration, the model algorithm checks whether the upper layer in each pair of two adjacent layers has a higher density than the lower layer. If this occurs, then the two layers are assumed to mix completely, which results in uniform temperature:

(37) T j _new = T j + 1_new = T j Δ z j + T j + 1 Δ z j + 1 / Δ z j + Δ z j + 1 ,

where Δzj and Δzj+1 represent the thickness of the jth (upper) and (j+1)th (lower) layers, respectively; Tj and Tj+1 are the water temperatures in these layers before convective mixing, respectively; Tj_new and Tj+1_new are the temperatures in these layers after convective mixing, respectively.

3.3 Model setup

The model code is written in MATLAB programming language. Equation (1) is discretized using the backward Euler scheme:

(38) c p ρ j Δ t T j n + 1 - T j n = 1 Δ z j k m + k t , j + 1 / 2 T j + 1 n + 1 - T j n + 1 z j + 1 - z j - k m + k t , j - 1 / 2 T j n + 1 - T j - 1 n + 1 z j - z j - 1 - 1 Δ z j Φ j + 1 / 2 n - Φ j - 1 / 2 n ,

where the subscript denotes the layer or the boundary between two layers, and the superscript denotes the time increment. Notice that the convective term from Eq. (1) is omitted in Eq. (38), since the algorithm employs convective mixing in a separate procedure after the integration step only if density inversion is detected in the water column, as explained in Sect. 3.2 and shown in Fig. 4. After the stability check, the algorithm performs a step which limits the temperature minimum to 0 C. Namely, as the model does not include an ice formation module, this step roughly assures no unreasonably low temperatures appear (Fig. 4). Equation (38) can be rearranged as follows:

(39) T j - 1 n + 1 - k m - k t , j - 1 / 2 z j - z j - 1 + T j n + 1 Δ z j c p ρ j Δ t + k m + k t , j + 1 / 2 z j + 1 - z j + k m + k t , j - 1 / 2 z j - z j - 1 + T j + 1 n + 1 - k m - k t , j + 1 / 2 z j + 1 - z j = Δ z j c p ρ j Δ t T j n + Φ j + 1 / 2 n - Φ j - 1 / 2 n .

Equation (39) can be written in matrix form as follows:

(40) M T n + 1 = A T n + B .

Then, the solution for Tn+1 is obtained as follows:

(41) T n + 1 = M - 1 A T n + M - 1 B .

The implicit Euler scheme is unconditionally stable and thus does not have an upper limit for the time increment. Considering the time resolution of the available input data, the model was run with a time step of one hour (runs with finer time steps were attempted; however, the performance improvements were not significant). The vertical resolution in the model corresponds to the measuring depths and decreases with lake depth. The depths of the integration points were consistent with the sensors' depths, while the boundaries of the layers were set halfway between each pair of consecutive points (Fig. 4b). The layer thicknesses ranged from 0.35 (surface layer) to 16 m (bottom layer). An overview of the model workflow is given in Fig. 4a.

Figure 4Model configuration. Panel (a) shows a schematic overview of the model workflow. The input consists of the initial time and date, the initial water temperature profile Tinitial, meteorological data (wind speed, air temperature, UVB radiation, relative humidity, precipitation, and atmospheric pressure), climatological data (mean annual air temperature, mean annual temperature range between the daily air temperature maximum and minimum, and mean monthly cloud cover), and location data (location latitude, longitude, and altitude above sea level). Panel (b) shows the layer setup. Points from 1 to J indicate the integration points where water temperatures are calculated, and zj is the depth of the jth point. The horizontal lines indicate boundaries between layers, Φj±1/2 are the heat fluxes across the layer boundaries, and Φ0 is the net surface heat flux.


4 Measures of the model performance

First, a sensitivity analysis was performed to assess the dependence of the model performance on the simulation length. A simulation run was initiated in every hour of the periods with available data, and each was run for up to 30 d. Measured water temperature profiles were used for simulation initialization. Predicted water temperatures and vertical temperature gradients obtained in each simulation after a certain amount of simulation time (from 1 to 30 d) were compared with the corresponding observed values. The results of this analysis are to show the model's ability to provide quality short-term prognoses and the rate of the result deterioration with the increase of the simulation length.

The model performance for each simulation length was evaluated by common bivariate measures. The mean bias error (MBE) is used to assess the tendency of the model to over- or underpredict the temperature. The mean absolute error (MAE) and the root mean square error (RMSE) both provide information about the error central tendency. However, RMSE also accounts for the distribution of the error and becomes larger as the error variability increases. RMSE places more weight on large errors, which makes it more sensitive to outliers. Due to all of the above, it has been argued that MAE is a more natural measure of average error than RMSE (Willmott and Matsuura, 2005). The maximum absolute error (MaxAE) is not a measure of systematic error, but it was calculated as a measure showing the most extreme outlier. The above measures are calculated from the following expressions:


where O and P correspond to the observed and predicted values, respectively, while n is the number of corresponding pairs of these values.

The index of agreement values were calculated using three different formulas proposed by Willmott et al. (2012), namely, the original (IAorig), modified (IAmod), and refined (IAref) index of agreement:



The IA represents a measure of the relative covariability of the observed and predicted values with respect to the observed mean. The original IA (Eq. 46) uses the square of the difference between predicted and observed values, which is why it overestimates the influence of large errors, similar to the RMSE, which is why the square is replaced with an absolute value in the modified version (Eq. 47); thus, IAmod is less sensitive to outliers than IAorig. IAmod approaches 1 (perfect agreement) more slowly than IAorig, which means that IAmod is more conservative and allows for finer comparisons of different models with relatively good performance. In IAref (Eq. 48), the prediction variability in the denominator is replaced with the observation variability. IAorig and IAmod range from 0 to 1, where a value of 0 means that the prediction and observation variabilities are out of phase, while a value of 1 means perfect fit. IAref ranges from −1 to 1 and has a well-defined lower boundary (Eq. 48b), which allows for a better comparison of models with poor performance. However, it should be stressed that IAref approaching the value of −1 does not necessarily indicate poor model performance, because it can also be a result of low observation variability.

The second goal of this study was to examine the ability of the model to predict the springtime onset of lake stratification assuming that there are no measured water temperature data available. For this purpose, a simulation initiated with approximately constant water temperature throughout the entire lake column, which is characteristic of the period when a lake is mixed, was run for the entire year, starting from 1 January 2019. Although accurate results were not expected for the yearlong simulation, the goal of this analysis was to evaluate the extent to which the model can provide relevant information regarding the stratification and/or thermocline depth. Such an approach is particularly appealing for lakes that are completely mixed during the winter, since it does not require measurement of the water temperature profile to initiate the simulation.

5 Results and discussion

Based on sporadic observations of the vertical temperature profiles in the Plitvice Lakes, previous studies suggest that Lake 1 and Lake 12 are dimictic (Klaić et al., 2018). Dimictic lakes are covered by ice during winter; they mix in spring and fall; and they are stratified in summer. The continuous observation data of the vertical temperature profiles in Lake 1 and Lake 12, shown in Fig. 2a and c, clearly illustrate for the first time that, during the field campaign, both lakes behaved as warm, monomictic lakes. Specifically, they were mixed during winter but stratified at other times. Furthermore, monomictic lakes (which are frequently found in temperate and tropical latitudes) typically do not freeze, and the two studied lakes did not freeze during the entire field campaign, since the wintertime temperatures in the top lake layers were above 0 C (Fig. 2b and d).

As the main driver of the lake temperature profile is the surface heat flux, it is interesting to first analyze its terms. Figure 5 shows the modeled mean diurnal variation in the total heat flux and the heat flux terms for a typical winter (a) and summer month (b). The solar heat flux is an order of magnitude higher than the other components of the total heat flux, which indicates that it is one of the main factors affecting the lake water temperature. Next in magnitude is the net longwave radiation, followed by the latent heat flux. The last two components are negative and are responsible for the negative heat flux, or cooling, at night.

Figure 5Modeled mean diurnal variations in the heat flux at the surface of Lake 12 for January (a) and July 2019 (b).


The observed and predicted water temperatures for various simulation lengths for 2019 are shown in Fig. 6 (Lake 1 – note that lake temperature measurements started in July) and Fig. 7 (Lake 12). The model performed reasonably well. Namely, the onset of the stratification period (Fig. 7) and both the vertical temperature profile and deepening of the thermocline over time were well captured (Figs. 6 and 7). Simulation results for Lake 12 reproduce the observed data more closely, while for Lake 1, higher discrepancies between simulated and observed data are present, especially for simulation lengths above 10 d. For Lake 1, the position of the maximum temperature gradient in the metalimnion between 12 and 16 m depth was captured, even in the 30 d simulations (Fig. A3), but the temperatures in the epilimnion are significantly overestimated in the stratification period (August) in the longer runs (Figs. 6 and A1).

For Lake 12, the difference between the predicted and observed position of the maximum temperature gradient is within 2 m, even for the 30 d simulations, but generally it is lower. Temperature overprediction is noticed in the epilimnion, especially towards the end of the year, for the simulation lengths above 10 d. The stratification began on 21 March, and in the 30 d simulations it was predicted on 23 March; the convective overturn began on 6 September, while in the 30 d simulation it was predicted on 10 September.

Figure 8 shows a closer view of the observed vs. predicted temperatures at depths of 0.2, 5, 15, and 27 m for the period between 6 July 2019 and 31 December 2019. This period was chosen because it is the longest period in which all necessary data (both meteorological and water temperatures for both lakes) were available. Additionally, observed vs. predicted temperature gradients and the prediction errors for both temperature profiles and temperature gradients for the same sample period are presented in Appendix A. As expected, the departure of the predicted from the observed quantities increases with the length of the simulation period. However, even the longest simulation runs (30 d) produced qualitatively acceptable results. Departures of the predicted hourly mean temperatures were mainly ±2 and ±1C for Lake 1 and Lake 12, respectively, except in the thermocline region, where they were mainly ±4 and ±2C for Lake 1 and Lake 12, respectively. The temporal temperature variations at various depths were satisfactorily simulated (Figs. 8, A1, and A2). Furthermore, thermocline depths and their deepening in time were well captured by the model (Figs. A3 and A4). However, the results also suggest that the lake temperatures are somewhat overpredicted in the epilimnion and, at times, may be slightly underpredicted in the hypolimnion (Figs. 8, A1 and A2).

Although the model satisfactorily reproduced the temporal variations in the lake temperatures at the hourly scale, it was not able to reproduce the internal seiches that were previously documented for both lakes (Klaić et al., 2020a, b). This finding is not surprising, since the present model is based solely on the energy budget approach; thus, except for vertical mixing of the two adjacent layers under unstable stratification, it does not account for any hydrodynamic behavior.

Figures 9 and 10 show the calculated model performance measures for both lakes. The model overestimates the water temperature in the epilimnion, especially near the surface and the thermocline region, with an MBE from 0.3 and <0.1C for 1 d simulations and up to 2.6 and 1.2 C (at 5 m depth) for 30 d simulations in Lake 1 and Lake 12, respectively (Figs. 9a and 10a). The MAE in the epilimnion in Lake 1 starts from <0.4C for 1 d simulations and increases relatively steadily to 2.6 C for 30 d simulations (Figs. 9b). In Lake 12, it also starts from <0.4C for 1 d simulations and slowly increases to 1.2 C as the simulation length reaches 30 d (Fig. 10b).

A couple of factors could lead to overestimated temperatures in the upper lake layers. The first is the underestimation of turbulent mixing and turbulent heat transfer, especially in periods of high winds. As seen from Figs. A1 and A2, this overestimation of the uppermost part of the lake is more pronounced for Lake 1 than for Lake 12. As argued in Sect. 2.2.2, measuring site M (where the data used for the atmospheric forcing of the model are measured) is more representative for Lake 1 than for Lake 12. Accordingly, due to its higher altitude and less sheltered position, Lake 1 is more likely exposed to winds stronger than those measured at site M, and thus, both the turbulent mixing and the consequent heat transfer are likely to be stronger than modeled.

The second possible reason is the overestimation of the shortwave radiation extinction coefficient. This coefficient depends on the amount of dissolved organics and particulates in the lake water and can thus be calibrated to reproduce the lake physical properties more closely. We did not proceed with extinction coefficient calibration, because our goal was to investigate the model performance and its general applicability without location-specific fitting.

Figure 6Observed (a) and predicted (b–f) water temperatures of Lake 1 for different simulation lengths in the period between 6 July 2019 and 31 December 2019. Periods with missing data are seen as white vertical stripes.


Figure 7Observed (a) and predicted (b–f) water temperatures of Lake 12 for different simulation lengths for 2019. Periods with missing data are seen as white vertical stripes.


Figure 8Observed and predicted water temperatures at depths of 0.2, 5, 15 and 27 m for different simulation lengths in (a) Lake 1 and (b) Lake 12 in the period between 6 July 2019 and 31 December 2019.


Also, it is possible that the surface heat flux has been overestimated, as the simplified approach used for its calculation is characterized by limited reliability. Finally, it should be pointed out that the influence of the tributary was not considered, and in case of Plitvice Lakes, it may be non-negligible.

In the hypolimnion, the values of the MBE, MAE, RMSE, and MaxAE remain particularly low, especially for the deepest layers, for both lakes regardless of the simulation length (Figs. 9 and 10). These low values are a result of the low temperature variability in the deep lake layers, which is not taken into account in the formulation of these measures. In Lake 12, the MBE in the hypolimnion stays below 0.5 C, and in Lake 1, it stays below 1 C. On the other hand, regardless of the formulation (original, modified, or reference), the IA takes the temperature variability into account and therefore decreases with the increase of the simulation length, even in the deep layers, indicating poorer performance as the simulations get longer.

Further inspection of the results for temperature and temperature gradients in Lake 12 (Figs. A2 and A6) shows that the temperature prediction in the metalimnetic layer (thermocline region), where the temperature gradients were the highest, was rather challenging. The model performed relatively poorly in this region, which is particularly noticeable for longer simulation periods.

As seen from Figs. 9 and 10, MaxAE did not increase significantly with increasing simulation length for either of the two lakes. As expected, MaxAE was highest near the surface, and the maximum values for both lakes were relatively high (6.7 C at 3 m depth and 4.9 C in the surface layer for Lake 1 and Lake 12, respectively).

Figure 9Model performance parameters for Lake 1 (calculated for all the periods with necessary data available: 6 July–31 December 2019 and 2 July–30 September 2020).


Figure 10Model performance parameters for Lake 12 (calculated for all the periods with necessary data available: 7 July–4 November 2018; 1 January–31 December 2019; 2 July–30 September 2020).


Figure 11 shows the monthly means of the observed and modeled vertical temperature profiles for Lake 12. The results for Lake 1 are not shown, because the necessary observation data were not available throughout any complete year during the observational campaign. The model successfully reproduced the annual variation in the temperature profile throughout 2019, including the stratification onset and its termination and the thermocline deepening over time. The model underpredicts the surface temperature in January and February (Fig. 11a and b). As the heating starts in spring and continues through summer, the model tends to slightly overpredict the temperatures (Fig. 11e–h). The difference between the predicted and observed surface temperature in the summer months stays below 1 C for the longest simulation period. This finding is consistent with the above presented discussion of the model performance measures. In August (Fig. 11h), although the stratification was still strong, the effects of convective mixing during the night started to affect the monthly mean. In the following months (Fig. 11i–l), the mixing depth grew and reached a maximum depth of approximately 20 m in December (Fig. 11l), while the lake stratification was much weaker than that in previous months. In these months of significant convective overturn, higher deviation of the predicted from the observed epilimnion temperature is observed. It becomes more significant with simulation length and reaches approximately 2 C for the longest simulation period.

The second goal of this study was to examine the ability of the model to predict the onset and termination of stratification and the deepening of the thermocline by yearlong simulation. Because all necessary data for the entire year were only available for 2019, the first day of the yearlong simulation was set to 1 January 2019. For Lake 12, the simulation was initiated with a nearly constant water temperature profile ( 4 C) that was observed for 1 January because these data were available, although a constant temperature of  4 C, which was generally observed over the entire water column (which is typical for the wintertime period, when a lake is mixed), can be used instead.

Figure 11Annual variation in the vertical profile of the water temperature. Panels (a) to (l) show the monthly means of the observed and predicted values in Lake 12 for 2019.


Figure 12Observed (a) and predicted (b) water temperature for Lake 12 in 2019. The predicted temperature is obtained by a single yearlong simulation run initiated on 1 January 2019.


Figure 13Observed (a) and predicted (b) water temperature gradient for Lake 12 in 2019. The predicted temperature gradient is obtained by a single yearlong simulation run initiated on 1 January 2019.


Figure 12 shows the contour diagrams of the observed and predicted water temperatures for Lake 12. Such results for Lake 1 are not shown, because they were almost identical to those obtained for Lake 12. Namely, meteorological forcing drove temperature changes. If the same forcing was used for both lakes, then the only other factor that can introduce a difference in the results was the initial vertical profile, which was very similar for both lakes. As previously pointed out in the discussion of model performance, the model generally overpredicted the temperatures of the upper layers, especially by the end of the year. However, the onset and termination of the stratification period were well predicted, with the onset being captured somewhat better than the termination. The first noticeable temperature rise and the early beginning of stratification appear on 21 March in the observed data and 18 March in the predictions (Fig. 12). Significant temperature gradients exceeding 2 C m−1 appear on 12 June; however, the maximum gradient appears at a depth of around 2.5 m in the predictions and at depth of 5 m (Fig. 13) in the observations. The thermocline depth increases during the summer, and the maximum temperature gradient appears on 21 September at a depth of 12 m and equals 2.5 C m−1. On the same date, the maximum predicted temperature gradient appears at the same depth but equals only 1.3 C m−1. Actually, Fig. 13 shows that, while the model accurately predicts the upper limit of the metalimnion, it consistently overestimates its thickness, which consequently leads to underprediction of the temperature gradient in it. With temperature gradients below 0.5 C m−1, 6 December may be considered as the point of complete end of stratification. Although the predicted mixing depth is in agreement with the measured data, the overestimated epilimnion temperature consequently leads to temperature gradients of around 0.6 C m−1, which persist until the end of the simulation.

Figure 14Monthly means of the observed (blue) and predicted (red) water temperature vertical profiles for Lake 12 in 2019. Predicted temperatures are obtained by a single yearlong simulation run initiated on 1 January 2019.


This is more clearly presented in Fig. 14, where the observed and modeled monthly profiles are shown. Here, it can also be seen how the model generally overpredicted the monthly mean lake temperatures. The discrepancies between the modeled and observed profiles were largest in the mixed layer during the fall and winter convective overturn. Nevertheless, the mixing depth was well captured. It is concluded that the modeled results satisfactorily reproduced the monthly mean profiles and their annual variation, except after the convective overturn, when higher temperature overestimation is observed.

6 Comparison with other models

To compare the performance of the proposed model with the performance of more complex models, we applied 1-D General Ocean Turbulent Model (GOTM;, last access: 30 August 2021) version 4.1.0 and Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM; Zhang et al., 2016) version 5.3 for Lake 12 for a one-year period, starting from 1 January 2019. GOTM model is a one-dimensional water column model designed for hydrodynamic, thermodynamic, and biogeochemical studies of lakes and enclosed or semi-enclosed marine water bodies. It simulates vertical transport of momentum, heat, and salt (Burchard et al., 1999). The model, which can be used as a standalone or coupled with other models, has several turbulence closure options. So far, the GOTM model has been applied in a number of oceanographic (e.g., Bruggeman and Bolding, 2014; Burchard et al., 2014; Li et al., 2021) and limnetic studies (e.g., Ciglenečki et al., 2015; Andersen et al., 2021; Nielsen et al., 2021). SCHISM is a three-dimensional (3-D) unstructured-grid model. It employs hydrostatic approximation and solves the Reynolds-averaged momentum as well as the continuity and the transport of salt and heat equations. Due to its unstructured grid, it is suitable for basins with complicated geometry. It has been widely used in hydrodynamic studies of rivers, coastal waters, seas and oceans (e.g., Jacob et al., 2016; Bubalo et al., 2018; Zhang et al., 2020; Burić et al., 2021), and lakes (e.g., Frishfelds et al., 2021). More details on both models and parameterizations employed in the present study are given in Appendix B.

As previously presented, the meteorological forcing for the SIMO simulation was modeled using solely measured data. Apart from the measured air temperature and wind data (GOTM simulation) and measured air temperature (SCHISM simulation), meteorological forcing was modeled with the atmospheric Weather Research and Forecasting (WRF) model (Skamarock and Klemp, 2008). In both GOTM and SCHISM simulations, freshwater was assumed. Also, due to consistency, in both model runs, the same k-ε turbulence closure scheme of Rodi (1984) was employed. Finally, both models were initialized with the lake temperatures observed at 1 January 2019 (same as SIMO).

The comparison of the water temperature at 0.2 m depth, as predicted by the three models (SIMO, GOTM, and SCHISM) and the observed values are shown in Fig. 15. SIMO generally outperforms both models until the middle of September, when it starts to consistently overpredict the temperature, while GOTM shows quite low error. On the other hand, SCHISM tends to underestimate the temperature by approximately 5 C from the beginning of June and almost until the end of the year. The performance measures for the three models are summarized in Table 3, which also shows that SIMO outperforms the other two models considering, all measures except the bias, which is 0.5 C lower for GOTM.

Figure 15Comparison of the near-surface water temperature for SIMO, GOTM, and SCHISM for the period 1 January–27 December 2019.


Table 3Comparison of performance measures for SIMO, GOTM, and SCHISM for the period 1 January–27 December 2019.

Download Print Version | Download XLSX

Numerous lake modeling studies report quantitative performance measures. However, the comparison of model performance with other models is not always straightforward, as there is no common systematic approach. Namely, different studies report different performance measures; sometimes the calculation methods, the observation periods, and the measurement frequency and depths are not clearly stated or measurements are too rare to represent short-term variations. Furthermore, no studies calculating the performance measures in relation to the simulation period using only the end results, as done here, were found. However, quite a few studies report on single, longer simulations. Some of these results are summarized in Table 4. For a yearlong simulation of the water temperature in a small dimictic lake, Martynov et al. (2010) reported a surface temperature RMSE of 1.8 C for an eddy diffusivity model (Hostetler model) and 3.2 C for a two-layer model (FLake). Bruce et al. (2018) ran a two-year simulation for 32 different lakes using the GLM model, and the calculated RMSEs for the entire vertical profile, epilimnion, and hypolimnion were 1.34, 1.62, and 1.31 C, respectively. MacKay (2012) ran a bulk mixed model simulation for approximately a month and a half and reported a surface temperature MBE <1C. Read et al. (2014) ran a 30-year simulation (restarted annually) for 434 temperate lakes and reported a RMSE of 2.78, 1.74, and 3.33 C for the entire vertical profile, epilimnion, and hypolimnion, respectively. Moore et al. (2021) ran four different models for a temperate monomictic lake and reported RMSE values from 0.8 to 2.96 C for the runs before the model parameter calibration and 0.61 to 1.17 C after it. The reported absolute MBE values ranged from 0.34 to 1.75 C for the runs before the model parameter calibration and 0.1 to 0.55 C after it.

Table 4Comparison of SIMO performance with other models.

Download Print Version | Download XLSX

The yearlong simulation in this study resulted in a surface temperature RMSE of 1.51 C; in the hypolimnion, the RMSE was lowest in the deepest layer at 0.33 C. The RMSE was the highest in the thermocline region, where it reached a maximum of 2.8 C at 17 m depth. The RMSE for the entire profile was 1.91 C. The surface temperature MBE was 0.88 C. The maximum MBE was again in the thermocline region and equaled 2.28 C. This systematic overprediction can also be noticed in Fig. 14. Considering the lake surface temperature and entire vertical profile as well as the epilimnion and hypolimnion temperature, the model performance for the yearlong simulation was satisfactory, since it was comparable with the performances of other models (Table 4). The performance for the thermocline region was somewhat poorer, but performance in that region was not specifically reported in the reviewed literature.

7 Summary and conclusions

The aim of this study was to offer a simple 1-D energy budget model for the prediction of the vertical temperature profile in a small, warm, monomictic lake that is forced by a reduced number of input meteorological variables. Specifically, these include meteorological variables that are routinely measured at meteorological stations (i.e., the air temperature, relative humidity, atmospheric pressure, wind speed, and precipitation) as well as UVB radiation data and climatological yearly mean temperature data. In addition, an observed vertical profile of the lake temperature was used as an initial condition.

The main challenge was to calculate the net heat flux on the lake surface and to determine its components (i.e., the shortwave and longwave radiation, sensible and latent heat flux, and precipitation heat flux) from the available data. The model performance was evaluated using lake temperatures measured continuously during an observational campaign in two lakes of Plitvice Lakes, Croatia: Lake 1 (Prošće Lake) and Lake 12 (Kozjak Lake). The necessary meteorological data were provided by a single meteorological station located approximately 2 and 1.6 km from the lake temperature measuring points for Lake 1 and Lake 12, respectively. Except being further away from the meteorological station, Lake 1 has an altitude approximately 100 m higher than Lake 12, is surrounded by more complex orography, and is very likely exposed to stronger winds and lower air temperatures than those used as meteorological input data. Accordingly, the model performance was somewhat poorer for Lake 1, which indicates the importance of the microlocation-specific input meteorological data, as the meteorological forcing is the main driver of the temperature profile evolution. In addition, the influence of the tributary water that inflows into Lake 1, which was not taken into account in the present model, could also contribute to higher differences between the modeled and measured temperatures in comparison to Lake 12.

Generally, epilimnion temperature was somewhat overestimated, especially with the onset of the convective overturn. The upper limit of the metalimnion was well captured, while its thickness was overestimated, leading to underestimation in the maximal temperature gradient. However, the model satisfactorily estimated the stratification and overturn dynamics. There are several possible causes of departures of modeled from measured temperatures. One of them is the underestimation of the turbulent heat transfer in the epilimnion, especially in periods of high winds. In addition, the model cannot simulate internal seiches and possible water exchange between the warmer epilimnion and colder hypolimnion. Other probable causes are the use of an inappropriate light extinction coefficient value and the limited reliability of the surface heat flux. However, considering all model simplifications, we conclude that the model performed reasonably well.

The sensitivity analysis of the model performance to the simulation length showed that, when using appropriate meteorological forcing (as is the case of Lake 12), the model performance, especially in the epilimnion, steadily deteriorated up to a simulation length of approximately 15 d; however, a further increase in the simulation length up to 30 d had little effect on the model performance parameters. This proves the model can be used for obtaining reasonable lake water temperature prognosis for at least 30 d-long periods.

Despite the model's shortcomings, the yearlong simulation showed that the model is able to predict the onset of stratification and convective overturn relatively precisely, and the values of the model performance measures were comparable to those reported for other models. Thus, for a certain lake with no water temperature measurement data available, a yearlong simulation such as this would provide an assessment of lake stratification establishment, which can be useful for various studies dealing with lake biology, geochemistry, sedimentology, etc.

To further corroborate the general applicability of the present model, it should be applied to a larger number of different monomictic lakes. Nevertheless, in the present study, no lake-specific parameter tuning was performed. Thus, we expect similar model performance for other lakes if adequate meteorological forcing is employed.

Appendix A

Figure A1Error in the predicted water temperature (PiOi) for Lake 1 for different simulation lengths for the period between 6 July and 31 December 2019. Panel (a) is omitted so that the panels' positions for different simulation lengths correspond to those in Fig. 6.


Figure A2Error in the predicted water temperature (PiOi) for Lake 12 for different simulation lengths in the period between 6 July and 31 December 2019. Panel (a) is omitted so that the panels' positions for different simulation lengths correspond to those in Fig. 7.


Figure A3Observed (a) and predicted (b–f) vertical gradients of water temperature for Lake 1 for different simulation lengths in the period between 6 July and 31 December 2019.


Figure A4Observed (a) and predicted (b–f) vertical gradients of water temperature in Lake 12 for different simulation lengths in the period between 6 July and 31 December 2019.


Figure A5Error in the predicted vertical gradient of water temperature (PiOi) in Lake 1 for different simulation lengths in the period between 6 July and 31 December 2019. Panel (a) is omitted so that the panels' positions for different simulation lengths correspond to those in Fig. 3.


Figure A6Error in the predicted vertical gradient of water temperature (PiOi) in Lake 12 for different simulation lengths in the period between 6 July and 31 December 2019. Panel (a) is omitted so that the panels' positions for different simulation lengths correspond to those in Fig. 5.


Appendix B

B1 Description of the SCHISM and GOTM model parametrization

The hydrodynamic model system SCHISM (Semi-implicit Cross-scale Hydroscience Integrated System Model; Zhang et al., 2016) was designed for simulations of the 3-D baroclinic and barotropic circulation at different spatio-temporal scales. The calculation is conducted on the points of a horizontal unstructured triangular grid, which is one of the most important features of the model, because the use of such a grid ensures a high spatial resolution. In the calculations, the model uses an efficient and accurate semi-implicit method on finite elements or volumes with the Euler–Lagrange algorithm to solve the Reynolds-averaged Navier–Stokes equations (in hydrostatic and non-hydrostatic form) in order to more realistically described a wide spectrum of physical and biological processes as well as atmospheric and river forcing. The equations are simplified by considering the hydrostatic and Boussinesq approximations.

The horizontal grid covers the entire area of 16 cascade lakes, and it is composed of 17 472 elements whose surfaces range from 1 to 200 m2. In the vertical, a hybrid SZ grid was used, whereby the hybrid Z layers are fixed at a certain depth, located below the S coordinates that follow the terrain (Song and Haidvogel, 1994) according to the prescribed distance. The Plitvice Lakes are shallow enough that it is not necessary to define Z layers, and 30 sigma levels were used in the vertical discretization. During model calibration, i.e., when adjusting various parameters, it turned out that the model gives the best results in simulations with a time step of 10 s. Bottom friction in the model is approximated by the quadratic law of friction, defined by the assigned coefficient of friction with the adopted standard value of 0.003. As the Plitvice Lakes are extremely transparent and clean, in order to simulate a realistic lake character, Jerlov I was taken as the type of water, whose extinction coefficients correspond to the clear water. For the lake albedo, the theoretical values of 6 %, which are usual for the ocean, were applied. Vertical mixing in the model is imposed through the k-ε scheme with the Kantha–Clayson stability function. TVD (total variation diminishing) scheme was used in the advective terms of the transport equation. TVD is a slower scheme but better displays sharp temperature gradients. A baroclinic mode was also included, by which the contribution of temperature to the density of the medium is included in the equations of motion.

GOTM is a 1-D water column model for simulating the most important hydrodynamic and thermodynamic processes related to vertical mixing in natural waters. The GOTM model is suitable for simulating and predicting the stratification and vertical temperature profile of closed systems, such as the Plitvice Lakes. It is configured in such a way that it can be connected to 3-D circulation models, such as SCHISM, and used as a module to calculate vertical turbulent mixing. The core of the model calculates solutions for one-dimensional versions of the momentum, salt, and heat transport equations. The strength of GOTM is in the large number of tested turbulence models implemented in the code. Calculations are made at only one point along the entire vertical where any number of layers can be placed. Its advantage is in its faster performance and better formulation of the heat flow between the atmosphere and water.

Model parameters such as water type and turbulence scheme in GOTM are the same as in the SCHISM model. Jerlov I (clear water) was taken as the type of water, and the k-ε scheme was used as the turbulent mixing scheme. The number of vertical layers at point K1 (Lake 12) was set to 91, because it is at a greater depth, while 60 layers were taken for point P1 (Lake 1).

Apart from the measured air temperature and wind data (GOTM simulation) and measured air temperature (SCHISM simulation), both models use time series of atmospheric variables from the WRF model and heat fluxes on the surface of the lake, which are the main driver of the physical processes that cause thermal stratification and vertical mixing in the lake. The models are forced by atmospheric input on an hourly basis, with SCHISM additionally having an hourly loop that simulates the exchange of heat, mass, and momentum between the lake and the atmosphere and the consequent heating and mixing processes that occur in the lake.

Code and data availability

The SIMO v1.0 code is published under Creative Commons Attribution 4.0 International license and it is available at (Šarović and Klaić, 2021).

Lake water temperature data are available on request for research purposes by contacting Zvjezdana B. Klaić ( Authors are not authorized to publicly share meteorological data. To access these data, requests should be sent to the Croatian Meteorological and Hydrological Service.

Author contributions

KŠ designed the SIMO model, wrote the model code, ran the simulations, and performed the result post processing. ZBK organized the lake temperature experiment and contributed to the model development and the result evaluation and interpretation. Both KŠ and ZBK contributed to the discussion and writing of the paper. MB produced composite pictures of the lake bathymetries and the digital orthophoto images, ran the GOTM and SCHISM model simulations, performed corresponding post processing, and contributed to writing of Sects. 6 and 7.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This study was performed under the project “Hydrodynamic Modeling of Plitvice Lakes System” founded by the Plitvice Lakes National Park, Croatia. We also appreciate the technical support of the PLNP during equipment installation and lake temperature data acquisition. Meteorological data were provided by the Croatian Meteorological and Hydrological Service. A 2019 bathymetry experiment was organized by the PLNP. We are grateful to Goran Gašparac of Croatia Control Ltd., Velika Gorica for performing WRF model simulations. Finally, we gratefully acknowledge the useful comments of two anonymous reviewers.

Financial support

This study was supported by the project “Hydrodynamic Modeling of Plitvice Lakes System” (PLNP; grant no. 7989/16).

Review statement

This paper was edited by Bethanna Jackson and reviewed by two anonymous referees.


Andersen, T. K., Bolding, K., Nielsen, A., Bruggeman, J., Jeppesen, E., and Trolle, D.: How morphology shapes the parameter sensitivity of lake ecosystem models, Environ. Model. Softw., 136, 104945,, 2021. 

Bahr, A., Evans, C., Martinoli, A., Huwald, H., Higgins, C., and Parlange, M.: Measuring sensible heat flux with high spatial density, in 2012 IEEE Sensors Applications Symposium Proceedings, SAS 2012, Brescia, Italy, 7–9 February 2012, 255–260,, 2012. 

Bell, V. A., George, D. G., Moore, R. J., and Parker, J.: Using a 1-D mixing model to simulate the vertical flux of heat and oxygen in a lake subject to episodic mixing, Ecol. Modell., 190, 41–54,, 2006. 

Benson, B. B. and Krause Jr., D.: The concentration and isotopic fractionation of gases dissolved in freshwater in equilibrium with the atmosphere. 1. Oxigen, Limnol. Oceanogr., 25, 662–671,, 1980. 

Bristow, K. L. and Campbell, G. S.: On the relationship between incoming solar radiation and daily maximum and minimum temperature, Agric. For. Meteorol., 31, 159–166,, 1984. 

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,, 2018. 

Bruggeman, J. and Bolding, K.: A general framework for aquatic biogeochemical models, Environ. Model. Softw., 61, 249–265,, 2014. 

Brunel, J. P.: Estimation of sensible heat flux from measurements of surface radiative temperature and air temperature at two meters: Application to determine actual evaporation rate, Agric. For. Meteorol., 46, 179–191,, 1989. 

Brutsaert, W.: On a Derivable Formula for Long-Wave Radiation From Clear Skies, Water Resour. Res., 11, 742–744,, 1975. 

Bryan, A. M., Steiner, A. L., and Posselt, D. J.: Regional modeling of surface-atmosphere interactions and their impact on Great Lakes hydroclimate, J. Geophys. Res.-Atmos., 120, 1044–1064,, 2015. 

Bubalo, M., Janeković, I., and Orlić, M.: Chrystal and Proudman resonances simulated with three numerical models, Ocean Dynam., 68, 97–507,, 2018. 

Burchard, H., Bolding, K., and Villarreal, M. R.: GOTM, a general ocean turbulence model. Theory, implementation and test cases, European Commission, Space Applications Institute, 103 pp., (last access: 30 August 2021), 1999. 

Burchard, H., Schulz, E., and Schuttelaars, H. M.: Impact of estuarine convergence on residual circulation in tidally energetic estuaries and inlets, Geophys. Res. Lett., 41, 913–919,, 2014. 

Burić, M., Grgurić, S., Mikulčić, H., and Wang, X.: A numerical investigation of tidal current energy resource potential in a sea strait, Energy, 234, 121241,, 2021. 

Chen, C. A. and Millero, F. J.: Thermodynamic Properties for Natural Waters Covering Only the Limnological Range, Limnol. Oceanogr., 31, 657–662,, 1986. 

Ciglenečki, I., Janeković, I., Marguš, M., Bura-Nakić, E., Carić, M., Ljubešić, Z., Batistić, M., Hrustić, E., Dupčić, J., and Garić, R.: Impacts of extreme weather events on highly eutrophic marine ecosystem (Rogoznica Lake, Adriatic coast), Cont. Shelf Res., 108, 144–155,, 2015. 

Crawford, T. M. and Duchon, C. E.: An improved parameterization for estimating effective atmospheric emissivity for use in calculating daytime downwelling longwave radiation, J. Appl. Meteorol., 38, 474–480,<0474:AIPFEE>2.0.CO;2, 1999. 

Forcat, F., Roget, E., Figueroa, M., and Sánchez, X.: Earth rotation effects on the internal wave field in a stratified small lake: Numerical simulations, Limnetica, 30, 27–42,, 2011. 

Frishfelds, V., Sennikovs, J., Bethers, U., Murawski, J., and Timuhins, A.: Modeling transit flow through port gates and connecting channel in Baltic Sea-Liepaja Port-Liepaja Lake system, Front. Mar. Sci., 8, 657721,, 2021. 

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., 107, 3230,, 2002. 

Henderson-Sellers, B.: New formulation of eddy diffusion thermocline models, Appl. Math. Model., 9, 441–446,, 1985. 

Henderson-Sellers, B.: Calculating the surface energy balance for lake and reservoir modeling: A review, Rev. Geophys., 24, 625–649,, 1986. 

Hipsey, M. R., Bruce, L. C., Boon, C., Busch, B., Carey, C. C., Hamilton, D. P., Hanson, P. C., Read, J. S., de Sousa, E., Weber, M., and Winslow, L. A.: A General Lake Model (GLM 3.0) for linking with high-frequency sensor data from the Global Lake Ecological Observatory Network (GLEON), Geosci. Model Dev., 12, 473–523,, 2019. 

Hostetler, S. W. and Bartlein, P. J.: Simulation of lake evaporation with application to modeling lake level variations of Harney-Malheur Lake, Oregon, Water Resour. Res., 26, 2603–2612,, 1990. 

Jacob, B., Stanev, E. V., and Zhang, Y. J.: Local and remote response of the North Sea dynamics to morphodynamic changes in theWadden Sea, Ocean Dynam., 66, 671–690,, 2016. 

Jassby, A. and Powell, T.: Vertical patterns of eddy diffusion during stratification in Castle Lake, California, Limnol. Oceanogr., 20, 530–543,, 1975. 

Klaić, Z. B. and Kvakić, M.: Modeling the impacts of the man-made lake on the meteorological conditions of the surrounding areas, J. Appl. Meteorol. Climatol, 53, 1121–1142,, 2014. 

Klaić, Z. B., Rubinić, J., and Kapelj, S.: Review of research on Plitvice Lakes, Croatia in the fields of meteorology, climatology, hydrology, hydrogeochemistry and physical limnology, Geofizika, 35, 189–278,, 2018. 

Klaić, Z. B., Babić, K., and Mareković, T.: Internal seiches in a karstic mesotrophic lake (Prošće, Plitvice Lakes, Croatia), Geofizika, 37, 157–179,, 2020a. 

Klaić, Z. B., Babić, K., and Orlić, M.: Evolution and dynamics of the vertical temperature profile in an oligotrophic lake, Hydrol. Earth Syst. Sci., 24, 3399–3416,, 2020b. 

Kristovich, D. R., Clark, R. D., Frame, J., Geerts, B., Knupp, K. R., Kosiba, K. A., Laird, N. F., Metz, N. D., Minder, J. R., Sikora, T. D., Steenburgh, W. J., Steiger, S. M., Wurman, J., and Young, G. S.: The Ontario winter lake-effect systems field campaign: Scientific and educational adventures to further our knowledge and prediction of lake-effect storms, B. Am. Meteorol. Soc., 98, 315–332,, 2017. 

Krumgalz, B. S.: Temperature Dependence of Mineral Solubility in Water. Part 3. Alkaline and Alkaline Earth Sulfates, J. Phys. Chem. Ref. Data, 47, 023101,, 2018. 

Kudish, A. I. and Evseev, E.: Statistical relationships between solar UVB and UVA radiation and global radiation measurements at two sites in Israel, Int. J. Climatol., 20, 759–770,<759::AID-JOC494>3.0.CO;2-K, 2000. 

Kudish, A. I., Lyubansky, V., Evseev, E. G., and Ianetz, A.: Inter-comparison of the solar UVB, UVA and global radiation clearness and UV indices for Beer Sheva and Neve Zohar (Dead Sea), Israel, Energy, 30, 1623–1641,, 2005. 

Ladwig, R., Hanson, P. C., Dugan, H. A., Carey, C. C., Zhang, Y., Shu, L., Duffy, C. J., and Cobourn, K. M.: Lake thermal structure drives interannual variability in summer anoxia dynamics in a eutrophic lake over 37 years, Hydrol. Earth Syst. Sci., 25, 1009–1032,, 2021. 

Li, Q., Bruggeman, J., Burchard, H., Klingbeil, K., Umlauf, L., and Bolding, K.: Integrating CVMix into GOTM (v6.0): a consistent framework for testing, comparing, and applying ocean mixing schemes, Geosci. Model Dev., 14, 4261–4282,, 2021. 

Liston, G. E. and Hall, D. K.: An energy-balance model of lake-ice evolution, J. Glaciol., 41, 373–382,, 1995. 

MacKay, M. D.: A process-oriented small lake scheme for coupled climate modeling applications, J. Hydrometeorol., 13, 1911–1924,, 2012. 

MacKay, M. D., Verseghy, D. L., Fortin, V., and Rennie, M. D.: Wintertime simulations of a boreal lake with the Canadian Small Lake Model, J. Hydrometeorol., 18, 2143–2160,, 2017. 

Martynov, A., Sushama, L., and Laprise, R.: Simulation of temperate freezing lakes by one-dimensional lake models: performance assessment for interactive coupling with regional climate models, Boreal Env. Res, 15, 143–164, 2010. 

Mironov, D., Heise, E., Kourzeneva, E., Ritter, B., Schneider, N., and Terzhevik, A.: Implementation of the lake parameterisation scheme FLake into the numerical weather prediction model COSMO, Boreal Environ. Res., 15, 218–230, 2010. 

Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surfacelayer of the atmosphere, Tr. Akad. Nauk SSSR Geophiz. Inst., 151, 163–187, 1954. 

Moore, T. N., Mesman, J. P., Ladwig, R., Feldbauer, J., Olsson, F., Pilla, R. M., Shatwell, T., Venkiteswaran, J. J., Delany, A. D., Dugan, H., Rose, K. C., and Read, J. S.: LakeEnsemblR: An R package that facilitates ensemble modelling of lakes, Environ. Model. Softw., 143, 105101,, 2021. 

National Park Plitvička Jezera (NPPL): Historical overview,, last access: 24 January 2021. 

Nielsen, A., Hu, F. R. S., Schnedler-Meyer, N. A., Bolding, K., Andersen, T. K., and Trolle, D.: Introducing QWET – A QGIS-plugin for application, evaluation and experimentation with the WET model Environmental Modelling and Software, Environ. Modell. Softw., 135, 104886,, 2021. 

NIWA: Latent and sensible heat fluxes from lake water surfaces,, last access: 28 January 2021. 

Pashiardis, S., Kalogirou, S., and Pelengaris, A.: Statistical Analysis and Inter- Comparison of Solar UVB and Global Radiation for Athalassa and Larnaca, Cyprus, SM J. Biometrics Biostat., 2, 1012,, 2017. 

Podstawczynska, A.: UV and global solar radiation in Łódź, Central Poland, Int. J. Climatol., 30, 1–10,, 2009. 

Pokhrel, R. P. and Bhattarai, B. K.: Relation between Global Solar Radiation and Solar Ultraviolet Radiation in Different Parts of Nepal, J. Inst. Eng., 8, 169–175,, 2012. 

Quay, P. D., Broecker, W. S., Hesslein, R. H., and Schindler, D. W.: Vertical diffusion rates determined by tritium tracer experiments in the thermocline and hypolimnion of two lakes, Limnol. Ocean., 25, 201–218, 1980. 

Rasconi, S., Winter, K., and Kainz, M. J.: Temperature increase and fluctuation induce phytoplankton biodiversity loss – Evidence from a multi-seasonal mesocosm experiment, Ecol. Evol., 7, 2936–2946,, 2017. 

Read, J. S., Winslow, L. A., Hansen, G. J. A., Van Den Hoek, J., Hanson, P. C., Bruce, L. C., and Markfort, C. D.: Simulating 2368 temperate lakes reveals weak coherence in stratification phenology, Ecol. Modell., 291, 142–150,, 2014. 

Rodi, W.: Turbulence models and their application in hydraulics – A state of the art review, State-of-the-Art Paper/International Association for Hydraulic Research, Delft, The Netherlands, OCLC Number/Unique Identifier: 1069725439, 104 pp., 1984. 

Šarović, K. and Klaić, Z.: SIMO v1.0: Simplified model of the vertical temperature profile in a small warm monomictic lake (v1.0), Zenodo [code and data set],, 2021. 

Skamarock, W. C. and Klemp, J. B: A time-split nonhydrostatic atmospheric model for weather research and forecasting applications. J. Comput. Phys. 227, 3465–3485,, 2008. 

Song, Y. and Haidvogel, D.: A semi-implicit ocean circulation model using a generalized topography-following coordinate system, J. Comput. Phys., 115, 228–244,, 1994. 

Stefan, H. G., Fang, X., and Hondzo, M.: Simulated climate change effects on year-round water temperatures in temperate zone lakes, Climatic Change, 40, 547–576,, 1998. 

Stepanenko, V. M., Machulskaya, E. E., Glagolev, M. V., and Lykossov, V. N.: Numerical modeling of methane emissions from lakes in the Permafrost Zone, Izv. – Atmos. Ocean Phys., 47, 252–264,, 2011. 

Stepanenko, V. M., Martynov, A., Jöhnk, K. D., Subin, Z. M., Perroud, M., Fang, X., Beyrich, F., Mironov, D., and Goyette, S.: A one-dimensional model intercomparison study of thermal regime of a shallow, turbid midlatitude lake, Geosci. Model Dev., 6, 1337–1352,, 2013. 

Stepanenko, V., Jöhnk, K. D., Machulskaya, E., Perroud, M., Subin, Z., Nordbo, A., Mammarella, I., and Mironov, D.: Simulation of surface energy fluxes and stratification of a small boreal lake by a set of one-dimensional models, Tellus A: Dynamic Meteorology and Oceanography, 66, p. 21389,, 2014. 

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,, 2016. 

Sun, S., Yan, J., Xia, N., and Sun, C.: Development of a Model for Water and Heat Exchange Between the Atmosphere and a Water Body, Adv. Atmos. Sci., 24, 927–938,, 2007. 

Thiery, W., Martynov, A., Darchambeau, F., Descy, J.-P., Plisnier, P.-D., Sushama, L., and van Lipzig, N. P. M.: Understanding the performance of the FLake model over two African Great Lakes, Geosci. Model Dev., 7, 317–337,, 2014. 

Vachon, D., Langenegger, T., Donis, D., and McGinnis, D. F.: Influence of water column stratification and mixing patterns on the fate of methane produced in deep sediments of a small eutrophic lake, Limnol. Oceanogr., 64, 2114–2128,, 2019. 

Verburg, P. and Antenucci, J. P.: Persistent unstable atmospheric boundary layer enhances sensible and latent heat loss in a tropical great lake: Lake Tanganyika, J. Geophys. Res.-Atmos., 115, D11109,, 2010. 

Råman Vinnå, L., Medhaug, I., Schmid, M., and Bouffard, D.: The vulnerability of lakes to climate change along an altitudinal gradient, Commun. Earth Environ., 2, 35,, 2021. 

Wald, L.: Basics in Solar radiation at Earth's surface, Lecture Notes, 1st edn., MINES ParisTech, PSL Research University O.I.E. – Observation, Impacts, Energy Center Sophia Antipolis, France, hal-02164311, (last access: 24 January 2021), 2019. 

Wang, K. and Dickinson, R. E.: Global Atmospheric Downward Longwave Radiation at the Surface From Ground-Based Observations, Satellite Retrievals, and Reanalyses, Rev. Geophys., 51, 150–185,, 2013. 

Willmott, C. J. and Matsuura, K.: Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance, Clim. Res., 30, 79–82,, 2005. 

Willmott, C. J., Robeson, S. M., and Matsuura, K.: A refined index of model performance, Int. J. Climatol., 32, 2088–2094,, 2012. 

Winslow, J. C., Hunt, E. R., and Piper, S. C.: A globally applicable model of daily solar irradiance estimated from air temperature and precipitation data, Ecol. Modell., 143, 227–243,, 2001. 

Wu, Y., Huang, A., Yang, B., Dong, G., Wen, L., Lazhu, Zhang, Z., Fu, Z., Zhu, X., Zhang, X., and Cai, S.: Numerical study on the climatic effect of the lake clusters over Tibetan Plateau in summer, Clim. Dynam., 53, 5215–5236,, 2019. 

Wu, Y., Huang, A., Lazhu, Yang, X., Qiu, B., Wen, L., Zhang, Z., Fu, Z., Zhu, X., Zhang, X., Cai, S., and Tang, Y.: Improvements of the coupled WRF-Lake model over Lake Nam Co, Central Tibetan Plateau, Clim. Dynam., 55, 2703–2724,, 2020. 

Zhang, Y. J., Ye, F., Stanev, E.V., and Grashorn, S.: Seamless cross-scale modeling with SCHISM, Ocean Model., 102, 64–81,, 2016.  

Zhang, Y. J., Ye, F., Yu, H., Sun, W., Moghimi, S., Myers, E., Nunez, K., Zhang, R., Wang, H., Roland, A., Du, J., and Liu, Z.: Simulating compound flooding events in a hurricane, Ocean Dynam., 70, 621–640,, 2020. 

Short summary
We develop a simple 1-D model for the prediction of the vertical temperature profiles in small, warm lakes. The model uses routinely measured meteorological variables as well as UVB radiation and yearly mean temperature data. It can be used for the assessment of the onset and duration of lake stratification periods when water temperature data are unavailable, which can be useful for various lake studies performed in other scientific fields, such as biology, geochemistry, and sedimentology.