Consequences of the Aral Sea restoration for its present physical state: temperature, mixing, and oxygen regime

The restoration of the North Aral Sea was an unprecedented effort to save a large water basin by construction of a dam that separates it from the rest of the desiccating Aral Sea area. As a result, the lake volume has stabilized at 27.5 km³, the area has increased from 2800 km² (2006) to 3400 km² (2020), and the salinity has dropped from 18 to 10 g kg⁻¹. The consequences of this unique experiment include highly dynamic changes in the thermal conditions, seasonal stratification, ice regime, and dissolved oxygen content and remain not fully quantified to date. We analyze the current state of the North Aral Sea with regard to stabilization of its long-term dynamics. We further consider the possible future projections in view of the global change effects on the regional hydrological regime and potential water management measures. Using data from a series of expeditions to the North Aral Sea in 2016–2019 and year-long continuous monitoring of the thermal and oxygen regime by an autonomous mooring station, we present the first comprehensive analysis of the North Aral system behavior on seasonal to interannual scales after its restoration. We demonstrate that the present seasonal mixing regime is distinguished by relatively weak summer thermal stratification occupying about 7 % of the lake volume. Salinity does not contribute to the summer density stratification, but a stable salinity stratification can develop during ice melt in late winter. On the background of weak thermal stratification, highly energetic internal waves with periods of ∼4.5 d dominate the near-bottom dynamics and facilitate mixing at the lake bottom. As a result, the bulk of the water column remains well-saturated with oxygen throughout the year. However, low oxygen conditions may develop in the deepest part of the lake in midsummer. In summary, the mixing regime of the restarted lake favors vertical transport of dissolved matter and water–sediment mass exchange, ensuring oxygenation of deep waters and supply of nutrients to the upper water column. While the North Aral Sea is restored to a well-mixed state similar to that before its desiccation started, its seasonal mixing regime is currently in unstable equilibrium, wobbling between polymictic and dimictic conditions. The fragility of this seasonal pattern is demonstrated by modeling results: slight changes in the water level or transparency may turn the Aral Sea to a steadily dimictic or polymictic state.

1 Introduction

Drying of lakes as a result of water resources mismanagement and climate change is a growing environmental concern around the world. Large endorheic lakes are particularly vulnerable to desiccation. Many of these lakes (Wang et al., 2018), such as the Aral Sea, Lake Chad, Lake Urmia, Great Salt Lake, and the Dead Sea, are experiencing acute shrinkage of their water volume, causing large-scale environmental, ecological, and socioeconomic consequences. As the problem of drying lakes grows globally, restoration efforts are becoming increasingly important. Various restoration measures have been proposed, such as regulating water withdrawal, water saving techniques in agriculture and industry, and mitigation of climate change through global cooperation and policy-making (Yapiyev et al., 2017). A particularly effective approach consists of lake water replenishment by managing the major inflows. Some notable examples include a proposal involving transfer of water from the Congo River Basin to Lake Chad (Adeniran and Daniell, 2021; Kitoto, 2021), which has shrunk by as much as 90 % since the 1960s due to overuse of the water for irrigation and the effects of climate change; the Dead Sea preservation project (Asmar and Ergenzinger, 2002; Asmar, 2003), which considers building a canal from the Red Sea to the Dead Sea to stabilize water levels; and the Lake Urmia restoration program (Danesh-Yazdi and Ataie-Ashtiani, 2019; Parsinejad et al., 2022), which includes diverting water from nearby rivers and a plan to release water from upstream dams when levels are critically low. Such “megaprojects” have multiple, barely predictable consequences for lake ecosystems and the regional water budget. In this regard, the outcomes of the Aral Sea restoration project are particularly insightful for current and future projects, serving as an example of a successful large-scale experiment with multifaceted consequences for the lake and its environment.

Once the fourth largest lake in the world, the Aral Sea has been shrinking since the 1960s after the rivers that fed it were diverted by Soviet irrigation projects. Rapid desiccation of the Aral Sea draws continuous attention from researchers as an example of fast anthropogenically driven change in a large aquatic ecosystem at unprecedented spatial scales. The formerly brackish lake was reduced to several (semi-)isolated water bodies with significantly divergent hydrological and biogeochemical conditions (Izhitskiy et al., 2016). By the 2000s, due to the active use of the two main tributaries, Amu Darya and Syr Darya, the lake had shrunk within 40 years to 10 % of its original surface area (Zavialov, 2007), resulting in severe ecological consequences, loss of fisheries, and local climate extremes. At the same time, the water salinity in isolated parts of the Aral Sea increased from 11 g kg⁻¹ to more than 200 g kg⁻¹ (Andrulionis et al., 2021, 2022; Izhitskiy et al., 2021).

As a countermeasure to prevent further desiccation, the 12 km long Kokaral Dam was constructed in 2005, separating the northern part of the Aral Sea from the rest of the former lake basin (Fig. 1), and was widely recognized as an exceptional success in large-scale water management and restoration (Micklin, 2014): the dam confined the discharge of the Syr Darya River to the North Aral Sea, increased its water level at 42.5 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (compared to about 40.5 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ before dam construction), and stabilized the volume at 27.5 km³. The lake area increased from 2800 km² in 2006 to 3400 km² in 2020, and the salinity dropped from 18 g kg⁻¹ in August 2002 (Friedrich and Oberhänsli, 2004; Friedrich, 2009) to 11 g kg⁻¹ in October 2014 (Izhitskiy et al., 2016), the latter value being comparable to that before the desiccation (Bortnik and Chistyaeva, 1990). Several studies have described significant changes in the biological characteristics of the lake during the period after its isolation (Plotnikov et al., 2016; Massakbayeva et al., 2020). Water level and salinity fluctuations led to diversity changes in plankton and zoobenthos communities of the North Aral Sea before (Aladin et al., 2005) and after the Kokaral Dam construction (Krupa et al., 2019; Klimaszyk et al., 2022) along with changes in ichthyofauna and fisheries (Ermakhanov et al., 2012). While rapid changes in fish and zooplankton communities were documented during the 15 years of restoration, little was known about the seasonal mixing regime of the restored Aral Sea, including the thermal and salinity structure of the water column.

Figure 1 Residual basins of the Aral Sea by MODIS Terra from 5 May 2018: 1 – the North (small) Aral Sea, 2 – the western South (large) Aral Sea, 3 – the eastern South (large) Aral Sea; approximate bathymetry map of the North Aral Sea with the location of measurement sites M2 (2016) and M1 (2018 and 2019).

The seasonal mixing regime is a fundamental lake characteristic, which governs the distribution of nutrients, oxygen, and organisms throughout the water column and determines the ecological balance and overall health of a lake (Hutchinson, 1957). Thereby, seasonal mixing regulates the primary production and sustains biodiversity of aquatic life. A crucial role of seasonal mixing consists of oxygenating the lake water column (Golosov et al., 2007; Valerio et al., 2019; Pilla et al., 2023). Oxygen from the surface gets mixed down to deeper waters, providing essential life support for aerobic organisms, including fish and bacteria. It also helps to release gases generated at the bottom of the lake, such as hydrogen sulfide and methane, and prevents toxic gas accumulation. In continental and cold arid climates, where the climatic mean air temperature crosses the freezing point of freshwater twice a year, two major seasonal mixing regimes are common: the dimictic regime, with two events of full mixing (spring and fall “overturns”) separating the winter and summer stagnation periods, and the polymictic regime, in which seasonal stagnation periods are prevented by wind and convective mixing regularly reaching the lake bottom (Kirillin and Shatwell, 2016).

The far-reaching implications of seasonal mixing disruption for the health and biodiversity of lake ecosystems have been demonstrated by the fate of the southern part of the Aral Sea, which was not affected by the restoration measures (Izhitskaya et al., 2019; Izhitskiy et al., 2021): after seasonal mixing was canceled by strong vertical salinity gradients, the residual Aral waters partially turned into anoxic environments with extremely low biodiversity and high rates of methane production in deep waters.

