Interactive comment on “ The internal seiche field in the changing South Aral Sea ( 2006 – 2013 )

The authors use two datasets from the largest residual basin of the former Aral Sea coupled with numerical modeling to analyse the internal wave weather in this recently formed natural water body. The basin-scale internal oscillations are known to be the major mediator of the energy transport to the bulk of the water column in enclosed stratified basins. In this regard, the study makes an important contribution into the fundamental understanding of the vertical energy transport in large strongly stratified lakes. The reported findings are especially relevant to dynamics of extreme aquatic environments, where strong vertical stratification is created by salinity gradients. The analysis is based on observational data separated by several years, allowing to trace the changes in the internal wave dynamics caused by continuous decrease of the water level and strengthening of the density stratification. The example of the Aral Sea is both unique and extendable on other saline lakes and seas subject to dessication/temporary


Introduction
The growing knowledge about surface and internal standing waves in lakes (seiches), beginning with the first documented observations of surface oscillations in Lake Michigan in the 17th century, has been exhaustively reviewed by Hutter et al. (2011).The modern era of internal seiche study began with the Defant-Mortimer model (Mortimer, 1979), which can be extended to high vertical modes, that is, when a lake responds as a multilayer system.Although it is now accepted that high vertical modes are often excited, observations of such modes were sparse until the end of the 20th century (Heaps, 1961;LaZerte, 1980;Csanady, 1982;Hutter et al., 1983).Currently, the importance of internal wave fields in redistributing wind energy within lakes is well known (Wüest et al., 2000;Stocker and Imberger, 2003;Shimizu and Imberger, 2008) and different authors have focused on its impact on mixing (Stevens, 1999;Planella et al., 2011;Bernhardt and Kirillin, 2013), sediment resuspension (Bogucki and Redekopp, 2008) and sediment and phytoplankton transportation (Ji and Jin, 2006;Rolland et al., 2013;Vidal et al., 2014).While studies of internal waves in small and large lakes are now quite common, for most of them, the characteristic internal seiche field has not been described in detail.Large lakes where internal seiches have been studied, either theoretically or experimentally to varying degrees, include the Kinneret (Ou and Bennet, 1979;Antenucci et al., 2000), Lake Champlain (Prigo et al., 1996), Lake Geneva (Lemmin et al., 2005), Lake Biwa (Saggio and Imberger, 1998) and Lake Michigan (Mortimer, 2004).
This paper presents the first study of the internal seiche field in the Aral Sea, which is located in Kazakhstan and Uzbekistan (central Asia).At present, it covers approximately 10 % of the surface area it covered in the 1960s and holds less than 10 % of its former volume, which is now mainly contained within two separate water bodies: the North and South Aral seas (also known as the Small Sea and the Large Sea).In 2009, the shallow western lobe of the South Aral Sea dried up completely, and after partial replenishment, it has dried up again (UNEP-GEAS, 2014).However, what remains of the South Aral Sea, a shadow of its former self, is still a large lake about 150 km long and 25 km wide, on average.
The Aral Sea's desiccation has had an enormous human, ecological and climatic impact (Micklin, 2007;Arashkevich et al., 2009;Zhitina, 2011;Rubinshtein et al., 2014) and has become a paradigm of the large number of lakes all over the world whose water levels have dropped (Bai et al., 2011;Yildirim et al., 2011;Lauwaet et al., 2012;Tourian et al., 2015).The ongoing changes in the Aral Sea are being closely followed by the scientific community (Singh et al., 2012;Schettler et al., 2013;Shi et al., 2014).Since 2002, the evolution of its stratification has been documented (Zavialov, 2005, Zavialov et al., 2009;Izhitskiy et al., 2014a, b), the predominant circulation has been described (Izhitskiy et al., 2014a) and the role played by geomorphological processes associated with hydrodynamics in the evolution of the lake has been discussed (Roget et al., 2009).However, due to a lack of data, only very preliminary works exist on internal waves.In 2013, during a field campaign that was not focused on the study of internal waves, continuous temperature data within the water column were recorded for the first time, making it possible to use the Princeton Ocean Model (POM) to assess the numerical results for the internal seiches.The POM, which is widely known, has been used to make marine forecasts for the US Great Lakes (NOAA-GLERL, 2016), to study internal waves (Rueda et al., 2003;Munroe and Lamb, 2005;Babu and Rao, 2011) and to study the Aral Sea (Roget et al., 2009;Izhitskiy et al., 2014a).
Here, we present a numerical study of the internal seiches for the conditions in autumn 2006 and 2013, and compare them with the measured data.More precisely, we analyze the modes of oscillation shorter than the inertial period in the region, which is about 17 h, and longer than 3 h.This range has been established by considering the length of the measured data series and the uncertainties in the simulation setup.This is discussed in the material and methods section, which has been organized into different subsections presenting a description of the measurements taken that are relevant to this study, the numerical simulations run for the two campaigns and the methods of analysis.
The results section is also organized into subsections.In the first subsection, the spectral analysis of the field data is compared with the numerical results predicted by the model at the same stations.Subsequent subsections analyze the structures of the modes identified in the spectral analysis of the field data in 2006 and 2013.The discussion section is organized into two subsections concerning the horizontal and vertical structure of the internal seiches; in both cases, the features that favor the excitation of different modes, the differences between the two campaigns and the possible relevance of these variations on the lake mixing dynamics are discussed.2 Materials and methods

