The world's largest heliothermal lake newly formed in the Aral Sea basin

Two short-term intense field surveys were carried out in Chernyshev on 13–17 October 2015 and 24–28 June 2016 including in situ measurements and water sampling at stations forming a cross-section across the deepest part of Chernyshev (figure 1). The in situ measurements at each station included vertical profiling of major physical–chemical characteristics and water sampling. Profiles of conductivity, temperature, depth, fluorescence and turbidity were taken in October 2015 with a SBE19plus profiler (SeaBird Scientific, USA), and in June 2016 with a Rinko CTD-Profiler (JFE Advantech, Japan), measuring the same characteristics plus dissolved oxygen (DO) content. Measurements of conductivity and fluorescence had been later discarded since the sensors malfunctioned in the hypersaline waters. Additionally, water transparency was measured using the standard Secchi disk. Water samples were taken by a 5 l Niskin Hydro-Bios bottle and analyzed later in laboratory for salinity by the gravimetric method [17].

Zoom In Zoom Out Reset image size

Continuous measurements of temperature and DO were conducted at two mooring stations (figure 1). Station 1 was installed at 14 m total water depth and was equipped with 12 temperature loggers TD-1060 (RBR Ltd, Canada) at depths of 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 14 m. Six DO loggers D-Opto (Zebra Tech, New Zealand) were added at 2, 4, 8, 10, and 13 m depths. Station 2 was set at 10 m water depth and was equipped with ten TD-1060 loggers at water depths of 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 m. The sampling period was set to 20 s for temperature and 1 h for DO measurements. Stations were deployed on 13 October 2015 and recovered on 26 June 2016. Depth-resolved measurements were however available for the period from 13 October to 21 November 2015 only. Later, when the surface water temperatures approached ∼0 °C, the dissolved salts started to crystallize on the upper part of the moorings and caused their eventual sinking to the bottom (supplementary figure S3). Afterwards, temperature and DO measurements were available only from the bottom, showing a constant temperature of 15.2 °C and zero oxygen content over the winter.

Damping of the vertical heat transport across the water column is the critical process leading to the meromixis formation. The time-depth temperature data collected between 13 October 2015 and 21 November 2015 allow direct estimation of the vertical mixing rates. To calculate the net vertical heat flux Q_TOTAL (z) across the lake water column and the vertical turbulent diffusion coefficient K_Z we applied the inverse solution of the heat transport equation,

$$\begin{equation}\begin{array}{*{20}{c}} {\dfrac{{\partial T\left( {t,z} \right)}}{{\partial t}} = - \dfrac{{\partial {Q_{{\text{total}}}}\left( {t,z} \right)}}{{\partial z}}\,} \end{array}\end{equation} \tag{ 1 }$$

where T (t, z) is the temperature [K] at time t [s] and depth z [m], Q_total (t, z) = Q_turb (t, z) + Q_I (t, z) + Q_c (t, z) is the bulk kinematic vertical heat flux [K m s⁻¹], whose three major components are the turbulent exchange Q_turb (z, t), the short-wave solar radiation Q_I (t, z), and the conductive (molecular) flux Q_c (t, z). Integration of equation (1) from the bottom z = H to some depth z yields the expression for the vertical profile of the heat flux

$$\begin{equation}{Q_{{\text{total}}}}\left( z \right) = - \int\limits_H^z {\frac{{\partial T\left( {t,z} \right)}}{{\partial t}}} {\text{d}}\zeta + {Q_{{\text{total}}}}\left( H \right)\end{equation} \tag{ 2 }$$

where the water-sediment heat flux Q_total (H) can be assumed negligible, since ∂T/∂t|_H ≈ 0 и ∂T/∂z|_H ≈ 0 from the observations. Assuming further that no solar radiation Q_I (t, z) penetrates below the layer of strong turbidity and salinity gradients z = h_SML (figure 2(a) and supplementary video 1), where h_SML is the depth of the surface mixed layer, Q_total (t, z) ≈ Q_turb (t, z) + Q_c (t, z) = [K_z (z) + κ] ∂T (t, z)∂z⁻¹, and one can calculate the vertical diffusion coefficient K_Z from

$$\begin{equation}{K_z}\left( z \right) + \kappa = \frac{{\int_H^z {\frac{{\partial T\left( {t,z} \right)}}{{\partial t}}} {\text{d}}\zeta }}{{{{\left( {\frac{{\partial T\left( {t,z} \right)}}{{\partial z}}} \right)}_z}}}{\text{ at }}{h_{{\text{SML}}}} < z < H\end{equation} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