Before desiccation, the Aral Sea was a typical dimictic lake (Zavialov, 2007). The North Aral Sea was ice-covered from late December or early January for 120–140 d. Summer was characterized by strong thermal stratification with surface temperatures >20 °C and bottom temperatures slightly above the maximum density value of ∼1.6 °C. Salinity varied seasonally between 9.0–10.8 g kg⁻¹ and was nearly evenly distributed vertically by spring and fall overturns.

Almost uniform vertical temperature profiles were observed in late summer 2002 before the Kokaral Dam construction (Friedrich and Oberhänsli, 2004) and in autumn 2014 after the restoration (Izhitskiy et al., 2016), suggesting that the lake mixes down to the bottom at least once a year. However, before restoration, hypoxic conditions were reported near the lake bottom (Friedrich and Oberhänsli, 2004; Friedrich, 2009), indicating active oxygen consumption in the deep thermally stratified part of the water column during the summer stagnation period. The seasonal mixing and stratification regime of the North Aral Sea after restoration remained unknown, but it is crucial for understanding the consequences of the restoration measures for the lake ecosystem at all trophic levels. Located in the arid continental climate, the North Aral Sea undergoes strong seasonal variations in the heat exchange with the atmosphere and, as a result, in the water temperatures and density stratification. The lake surface remains isolated from the atmosphere by the ice cover during winter months (Kouraev et al., 2004), which leads to winter hypoxia development in lakes (Golosov et al., 2007). In summer, deep hypoxia can be accelerated by the strong surface heating and development of thermal stratification, preventing oxygen exchange between surface and bottom waters. The North Aral ecosystem remains one of the largest ecosystems in the arid region of Central Asia. Taking into account its vulnerability to water use and climatic changes in hydrological balance in the second largest endorheic lake basin (after the Caspian Sea), understanding the response of the lake to mitigation measures is crucial for planning future activities in the context of climate-driven water scarcity and growing irrigation.

Summarizing the information above, the following basic hypotheses can be formulated on the present and future mixing and oxygen dynamics of the North Aral Sea: (i) the lake may have re-established a dimictic pattern, with two mixing events in spring and fall, separated by periods of stagnation in winter and summer; (ii) hypoxic conditions could develop in the lake's deep areas during stagnation phases due to limited deep oxygen replenishment; and (iii) gradual shifts in water level may cause the lake to oscillate between polymictic and dimictic mixing regimes, depending on changes in lake depth. Below, we use data from several measurement campaigns, including year-round monitoring of temperature and near-bottom oxygen concentrations, as well as modeling, climate scenarios, and remote sensing data, aiming at

• quantifying the consequences of large-scale lake restoration measures for the annual thermal and mixing regime of the North Aral Sea, including potential stagnation;

• revealing the effects of the mixing regime on dissolved oxygen levels in deep waters;

• assessing possible future impacts of regional hydrological regime changes on lake mixing.

2 Materials and methods

2.1 Sampling sites and procedures

We performed our observational studies in the western part of the North Aral Sea, the Shevchenko Bay (Fig. 1), during three short-term intense field surveys that took place between 2016–2019. On 24 June 2016, direct measurements were carried out in the southern part of Shevchenko Bay at coordinates 46.41° N, 60.27° E and water depth of 13.1 m (site M2). On 24 September 2018 and 28 September 2019 direct measurements were carried out in the western part of Shevchenko Bay at coordinates 46.54° N, 60.06° E and water depth of 11.9 m (site M1). The bathymetry of the North Aral Sea can be estimated only approximately and varies within the seasonal hydrological cycle. The last bathymetric map of the Aral Sea was created in the 1960s (Bortnik and Chistyaeva, 1990). Extrapolating these data to the present satellite altimetry and using selected depth measurements from our surveys, the mean depth of the North Aral Sea varies in the range of 6–8 m with maximum depths around 15–18 m. The mean depth of the Shevchenko Bay is 7–8 m, with maximum depths reaching up to 15–16 m in its central part.

During each of the surveys, in situ measurements included vertical profiling of the main physical and chemical characteristics and water sampling. Profiles of conductivity, temperature, depth, fluorescence, and dissolved oxygen were taken in June 2016 and September 2018 with a Rinko conductivity–temperature–depth (CTD) profiler (JFE Advantech, Japan). In September 2019, profiles of conductivity, temperature, and depth were taken with an RBRConcerto CTD profiler (RBR Ltd., Canada). Measurements of water transparency were conducted using a standard limnological Secchi disk. The water samples were taken with a 5 l Niskin Hydro-Bios bottle and analyzed later in the laboratory for salinity and ionic compositions using the methods described by Andrulionis et al. (2022).

2.2 Long-term autonomous monitoring

An autonomous moored chain was installed in Shevchenko Bay at the coordinates 46.54° N, 60.06° E and a water depth of 11.8 m, just near the sampling site M1. The chain was deployed on 24 September 2018 and recovered on 28 September 2019. The moored chain was equipped with nine temperature sensors, TR-1060 (RBR Canada, accuracy 0.002 K), distributed between the surface and the bottom of the water column at depths of 2.22, 2.77, 3.42, 4.82, 5.82, 6.92, 8.92, 10.02, and 11.77 m. In addition, a dissolved oxygen (DO) logger, D-Opto (ZebraTech, New Zealand, accuracy 1 % of measurement or 0.02 mg L⁻¹, whichever is greater), was deployed at 11.77 m water depth. In situ measurements were performed with sampling rates of 30 s for water temperature and 1 h for dissolved oxygen concentration.

2.3 Framework of stratification analysis

We analyzed thermal stratification in terms of the buoyancy frequency,

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }N\left(z\right)=\sqrt{{\displaystyle \frac{g}{\mathit{\rho }}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial z}}},\end{array}$$

where g=9.81 m s⁻² is the gravity acceleration, and ρ is the water density, calculated as a function of temperature and salinity. ρ₀ is the reference water density in the Boussinesq sense, which can be taken as the mean density constant in time and depth across the water column.

For estimation of vertical density gradients and N(z) we used the TEOS-10 toolbox (McDougall and Barker, 2011). The salt composition of the North Aral water (Andrulionis et al., 2022) (see also Table 1 in the Results section) is slightly different from that of the seawater (Millero et al., 2008); however, the water column is well-mixed with respect to salinity, making the contribution of salinity to vertical density stratification nearly negligible.

Table 1 Salinity and major ions in the North Aral Sea waters in g kg⁻¹. Data for August 2002 as pre-Kokaral conditions are given from Friedrich and Oberhänsli (2004) for comparison. Values for September 2019 are taken from Andrulionis et al. (2022).

Density stratification estimates were further applied to evaluate characteristics of internal gravity waves. For a water column of depth H, the buoyancy frequency determines the celerity (phase speed) of a long internal wave,

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{C}_m={\displaystyle \frac{\mathrm{1}}{m\mathit{\pi }}}\underset{\mathrm{0}}{\overset{H}{\int }}N\left(z\right)\mathrm{d}z,\end{array}$$

where $m=\mathrm{1},\mathrm{2},\mathrm{3}\mathrm{\dots }$ is the vertical “natural wave mode”, and m=1 corresponds to the gravest (slowest) wave C₁ (Gill, 1982). For the two-layered vertical stratification with a sharp density gradient separating a lighter homogeneous layer (epilimnion) from the denser hypolimnion beneath, typical for seasonally stratified enclosed lakes, the expression for the first-order internal wave speed C₁ reduces to

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{C}_{\mathrm{1}}=\sqrt{\left(g{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}}}\right)\left({\displaystyle \frac{h_{\mathrm{1}}{h}_{\mathrm{2}}}{h_{\mathrm{1}}+h_{\mathrm{2}}}}\right)}=\sqrt{g^{\prime }{h}_{\text{eq}}}.\end{array}$$

