**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ć

^{1},

^{2},

^{1}

**Kristina Šarović et al.**Kristina Šarović Melita Burić and Zvjezdana B. Klaić

^{1},

^{2},

^{1}

^{1}Department of Geophysics, Faculty of Science, University of Zagreb, Zagreb, 10000, Croatia^{2}Gekom Geophys & Ecol Modelling Ltd, Zagreb, 10000, Croatia

^{1}Department of Geophysics, Faculty of Science, University of Zagreb, Zagreb, 10000, Croatia^{2}Gekom Geophys & Ecol Modelling Ltd, Zagreb, 10000, Croatia

**Correspondence**: Kristina Šarović (kristina.sarovic@gmail.com)

**Correspondence**: Kristina Šarović (kristina.sarovic@gmail.com)

Received: 12 Apr 2021 – Discussion started: 16 Jun 2021 – Revised: 03 Aug 2022 – Accepted: 19 Sep 2022 – Published: 18 Nov 2022

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.

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 CO_{2}) 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.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).

## 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.20 ^{∘}C for
temperatures between 0 and 70 ^{∘}C and ±0.25 ^{∘}C
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.

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.

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

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:

where *c*_{p} 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), *k*_{m} and *k*_{t} are the
molecular and turbulent thermal conductivity (W m^{−1} K^{−1}), *Φ* is
the heat flux (W m^{−2}), and *M*_{conv} is the convective mixing term (W m^{−3}).

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

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:

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, *Pr*_{0}=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):

where *φ* is the latitude and *U*_{2} is the wind speed at 2 m above
the water surface (m s^{−1}), and

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

The wind speed *U*_{2} is determined from the logarithmic formula:

where *z*_{0} is the roughness length (m). The air shear velocity *u*^{*} and the
roughness length *z*_{0} 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 (*S*_{n}), net longwave radiation (*L*_{n}), sensible heat flux
(*H*_{s}), latent heat flux (*H*_{l}), and heat flux brought by precipitation
(*H*_{p}). The surface boundary condition can be written as follows:

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 $\mathrm{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 $\mathrm{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

where *S*_{surf} 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, rh${}_{{T}_{\mathrm{max}}}$ is the relative humidity
at the moment when daily maximum air temperature (*T*_{max}) is reached, and
*S*_{top} 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):

where *S*_{0}=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
(*c*_{elev}):

To calculate *τ*_{0}, *τ*_{a}, *τ*_{v}, *c*_{elev}, *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:

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:

where *T*_{mean} 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 (*c*_{elev}) is calculated as
follows:

where ${z}_{\mathrm{a}.\mathrm{s}.\mathrm{l}.}$ is the site elevation (m).

From Eq. (11), *τ*_{cf} *S*_{top} is the
maximum cloud-free value of *S*_{surf}. The effect of cloudiness is
indirectly taken into account by introducing the factor *D*(1−*β*rh${}_{{T}_{\mathrm{max}}}$). This is based on the finding that the solar irradiation
from sunrise, when minimum humidity is expected (rh${}_{{T}_{\mathrm{min}}}\approx \mathrm{1}$),
until the maximum daily air temperature (and minimum humidity
rh${}_{{T}_{\mathrm{max}}}$) is reached, is proportional to the decline of the relative
humidity, ${S}_{\text{surf\_}{T}_{\mathrm{max}}}\propto $ (1−*β*rh${}_{{T}_{\mathrm{max}}}$). The factor $D={S}_{\mathrm{surf}}/{S}_{\text{surf\_}{T}_{\mathrm{max}}}$ 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:

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

where Δ*T*_{m} 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:

where *S*_{surf} and UVB_{day} are the daily values (J m^{−2} d^{−1})
and *S*(*h*) and UVB(*h*) are the mean values (W m^{−2}) for the *h*^{th} 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:

where *S*_{n}(*z*) is the net shortwave radiation at water depth *z* (W m^{−2}),
*α*=0.06 is the water surface albedo, and *K*_{1}, *K*_{2}, and *K*_{3}
are empirical constants. *K*_{1} corresponds to the light extinction
coefficient *λ*_{e}=0.1 (value of 0.1 is appropriate for clear
oligotrophic lakes). *K*_{2} is calculated as

where *β*=0.4 accounts for the absorption in the surface layer, and
*z*_{A}=0.6 m is the depth of the surface absorption layer, where the
exponential decay starts. The third parameter, *K*_{3}=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 (${L}_{\mathrm{a}}^{\downarrow}$) and the outgoing upwards radiation from the lake surface (${L}_{\mathrm{s}}^{\uparrow}$). As direct measurement data of longwave radiation by pyrgeometers are not routinely available, longwave radiation may be calculated using the following formula:

where *r* is the water reflectivity for longwave radiation, *ε* and
*ε*_{a} are the emissivities of the lake surface and the
atmosphere, respectively; *T*_{s} is the water surface temperature
(^{∘}C), *T*_{a} 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:

where *e*_{a} is the water vapor pressure (hPa), which is related to the
relative humidity (rh) and saturation vapor pressure (*e*_{s}):

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

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:

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

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:

### 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 *C*_{H} and *C*_{E} are the transfer coefficients for sensible and
latent heat flux, respectively; *c*_{a}=1005 J kg^{−1} K^{−1}
is the specific heat of air, *L*_{v}≈2500 kJ kg^{−1} is the
latent heat of evaporation, *ρ*_{a} is the air density (kg m^{−3}),
and *q*_{s} and *q*_{a} 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 *C*_{D} is the drag coefficient, *h* is the height above ground (m),
*z*_{0} and *z*_{E} 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 $\mathit{\zeta}=h/L$, where *L* is the Monin–Obukhov length:

where *T*_{v} is the virtual temperature. Obviously, *L* depends on *H*_{s} and
*H*_{l}, while *H*_{s} and *H*_{l} depend on the stability of the atmosphere.
Therefore, to calculate *H*_{s} and *H*_{l}, an iterative procedure has to be
used. The procedure is initiated by assuming neutral conditions (${\mathit{\psi}}_{\mathrm{M}}={\mathit{\psi}}_{\mathrm{E}}=\mathrm{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

where *T*_{1} and *T*_{1+p} represent the water temperature of the first
layer before and after the precipitation has been introduced
(^{∘}C), *T*_{prec} is the precipitation temperature (^{∘}C),
Δ*z*_{1} is the thickness of the first layer (m), and *P* is the hourly
precipitation (mm h^{−1}). The heat flux brought in by precipitation
*H*_{p} (W m^{−2}) can then be calculated as

Since *T*_{prec} 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:

where Δ*z*_{j} and Δ*z*_{j+1} represent the thickness of
the *j*th (upper) and (*j*+1)th (lower) layers, respectively; *T*_{j} and
*T*_{j+1} are the water temperatures in these layers before convective
mixing, respectively; *T*_{j_new} and
*T*_{j+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:

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:

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

Then, the solution for *T*^{n+1} is obtained as follows:

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.

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 (IA_{orig}),
modified (IA_{mod}), and refined (IA_{ref}) 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, IA_{mod} is less sensitive to
outliers than IA_{orig}. IA_{mod} approaches 1 (perfect agreement)
more slowly than IA_{orig}, which means that IA_{mod} is more
conservative and allows for finer comparisons of different models with
relatively good performance. In IA_{ref} (Eq. 48), the prediction variability in the denominator
is replaced with the observation variability. IA_{orig} and IA_{mod}
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. IA_{ref} 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 IA_{ref}
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.

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.

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 $\le \pm \mathrm{2}$ and $\le \pm \mathrm{1}$ ^{∘}C for Lake 1 and Lake 12, respectively, except in the
thermocline region, where they were mainly $\le \pm \mathrm{4}$ and
$\le \pm \mathrm{2}$ ^{∘}C 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.1 ^{∘}C 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.4 ^{∘}C 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.4 ^{∘}C 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.

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

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.

To compare the performance of the proposed model with the performance of more complex models, we applied 1-D General Ocean Turbulent Model (GOTM; https://gotm.net/about/, 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.

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 *<*1 ^{∘}C.
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.

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.

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.

## 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 m^{2}. 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.

The SIMO v1.0 code is published under Creative Commons Attribution 4.0 International license and it is available at https://doi.org/10.5281/zenodo.6367810 (Šarović and Klaić, 2021).

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

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.

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.

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

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, https://doi.org/10.1016/j.envsoft.2020.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, https://doi.org/10.1109/SAS.2012.6166293, 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, https://doi.org/10.1016/j.ecolmodel.2005.02.025, 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, https://doi.org/10.4319/lo.1980.25.4.0662, 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, https://doi.org/10.1016/0168-1923(84)90017-0, 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, https://doi.org/10.1016/j.envsoft.2017.11.016, 2018.

Bruggeman, J. and Bolding, K.: A general framework for aquatic biogeochemical models, Environ. Model. Softw., 61, 249–265, https://doi.org/10.1016/j.envsoft.2014.04.002, 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, https://doi.org/10.1016/0168-1923(89)90063-4, 1989.

Brutsaert, W.: On a Derivable Formula for Long-Wave Radiation From Clear Skies, Water Resour. Res., 11, 742–744, https://doi.org/10.1029/wr011i005p00742, 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, https://doi.org/10.1002/2014JD022316, 2015.

Bubalo, M., Janeković, I., and Orlić, M.: Chrystal and Proudman resonances simulated with three numerical models, Ocean Dynam., 68, 97–507, https://doi.org/10.1007/s10236-018-1146-8, 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., https://op.europa.eu/en/publication-detail/-/publication/5b512e12-367d-11ea-ba6e-01aa75ed71a1/language-en/format-PDF/source-272420379 (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, https://doi.org/10.1002/2013GL058494, 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, https://doi.org/10.1016/j.energy.2021.121241, 2021.

Chen, C. A. and Millero, F. J.: Thermodynamic Properties for Natural Waters Covering Only the Limnological Range, Limnol. Oceanogr., 31, 657–662, https://doi.org/10.4319/lo.1986.31.3.0657, 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, https://doi.org/10.1016/j.csr.2015.05.007, 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, https://doi.org/10.1175/1520-0450(1999)038<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, https://doi.org/10.23818/limn.30.04, 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, https://doi.org/10.3389/fmars.2021.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,
https://doi.org/10.1029/2001jc000954, 2002.

Henderson-Sellers, B.: New formulation of eddy diffusion thermocline models, Appl. Math. Model., 9, 441–446, https://doi.org/10.1016/0307-904X(85)90110-6, 1985.

Henderson-Sellers, B.: Calculating the surface energy balance for lake and reservoir modeling: A review, Rev. Geophys., 24, 625–649, https://doi.org/10.1029/RG024i003p00625, 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, https://doi.org/10.5194/gmd-12-473-2019, 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, https://doi.org/10.1029/WR026i010p02603, 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, https://doi.org/10.1007/s10236-016-0949-8, 2016.

Jassby, A. and Powell, T.: Vertical patterns of eddy diffusion during stratification in Castle Lake, California, Limnol. Oceanogr., 20, 530–543, https://doi.org/10.4319/lo.1975.20.4.0530, 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, https://doi.org/10.1175/JAMC-D-13-0163.1, 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, https://doi.org/10.15233/gfz.2018.35.9, 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, https://doi.org/10.15233/gfz.2020.37.11, 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, https://doi.org/10.5194/hess-24-3399-2020, 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, https://doi.org/10.1175/BAMS-D-15-00034.1, 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, https://doi.org/10.1063/1.5031951, 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, https://doi.org/10.1002/1097-0088(20000615)20:7<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, https://doi.org/10.1016/j.energy.2004.04.033, 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, https://doi.org/10.5194/hess-25-1009-2021, 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, https://doi.org/10.5194/gmd-14-4261-2021, 2021.

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

MacKay, M. D.: A process-oriented small lake scheme for coupled climate modeling applications, J. Hydrometeorol., 13, 1911–1924, https://doi.org/10.1175/JHM-D-11-0116.1, 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, https://doi.org/10.1175/JHM-D-16-0268.1, 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, https://doi.org/10.1016/j.envsoft.2021.105101, 2021.

National Park Plitvička Jezera (NPPL): Historical overview, https://np-plitvicka-jezera.hr/en/scientific-research/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, https://doi.org/10.1016/j.envsoft.2020.104886, 2021.

NIWA: Latent and sensible heat fluxes from lake water surfaces, https://niwa.co.nz/our-services/software/heat-fluxes-from-lakes/, 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, https://doi.org/10.36876/smjbb.1006, 2017.

Podstawczynska, A.: UV and global solar radiation in Łódź, Central Poland, Int. J. Climatol., 30, 1–10, https://doi.org/10.1002/joc.1864, 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, https://doi.org/10.3126/jie.v8i3.5942, 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, https://doi.org/10.1002/ece3.2889, 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, https://doi.org/10.1016/j.ecolmodel.2014.07.029, 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], https://doi.org/10.5281/zenodo.6367810, 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, https://doi.org/10.1016/j.jcp.2007.01.037, 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, https://doi.org/10.1006/jcph.1994.1189, 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, https://doi.org/10.1023/A:1005371600527, 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, https://doi.org/10.1134/S0001433811020113, 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, https://doi.org/10.5194/gmd-6-1337-2013, 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, https://doi.org/10.3402/tellusa.v66.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, https://doi.org/10.5194/gmd-9-1977-2016, 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, https://doi.org/10.1007/s00376-007-0927-7, 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, https://doi.org/10.5194/gmd-7-317-2014, 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, https://doi.org/10.1002/lno.11172, 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, https://doi.org/10.1029/2009JD012839, 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, https://doi.org/10.1038/s43247-021-00106-w, 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, https://hal-mines-paristech.archives-ouvertes.fr/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, https://doi.org/10.1002/rog.20009, 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, https://doi.org/10.3354/cr030079, 2005.

Willmott, C. J., Robeson, S. M., and Matsuura, K.: A refined index of model performance, Int. J. Climatol., 32, 2088–2094, https://doi.org/10.1002/joc.2419, 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, https://doi.org/10.1016/S0304-3800(01)00341-6, 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, https://doi.org/10.1007/s00382-019-04856-4, 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, https://doi.org/10.1007/s00382-020-05402-3, 2020.

Zhang, Y. J., Ye, F., Stanev, E.V., and Grashorn, S.: Seamless cross-scale modeling with SCHISM, Ocean Model., 102, 64–81, https://doi.org/10.1016/j.ocemod.2016.05.002, 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, https://doi.org/10.1007/s10236-020-01351-x, 2020.

- Abstract
- Introduction
- Study area and measurements
- Model description and governing equations
- Measures of the model performance
- Results and discussion
- Comparison with other models
- Summary and conclusions
- Appendix A
- Appendix B
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Study area and measurements
- Model description and governing equations
- Measures of the model performance
- Results and discussion
- Comparison with other models
- Summary and conclusions
- Appendix A
- Appendix B
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References