where κ = Q_c (∂T/∂z)⁻¹ = O (10⁻⁷) m²s⁻¹ is the molecular heat conduction coefficient. The net flux at the lake surface Q_s was estimated from extrapolation of the flux profile in the surface mixed layer (SML). For a mixed layer with homogeneous vertical temperature distribution (0 ⩽ z ⩽ h_SML), a linear vertical flux distribution follows directly from equation (1) and is supported by our calculations (figure 3(b)):

$$\begin{align} {Q_{{\text{total}}}}\left( z \right)& = { }{Q_{\text{s}}}\left[ {1-\frac{z}{{{h_{{\text{SML}}}}}}} \right] \nonumber \\ & \quad + {Q_{{\text{total}}}}\left( {{h_{{\text{SML}}}}} \right)\frac{z}{{{h_{{\text{SML}}}}}}{\text{at }}0 < z < {h_{{\text{SML}}}}{ } .\end{align} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Herewith, the surface flux Q_s can be calculated by linear extrapolation of equation (2) to the water surface as Q_s = Q_total (z₁) − (Q_total (z₂) − Q_total (z₁)) z₁/(z₂ − z₁), where z₁ and z₂ > z₁ are any two depths within the SML, at which Q_total is calculated from equation (2). We used the two uppermost calculation depths z₁ = 2.5 and z₂ = 3.5 m for the data from the mooring Station 1.

Short-term diurnal and sub-diurnal temperature variations have correspondingly smaller spatial scales, for which the 1 m vertical resolution of measurements can appear too coarse, resulting in unrealistic flux fluctuations. To filter out diurnal variations from the seasonal heat budget, the original time series of temperature were averaged over 48 h prior the calculations. As an independent cross-validation, the estimates of the surface heat flux were compared with the sum of sensible and latent heat fluxes from ERA-Interim reanalysis [18] and with the surface fluxes computed by the bulk aerodynamic approach [19] using the measured surface water temperatures as input. Thereby, the backward lake effect on the air–land interaction was indirectly estimated by comparison of the ERA-based results (no lake effect) vs measurement-based fluxes.

To estimate the seasonal variability in the heat budget, we applied a one-dimensional (1D) two-layered heat budget and temperature model. The model is aimed at (a) revealing the major physical mechanisms producing the heliothermal behavior of the lake, (b) reconstruction the seasonal variability in the heat budget, (c) estimation of the subsurface temperature maximum lifetime and the range of the temperature variations within it, and (d) estimation of the heat accumulation rate in the deep monimolimnion. The model is based on the heat conservation only, without any artificial fitting to the observed temperatures and is validated by comparison of the modeling results against the observed temperatures. Thereby, the model provided the opportunity to quantify the seasonal balance of the heat sources and sinks, which leads to the formation of the heliothermal regime. Based on the empirical evidence, the model assumptions were formulated as follows:

• The vertical position of the chemocline remains constant.
• The vertical turbulent flux across the chemocline is negligible.
• The temperature profile in the SML (mixolimnion) is homogeneous.
• No solar radiation penetrates beneath the chemocline.

Using K_Z from equation (3) as a seasonal-mean estimation, the accumulation of heat in the deep monimolimnion (h_SML < z < H) is estimated as

$$\begin{equation}\begin{array}{*{20}{c}} {\frac{{\partial T\left( {t,\,z} \right)}}{{\partial t}} = {K_Z}\frac{{{\partial ^2}T\left( {t,\,z} \right)}}{{\partial {z^2}}}\,} \end{array}\end{equation} \tag{ 5 }$$

with the boundary conditions at the bottom, z = H,

$$\begin{equation}\begin{array}{*{20}{c}} {{K_Z}\frac{{\partial T\left( {t,{ }z} \right)}}{{\partial z}} = 0{ }} \end{array}.\end{equation} \tag{ 6 }$$

At the top of the monimolimnion, the temperature jump produced by absorption of solar radiation is balanced by the nearly conductive heat flux,

$$\begin{equation}\begin{array}{*{20}{c}} {{K_Z}\left( {{h_{{\text{SML}}}}} \right)\frac{{\partial T\left( {t,{ }z} \right)}}{{\partial z}} = - {Q_{\text{I}}}\exp \left( { - \gamma {h_{{\text{SML}}}}} \right){ }} \end{array}\end{equation} \tag{ 7 }$$

where γ is the extinction coefficient in the one-band exponential radiation decay, K_Z (h_SML) ≈ κ= 10⁻⁷ m²s⁻¹. In the SML, temperature is modeled by integration of equation (1) from the surface to z = h_SML