Here, g^′ is the “reduced gravity” expressing the density jump across the infinitesimal boundary between the two layers (ρ₂>ρ₁), and h_eq is the “equivalent depth”. The ratio of C₁ to the Coriolis acceleration f represents the first-order baroclinic Rossby radius, a fundamental length scale expressing the effect of the Earth's rotation on the internal wave characteristics. Here, $f=\mathrm{2}\mathit{\omega }\sin \mathit{\varphi }\approx {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ is the Coriolis frequency, $\mathit{\omega }=\mathrm{7.29}\times {\mathrm{10}}^{-\mathrm{5}}$ rad s⁻¹ is the angular speed of the Earth's rotation, and ϕ is the geographic latitude (∼46.5° N for the North Aral Sea). The dimensionless ratio of the baroclinic Rossby radius to the lateral scale of the lake L is the Burger number (Gill, 1982):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }S={\displaystyle \frac{C_{\mathrm{1}}}{fL}}.\end{array}$$

As a characteristic horizontal scale L of the basin-scale motions, the average radius of a circle approximating the Shevchenko Bay (Fig. 1) was chosen. The reason for this choice is isolation of the bay from the eastern part of the lake constraining the large-scale water motions within its nearly round-shaped basin. The model of basin-scale internal waves in a large circular lake (Csanady, 1967; Lamb, 1932) was applied in the further analysis of these motions. According to the model solution of the water motion equations, the internal waves are characterized by a set of frequencies σ, determined by two sets of eigenvalue problems.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(\sqrt{{\mathit{\varsigma }}^{\mathrm{2}}-\mathrm{1}}/S\right)J_{n-\mathrm{1}}\left(\sqrt{{\mathit{\varsigma }}^{\mathrm{2}}-\mathrm{1}}/S\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }+n(\mathrm{1}/\mathit{\varsigma }-\mathrm{1})J_n\left(\sqrt{{\mathit{\varsigma }}^{\mathrm{2}}-\mathrm{1}}/S\right)=\mathrm{0}\text{ at }\mathit{\varsigma }>\mathrm{1}{\displaystyle }& {\displaystyle }\left(\sqrt{\mathrm{1}-{\mathit{\varsigma }}^{\mathrm{2}}}/S\right)I_{n-\mathrm{1}}\left(\sqrt{\mathrm{1}-{\mathit{\varsigma }}^{\mathrm{2}}}/S\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }+n(\mathrm{1}/\mathit{\varsigma }-\mathrm{1})I_n\left(\sqrt{\mathrm{1}-{\mathit{\varsigma }}^{\mathrm{2}}}/S\right)=\mathrm{0}\text{ at }-\mathrm{1}<\mathit{\varsigma }<\mathrm{0}\end{array}$$

Here, $n=\mathrm{1},\mathrm{2},\mathrm{\dots }$ is the radial wave mode, S is the Burger number (Eq. 4), J_n and I_n are the Bessel functions and the modified Bessel functions of the first kind of order n, respectively, and $\mathit{\varsigma }=\mathit{\sigma }/f$ is the dimensionless wave frequency. The solution of Eqs. (5)–(6) gives a pair of frequencies, one corresponding to a cyclonic (ς<0) Kelvin wave with amplitudes growing towards the lake shores due to the exponential form of I_n and the other corresponding to an anticyclonic (ς>1) Poincaré wave with amplitudes varying periodically across the lake due to the periodic form of J_n. For a real lake, the wave frequencies can be estimated from data on vertical density stratification and the lake size and may be identified in the spectra of density (temperature) oscillations.

To estimate the stratification stability on lake-wide scales we utilized the Schmidt stability parameter,

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\text{St}={\mathit{\rho }}_{\mathrm{0}}g\overline{h}(Z_{\mathrm{v}}-Z_{\mathrm{g}}),\left[\mathrm{kg}{\mathrm{s}}^{-\mathrm{2}}\right],\end{array}$$

where

$$Z_{\mathrm{v}}=V^{-\mathrm{1}}\underset{\mathrm{0}}{\overset{H}{\int }}A\left(z\right)z\mathrm{d}z\text{ and }{Z}_{\mathrm{g}}=M^{-\mathrm{1}}\underset{\mathrm{0}}{\overset{H}{\int }}A\left(z\right)\mathit{\rho }\left(z\right)z\mathrm{d}z$$

are the heights of the lake center of volume and that of gravity, respectively; $M={\int }_{\mathrm{0}}^H\mathit{\rho }\left(z\right)A\left(z\right)\mathrm{d}z$ is the lake mass; $V={\int }_{\mathrm{0}}^{H}A\left(z\right)\mathrm{d}z$ is the lake volume; ρ₀ is the mean lake density; and $\overline{h}$ is the mean lake depth. Following from the definition, St is a measure of the potential energy excess created by vertical density stratification. Accordingly, St can be directly compared to the kinetic energy input from wind and (or) buoyancy loss at the lake surface due to heat exchange with the atmosphere. In particular, the ratio of the surface wind stress τ to $\text{St}/h_{\text{eq}}$ is an approximate measure of the isopycnal slope due to the wind energy input at the lake surface and may be used to estimate the probability of lake overturn by wind mixing.

The spectra of temperature oscillations were estimated by the Welch method: splitting of the original temperature time series into several 50 % overlapping segments with subsequent tapering of the segments using a Hamming window to reduce the cut-off effects of the limited data length, calculation of the discrete Fourier transform, and averaging the resulting power spectral densities over all segments. The segment length was chosen as 40 d, being significantly shorter than the period of seasonal variations but exceeding the typical synoptic timescales of 5–10 d.

2.4 Modeling

We applied the lake model FLake (Mironov et al., 2010; Kirillin et al., 2011) for simulation of the lake temperature cycle, water mixing, and seasonal ice cover. FLake is a one-dimensional model using a multilayer parametric representation of the vertical temperature structure in the ice–water–sediment system combining integration of the heat budget in each layer with semi-empirical knowledge on the vertical thermal structure. The upper water layer is modeled following the mixed boundary layer approach (Kraus and Turner, 1967; Kraus, 1972), while temperature profiles in the stably stratified part of the water column (hypolimnion), within the ice cover, and in the thermally active upper layer of lake sediments are parameterized using self-similar time-dependent functions of the vertical coordinate. The hypothesis of self-similarity (Kitaigorodski and Miropolski, 1970; Barenblatt, 1978; Zilitinkevich et al., 1979) applied to vertical temperature distributions in layers with quasi-homogeneous physical conditions reduces the modeling approach to solving a set of ordinary differential equations in the time domain instead of discretizing the original partial differential equations along the vertical coordinate. The solution of the initial value problem for the set of three ordinary differential equations is performed using the Newton method with a constant time step. This approach distinguishes FLake from other lake temperature models and provides it with high computational efficiency (Kirillin et al., 2011; Thiery et al., 2014b), allowing flexible applications in long-term scenarios of lake response to climate change or to the anthropogenic change in external fluxes. FLake has proven its robustness and reliability for simulating long-term temperature and mixing conditions in various freshwater lakes and reservoirs (Kirillin et al., 2011, 2017; Thiery et al., 2014a; Su et al., 2019; Almeida et al., 2022).