Site and measurements
The field data analyzed to support the numerical simulations were recorded during two field campaigns (27-30 September 2006 and 29 October-3 November 2013) described in detail by Zavialov et al. (2008) and Izhitskiy et al. (2014b), respectively.Figure 1 presents the shorelines of the South Aral in 2006 and 2013 together with the bathymetric levels for 2013, when the sea surface level decreased by 3.2 m.The shoreline has mainly receded along the eastern shore, where the lake is shallower.Furthermore, the neck connecting the northern part of the South Aral Sea (Chernyshev Bay) to the main body of the lake has narrowed considerably.The measuring stations in 2006 and 2013 are also indicated in Fig. 1.The choice of sampling sites was limited due to difficult access to the sea and resulting logistical restrictions.
In 2006, a Nortek Aquadopp acoustic Doppler current meter was deployed for 3 days at a depth of 39 m at station A2, in the deeper part of the lake.The uncertainty of the measurements was about 1 % of the measured value ±0.5 m s −1 .In 2013, bottom velocities were recorded with a SeaHorse tilt current meter (TCM) (Sheremet, 2010) at a depth of 25 m at station W1, close to the western shore.The TCM was deployed in the high current mode (accuracy range from 0 to 80 cm s −1 ) with a velocity resolution of 0.1 cm s −1 and a direction resolution of 0.1  numerical results.The resolution and accuracy of the temperature sensors was 0.032 and ±0.1 • C, respectively.All the data series recorded in 2013 and analyzed here are from a 2-day period.During the campaigns, other moorings were deployed but only during shorter periods, so they had to be disregarded.Wind speed was recorded with a portable automatic meteorological station installed on the western shore of the lake near station W1.
In 2006, and again in 2013, all data were recorded continuously at a sampling rate of 30 s, but for this study they have been averaged every 10 min.At the beginning of both campaigns, a conductivity-temperature-depth (CTD) profile was taken at station A2 and water samples were collected at different depths.For each campaign, salinity values were obtained from samples processed in the laboratory using the standard dry rest method.Next, a linear regression was calculated for each year between the sample salinity values and the corresponding salinity obtained from the CTD data using the UNESCO formula.These regressions were used to obtain a salinity profile from the CTD data before density profiles were calculated using the Aral Sea density formula proposed by Gertman and Zavialov (2011).
The characteristic temperature, salinity and density profiles for both campaigns are presented in Fig. 2. The stratification was weak in 2006 and strong in 2013, but this is merely circumstantial; it depends solely on the time when the campaigns were carried out.In general, this hypersaline lake becomes saltier but the maximum annual density gradient decreases (Izhitskiy et al., 2014a).The elongated horizontal dots on the 2013 profile (Fig. 2c) indicate the depths at which the temperature sensors in the chain were placed.

