Mid-depth temperature maximum in an estuarine lake

The mid-depth temperature maximum (TeM) was measured in an estuarine Bol’shoi Vilyui Lake (Kamchatka peninsula, Russia) in summer 2015. We applied 1D k–ε model LAKE to the case, and found it successfully simulating the phenomenon. We argue that the main prerequisite for mid-depth TeM development is a salinity increase below the freshwater mixed layer, sharp enough in order to increase the temperature with depth not to cause convective mixing and double diffusion there. Given that this condition is satisfied, the TeM magnitude is controlled by physical factors which we identified as: radiation absorption below the mixed layer, mixed-layer temperature dynamics, vertical heat conduction and water-sediments heat exchange. In addition to these, we formulate the mechanism of temperature maximum ‘pumping’, resulting from the phase shift between diurnal cycles of mixed-layer depth and temperature maximum magnitude. Based on the LAKE model results we quantify the contribution of the above listed mechanisms and find their individual significance highly sensitive to water turbidity. Relying on physical mechanisms identified we define environmental conditions favouring the summertime TeM development in salinity-stratified lakes as: small-mixed layer depth (roughly, ~< 2 m), transparent water, daytime maximum of wind and cloudless weather. We exemplify the effect of mixed-layer depth on TeM by a set of selected lakes.