• Salinity extension. To account for the effects of dissolved salts on vertical mixing and ice formation in the brackish waters of the North Aral Sea, FLake was extended to include salinity effects on the thermal conditions. The contribution of salinity to the vertical density stratification in brackish lakes is typically minor compared to that of the vertical temperature gradients: since the temperature of the maximum density T_md of brackish waters is above the freezing temperature, the lakes are effectively vertically mixed with respect to salinity by convection when surface temperatures cross T_md due to seasonal surface cooling. The salt fluxes by freshwater discharge into the lake, evaporation, and cryoconcentration (salt exclusion due to ice formation) are usually too low to produce significant stratification, although freshwater release from melting ice may affect the upward heat transport from water to the ice base and decelerate the ice thaw. Therefore, salinity was taken as constant in time and depth, i.e., was assumed to be evenly distributed vertically, and the salt fluxes from (to) the atmosphere, ice, sediment, and tributaries were neglected. However, salinity significantly affects the temperature and ice regime by changing the maximum density temperature, thermal expansion coefficient, and freezing point. In the current version of FLake, the three properties were parameterized as functions of water salinity. The assumption allowed the major salinity effects to be incorporated into the bulk of brackish lakes without increasing the model complexity.

  To parameterize salinity effects on water density, we adopted the quadratic equation of state,

  $$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\mathrm{w}}={\mathit{\rho }}_r\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_T{\left(T-T_{\text{md}}\right)}^{\mathrm{2}}\right],\end{array}$$

  where ρ_w is the water density, ρ_r is the maximum density point at temperature T_md, and a_T is an empirical coefficient. Equation (8) is the simplest equation of state that accounts for the fact that the temperature of the maximum density of the water exceeds its freezing point, T_f0=273.15 K (Farmer and Carmack, 1981). Equation (8) enters the model via the buoyancy parameter

  $$-{\displaystyle \frac{g}{{\mathit{\rho }}_r}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial T}}=g{a}_T(T-T_{\text{md}}).$$

  To incorporate the salinity effects on water density, the temperature of maximum density T_md was represented by a function of salinity S [g kg⁻¹ or ‰] following Caldwell (1978):

  $$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{T}_{\text{md}}\left(S\right)=T_{\text{md0}}-\mathrm{0.2229}S,\end{array}$$

  where T_md0=277.13 K is the maximum density temperature of freshwater. The salinity effect on the thermal expansion coefficient α_T was parameterized by simplified polynomials of Bryden (1973) with coefficients adjusted according to Caldwell (1978),

  $$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{\alpha }}_T(T,S)& =[-\mathrm{56.6537}+\mathrm{3.3858}S-\mathrm{0.81218}\\ & \times {\mathrm{10}}^{-\mathrm{2}}{S}^{\mathrm{2}}+\mathrm{14.2407}(T-T_{\mathrm{f}\mathrm{0}})]{\mathrm{10}}^{-\mathrm{6}}.\end{array}\end{array}$$

  The melting temperature was calculated from water salinity after Feistel and Hagen (1998):

  $$\begin{array}{rl}{T}_{\mathrm{f}}\left(S\right)& =T_{\mathrm{f}\mathrm{0}}\\ & +S(-\mathrm{57.5}+\mathrm{1.710523}\sqrt{S}-\mathrm{0.2154996}S){\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

• Model setup. Flake was forced with meteorology from the European Centre for Medium-Range Weather Forecasts reanalysis v5 (ERA5) (Hersbach et al., 2020), including shortwave and longwave radiation, 10 m wind speed, air temperature, and humidity. The lake-specific parameters were modeled on Shevchenko Bay, with salinity of 11 ‰, characteristic fetch of 53 km, and light extinction coefficient of 0.5 m⁻¹ corresponding to a Secchi depth of 3.5 m using the equation of Poole and Atkins (1929). The simulation period was from 1 September 2018 to 31 December 2020, with a 9-month spin-up period and 1 h time step. The initial profile on 1 January was set with a surface and bottom temperature of 0 and 4 °C. Stratification was inferred when the surface-bottom temperature difference exceeded 1 °C, and the length of the longest uninterrupted period was considered its seasonal duration. Model performance was assessed with root mean square error (RMSE) of surface and bottom temperatures as

  $$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(o_i-m_i{)}^{\mathrm{2}}}{n}}},\end{array}$$

  where o_i and m_i are observed and modeled values, respectively. The bathymetry of the North Aral Sea after restoration is not exactly known, and the lake experiences seasonal variations of the water level. Therefore, several model runs were performed with the mean lake depth varying between 6–14 m. To examine potential effects of changes in water level and water quality on the mixing regime, additional series of sensitivity model runs were performed with both mean depth and water transparency varying in the ranges 6–14 m and 0.3–0.9 m⁻¹, respectively.

3 Results

3.1 Salt composition and vertical distributions of major water characteristics by results of 2016–2019 surveys

According to the observational results, the North Aral Sea is currently a brackish lake without pronounced vertical haline stratification in the water column. Salinity of surface waters ranged from 10.71 g kg⁻¹ in September 2018 to 10.68 g kg⁻¹ in September 2019. The differences between surface and bottom values were less than 0.1 g kg⁻¹ in 2018 and about 0.2 g kg⁻¹ in 2019. Note that these salinity values refer to the western part of the North Aral Sea, which is the farthest from the Syr Daria delta. Salinity is expected to decrease eastwards: in 2015, salinity values as low as 9.94 g kg⁻¹ were reported in the mouth area of the Syr Darya River (Andrulionis et al., 2022). Table 1 summarizes the ionic composition of the North Aral Sea before and after the construction of the Kokaral Dam. The data demonstrate a decrease in salinity by almost 7 g kg⁻¹ relative to the pre-Kokaral conditions. The freshening of lake water was accompanied by changes in the ionic composition, mainly associated with an increase in the ${\mathrm{SO}}_{\mathrm{4}}/\mathrm{Cl}$ ratio (Andrulionis et al., 2022). Throughout the period 2016–2019 the vertical distribution of temperature, conductivity, DO, and chlorophyll a was nearly uniform in the upper 12 m of the water column (Fig. 2). In summer (see data from June 2016, Fig. 2a), a stable temperature stratification ∼1 K m⁻¹ develops below 12 m water depth, accompanied by a slight (from 13.78 to 13.83 mS cm⁻¹) downward increase in electrical conductivity, which roughly (±30 %) corresponds to the salinity gradient of 0.02 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ (see, e.g., Pawlowicz, 2008). Herewith, the density ratio (McDougall, 1987) is about 15; i.e., the temperature contribution to stability is 15 times higher than that of salinity. The vertical conductivity gradient in this layer can therefore be interpreted as a consequence of the temperature-driven stratification and subsequent accumulation of dissolved matter near the sediment surface, which does not significantly affect the density stratification. The downward DO decrease of ∼1 $\mathrm{mg}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ within 1 m of the near-bottom stratified layer indicated suppression of the dissolved oxygen supply from the mixed layer by thermal stratification. Fall temperatures, DO concentrations, and conductivities do not reveal any strong vertical gradients (0.025 K m⁻¹ in 2018 and 0.030 K m⁻¹ in 2019, Fig. 2b and c). The slight increase in temperature values near the surface apparently resulted from the diurnal cycle of surface heating. It is noteworthy that the mean values of electrical conductivity continuously decreased from 13.5 in 2016 to 11.85 in 2018 to 11.3 in 2019. While the records might be affected by seasonal and interannual variability of the freshwater budget, they suggest ongoing freshening of the Aral water. The recorded chlorophyll a values of ∼2 µg L⁻¹ are characteristic of oligo-mesotrophic conditions (Carlson, 1977). The Chl a values from fluorescence measurements were not calibrated against standard laboratory methods and may deviate from the real values. Nevertheless, the values are 3–10 times below the typical eutrophic Chl a concentrations, suggesting a mesotrophic state of the lake.

Figure 2 Vertical profiles of temperature (red), electrical conductivity compensated for 25 °C (black), dissolved oxygen (blue), and fluorescence (green) observed in the North Aral Sea in June 2016 (a), September 2018 (b), and September 2019 (c).

3.2 Seasonal thermal regime from continuous mooring observations

According to the autonomous measurements, the temperature cycle of the North Aral Sea is characteristic of holomictic lakes (those entirely mixing at least once during the annual cycle) and is distinguished by two periods of vertical thermal stratification separated by periods of nearly homothermal conditions (Fig. 3a). The period of winter stratification coincides with the period of ice cover and is characterized by a downward temperature increase from the ice base to the sediment surface. In summer, the water column is generally stratified with a temperature drop of up to ∼10 °C from the surface to bottom, whereas the stratification at the measurement site was periodically interrupted by short-term full mixing events. The overall seasonal stratification pattern of the North Aral Sea is thereby between cold monomicitic and dimictic according to the classification of Hutchinson (Hutchinson, 1957). Four periods are distinguishable in the seasonal cycle characterized by specific mixing and stratification conditions.

• Autumn–winter cooling (the fall/autumn overturn, FA in Fig. 3a). This is when convection by surface cooling effectively mixes the entire water column. The period starts in early September and lasts around 3 months until ice cover formation. In our observations, the lake water column stayed thermally homogeneous from the start of the observational period on 21 September 2018 to the satellite-derived date of ice-on on 9 December 2018 (Fig. 4) and from 1 September 2019 to the end of the observational period; it was characterized by a uniform decrease in temperature over time across the entire water column from about 18 to −0.44 °C (which corresponds to the freezing point at ≈8 ‰ for seawater salt composition). During this cooling phase, the daily mean surface-bottom temperature difference ΔT had a median value of 0.00 °C and did not exceed 0.2 °C on 80 % of days, interrupted by one to two day-long events of stratification with ΔT≤1 °C.