Numerical simulations
Numerical simulations for 2006 and 2013 were performed using the POM, which is based on hydrostatic threedimensional primitive equations with the Mellor-Yamada closure scheme and terrain-following sigma coordinates (POM, 2016).The density formula for the current Aral Sea obtained by Gertman and Zavialov (2011) was included in the model.The simulation was initialized with the temperature and salinity profiles measured in the field.The Aral Sea bathymetry used has a 967 m north-south and 538 m eastwest resolution (Roget et al., 2009).For the vertical coordinates, 12 equidistant sigma levels were considered in 2006 and 17 in 2013.The model has a free surface and a split time step.For the simulations presented here, the time steps for the external and internal modes were 10 and 500 s, respectively, both of which fulfill the Courant-Friedrichs-Lewy criterion.However, given the uncertainties introduced into the simulations by using equal initial stratification conditions over the whole lake, and considering the limitation of the measured data series for comparison, the output of the model was fixed for every hour.
Wind stress on the surface was calculated as τ = ρ a C D V 2 10 with the drag coefficient, C D , as defined by Hasselmann (1988).A similar approach was used for the bottom stress, where the drag coefficient was calculated to fit the velocities at the first grid point nearest the bottom boundary with the logarithmic law of the wall taking a bottom roughness of 0.01 m into account.For all simulations, the heat flux at all boundaries was set to zero and no mass input or output was considered.Zero normal velocities were used as the lateral boundary conditions.
In 2006, the model was forced by hourly wind data downloaded from the coupled global NCEP Climate Forecast System Reanalysis.Data were interpolated from a Gaussian to a regular geographic grid with a 0.31250 • step and later to the grid model itself.The model ran for 10 days, but the results were only analyzed after the fifth day, coinciding with the experimental field campaign.During this time, there was a mild northeasterly wind, which is the prevailing wind direction in the area.For the 2013 simulation, the model was forced over the first 5 days by a constant and spatially uniform easterly wind of 3 m s −1 , which coincides with the mean speed and the dominant direction prior to and during the campaign, and afterwards by the wind measured in situ.
The internal seiches are, by definition, standing waves and they depend only on bathymetry and stratification.In fact, standard models for the study of the internal seiches find the characteristic roots (frequencies) and characteristic vectors (structure of the vertical displacements) of the momentum balance equations by considering the appropriate boundary conditions for stationarity but with no external forcing (Salvadé et al., 1988;Guyennon et al., 2014).However, the seiches excited from all the possible stationary waves do depend on the pattern of the wind that has generated them (Fricker and Nepf, 2000).That is why we have forced the POM model with the characteristic wind prior to and during the campaigns.For this work, a 5-day spin-up period of the model was found to be enough for the standing waves to develop, although accurately settling the circulation of the lake might require a larger spin-up period (Korotenko et al., 2010).Lorrai et al. (2011) also used a spin-up period of 5 days of the 3-D hydrostatic Boussinesq model to study the seiching dynamics of a medium-size lake.Note that if the POM model had not reproduced the modes studied in this paper, the surface elevation filtered around the periods of the observed oscillations would not present a coherent structure over the lake that is repeated for every period, as presented in the results section.

Methods of analysis
The horizontal structure of the internal seiche field is analyzed based on surface level deviations from the equilibrium in the numerical results.This method assumes that internal oscillations can be traced by the corresponding oscillations on the lake surface with a characteristic amplitude, ζ s , which can be approximately related to the amplitude of the displacement of internal the density interface (halocline), ζ h , as ζ s = ((ρ h − ρ e ) /ρ h ) ζ h , where ρ e and ρ h are the densities at the surface (epilimnion) and bottom (hypolimnion) layers (Hellström, 1941).In our study, which considers the density stratification during the 2006 and 2013 campaigns, the amplitudes of the internal waves at the halocline induced by the corresponding wind fields are, respectively, about 1000 and 100 times larger than at the surface.
To study the modes of the internal seiches identified in the spectral analysis of the field data, the surface elevation predicted by the model is band-pass filtered around the corresponding frequency (or period).When the oscillation is re-produced by the model, the filtered surface elevation shows the horizontal structure of the standing waves, which is repeated after a complete oscillation.The method was initially proposed by Mortimer (1963) and has subsequently been used by different authors, including Caloi et al. (1986), Sirkes (1987), Lemmin et al. (2005) and Forcat et al. (2011).This approach avoids the interpolation of the numerical results in the water column.Note that the periods of the internal seiches identified at the surface are in a range completely different from the periods of the surface seiches (Kirillin et al., 2015).In general, the mode number of the horizontal oscillation corresponds to the number of lines where the surface elevation is zero (nodal lines).
In this study, the numerical results are compared with the spectra of the vertical position of the isotherms (in 2013) and of the bottom velocities (in 2006 and 2013) measured in the lake.Velocity components along and across the lake were computed considering its 14 • tilt with respect to the northsouth direction.Limitations resulting from the length of the field data forced us to focus this study on short modes.
For the 2013 campaign, the density and the temperature profiles were monotonous in depth.Therefore, the vertical displacements of the isopycnals (depths of equal density) were obtained from the position of the isotherms (depths of equal temperature), which were calculated by linear interpolation of the temperature values measured at fixed depths.To analyze the measured and numerical data, spectral, filter and correlation analyses were carried out with standard MATLAB functions.
The vertical structure of the identified modes of the internal seiches is analyzed from the numerical results by bandpass filtering the horizontal velocity around the corresponding frequency (period).This is done in all grid depths at a fixed point.If the oscillation is excited, the velocity sign within each layer is the same and changes at the interfaces between them.Note that if the velocities are band passed around a period which does not correspond to a mode reproduced by the model, no coherent structure regarding the velocity sign (direction) is observed in the water column.In general, the lake behaves as a multilayer system, so there can be multiple interfaces whose number is one less than the number of layers.The number of interfaces defines the order of the vertical mode.