Introduction
A remarkable feature of estuaries and some saline lakes is the presence of a freshwater layer on the top of a salty water mass. In some cases it leads to development of a mid-depth or bottom temperature maximum (TeM), found in natural lakes (Kalecsinsky 1902, Egorov 1991, Wetzel 2001) (called 'heliotermic'), some of the icecovered lakes (Spigel andPriscu 1996, Kirillin andTerzhevik 2011), and artificial 'solar ponds' (Velmurugan and Srithar 2008). The common feature of all these water bodies is that the stratification below the freshwater mixed layer is controlled by salinity gradient, rather than by temperature gradient. Many of such lakes are meromictic lakes, i.e. the lower stable water layer persists permanently, unmixed with the water masses circulating above .
Although at least 177 meromictic lakes worldwide are reported in literature (Hall and Northcote 2012), the TeM phenomenon has been observed on a small subset of these water objects. Hence, the question arises, what are the conditions necessary for TeM development, and what are the factors causing its magnitude?
The temperature maximum in metalimnion may pose a pronounced effect on a lake ecosystem. The pycnocline depth defines to a large extent the position of chorophyll concentration maximum (Cantin et al 2011), that is usually located inside this layer. This maximum often causes a major contribution to the whole lake's primary production (e.g. Camacho (2006) and references therein). Concomitantly, the phytoplankton primary production is exponentially dependent on water temperature (Eppley 1972, Quiros 1988, Kishi et al 2007. Therefore, metalimnetic temperature is Figure 1. The Estuary of Bol'shoi Vilyui River and its hydrological zones: I-a-riverine water masses; I-b-marine water masses; II-mixing zone; III-secondary water masses (secondary water masses are formely riverine of marine masses, which thermodynamic properties have been transformed); 1-river discharge direction (two-way arrows denote tidal currents); 2-the measurement raft location in 2015; 3-a line along which a vertical cross-section of the water body is depicted in a lower inset; 2.3 km... etc.-the distance from the estuary mouth; 6.6 km is the location of Vilyuchinsk salmon hatchery, where on-shore measurements were carried out; on the lower inset-vertical cross-section of the estuary and Bol'shoi Vilyui lake, made along a line 3 (distance from the estuary mouth is provided at horizontal axis, the water level is given for a high tidal phase).
important for the primary production, and for the lake ecosystem in general. In turn, ecosystem productivity causes the trophic status of the lake, which is correlated to greenhouse gas production (Casper 1992), important actors in the climate system. So far, mathematical models simulating the deepened temperature maximum in lakes stratified by salinity have been developed either specifically for ice-cover conditions (Shirtcliffe 1964, Kirillin andTerzhevik 2011) or designed to reproduce artificial solar ponds (Rabl and Nielsen 1975, Atkinson and Harleman 1983, Bernad et al 2013. The bulk models such as presented in Egorov and Zilitinkevich (1999) do not reproduce vertical temperature profile explicitly.
Here we use a more general model LAKE, that is a representative of k-family broadly applied either for freshwater lakes (Jöhnk and Umlauf 2001, Perroud et al 2009, Stepanenko et al 2014 or for marine boundary layers (Burchard and Petersen 1999). It is a one-dimensional model, explicitly reproducing the vertical distribution of velocity, temperature, salinity and other thermodynamic and chemical quantities. The model is validated using observation data collected on the estuarine Bol'shoi Vilyui Lake (Kamchatka peninsula, Russia) in July 2015. The small magnitude of tides in this lake allows one to neglect their dynamic and advective effect at timescales considered. Based on the model simulations, we elucidate the contribution of different physical mechanisms to the temperature maximum development, and discuss implications of this analysis for other lakes, exhibiting TeM phenomenon.

The site and measurements
Bol'shoi Vilyui Lake (52.825 • N, 158.55 • E) is located at the Eastern coast of Kamchatka peninsula (Russia) (figure 1). Together with a smaller Malyi Vilyui Lake it forms an estuary, connecting to the ocean through a small channel. The area of Bol'shoi Vilyui is 4.3 km 2 , the mean depth is 3 m. The lake depth increases from 1.5-2 m at its southern part to 6.5-7 m in the northwestern part (figure 1). The estuarine channel is of ca. 2.5 km length, 50 to 200-300 m width, and with depth of 1 m at the low tidal phase. The tide magnitude is 1-2 m at the ocean coast. In the lake, this magnitude is significantly attenuated, down to 0.01-0.03 m in the north-western part. The major inlet to Bol'shoi Vilyui Lake is Bol'shoi Vilyui River, where the mineralization does not exceed 20-40 mg l −1 . Due to sources of both freshwater and seawater, the lake body consists of mainly two water masses: the surface freshwater and the saline one below. This causes the existence of sharp pycnocline, its top locating at ∼1.5 m depth in summer (Gorin 2013). While the freshwater is supplied to the lake continuously, it is suggested that the marine water enters the lake trough mainly during the winter storms (Gorin 2013). The salinity of seawater varies from 25‰ in summer to 32‰ in winter. The climate of the region is temperate maritime, summers are cloudy, with weak rain, drizzle and fog often observed. Breezes are a characteristic feature of the local atmospheric circulation.
High-frequency meteorological measurements were performed in 2015 at two locations-on a raft at the lake and on a shore of the lake (figure 1), for a period from 15-28 July. The raft hosted eddy covariance system for heat, momentum and methane fluxes, and a termistor chain in water. On-shore station included sensors for conventional meteorological variables, shortwave and longwave radiation fluxes (see instrumentation summary in table S1 available at stacks.iop.org/ERL/13/035006/mmedia). Besides these observations we also conducted CTD soundings near the raft and whole-lake CTD surveys.
On-raft sonic anemometer data were synchronized with those of gas analyzer to compute eddy covariance fluxes of sensible heat (and not of latent heat as we had no high-frequency humidity sensor), momentum and methane. Eddy covariance fluxes were obtained involving spectral corrections (Moncrieff et al 2004), compensation for density fluctuations (Webb et al 1980), sonic temperature correction for humidity (Van Dijk et al 2004), time lag compensation, double rotation for tilt correction, block averaging, and statistical tests (Vickers and Mahrt 1997). The running averaging interval was chosen 20 min with 10 m in time step of resulting time series of turbulent fluxes. The application of the Kormann-Meixner footprint model (as described in Wilson (2015)) indicated that 80%-cumulative footprint area was confined within the lake surface during 88.5% of observation timespan, i.e. most time the contribution of lake surface to measured turbulent fluxes was dominant.
For a rough estimate of a water transparency, we dipped a white object into the water until it was visible, measuring the respective depth-similar to Secchi disk methodology.

The lake model and simulation setup
The LAKE model is the horizontally averaged, 1D model for heat, momentum, salinity and greenhouse gas dynamics in an enclosed water body (Stepanenko et al 2011, Stepanenko et al 2016, with 40 numerical layers used in this study. It uses k-turbulence closure with stability-dependent coefficients after (Canuto et al 2001) for vertical turbulent viscosity/diffusivity. Water density is given as a function of temperature and salinity according to McCutcheon et al (1993). Sensible, latent heat fluxes, and momentum flux from a lake to the atmosphere are calculated via Monin-Obukhov similarity theory using stability corrections from Paulson (1970), Businger et al (1971). The momentum flux from the atmosphere is distributed between wave development and current acceleration according to Stepanenko et al (2014).
Upwelling shortwave and longwave radiation fluxes at the water surface are computed using albedo (6%) and Stefan-Boltzmann law. Solar (shortwave) radiation S penetrated to a depth z in a lake is given by Beer-Lambert law with extinction coefficient , m −1 , allowed to be z-dependent. In the pycnocline, we chose to increase linearly towards twice of its value in the mixed-layer, i.e. 2 ML 8 , at the bottom, because the lake water was observed to be more turbid there.
For the first estimate of ML , we applied Poole and Atkins formula (Poole and Atkins 1929) to our water transparency measurements, obtaining ML = 0.57 m −1 . Then, the best match between modeled temperature time series and the measured ones was attained with a slightly different value, ML = 0.5 m −1 . Extinction coefficient was the only parameter calibrated in the model.
Bottom sediments are implemented as a set of vertical columns, originating at different lake depths (Stepanenko et al 2016), enabling lateral soil-water heat exchange to be taken into account. The parameterization of surface seiches (Stepanenko et al 2016) is switched off, as along with the wind forcing, there is dynamical forcing by tides, that should interact with seiches in intricate way, not included in the model. The model involves advection terms representing the tributaries and effluents, given their discharges and respective properties (temperature, salinity, velocity). The LAKE model has been validated versus other models and field experimental data in LakeMIP project and other intercomparisons (Stepanenko et al 2013, Stepanenko et al 2014, Thiery et al 2014, Heiskanen et al 2015. For Bol'shoi Vilyui Lake simulation we chose the following settings for the model: the lake depth 7 m, the wind fetch 1.5 km, the total discharge of rivers and brooks into the lake's mixed layer 2.8 m 3 s −1 (Gorin 2013) (the same discharge is assumed from the lake's mixed layer to the ocean), and the temperature of rivers 6 • C. We neglected advection of temperature and salinity by tides, because (i) in the north-western part of the lake the tidal variation of water level is very small (∼1 cm) and (ii) the volume of marine water entering Bol'shoi Vilyui during tidal cycle is highly uncertain.
To set roughness parameter for momentum z 0 in the model, we retrieved z 0 from eddy covariance measurements, using logarithmic law corrected for stratification (Repina et al 2012). This parameter demonstrated correlation with the wind speed, however, existent formulations (including Charnock formula and its extensions) could not satisfactorily reproduce this dependence at wind speed < 4 m s −1 . Therefore, in the model, we chose a somewhat practical approach, i.e. calculated z 0 as a 4th order polynomial of wind speed, approximating our empirical z 0 -wind curve.
The initial conditions for water temperature and salinity in the model were taken from CastAway-CTD sounding conducted from the raft. For model boundary conditions, air temperature and wind speed were provided from the on-raft USA Gill measurements, humidity and wind direction-from the onshore Davis Instruments station data, downward shortwave and longwave radiation-from KippandZonen sensors recordings, all data with 10 min time step. Although wind direction data were evidently affected by the onshore bluff topography, it is wind speed being the most important external variable for vertical mixing in the lake model, and it was measured at the raft, where the surface layer flow is a fair representative for the most of the lake area.

Heat balance of the temperature maximum: quantification of physical mechanisms
The temperature maximum (TeM) develops in the LAKE model and in the observation data below the mixed layer (ML) depth z ML (see section 3), where the salinity gradient causes extremely stable stratification (section 3.3.1) and where, in the LAKE model, turbulent conductivity falls to its minimal value k min . This makes grounds to seek temperature maximum in the pycnocline as an analytic solution of 1D unsteady molecular heat conductance equation with radiation source according to Beer-Lambert law (see text S1).
Remarkably, such a simple approach allows to get TeM in the pycnocline qualitatively and quantitatively similar to that observed in Bolshoi Vilyui case. This implies that radiation absorption and heat conductance comprise 'minimal set' of processes creating TeM. However, more realistic and general consideration of the problem would introduce other factors of TeM such as the lateral heat exchange of the water body with sediments, lake morphometry effects, variability of eddy conductivity with depth, variability of mixedlayer temperature and mixed-layer depth, and diurnal cycle of radiation flux. In the LAKE model, all these factors are taken into account. To quantify them, we consider the heat balance of the water volume confined in the temperature maximum interval (TMI)-a depth interval (z ML , Z B ), where the water temperature exceeds the mixed-layer temperature, T ML (figure 2). The heat balance of TMI is given as (see text S2.1): wherẽ′ is a mean value of T ′ =T-T ML > 0 in the TMI (hereafter referred to as Mean temperature in the temperature maximum interval, MTTM, geometrically, the red area at figure 2 divided by TMI thickness, z -z ML ), F is the change of mean temperature in TMI due to the change of water volume bounded in TMI, F -'dissipation' of MTTM by heat conduction, F is a heating in TMI by shortwave radiation, F is a sink of heat from TMI to sediments that are in contact with it, and ML= − dT ML dt . The term F ML reflects the fact that the excess of TeM over the mixed-layer temperature may be increased solely by the mixed-layer cooling.
Equation (1) may be used to calculate contributions to MTTM of five physical mechanisms mentioned above, averaged over the LAKE model integration time, specifically (see text S2.2): wherẽ′ is the time-averaged MTTM, subscripts in , = , , , , ML stand for the sameprocesses as in , = , , , , ML (see equation 1), and is a residual arising from finite-difference discretization of the model equations. The latter was negligible compared to other terms and hence is not discussed in this paper.

Weather conditions and lake stratification
The fair weather was dominating during the expedition with very small rain events, low cloudiness (untypical for local climate) and high shortwave downward radiation ( figure 3(a), lower panel). Breeze circulation caused calm nights and maximum wind speed during daytime ( figure 3(a), upper panel). The prevailing wind directions were S and SE (sea breeze), as well as NE and NW. High insolation and diurnal cycle of wind created favorable conditions for TeM development (see below). Figure 3(b) shows temperature and salinity profiles for 17.07.2015 used as initial conditions for the LAKE model. The abrupt salinity change from 1-2‰ at 1.5 m depth to 20‰ at 3 m causes very stable stratification and TeM emergence, the latter is, however, very small at this date. Salinity vertical distribution did not change significantly during the coming 12 d of measurements, whereas temperature maximum enhanced considerably. Remarkably, very similar temperature and    Figure 4, left panel, displays the temperature distribution in time-depth space in the model during the simulation period. The warming trend at all depths is conspicuous. Temperature profile consists of the mixed layer of ca. 2 m depth, experiencing clear diurnal cycle, bounded by thermocline below. The decrease of temperature down to 5 • C at the bottom is accompanied by increase of salinity up to 29-30‰ (right panel of figure 4). Due to sharp salinity increase and negligible velocity in the pycnocline, the modeled Richardson number is orders of magnitude above critical value (0.25) in this layer (not shown). Contrary to most freshwater bodies, the temperature at almost all instants 1 5 -0 7 -1 8 1 5 -0 7 -1 9 1 5 -0 7 -2 0 1 5 -0 7 -2 1 1 5 -0 7 -2 2 1 5 -0 7 -2 3 1 5 -0 7 -2 4 1 5 -0 7 -2 5 1 5 -0 7 -2 6 1 5 -0 7 -2 7 1 5 -0 7 -2 8   reaches its maximum below the ML, at a depth of ≈ 2 m, most prominentlyduring the last three days of our field measurements (26-28 July). This temperature pattern is well consistent with measurements ( figure 5). In the mixed layer, where sensors were located at depths 0. The modelled temperature in the pycnocline exhibits much smoother temporal variability than that in the mixed layer, with clear diurnal cycle and warming trend. Both trend and cycle are mostly caused by absorption of the shortwave radiation: increasing turbidity in the model up to ML = 1 m −1 eliminates the trend and significantly reduces the magnitude of Figure 6. The components of TMI heat balance: -'dissipation' of TeM due to heat conductance, -heating in TMI by shortwave radiation, F S is a sink of heat from TMI to sediments;̃′ is a mean excess of temperature in TMI over the mixed-layer temperature T ML . The term is not shown as it is zero almost all times except a small number or model timesteps, where it has huge values compared to other terms and thus is not convenient to depict at the same plot. diurnal cycle (at 3.05 m, from ≈ 0.2 • C to ≈ 0.1 • C, not shown). The role of vertical heat conduction in the diurnal variability and trend of temperature in pycnocline is of minor importance, as under very stable stratification, the LAKE model produces negligible eddy conductivity under the ML (figure 4, right panel).

Salinity-related effects
In the upper part of TMI both temperature and salinity increase with depth. As these thermodynamic variables act on density oppositely, a density ratio is often considered to establish whether stratification is stable or unstable: where = ∕ , = − ∕ are defined in LAKE model according to the equation of state from McCutcheon et al (1993). Stratification is stable when R > 1, but in the interval 1 < R < 10 the diffusive regime of double-diffusive convection is anticipated (Kelley et al 2003, Toffolon et al 2015. Double diffusion leads to effective vertical mixing, potent to hinder the TeM development. However, in Bolshoi Vilyui simulations, R ≫ 10 in TMI almost all time of model integration (figure S7), so that the stratification is enough strong and double diffusion is highly unlikely to cause excessive mixing in the upper part of temperature maximum depth interval.
Salinity profile changes in LAKE model very slowly throughout the period considered, due to the absence of significant sinks and sources of salinity and the fact that only molecular diffusion acts in vertical below the mixed layer.

The heat balance of the temperature maximum in the LAKE model
In the following discussion we will mostly refer to MTTM as a measure of TeM 'intensity', because we have equations (1) and (2) at our disposal to analyze the origins of its value. All terms of TMI heat balance are non-negligible in the LAKE model experiment ( figure 6). Heat conduction F is always negative, i.e. it reduces temperature maximum, with higher absolute nocturnal values, due to larger MTTM and consequently larger temperature gradients in the TMI. Soil-heat-exchange term F is negative as well, but exhibits small temporal variability. Radiation term F is positive, as it increases TeM, or zero (at nighttime). Mixed-layer heating (F ML ) during daytime reduces MTTM, while radiation heating (F ) about halfcompensates this effect (see blue, yellow and black curves at figure 6). On the contrary, during night, the decrease of T ML obviously enhances MTTM.

Water turbidity effects
In this section we examine the sensitivity of MTTM in the LAKE model to water turbidity ML , a crucial parameter controlling TeM development.
First, consider NOMORPH set of numerical experiments, where no lake morphometry effects are included. Figure 7 shows MTTM and its 'components' . An important feature of simulations is that in highly turbid water ( ML = 0.7, 0.8, 0.9 m −1 ), ML temperature effect on MTTM magnitude becomes more important than the radiation absorption ( ML > ). This is due to that TeM occurs mostly at nights in these numerical experiments, when it is caused by the ML cooling down. At low turbidity ( ML = 0.1 m −1 ) we observe ML ≈ ; in this case radiation effect on TeM drops down due to less radiation absorption under the ML, while the ML itself gets less heated from radiation as well, with corresponding decreasing the ML temperature and increasing ML .
An interesting effect is represented by a volumeeffect term , which is positive and one of the largest terms. The thickness (∼volume) of TMI, Δz TMI =zz ML , is actually linked to ML depth z ML , since TMI and ML share a common boundary (figure 2). In the model, the lower boundary of TMI, z , changes much less, than z ML , so that dΔz TMI /dt ≈−dz ML /dt. Mixed-layer thickness undergoes diurnal cycle: during daytime, the ML deepens by enhanced breeze wind forcing, and retreats in evening; therefore, TMI thickness, increases in evenings and decreases in mornings. Hence, F ∼−dΔz TMI /dt * MTTM changes sign in a diurnal course, and one could anticipate that averaging of F over the cycle may yield ≈0. However, the changes of Δz TMI in the morning and in the evening take place at different MTTM values, specifically: contraction of TMI during mornings happens at higher MTTM (because MTTM mostly increases during night) than the TMI expansion in evenings (figure S8). Hence averaging of F over the diurnal cycle provides positive value. In respect to MTTM, corresponding warming effect resembles pumping, because the changes of TMI are alternating contractions and expansions, and at each new cycle the MTTM becomes greater than during the previous one. In other words, the above described mechanism happens due to a ∼ /2 phase shift between diurnal cycles of TeM magnitude and ML depth. Now, examine the lake morphometry effect on temperature maximum by comparing the results from NOMORPH model experiments to these from RAD series ( figure 7). First, water-sediment heat flux term appears when morphometry is included, increasing in magnitude with increasing MTTM. Despite in-soil heat sink from the TMI, MTTM is much higher when morphometry is taken into account, especially at lower ML . The reason is that the radiation forcing term is augmented when the area of the top surface of the volume confined in TMI is much larger than the bottom area of TMI volume, so that the radiation influx to TMI increases compared to radiation loss from TMI (see the form of F in text S2.1). In NOMORPH model runs the top and bottom areas of TMI volume are the same. In contrast, at high ML values, radiation term is of minor importance for TeM development, as discussed above, and MTTM values almost coincide in RAD and NOMORPH simulations.

Implications for other salinity-stratified lakes
Given the analysis above, we can formulate conditions favouring mid-depth temperature maximum development in salinity-stratified lakes: 1. small mixed-layer depth: it allows more radiation to penetrate into pycnonline and causes larger windinduced diurnal variation of ML thickness (leading to larger pumping effect); 2. the water in the mixed-layer is not turbid; The effect of conditions (i) and (ii) may be illustrated by comparing characteristics of salinitystratified lakes featuring TeM in summertime, to those lacking this phenomenon (table 1). Although being limited in number of lakes, the table data corroborates our statement that TeM developes under shallow mixed layers and when the water in the mixed layer is not turbid. On the contrary, meromictic lakes with deep chemocline position tend to lack temperature maximum. Shallow and relatively transparent mixed layer, however, does not ensure TeM emergence, if the salinity-driven stratification below is not strong, so that double diffusion occur (R < 10). This takes place, e.g. in Waldsee Lake (Boehrer et al 2009). Daytime wind speed maximum in surface layer is a typical feature observed over land, though its impact should be reduced for lakes bounded by bluff topography (forest) (Markfort et al 2010). Finally, the last condition (v) can be exemplified by our Bolshoi Vilyui Lake case. Before the expedition period, cloudy weather had been prevailing, so that at the beginning of the period TeM was negligible ( figure 3(b)).
Note that in late summer and autumn in midlatudinal salinity-stratified lakes, TeM typically arises from ML cooling (term F ML in (equation 1)) induced by sensible, latent heat and longwave radiation losses to the atmosphere. This takes place, e.g. in the meromictic lakes of the White Sea coast, where the deep temperature may exceed that at the surface by 5 • C-6 • C in August-September (personal communication with E D Krasnova), or in Fayetteville Green Lake (Brunskill and Ludlam 1969). We also suggest that 12 • C jump of temperature in one of residual basins of Aral Sea reported in (Izhitskiy et al 2016) is of the same nature. However, condition (iii) holds in this case as well.

Conclusions
The 1D k-LAKE model demonstrated a good skill of simulating the thermal regime and sensible heat flux to the atmosphere of an estuarine Bol'shoi Vilyui Lake on Kamchatka peninsula (Russia) at diurnal and weekly time scales, for 11 days in July 2015. The remarkable feature of the lake's thermal regime-middepth temperature maximum (TeM)-was observed and captured by the LAKE model. This feature is likely to be typical for at least summertime at this lake, as was observed in July 2004 as well.
The main prerequisite for mid-depth temperature maximum development in the lake is the salinity increase below the freshwater mixed layer, sharp enough to prevent convective mixing and double diffusion. Given that this condition is satisfied, there are five factors controlling the magnitude of the maximum: 1. radiation absorption below the mixed layer, 2. mixed-layer temperature variability, 3. temperature maximum 'pumping', resulting from the phase shift between the diurnal cycles of mixed-layer depth and the temperature maximum magnitude, 4. heat conduction, 5. water-sediment heat exchange.
Given the mechanisms above, we formulate the conditions favouring emergence of temperature maximum in lakes: 1. small mixed-layer depth; the rough condition on depth, depending on water turbidity, is ∼< 2 m; 2. the water in the mixed-layer is not turbid; 3. strong salinity gradient in chemocline enough to prevent double diffusion; 4. diurnal wind speed cycle with daytime maximum; 5. cloudless weather.
Based on literature data, we provide a set of eight salinity-stratified lakes, where four lakes exhibit temperature maximum and four lakes lack this phenomenon. These two groups are sharply differing in mixed-layer depth, empirically supporting the importance of condition (i). The empirical assessment of other conditions would demand a separate study.
We are not aware of any overview providing estimates of the global distribution of water bodies exhibiting TeM. This study provides theoretical grounds to determine whether a lake is likely to have TeM, based on such general lake characteristics as the freshwater layer depth, water turbidity and salinity gradient.

Acknowledgments
The LAKE model development and analysis of simulation results were supported by the Russian Science Foundation, grant 17-17-01210. Analysis of observation data was supported by Russian Foundation of Basic Research, grant 17-05-01224. The authors are indebted to a team of Vilyuchinsk salmon hatchery for warm reception and assistance in multiple issues of everyday life on Bol'shoi Vilyui. We express our gratitude to O B Tepnin from Kamchatka Research Institute of Fisheries and Oceanography for instrumentation he kindly provided in our disposal. The data used in this paper are listed in the references.