• Ice-covered period and winter stratification (WS in Fig. 3a). The ice-covered period lasted from early December 2018 to 31 March 2019 (Fig. 4). First ice formed in shallow areas of the lake on 25–27 November 2018. On 6 December 2018 the ice cover started to form in the area of the autonomous measurement station that was reflected in a fast water temperature drop from ∼0.7 to −0.44 °C within 3 h over the entire water column. On 11 December 2018, a stable ice cover formed, followed by an immediate increase in the near-bottom water temperatures and development of winter “inverse” stratification, characterized by a downward temperature increase. The stratification had a two-layer structure with a thermally homogeneous mixed layer at temperatures about −0.35 °C occupying the upper 10 m of the water column and ∼1.4 K m⁻¹ downward temperature increase down to the water depth of 12 m, where values of about 1.8 °C were observed, which is about the maximum density temperature T_md as determined from the full mixing event in spring (see below). The lower stratified part of the water column revealed quasi-periodic temperature oscillations with periods of 3–9 d (Fig. 4), which can be ascribed to under-ice basin-scale internal waves (see, e.g., Kirillin et al., 2009). The origin and characteristics of these waves are the subject of dedicated research and are not further discussed in the present study.

  The temperature and stratification regime changed remarkably in early March, when the snow cover melted and solar radiation started to warm up the water under the ice cover. From that moment the temperature of the upper water column continuously increased from $<-\mathrm{0.3}$ to 2.0 °C within less than a month. At the moment of ice breakup on 31 March, the heat gain from solar heating (estimated from the mean increase in the heat content in the water column $Q_{\mathrm{N}}=C_p\mathit{\rho }H\mathrm{d}{T}_{\text{mean}}/\mathrm{d}t$) amounted to 80 W m⁻² (Fig. 3B). It is worth noting that temperatures at the uppermost measurement depth of 2.2 m under the surface increased faster than below and created a vertical thermal instability. The thermal instability did not initiate convective mixing (cf. late March data in Figs. 3 and 4) and persisted until the ice breakup. This fact can be interpreted as indirect evidence of the water salinity contribution to the vertical stratification under ice: the apparent source of salt stratification is freshening of the upper waters by the melt at the ice–water interface, which prevents the convection typical for thermally stratified freshwater lakes (Mironov, 2002; Kirillin et al., 2012).

• Spring mixing (SM in Fig. 3a). This lasted in our observations from the ice breakup on 31 March to mid-April 2019. During this period, solar and atmospheric heating at the surface causes the increase in surface temperatures from values close to the freezing point ($\sim -\mathrm{0.4}$ °C) to the temperature of maximum density and results in breaking of winter temperature inversion in the bottom layer. It is worth noting that at the moment of ice breakup the mean water column temperature (∼1.5 °C) is significantly higher than the freezing temperature because of the solar heating under ice in previous months. Melting of the remaining ice floes and heat release to the atmosphere produce short-term cooling and homogenization of the entire water column to ∼0.5 °C on 1 April. Vertical thermal stratification with a downward temperature decrease starts to develop around 9 April, when average water temperatures reach ∼1.8 °C. Accordingly, this temperature represents the temperature of maximum density for the salt composition and concentration of the North Aral Sea waters.

• Summer stratification (SS in Fig. 3a). This lasted from late April to the end of August. After the initial nearly linear vertical thermal stratification formed in late spring, the surface temperatures reached their maximum in late July with daily means of >25 °C and peak values of up to 27 °C. Temperatures in the bottom layer remained generally lower, with temperature differences across the water column of 5–10 °C. However, the thermal stratification at the measurement site was repeatedly interrupted by short-term events of full vertical mixing, indicating penetration of the surface waters down to the bottom (Fig. 5A). The periodic character of these events suggests their wave origin. A deeper insight into the character of wave motions beneath the mixed layer is provided by the spectral analysis of temperature oscillations combined with the theoretical model of basin-scale waves in a circular basin (Eqs. 5–6): the oscillations have periods of 4.5 d, which are slower than the inertial period 2πf⁻¹ (Fig. 5B) and are therefore strongly affected by the Earth's rotation. According to the theoretical model, the period of 4.5 d in a basin with L=16 km (characteristic radius of Shevchenko Bay) refers to the first-mode Kelvin wave with corresponding Burger number S≈0.12 (Eq. 4), Rossby radius $S\times L\approx \mathrm{2}$ km, and wave speed C₁=0.16 m s⁻¹ (Eq. 2). The same value of C follows from the two-layered vertical temperature distribution (Eq. 3) in a 12 m deep basin with a 10 m thick upper mixed layer and water temperatures of 24 and 12 °C in the upper and the lower layer, respectively. These values agree well with the vertical temperature profiles observed at the monitoring site (Fig. 5A).

  For a Kelvin wave with frequencies σ≪f the theory predicts that a counterpart Poincaré wave exists at frequencies σ≥f. A statistically significant peak at the second frequency from Eqs. (5)–(6) corresponding to a Poincaré wave exists in the temperature spectrum, although it is much less energetic than that for the Kelvin wave. On the one hand, this result is evidence that the observed temperature oscillations were produced by a Kelvin–Poincaré pair of basin-scale waves. On the other hand, it demonstrates that the contribution of the Poincaré wave to the oscillations was weak. The latter conclusion is confirmed by previous studies (Antenucci and Imberger, 2001; Kirillin et al., 2009): while Kelvin waves have their maximum energy at the lateral boundaries, with amplitudes decreasing exponentially towards the lake center, the amplitudes of Poincaré waves are concentrated in the horizontal-plane circular motions around the mid-part of the lake and weakly affect the vertical motions of isotherms.

  More insight into external forcing behind the temperature conditions in the North Aral Sea is provided by the rate of change of the bulk heat content in the lake, $Q_{\mathrm{N}}={\mathit{\rho }}_{\mathrm{0}}{C}_p{\int }_{\mathrm{0}}^H\partial T/\partial t\mathrm{d}z$, where C_p is the water heat capacity. The value represents a measure of the net heat flux to the water column (Fig. 3b) and demonstrates strong seasonal flux variations reflected in the thermal stratification. During the autumn cooling the lake lost energy at ∼100 W m⁻² (0.3 K d⁻¹ lake mean cooling rate), with peaks as low as −300 W m⁻². The average heating rates in spring amounted to ∼100 W m⁻², while summer fluxes experienced strong variations with occasional negative values, which can be ascribed to cold water advection by basin-scale internal waves rather than to surface cooling events.

  Apart from the summer stratification period, the water column remained well-mixed vertically with the Schmidt stability (St) close to zero. In summer, St followed the internal-wave-driven temperature fluctuations at the measurement site (Fig. 6), varying between 0 and >100 kg s⁻². The maximum values peaked at St≈130 kg s⁻², which is about 2 times lower than those typical for seasonally stratified dimictic lakes (Read et al., 2011). Hence, a typical surface wind stress of τ≈0.05 N m⁻², roughly corresponding to wind speeds of 4–6 m s⁻¹, would result in the isotherm slope $\mathit{\tau }{h}_{\text{eq}}{\text{St}}^{-\mathrm{1}}\approx \mathrm{1.5}$ m km⁻¹. Approximating the Shevchenko Bay by a cylindrical shape with radius L≈16 km and the height of the surface mixed layer h_mix≈10 m, the ratio $(\mathit{\tau }{h}_{\text{eq}}/\text{St})(h_{\text{mix}}/L{)}^{-\mathrm{1}}\approx \mathrm{2}>\mathrm{1}$, indicating that the slope would produce surfacing of the deep waters at the upwind shore and could potentially destroy vertical stratification completely. The persistence of stratification throughout the summer can be attributed to the irregular bathymetry of the North Aral Sea and the nonlinear slope of the isotherms influenced by the Coriolis effect, which transforms the wind energy to rotational wave motions, with the isotherm slope concentrated within a distance of the Rossby radius L_R≈2 km from the lateral boundaries (Lamb, 1932; Gill, 1982; Antenucci and Imberger, 2001).