$$\begin{equation}\begin{array}{*{20}{c}} {{h_{{\text{SML}}}}\frac{{\partial {T_{{\text{SML}}}}}}{{\partial t}} = {Q_{\text{s}}} - {Q_{\text{b}}} + {Q_{\text{I}}}\left( 0 \right)\left[ {1 - \exp \left( { - \gamma {h_{{\text{SML}}}}} \right)} \right]{ }} \\[6pt] \end{array}\end{equation} \tag{ 8 }$$

where Q_b is the heat flux at the SML base.

The surface heat flux Q_s is calculated from the air-lake heat balance [19] using the surface water temperatures at the previous time step and atmospheric characteristics from the ERA-Interim reanalysis. The chemocline is simulated in the model by a zero-thickness layer at z = h_SML, with the heat flux across it Q_b = χ(T (0) − T_SML). The heat transfer rate across the chemocline is parameterized as χ = 2.7 × 10⁻⁷ m s⁻¹ that corresponds to molecular heat flux across a 0.5 m thick layer. As initial condition, T (z) = const = T (H) = 15 °C (corresponding to the measured annually constant temperature of the deep monimolimnion) and T_SML= 10 °C (roughly corresponding to the 'equilibrium' water surface temperature approximated as the annual mean air temperature) were used, corresponding to some date t₀ at the initial stage of the summer heating period, when the mixolimnion stayed colder than monimolimnion and stability of the water column was controlled by salinity gradient. After several test model runs, t₀ was set to 1 May 2015, whereas the model was weakly sensitive to the variation of the date within ∼30 days interval. The radiation extinction coefficient γ was set to 0.4 m⁻¹, based on Secchi disc measurements and model sensitivity experiments. Equations (5)–(8) were solved numerically using a Python script.

The area and the mean depth of Chernyshev at the moment of measurements varied around ∼80 km² and 4 m, respectively, with the maximum depth of about 15 m (figure 1). Chernyshev remains connected to the Large Aral through a narrow straight in its south, and, intermittently, to ephemeral lakes on the dry bottom of the former Aral Sea in the east (figure 1). The current water budget is unknown, consisting presumably from episodic releases of water from the Small Aral, meltwater, precipitation-evaporation, and, possibly, groundwater discharge.

Salinity of the lake exhibited a three-layered structure (figure 2(b)): upper mixed layer, a layer with a salinity jump of ∼100 g kg⁻¹ between 4 and 6 m water depth, and a deep nearly homogeneous layer extending to the bottom (in terms of meromixis: mixolimnion, chemocline, and monimolimnion, respectively). In the monimolimnion, salinity slightly increased downwards from 233 to 244 g kg⁻¹ between 6 and 12 m depth. Light conditions varied strongly across the water column: the upper transparent layer had sharp lower boundary with turbidity peak at the chemocline (figure 2(a)). Underneath, turbidity decreased, but remained twice as high as those in the mixolimnion. A vivid qualitative picture of the transparency changes in the chemocline is provided by the underwater camera moving from water surface down to the bottom (supplementary video 1): the water column beneath the chemocline remains completely impermeable to the solar radiation beneath the sharp turbidity peak. The white (Secchi) disc measured transparency depth of 6.5 m, which value corresponded to the depth of the turbidity peak and likely underestimated the transparency of the upper layer.

Vertical profiling in the middle of the autumn cooling period (October 2015) and in the middle of the summer heating period (June 2016) provided an estimation of the range of seasonal temperature and oxygen variations (figure 2(a)): the most striking feature of the thermal regime is the formation of the mid-depth temperature inversion in October with the downward temperature increase from 12.5 ℃ at the water surface to the maximum of 30.7 ℃ at ∼5 m depth. This subsurface temperature maximum was followed by a gradual temperature decrease down to 8–9 m depth and a homogenous deep layer with temperature values about 15.2 ℃ extending to the bottom. According to the ERA Interim reanalysis data [18], the mean air temperatures for the period 10–15 October 2015 amounted at 0.0 ℃ with the minimum of −5.7 ℃. The inversion degraded in summer (figure 2(a)), when the upper mixolimnion was heated up to 29.4 ℃.

The mixolimnion was well-oxygenized, with DO concentrations close to 70%–80% of saturation (9.0 mg l⁻¹ in October and 5.5 mg l⁻¹ in June, note that the real absolute concentrations may differ due to effects of the high salinity on DO solubility). Under the chemocline, the environment was completely anoxic throughout the year (figure 2(b)).