Results
Here, we analyze the internal seiches in the South Aral during the 2006 and 2013 campaigns.As already mentioned, in 2006 the stratification was weak and there was a moderate north wind, while in 2013 the stratification was strong and there was a light northeasterly wind.Note, however, that the density step between 25 and 30 m of depth in 2006 (Fig. 2c) would be equivalent to a temperature step of more than 5  in a freshwater lake and would be considered an example of strong stratification (Boehrer and Schultze, 2008).Regarding the data measured in 2006, the spectra of the bottom velocity components along and across the main axis are presented in Fig. 4. As it can be seen, there are energy peaks at around 11, 7 and 5 h, which agree with the spectrum of the numerical results presented in Fig. 3a.The shape of the spectrum at the lower frequencies in the component along the lake is affected by the window of 128 points used to give more statistical relevance to the shorter modes.Figure 5 shows the spectral analysis of the field data recorded in 2013 based on the vertical displacements of the 6, 7 and 10 • C isotherms (Fig. 5a) and on the bottom velocities (Fig. 5b).In both figures, peaks of energy are observed at frequencies corresponding to periods of around 14, 7 and 4 h, in agreement with the spectra of the numerical results for 2013 (Fig. 3b).As expected, the confidence intervals, displayed on the left side of the plots, are relatively large due to the short length of the experimental data series.However, the results are still robust, especially if we consider that the spectral analysis of two completely different measurements shows the same results, which are also in agreement with the numerical simulations.Given these results, the structure of the oscillations of 11, 7 and 5 h in 2006, and of 14, 7 and 4 h in 2013, are analyzed in the following sections based on the outputs of the numerical simulations mentioned earlier in the material and methods section.For 2006, a longer mode of 36 h predicted by the simulations is also discussed in the following section.