Figure 3 (a) Annual temperature cycle in the North Aral Sea in 2018–2019. The thick black isoline is the boundary of the bottom stratified layer defined as $N^{\mathrm{2}}>{\mathrm{2}}^{-\mathrm{3}}$ s⁻². White vertical dashed lines mark the transitions of “lake seasons”: FA refers to fall (autumn) overturn, WS to winter stratification, SM to spring mixing, and SS to summer stratification. (b) Rate of change of the lake heat content (net heat flux into the lake Q_N).

Figure 4 Temperature records between ice-on and ice-off in winter 2018–2019 at four selected water depths. The blue bar on top of the panel marks the ice cover duration. The four MODIS satellite images (NASA Worldview) with corresponding vertical lines across the temperature panel illustrate the ice conditions at the times of lake surface freezing (I–II), snowmelt (III), and ice thaw (IV). The location of the mooring site is shown with the circle on the satellite images.

Figure 5 Summer temperature stratification and internal wave characteristics. (a) Time series of the (1) surface and (2) near-bottom temperatures in summer 2019. The bottom temperatures are complemented by the lowest temperature envelope (3) as a proxy for the temperature of cold deepwater masses ascending to the measurement site. (b) Spectral density of near-bottom (11.8 m depth) temperature oscillations. The red dashed line corresponds to the lower confidence boundary for the spectral peaks defined as the upper 95 % for the red noise signal (first-order autoregressive process) containing the same variance as the observed time series (Gilman et al., 1963). Black dotted lines show the Burger number dependence of the wave frequency S(σ) following Eqs. (5)–(6) and the inertial frequency f. Gray circles connected with gray lines mark the periods of the Kelvin wave corresponding to the spectral energy peak (4.5 d, large circle) and to the associated Poincaré wave (close to the inertial period of 16.5 h).

Figure 6 Annual course of thermal stratification in the North Aral Sea expressed in terms of Schmidt stability (Eq. 7).

3.3 Annual oxygen regime

Our year-long DO record provides first insights into seasonal DO dynamics in the North Aral Sea. It demonstrates that the lake waters are generally well-saturated with oxygen in the annual cycle: during most of the year, the near-bottom DO saturation was around 100 % (Fig. 7A). The absolute concentrations varied between ∼15 mg L⁻¹ in winter and ∼8 mg L⁻¹ in summer, following the seasonal variations in water temperature and the corresponding saturation levels. During autumn–winter cooling, DO concentrations increased from 10 to 15 mg L⁻¹, apparently caused by an increase in solubility due to water cooling (Fig. 7B) and absorption of oxygen from the atmosphere facilitated by strong convective mixing. In the ice-covered period, the DO saturation was close to 100 % and slightly increased after February, indicating possible primary production under ice. From early March, coinciding with the start of under-ice heating by solar radiation (“Winter II” period), the DO concentrations decreased from 15 mg L⁻¹ to the local minimum of 13 mg L⁻¹ (from 105 % to 90 % of saturation) around 15 March. The decrease might be caused by entrainment of undersaturated water from the deep parts of the lake. Quantification of the mechanism requires additional spatially resolved observations. Later, the mean concentrations stabilized at ∼14 mg L⁻¹, while the DO saturation gradually increased to 120 % in mid-May. In summer, the lake remained 100 % saturated with DO, while absolute concentrations decreased following the temperature-related decrease in solubility.

Figure 7 Annual cycle of (a) near-bottom dissolved oxygen and (b) near-bottom water temperature. Values of DO concentration in mg L⁻¹ are shown with a blue line, while values of DO saturation are shown with a brown line. Thick lines are 5 d running averages filtering out short-term DO oscillations. The gray bar indicates the ice-covered period, and the black vertical dashed lines mark the transitions of “lake seasons”: FA refers to fall/autumn overturn, WS to winter stratification, SM to spring mixing, and SS to summer stratification.

After the beginning of solar heating in March, oscillations with amplitudes of 2–4 mg L⁻¹ were present in the near-bottom DO records and persisted after the ice break, with a short interruption during the spring mixing on 15–20 April, resuming in summer. The drops in DO content coincided with the drops in near-bottom temperature (Fig. 8a), and the DO oscillations had the same periods of 4.5 d as the temperature oscillations associated with the Kelvin wave (cf. low-frequency spectral peaks in Figs. 8b and 5b). Hence, the observed DO fluctuations can be ascribed to wave-driven events of colder undersaturated water upwelling from deeper lake areas to the measurement location. These DO dynamics indicate the presence of a water mass with lower oxygen concentrations in the deep central part of the lake during both summer and winter stratification periods. Driven by basin-scale internal waves, these deep waters periodically flowed over the measurement site and caused DO saturation drops down to about 60 %, which, however, were shorter than the cold events (Fig. 8a), suggesting DO replenishment by local primary production or by DO entrainment from the upper saturated waters due to wave breaking. Local primary production is supported by the significant peak at 24 h period in spectra of near-bottom DO variations (Fig. 8b) and a weaker spectral peak at 12 h period, as an indicator of nighttime canceling of production. The up to 130 % oversaturation in May–June can also be attributed to primary production during the spring phytoplankton bloom as its most probable source. The appreciable effect of primary production on DO content at 12 m water depth is a remarkable feature indicating sufficient availability of light and nutrients over the major part of the North Aral bottom. The diurnal oscillations were, however, absent in late winter DO fluctuations (not shown), suggesting limited light availability at the 12 m water depth during this period.

Figure 8 Summer dynamics of the DO content at the bottom slope. (a) Time series of the DO saturation (left axis, black lines) and the near-bottom temperature fluctuations around the summer mean value (right axis, red line). The thick black line is the DO saturation filtered with a 24 h moving average to remove the daily primary production/respiration cycle. Thin dashed lines are the original DO saturation data. Coincident temperature and DO saturation drops indicate upwelling of hypoxic deep water at the mooring site. (b) Spectral density of near-bottom DO content oscillation in June–August. Vertical dashed lines mark the periods of significant energy-containing oscillations. The Kelvin wave period corresponds to the energy-containing peak of temperature oscillations in Fig. 5. The dotted line shows the lower confidence boundary for the spectral peaks defined as the upper 95 % for the red noise signal (first-order autoregressive process) containing the same variance as the observed time series (Gilman et al., 1963).

3.4 Modeled thermal regime

The results of the uncalibrated model FLake validation showed an excellent fit to the observed surface temperature over the whole measurement period (RMSE=0.80 °C). The prediction of the ice-covered period was equally good. The model predicted stable ice duration from 13 December 2018 to 3 April 2019, thus overestimating ice-on and ice-off dates by only 2–3 d and the total ice duration (110 d) by only 1 d. The modeled bottom temperature closely followed the observed bottom temperature (RMSE=2.60 °C) when the mean lake depth was set to 7 m (Fig. 9a), corresponding to the mean depth of the whole North Aral Sea according to satellite altimetry (about 7.3 m). However, the one-dimensional model expectedly failed to reproduce the observed bottom temperature fluctuations with periods of several days, produced by wave-driven advection of water along the bottom slope. As a result, the model predicted a dimictic regime of the bay with stable, continuous summer stratification and vertical temperature differences from 5 to 17 °C, in contrast to the intermittently interrupted stratification that was observed. On the other hand, if the water level was decreased by 1 m to a mean lake depth of 6 m, the model predicted polymictic behavior with the bottom temperature equal to the surface temperature throughout the year. The result suggests a high sensitivity of the mixing regime in the North Aral Sea to water level fluctuations. Sensitivity model runs with a stepwise increase in the lake depth up to 14 m (Fig. 9a) revealed a gradual decrease in the mean bottom temperature in summer from ∼15 to ∼5 °C and a slight increase in the summer stratification duration. We performed a sensitivity analysis to explore the combined effect of changes in water depth and another length scale, the light extinction (Secchi depth), on the thermal regime, where the variable lake depth simulates potential variations in the lake level due to water regulations at the lake catchment and at the dam, and extinction variations reflect the bulk effects of changes in trophic level (productivity) on vertical redistribution of solar radiation over the water column. The thermal regime appears to be very sensitive to both a change in water column depth and water transparency (Fig. 9b), with both polymictic and dimictic regimes possible. The model suggested that the North Aral Sea would be polymictic at lower depth and light extinction and dimictic at higher depth and extinction. Further, the model suggested that the North Aral Sea would switch from a polymictic to a dimictic regime at a lake depth around 6–9 m and corresponding extinction between 0.3–0.8 m⁻¹. With a mean depth of 7 m and extinction of 0.5 m⁻¹, the North Aral Sea is currently close to a transitional state in the mixing regime. The result also suggests that changes in the regional water budget, modifications of the Kokaral Dam that affect water level, or changes in water quality that affect water transparency may considerably affect the thermal regime. For instance, a 1 m increase in depth from 7 to 8 m would increase stratification duration by 40 d and the July vertical temperature difference by about 3 °C according to the model. A depth increase from 10 to 11 m, on the other hand, would increase stratification by 15 d and the July temperature difference by 2 °C.