The spatial structure of temperature and DO followed the three-layered meromictic pattern (figures 2(a) and (b)): the vertically homogeneous mixolimnion was well-saturated with oxygen and revealed strong seasonal temperature variations. In turn, beneath the ∼6.5 m depth, the oxygen was absent, and temperatures remained at the nearly constant value of 15.2 ℃ beneath ∼10 m depth. Radial transects demonstrated the persistence of the three-layered structure across the lake (figures 2(c) and (d)).

Reliable continuous temperature and DO records from autonomous moored stations were available from 13 October to 21 November 2015 only. Later, the submerged temperature-oxygen chains started to gradually sink due to crystallization of salt on buoys and sensors in cold surface waters and eventually drowned to the bottom around 31 December (see section 2 and supplementary information 1). Nevertheless, the mooring data provide a valuable insight into the heat storage and vertical heat transport: during the intense cooling of the mixolimnion, temperatures in the upper mixed layer decreased from 12.5 ℃ on 13 October to 2.6 ℃ on 21 November, while the intermediate temperature maximum decreased from 30.6 °C to 21.8 °C (figures S2(a) and (b)). Herewith, the temperatures in the chemocline in November remained ∼6.2 ℃ higher than at the bottom, ∼19.2 ℃ higher than at the surface, and >27 ℃ higher than the minimum air temperatures. The vertical position of the chemocline was constant, so that the autumn cooling affected only the upper water column, without entrainment of cold surface water into the monimolimnion. Correspondingly, anoxic conditions in the bottom layer along with the temperature values of 15.2 ℃ remained stable beneath the chemocline. Salt precipitation started after the mixolimnion temperatures dropped to ∼3 °C (supplementary figures S1–S3) causing gradual sinking of the moorings. The lowest temperatures in the mixolimnion of ∼0 °C were measured around 20–25 December; afterwards, all sensors settled on the bottom and continued to record during January–June 2016 constant temperatures of 15.2 °C and zero DO content (supplementary figure S1).

The surface heat fluxes (figure 3(a)) demonstrated a strong (100–200 W m⁻²) heat release to the atmosphere from the upper mixolimnion, contributed mostly by evaporation. The ERA-interim reanalysis showcased a higher amplitude of the flux variability than data-based calculations, indicating the damping effect of the lake on the air-land heat exchange not accounted for in the reanalysis calculations. The vertical profiles of the heat flux Q_total (z) and the vertical exchange coefficient K_z (z) averaged over the period of moored observations (figure 3(b)) showed strong mixing in the mixolimnion with linear decrease of Q_total (z) from the maximum at the water surface to nearly zero in the chemocline, characteristic of the well-mixed surface layer with homogeneous vertical temperature distribution (equation (1)). K_z in the chemocline dropped to the level of the purely conductive exchange, close to the molecular heat conductivity κ ≈ 10⁷ m²s⁻¹, showcasing effective cancellation of the vertical transport by the density (salinity) gradient in the chemocline. Herewith, the strong surface cooling (figure 3(a)) and cancelling of the heat flux in the chemocline at ∼6 m depth (figure 3(b)) can be identified as the ultimate mechanisms resulting in the intermediate temperature maximum formation (figure 3(c)). In the deeper monimolimnion, K_z varied around the value 10κ, characteristic of non-zero vertical turbulent mixing rates, though orders of magnitude smaller than in the surface layer. The mixing in the wind-isolated monimolimnion is apparently enhanced by internal waves, which have been shown to produce the same tenfold increase of the background mixing over the molecular level in stratified lakes isolated from the atmosphere by the ice cover [20, 21].

The model simulated maximum subsurface temperatures of up to 37.4 °C in mid-August (figure 3(c)) caused by accumulation of solar radiation in the chemocline. Maximum inverse temperature gradients at the base of the upper mixed layer reach 18.2 K m⁻¹ in late October, in a good agreement with temperature measurements (figure 3(d)). The summer temperatures also agreed qualitatively with the observations, with larger deviations in the upper mixed layer, where the model assumes vertical homogeneity, but the real temperature profiles are affected by short-term diurnal fluctuations. The heat storage in the deep monimolimnion (8–12 m depth from the surface) stays nearly constant; the bottom temperature increases at ∼0.8 °C per annum. The latter value implies continuous accumulation of heat by Chernyshev on inter-annual time scales, but should be treated with care, because 1D models are prone to underestimate the near-bottom mixing [22, 23].