The structure of the internal seiches in 2006
Based on the simulations, a standing oscillation of 36 h was identified in 2006.The horizontal structure of this mode can be observed from the bypass-filtered surface elevation around a frequency of 1/36 h −1 .The corresponding plots from a reference time (t = 0 h) and after 18 and 36 h, respectively, are presented in Fig. 6a, showing a complete oscillation.As observed in all three plots, the surface elevation shifts smoothly along the lake and, at about 45.38 • N, there is a line where the surface elevation is zero (nodal line).At one side of the nodal line the surface level is above the equilibrium level (positive values) and at the other side it is below it.The direction of the vertical displacements changes from one plot to the next when the time shift is half the oscillation period.Because this oscillation has a single nodal line it is a horizontal first mode.The mode of 36 h, however, was not observed in the measurements.This may be because the short length of the data series does not allow its resolution or because, as discussed by Imam et al. (2013), the first horizontal modes in multi-basin lakes can be strongly damped because of the drag in the constricted areas between the basins.
The structure of the mode of 11 h identified in the spectra of the field data in 2006 can be observed in Fig. 6b, where the simulated surface elevation bypass filtered at around 11 h is also presented for three different times covering the oscillation.As in Fig. 6a, there is also a single nodal line.This line is mainly transversal to the main axis of the lake and located around 45.40 • N. Accordingly, the 11 h oscillation is also a first horizontal mode that oscillates longitudinally; that is, the surface elevations at the north and south sides of the transversal nodal line have opposite signs.However, in the southern part of the lake, for the mode of 11 h, the lines of equal surface elevation are mainly longitudinal, and therefore the wave presents a transversal structure.The transversal character of this oscillation at the southern part of the lake is in accordance with the results shown in Fig. 4, where a peak at 11 h is observed in the transverse component of the velocity measured at station W1.
Analogously, the oscillation of about 7 h observed in the spectra of the bottom velocity in 2006 can be interpreted, in agreement with the POM results, as a second horizontal mode because the corresponding horizontal oscillation traced in the surface layer presents two nodal lines.Fig. 6c shows a northern nodal line located between 45.74 and 45.92 • N, close to the entrance to Chernyshev Bay, and a second one at about 45.05 • N. In this case, however, the nodal lines are less well defined than those in the figures of all the other analyzed modes.This may be the result of the output of the model being fixed every hour, but the simulated oscillation has a period of 7.5 h (Fig. 3a).
The oscillation of about 5 h identified in the velocity spectra for 2006 is also reproduced by the model as a horizontal standing wave, represented in Fig. 6d.In this case, the oscillation presents clear two nodal lines, so it is identified as another second horizontal mode.As observed, the nodal lines for the mode of 5 h are located close to those of the mode of 7 h.
Regarding the vertical structures of the modes observed in 2006, they are analyzed from the evolution in time of the filtered horizontal velocity along the lake, at all depth grids of the water column in the deeper part of the lake, as discussed in the methods section.Despite the low velocity values filtered around the corresponding period, the velocity sign changes only at a few fixed depths that coincide with the interfaces between the layers with different velocity directions.The location of these interfaces remains constant over time, showing the persistence of the mode.
Figure 7a presents the velocities filtered around 36 h.In the figure, a single discontinuous horizontal line at about 14 m depth shows the single depth where the horizontal velocity changes direction.So, for this mode, the lake responds as if formed by two layers moving in opposite directions.Accordingly, the mode of 36 h is a first vertical mode.Because at the same time it is also a first horizontal, this mode is known as the fundamental mode.
Figure 7b-d show the vertical structures of the modes of 11, 7 and 5 h based also on the filtered horizontal velocities.In all cases, two dashed horizontal lines indicate the interfaces between the layers: located about 7 and 22 m for the mode of 11 h, at 5 and 26 m for the mode of 7 h and at 12 and 27 m the mode of 5 h.Accordingly, for these oscillations, the lake responds as if formed by three layers moving in opposite directions and the modes are consistent with the excitation of second vertical modes.In summary, we can conclude that during the 2006 campaign the fundamental mode of internal seiches was 36 h and that the measured oscillation of 11 h was a first horizontalsecond vertical mode, while those of 7 and 5 h were both second horizontal-second vertical modes with different structure.