Figure 9 FLake model results and sensitivity of the thermal regime to lake depth and light extinction. (a) Modeled bottom temperatures at mean lake depths ranging between 6–14 m (colored lines). (b) Modeled mean temperature difference between the surface (T_s) and bottom (T_b) in July as a function of mean lake depth and light extinction. The black contour lines indicate the corresponding summer stratification duration in days. The blue area and the 50 d isoline correspond to the polymictic regime. The model runs from (a) correspond to a horizontal line at extinction coefficient 0.5 m⁻¹.

4 Discussion

Our results demonstrate significant changes in the physical and chemical conditions in the North Aral Sea after construction of the Kokaral Dam. Salinity of the North Aral waters has dropped significantly, by almost 7 g kg⁻¹, and is currently stabilized around ≤11 g kg⁻¹, which is close to the brackish conditions observed during the “natural” period before the shrinkage (Bortnik and Chistyaeva, 1990; Kosarev and Kostianoy, 2010). Apart from the manifold consequences for the biotic components of the lake ecosystem, which are out of the scope of our study, the salinity decrease implies several important consequences for the seasonal mixing and oxygen regime of the lake. These consequences are manifested primarily through the salinity effects on the freezing temperature T_f and the temperature of the maximum density T_md. Our data on vertical mixing allowed these values to be estimated in the North Aral Sea as $T_{\mathrm{f}}\approx -\mathrm{0.4}$ and T_md≈1.8 °C. These estimates correspond to the observed salinity of 10.5 g kg⁻¹ and slightly differ from the average values of −0.57 and 1.60 °C, respectively, given for the entire lake at the salinity of 10 g kg⁻¹ during pre-desiccation conditions (Bortnik and Chistyaeva, 1990). The observed difference can also be linked to a transformation of the salt composition, which occurred during the lake shrinkage and separation of the Aral's residual basins from each other.

In contrast to ocean waters and hypersaline lakes, brackish lakes have a temperature of maximum density above the freezing temperature, T_md>T_f, which causes strong convective mixing during the autumn cooling period and homogenization of the entire water column. As a result, the effect of salinity on the vertical density stratification is minor compared with that of temperature, which prevents the formation of permanent deepwater stagnation (meromixis) as observed in the southern parts of the Aral Sea (Izhitskiy et al., 2021). This fact also implies the presence of warmer waters near the bottom during winter stratification, which is expected to be common in large brackish lakes of arid climate, in contrast to the nearby hypersaline waters of the South Aral Sea (Izhitskiy et al., 2016). Hence, like in the past natural state, autumn overturn is one of the most important processes forming the present hydrological structure of the lake.

On the one hand, relatively low salinity ensures the formation of the seasonal ice cover in winter, which is completely absent in the hypersaline southern Aral Sea (Kouraev et al., 2004; Kouraev and Crétaux, 2010). On the other hand, brackish waters freeze at lower temperatures than freshwater: this is one of the factors explaining low bulk water temperatures during the ice-covered period. Another factor is a strong heat loss during the autumn cooling before the ice formation, typical for the continental (cold desert) climate of Central Asia. In this regard, the autumn and winter regime of the North Aral Sea is similar to that of freshwater lakes of the Tibetan Plateau, where intense surface cooling in autumn is accompanied by strong winds delaying ice cover formation (Kirillin et al., 2021).

Despite the weak vertical salinity gradients in the brackish North Aral Sea, these gradients may still affect vertical mixing at certain stages of the seasonal stratification. Convective mixing during the autumn–winter cooling occurs in two stages: the first stage begins with the onset of surface cooling, characterized by strong convective mixing, which continues until the mixed water layer reaches the temperature of maximum density. Afterwards, convection due to the surface heat release ceases. In freshwater lakes, the inverse stratification starts to form at this moment (Kirillin et al., 2021). After the water temperatures reach the freezing point and the ice cover starts to form, the second phase of convective mixing begins, driven by the density increase near the surface due to salt exclusion by ice formation (Pieters and Lawrence, 2009). This effect is more pronounced in brackish waters and may contribute to vertical mixing in early winter (cf. Stage I in Fig. 4). Vertical salinity gradients have an opposite effect in late winter with the onset of under-ice heating by solar radiation. In freshwater lakes, the resulting increase in the surface water temperatures produces convective mixing (Matthews and Heaney, 1987; Mironov, 2002), which is prevented in the brackish conditions by freshening of surface waters due to melting at the ice base, similarly to the oceanic ice boundary layer (McPhee, 1992). Our modeling projections, while reliably simulating the seasonal thermal stratification pattern, did not take these salinity gradients into account. Incorporating the salinity effects of vertical salt gradients may improve short-term model predictions of the vertical stratification in brackish lakes, especially in the ice-covered period and during the spring overturn, when the density dependency on temperature vanishes and the near-bottom salinity gradients may prevent convective mixing at the water–sediment interface (Mironov, 2002). It is worth noting that the period of under-ice heating in the North Aral Sea lasts for almost 2 months, which is significantly longer than in temperate and polar lakes (Mironov, 2002; Bouffard et al., 2019). The long duration of this phase of the ice-covered period (sometimes termed “Winter II”) is a distinguishing feature of low-latitude lakes, which gain a larger amount of radiation than typical temperate and (sub)arctic ice-covered lakes. In this regard, the North Aral Sea shares the features of the lakes of Tibet and the Mongolian Plateau, which are characterized by strong heating of the under-ice water column (Kirillin et al., 2021; Huang et al., 2022). Accordingly, the late winter is a crucial period in these lakes, with a warm and quiet environment under ice favorable for planktonic primary production and fish.

The cold arid climate is characterized by strong seasonal variations in the heat flux at the lake–atmosphere interface, resulting in a wide seasonal range of lake temperatures. Before the desiccation, the range of seasonal temperature fluctuations at the Aral surface did not exceed 23–25 °C (Bortnik and Chistyaeva, 1990), whereas now this range exceeds 27 °C. This is apparently due to the significant decrease in the lake volume and, as a result, the lower heat storage capacity of the reservoir. This fact also explains the more intense cooling and heating processes observed during the autumn and spring mixing period: the net rates of surface cooling in autumn and surface warming in summer exceed 200 W m⁻², which is several times higher than the values reported in temperate lakes (Henderson-Sellers, 1986).