The structure of the internal seiches in 2013
According to Mortimer's (1953) formula, the ratio of the periods of the fundamental internal seiche for two different stratifications should be equal to the inverse ratio of the square root of their normalized density difference at the halocline.Given that the normalized density differences in 2006 and 2013 were, respectively, 2/1077 and 12/1088 (Fig. 2c), the period of the fundamental mode in 2013 is expected to be about 2.4 times less than in 2006, when it was shown to be 36 h.So, it is reasonable to consider that the peak of around 14 h at station W1, observed in both the numerical and the experimental spectra, corresponds to the fundamental mode, which is also the more easily excited mode.However, the 14 h oscillation observed in the spectra of both the field and the numerical data at station W1 cannot be identified over the whole lake from the bypass-filtered surface elevation as was the case for the oscillation of 36 h in 2006 (Fig. 6a).Note from Fig. 3b that the peak at around 14 h at station W1 shows two maxima, which may indicate two oscillations with similar periods but different structures that cannot be isolated from each other because the output of the model was fixed for every hour.In fact, the location (45.4 • N, 58.6 • E) represented in Fig. 3b coincides with the nodal line of the mode of 36 h in 2016 (Fig. 6a) and therefore, if we assume that the oscillation of 14 h corresponds to the fundamental mode in 2013, no vertical displacements are expected at this location.Furthermore, because of the limited length of the field data series, the resolution might not be high enough to determine two close oscillations from the measured data, which in Fig. 5 also show a wide energy peak at around 14 h.All in all, it can be concluded that the 14 h peak observed in both the numerical and the field data spectra contains energy of the fundamental mode.
Figure 8a, b present the surface elevations bypass filtered around the 7 and 4 h periods identified in the spectra of the field data in 2013 (Fig. 5).In both cases, the filtered surface elevations are shown for three different times covering the corresponding periods.In Fig. 8a, the whole lake is observed to oscillate longitudinally with a period of 7 h and two nodal lines, indicating that it is a second horizontal mode.The first nodal line is located between 45.78 and 45.98 • N, at the entrance to Chernyshev Bay (see Fig. 1), and the second, which has a crosswise character, extends from 44.98 • N at the western shore to 45.38 • N at the eastern shore.As presented in Fig. 8b, the 4 h standing wave has three nodal lines located at around 44.78, 45.78 and 45.98 • N, at the entrance to Chernyshev Bay, corresponding to a third horizontal mode.
The vertical structures of the coexisting modes in 2013 are observed from the evolution in time of the filtered horizontal velocities presented in panels b and c of Fig. 9.In Fig. 9b, which corresponds to the 7 h mode, there are three horizontal lines at depths of about 6, 16 and 26 m, where the horizontal velocity is zero (nodal levels) and so where the interfaces between layers are located.Therefore, the mode of 7 h is a third vertical mode and the lake behaves like a four-layer system.Temperature data recorded at fixed depths within the water column support the idea that the 7 h mode is a third vertical mode.Considering the depths of the temperature sensors along the thermistor chain, we were able to obtain the vertical oscillation of the isotherms of 11 and 6.8 • C. Figure 9a shows the depths of these isotherms which, as observed in Fig. 9b, are located in the upper and lower parts of the third layer, which should alternately approach and retreat and therefore move out of phase because of the standing oscillation (Pérez-Losada et al., 2003).Unfortunately, the measuring station W1 was located relatively close to the horizontal nodal line so that the vertical displacements at this location are expected to be relatively small.However, even in this case, and considering the short length of the recorded data series, the correlation function between the 11 and 6.8 • C isotherms, presented in Fig. 10, shows a minimum when the time shift is zero, confirming that both isotherms oscillate out of phase.No temperature data were recorded within the other layers, but because the data available corroborate the location and the behavior of one of the predicted layers, and the period of the observed and numerical oscillations coincides, the excitation of the third vertical mode is indirectly corroborated.As for the vertical structure of the mode of 4 h, Fig. 9c shows that it corresponds to a second vertical mode.That is, the lake behaves as a three-layer system and, at two depths within the water column (about 15 and 25 m), the horizontal velocity is zero because it changes its direction.
All in all, we can conclude that the fundamental mode in 2013 was of 14 h, that the standing oscillation of 7 h corresponds to a second horizontal-third vertical mode and that of 4 h corresponds to a third horizontal-second vertical mode.

Horizontal structure of the internal seiches
In 2006, two first horizontal modes (including the fundamental mode) and two second horizontal modes were found to be excited, while in 2013 one first (the fundamental mode), one second and one third horizontal mode were found to be excited.So, in general, the horizontal modes in 2013 were higher.The development of nodal lines at locations where the transversal section is relatively small in comparison to the mean transversal area of the lake is well known (Roget et al., 1997;Imam et al., 2013).According to our results, the decrease of 3.2 m in the surface level between 2006 and 2013 seems to have favored the development of a nodal line at a seal near the entrance to Chernyshev Bay, where the depth was less than 5 m in 2013.Note that the development of a nodal line between Chernyshev Bay and the main body of the lake is not equivalent to the fact that Chernyshev Bay would have become detached from the South Aral Sea: while maximum horizontal velocities are found at the nodal lines, they are zero at the coast (Hutter et al., 2011).According to Izhitskiy et al. (2016), the separation of Chernyshev Bay from the rest of the basin is imminent, so the distribution of velocities affecting the exchange between different sub-basins (Umlauf and Lemmin, 2005) and the mixing activity (Vidal et al., 2013;Guyennon et al., 2014) in the actual Large Aral Sea could change relatively soon.
Although the predominance of the fundamental mode is generally expected, in 2006 the 36 h fundamental mode predicted by the model was not identified from the short series of data.Furthermore, its amplitude in the simulations was smaller than the amplitude of the other modes (see the scales in Fig. 6).This situation could be explained by considering the work of Imam et al. (2013), who showed that, in an elongated large lake, the first horizontal mode is strongly damped by the drag in the region between the sub-basins.The damping of the fundamental mode of 14 h in 2013 could also explain why, apparently, it could not have been separated from other modes.
During the 2006 and 2013 campaigns, which were carried out during two periods when there were mild northeastern and northern prevailing winds, all the observed modes have been identified as longitudinal modes; that is, they oscillate along the main axis of the lake which is orientated about 14   from the north-south direction.In 2006, however, the longitudinal mode of 11 h also had a transversal structure at the southern part of the lake and in 2013 the longitudinal mode of 7 h had a crosswise nodal line encompassing 0.4 • of latitude.The importance of the wind pattern in determining the standing waves that become excited from the spectrum of possibilities has been illustrated by different authors (Valerio et al., 2012) and the excitement of longitudinal seiches with a transversal structure at the wider and deeper lobe has also been documented in a small lake as the standard response to an oblique or transversal wind forcing (Roget et al., 1997).

Vertical structure of the internal seiches
Although only first (the fundamental mode) and second vertical modes were measured in 2006, a third vertical mode was excited in 2013 (Fig. 9b).Comparing the 2013 density profile in Figs.2c and 9b shows that the wide and strong halocline located between about 15 and 25 m of depth behaves as a layer by itself.Opposite to that, the main halocline in 2006, which extends from a depth of 25 to 28 m (Fig. 2c), coincides with the interface between the layers of the mode of 4 h (Fig. 9d).Although the density gradient of the main halocline in 2006 is apparently small in comparison with 2013, it is equivalent to a variation in temperature of 5 • C in 3 m, which in a freshwater lake would be considered as a sharp stratification.
The importance of the stratification to the development of high vertical modes has been discussed in detail by Pérez-Losada et al. (2003) based on the bulk Richardson number.Also, Boehrer et al. (2000) and Bastida et al. (2012) presented how variations in the wind produce changes in the order of the vertical modes.The results presented here show that, at least under the light breeze conditions prior to and during the campaigns, the deep quasi-homogeneous mixed upper layers (∼ 10 −2 kg m −4 ) sustain internal waves which increase the number of vertical modes.Otherwise, in 2006, when there is sharp pycnocline, the lower mode would be a first vertical, and in 2013, when there is a thicker pycnocline, a second vertical mode.
During the last decade, the water exchange between the eastern basin (which has vanished by now) and the western basin of the South Aral increased the salinity of the western lobe up to more than 120 ppt in the bottom layer and favored a layered structure of the water column (Zavialov, 2005) which is not always present in the density profile.The 2016 profiles (Fig. 2) show that the first 13 m of the upper mixed layer were the warmest and saltiest and that, from 13 to 25 m, the temperature and salinity decreased below the values of the bottom layer.In 2013, however, the water exchange between the western and eastern lobe almost ceased, and the upper and bottom layers looked relatively homogeneous.Recent measurements (Izhitskiy et al., 2016) suggested sporadic inflows through the remaining channel into Chernyshev Bay, which is saltier than the bottom layer in the central part of the eastern basin.Thus, saltier water inflowing from Chernyshev Bay into the main body of the Large Aral Sea can still influence its evolution.Therefore, the variation in the annual stratification of the Large Aral Sea is still expected and, as a result, the internal wave field could change.Several authors have shown the effect of changing modes on the spatial variability of turbulence.Lorrai et al. (2011) have shown how the mixing rates present variability at twice the seiching frequency.Sakai and Redekopp (2011) have also shown that, in a lake where super-inertial rotatory waves dominated, the superposition of higher modes enhanced particle transport.Accordingly, given the expected variation on the annual stratification of the Large Aral Sea, the effective basin-scale diffusivities could also be expected to change.

Conclusion
Our results show that the fundamental modes in 2006 and 2013 were of 36 and 14 h, respectively, and they might be dampened by the drag in the region between the main subbasins.In autumn 2006, when the lake was relatively weakly stratified (up 0.2 kg m −4 ), there was a second vertical-first horizontal mode of 11 h and two different second verticalsecond horizontal modes of 7 and 5 h.In autumn 2013, when the stratification was very strong (up to 3.2 kg m −4 ), there was a third vertical-second horizontal of 7 h and a second vertical-third horizontal mode of 4 h.So, only in 2013 was a third vertical mode excited, favored by the wide strong halocline that behaves as a layer by itself.In both years, a deep quasi-homogeneous mixed upper layer (∼ 0.01 kg m −4 ) was found to sustain internal waves, at least under light breeze conditions.
All the studied modes were longitudinal; however, favored by northeasterly or northern winds, the first horizontalsecond vertical mode of 11 h in 2006 had a transversal structure at the southern part of the sea, and the second horizontalthird vertical mode of 7 h in 2013 presented a crosswise nodal line.
This study points to the relevance of sea level variation to the excitation of higher horizontal modes.Shallower areas connecting different sub-basins can become critical for nodal line development, as is the case of the neck connecting the northern part (Chernyshev Bay) with the main body of the South Aral Sea.Considering that, in the current climate scenario, the water level in many lakes varies, it is important to remember that in stratified, shallow lakes, water level variations can induce great modifications in the dynamics of the lake due to a sudden change in the standing oscillations of horizontal velocities and vertical displacements, which could influence the distribution of energy mixing and phytoplankton exposure to sunlight.
This work is a first approach to the seiches in the presentday Aral Sea and it raises several aspects which could be of relevance to other stratified water bodies subject to desiccation.However, further investigation is needed to have a more holistic view of the internal wave field in the large Aral Sea, including larger periods, rotation effects and consideration of the heterogeneity caused by complex bathymetry.