The relative shallowness of the North Aral Sea determines another characteristic feature of its seasonal mixing regime: despite the strong surface heating, the summer stratification is relatively weak, as demonstrated by the low Schmidt stability. The latter suggests that moderate winds could completely destroy the summer stratification and could turn the lake polymictic, i.e., fully mixed vertically during most of the year (Hutchinson, 1957). However, our analysis of seasonal temperature time series in time and frequency domains reveals a 2–3 m thin near-bottom layer with temperatures decreasing downwards, which is persistent throughout the summer period. Due to the decrease in the surface area with depth, the relative volume of the layer is only about 7 % of the entire lake. As a result, the wind-induced basin-scale internal waves, which are inevitably produced in the stratified water column, move the entire cold near-bottom water along the bottom slope instead of mixing it with the overlying warm mixed layer (see the middle panel in Fig. 10). Thanks to the large horizontal dimensions of the lake and the low water column stability, the period of these waves is rather long, reaching several days, which significantly exceeds the inertial period at the lake latitude. Hence, the Coriolis effect captures the wind energy in the rotational motion instead of vertical mixing (Lamb, 1932; Gill, 1982). Thereby, the large horizontal dimensions and the weak thermal stratification prevent the lake from becoming polymictic. It should be noted that our observations are confined to the western part of the North Aral Sea, the Shevchenko Bay. However, taking into account the similar morphometries of the Shevchenko Bay and the central part of the Aral Sea, the same wave pattern should dominate the mixing in the bulk of the sea except shallow bays in the north and at the river estuary.

Figure 10 Sketch illustrating the current state of the vertical mixing regime in the North Aral Sea and possible scenarios of its transformation depending on the water level fluctuations. h_mix is the summer mean thickness of the upper mixed layer with a high level of turbulent mixing due to wind shear and thermal convection, H is the lake depth, and O₂ is the dissolved oxygen saturation level. The gray “stagnation” areas are occupied by the stratified hypolimnion, which currently oscillates with high amplitudes. A decreased water level potentially leads to complete vanishing of the hypolimnion. If the water level rises, the hypolimnion volume will increase and the wave amplitudes will decrease, additionally reducing bottom mixing.

Apart from preventing full vertical mixing, the basin-scale waves play a potentially important role in the water–sediment mass exchange by enhancing shear mixing in the bottom boundary layer and increasing the supply of dissolved oxygen to the sediment surface and transport of dissolved nutrients from the sediment to the water column. This suggestion is supported by our data on near-bottom DO concentrations: the waters at 12 m depth remain well-saturated with oxygen in summer, while short-term wave-driven drops in concentrations to ∼60 % of saturation indicate the existence of a low-oxygen volume in the deepest part of the lake, which, however, is small relative to the entire volume of the lake. The specific thermal and mixing conditions also determine the oxygen regime in winter: the DO saturation under ice is about 100 % throughout the ice-covered period, indicating no current risk of strong hypoxia in the North Aral Sea, although hypoxic conditions can still locally develop in the deepest parts of the lake. This fact distinguishes the North Aral Sea from the majority of temperate and boreal lakes, where winter hypoxia is a common phenomenon with negative consequences for lake ecosystems (Malm, 1998; Leppi et al., 2016; Steinsberger et al., 2020; Perga et al., 2023; Schwefel et al., 2023). A distinct feature of the under-ice oxygen regime in the North Aral Sea is the high absolute DO concentration, reaching ≥15 mg L⁻¹. These values are 15 %–25 % higher than typical concentrations in freshwater lakes at higher latitudes (Terzhevik et al., 2010; Bengtsson, 2011; Bouffard et al., 2013) and are due to the low under-ice temperatures of the brackish lake water and the corresponding high DO solubility. Therefore, brackish seasonally ice-covered lakes, which are common for the endorheic zone of Eurasia, are more effective in storing oxygen under ice and are less prone to the oxygen deficit in winter than temperate freshwater lakes.

Summarizing our results, we can conclude that the North Aral Sea is now in a transitional mixing regime between dimictic and polymictic. This state is intrinsically unstable and is prone to a quick regime shift. As demonstrated by the model experiments, such a shift can be triggered by variations of the water level or by a change in the water transparency due to, e.g., an increase in primary production. The sensitivity model runs showed that under the currently observed stratification, a change in depth of 2 m can shift the mixing regime and potentially alter stratification duration by up to 2 months. This magnitude of water level variations is not unrealistic: according to recent assessments of the variability of water balance conditions in the North Aral Sea, the discharge of the main inflow, the Syr Daria River, will undergo significant fluctuations in the future (Ayzel and Izhitskiy, 2018, 2019). Combined with the variability of evaporation and precipitation, this could lead to significant changes in the lake surface level. According to climate scenarios, the lake surface may drop by several meters by the end of the century (Izhitskiy and Ayzel, 2023). In this case, a change in the seasonal mixing regime to polymictic (Fig. 10) can be expected. While the risk of deep hypoxia would be minimal in this case, negative consequences for the ecosystem may emerge as a strong increase in primary production and eutrophication of the lake due to increase in the nutrient concentration.

On the other hand, the water authorities are considering further restoration measures including raising the Kokaral Dam crest by another 6 m (Izhitskiy and Ayzel, 2023). The resulting water level increase would enhance the vertical thermal stratification, eventually turning the seasonal mixing regime to dimictic with the summer stagnation period lasting for several months (Fig. 10). This scenario suggests development of deep hypoxia or even anoxia with long-lasting oxygen-free conditions near the lake bottom. These possible changes in mixing regime can have severe consequences for the biogeochemical cycle of the lake. Currently, concentrations of dissolved methane in the North Aral Sea exceed atmospheric equilibrium by a factor of 13 in the southern part of the lake and by a factor of 48 in the area near the Kokaral Dam (Izhitskaya et al., 2019). High inflows of riverine organic matter (Klimaszyk et al., 2022) and its reduction in anaerobic microenvironments and sediments in combination with a possible increase in stratification may lead to a strong increase in anoxic methane production near the lake bottom, similar to that currently observed in the western South Aral Sea and Lake Chernyshev (Izhitskaya et al., 2019).

Our measurements of Chl a concentrations suggest the North Aral Sea to be currently in a mesotrophic state. The mesotrophic character of the lake is also supported by our Secchi disk measurements of 3.5 m and by the trophic index estimates based on inorganic phosphorus and nitrogen by Klimaszyk et al. (2022). These authors also found a strong increase in the nutrient concentrations towards the river inlet, which indicates a high nutrient supply by the inflow. The trophic state provides a strong feedback on mixing conditions via the phytoplankton effect on water transparency (Shatwell et al., 2016). As demonstrated by our sensitivity model experiments, this feedback may hypothetically reveal itself in scenarios of both increased and decreased water level. While the water level drop forces polymictic conditions, it can intensify primary production, which may in turn reduce transparency and produce re-stratification of the water column with a shallow mixed layer depth (Kirillin and Shatwell, 2016). Such conditions are favorable for cyanobacteria blooms, deteriorating lake water quality (Bartosiewicz et al., 2019). On the other hand, a possible water level increase and subsequent transition to dimictic conditions can increase nutrient retention in the lake and provoke a transition to eutrophic conditions (Huisman et al., 1999) with a transparency decrease and a positive feedback on stratification. The latter can again reduce the mixed layer depth and increase the probability of plankton blooms. Hence, the nonlinear feedbacks between water level, water transparency, and stratification suggest manifold possible scenarios of lake ecosystem changes. The present mixing regime, which prevents strong deep anoxia and ensures a deep surface mixed layer, appears in this regard to be a favorable one compared to both dimictic and polymictic regimes.

5 Conclusions

Terminal lakes are sensitive indicators of anthropogenic activity, such as water use for irrigation, and of the climatically driven changes in the hydrological regime of endorheic areas. The area of the Aral Sea catchment is 3×10⁶ km² (Zavialov, 2007), which makes it the second largest endorheic catchment after the Caspian Sea. In this regard, the results of the Aral Sea conservation experiment represent a unique opportunity to assess the outcomes of large-scale lake manipulation for its thermal and mixing regime. The evident decrease in salinity and restoration of the biological communities including the fish population after the dam construction were acknowledged as indicators of the success of the conservation experiment. While the observed temperature and salinity regime of the North Aral Sea is now similar to the pre-desiccation state of the lake, our results demonstrated distinct mixing characteristics on seasonal and sub-seasonal timescales, which make the regime unstable. This unstable mixing regime favors, on the one hand, heat and mass exchange at the lake bottom but can be quickly disturbed by further discharge manipulation or by the consequences of global change.