Data availability
Experimental data are archived at the Laboratory of landocean interactions and the anthropogenic impact (P.P. Shirshov Institute of Oceanology -RAS) and available upon request.The simulations with the POM model implemented in the Aral Sea are available at the University of Girona.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Shoreline of the South Aral Sea in 2006 and 2013 and bathymetric levels in 2013.Locations of the measurement stations (station W1 in 2013 and station A2 in 2006) are indicated.

Figure 2 .
Figure 2. Characteristic profiles of (a) temperature, (b) salinity and (c) density during field expeditions in 2006 and 2013.Elongated dots on the 2013 density profile show the location of the temperature sensors.

Figure 3 .
Figure 3. Spectra of the simulated surface level (a) in 2006 at station A2 and (b) in 2013 at station W1 and at 45.4 • N, 58.6 • E.

3. 1
Figure3apresents the spectra of the simulated surface elevation in 2006 at measuring station A2.It shows two peaks of energy (at around 11 and 7 h) that cannot be completely separated, and a third peak at 5 h.The 95 % confidence interval shown in the figure indicates that the peaks are statistically significant.Figure3bpresents the equivalent spectra in 2013 at measuring station W1 and at another location (45.4 • N, 58.6 • E).Both spectra show peaks at 7 and 4 h, and at station W1 there is also a wide peak at around 14 h.All the peaks are statistically significant.Regarding the data measured in 2006, the spectra of the bottom velocity components along and across the main axis are presented in Fig.4.As it can be seen, there are energy peaks at around 11, 7 and 5 h, which agree with the spectrum of the numerical results presented in Fig.3a.The shape of the spectrum at the lower frequencies in the component along the lake is affected by the window of 128 points used to give more statistical relevance to the shorter modes.Figure5shows the spectral analysis of the field data recorded in 2013 based on the vertical displacements of the 6, 7 and 10 • C isotherms (Fig.5a) and on the bottom velocities (Fig.5b).In both figures, peaks of energy are observed at frequencies corresponding to periods of around 14, 7 and 4 h, in agreement with the spectra of the numerical results for 2013 (Fig.3b).As expected, the confidence intervals, displayed on the left side of the plots, are relatively large due to the short length of the experimental data series.However, the results are still robust, especially if we consider that the spectral analysis of two completely different measurements shows the same results, which are also in agreement with the numerical simulations.

Figure 4 .
Figure 4. Spectra of bottom velocities recorded at station A2 in 2006.

Figure 5 .
Figure 5. Spectra of (a) the vertical displacements of different isotherms and (b) bottom velocities measured at station W1 in 2013.

Figure 6 .
Figure 6.Band-pass-filtered surface elevations at three different times covering one oscillation of the internal seiches of (a) 36 h, (b) 11 h, (c) 7 h and (d) 5 h in 2006.

Figure 7 .
Figure 7. Filtered horizontal velocity along the main axis of the lake for internal seiches of (a) 36 h, (b) 11 h, (c) 7 h and (d) 5 h in 2006.The grey horizontal lines indicate the interfaces between layers.

Figure 8 .
Figure 8. Band-pass-filtered surface elevations at three different times covering one oscillation of the internal seiches of (a) 7 h and (b) 4 h in 2013. •

Figure 9 .
Figure 9. (a) Measured temperature profile in 2013 showing the depths at 6.8 and 11 • C. (b) Filtered horizontal velocity along the main axis of the lake for the internal seiche of 7 h and (c) filtered velocity for the internal seiche of 4 h in 2013.The dashed horizontal lines indicate the interfaces between layers.

Figure 10 .
Figure 10.Correlation function between the vertical displacements of the isotherms of 11 and 6.8 • C measured at station W1 in 2013.