Abstract and keywords
Abstract (English):
The paper presents the quantitative assessment of the hydrodynamic parameters of the convectively mixed layer (CML) arising in ice-covered boreal lakes as a result of non- homogeneous heating of the water column. The study is focused on revealing the features of CML spatial structure and calculating all six elements of the turbulent stress matrix. The main feature of the experimental technique is application of two rigidly coupled acoustic current profilers (ADCP) installed on ice and operating in the asynchronous measurement mode. In the case of down-looking axes of both devices, their positions relative to each other were chosen so that one or two pairs of beams intersected at a certain depth. Due to this configuration, it was possible to perform a rigorous calculation of all turbulent stresses based on the data obtained from all six beams. The intensities of pulsations were estimated along all three axes, complete with error analysis. A high degree of anisotropy of the pulsations and a periodic nature of its change over time were detected. An analysis of the dynamics of average velocities carried out for depths of up to 2.87 m with discreteness of 2.5 cm (for CML thickness of 3–6 m) revealed the presence of quasi-deterministic convective cells. In the horizontal plane, a systematic “drift” was found, due to the presence of large-scale geostrophic circulation. The presence of such a drift made it possible, as a first approximation, to convert the experimentally obtained Eulerian characteristics into Lagrangian ones and, accordingly, draw conclusions about the spatial structure of the cells. In particular, based on the analysis of progressive vector diagrams in vertical planes, the depth ranges were determined at which the zones of up- or downwelling prevailed. The distribution of these zones is of great importance in studying the spatial dynamics of plankton. KEYWORDS: Boreal lakes; under-ice radiation; convectively mixed layer; acoustic Doppler current profilers; quasi-deterministic structures; progressive vector diagrams; up- and downwelling zones; turbulent stresses.

Boreal lakes; under-ice radiation; convectively mixed layer; acoustic Doppler current profilers; quasi-deterministic structures; progressive vector diagrams; up- and downwelling zones; turbulent stresses.
Publication text (PDF): Read Download


Convective processes, along with tides and internal waves, play a decisive role in the mixing of water masses. At that, as compared to oceans, where the main source of convection is surface cool-

ing, such sources in lakes may have a different nature [Bouffard and Wu¨est, 2018]. In particular, the main driver of convection in ice-covered lakes is inhomogeneous heating of the water column. This phenomenon usually occurs during late winter or early spring stages, when, after the snow cover has melted, a significant part of the solar radiation penetrates the ice. Before the onset of radiative warming, the water column in ice-covered lake is stable and inversely stratified. Inhomogeneous radiation heating of the water column in a temperature range from 0C to the maximum density temperature 𝑇md (3.98C) gives rise to vertical anomalies of the water density, which leads to the sinking of warmer and heavier liquid volumes, and the development of convective movements. A convectively mixed layer (CML) is gradually formed, its lower boundary gets deeper over time, and the temperature approaches 𝑇md. The duration of this process varies from weeks in shallow lakes to several months in large ones [Austin, 2019; Jewson et al., 2009].

In a broader context, the problem of convection

developing under the action of internal heat sources arises in the study of a number of geophysical phenomena [Lepot et al., 2018] and can be attributed to one of the fundamental challenges in the study of turbulence in a stratified fluid. Interest in its solution is due, in particular, to the peculiarities of turbulent heat transfer and the search for a universal relation connecting the Nusselt (Nu) and Rayleigh (Ra) parameters. As shown in [Bouillaut et al., 2019; Lepot et al., 2018], it is just for this type of convection that the asymptotic variant of transport («ultimate regime») is realized, in which the exponent 𝛾 in the dependence Nu Ra𝛾 is close to the value ½. Another important feature is the abnormally large value of the mixing coefficient Γ which characterizes the distribution of the incoming energy between the potential and internal energy during the mixing process. For example, in ice-covered lakes the energy of the negative buoyancy flux is 65% converted into potential energy [Ulloa et al., 2018], while in the energy balance with other types of convection, the corresponding share is only 20% [Maffioli et al., 2016].

As to convection in ice-covered lakes, it was

largely due to the study of this phenomenon [Farmer, 1975; Jonas et al., 2003; Mironov et al., 2002] that the concept of the winter period as a passive interval of the annual thermal and biological cy-

cle of lakes was revised [Pernica et al., 2017]. The emergence and development of CML, as shown in a number of recent studies [Yang et al., 2017], significantly affects not only the transformation of the temperature profile, but also the distribution of oxygen and plankton. This period also largely governs the subsequent evolution of plankton, and also affects the characteristics of the entire annual temperature and oxygen cycles. In this regard, the tasks related to the study of the structure and dynamics of the CML in ice-covered lakes are still of high relevance and importance. The aspects of special interest are the rate of CML deepening, the features of its hydrodynamic structure, estimates of the coefficients of turbulent transfer, and the ratio of the upand downwelling zones.

As shown by Kelley [1997], even where the CML is modeled by a system of ordered cells, the distribution of plankton in the photic zone has a rather complex form, determined by the features of the mean velocity field. Moreover, the hydrodynamic parameters of the cells can also have a differentiated effect on different fractions of plankton [Palshin et al., 2019; Woodward et al., 2019]. In turn, reliable estimates of the coefficient of vertical turbulent diffusion are required to study the stability of plankton populations in a turbulent CML [McNeill, 1990].

For all the urgency of the problems associated with the study of the CML hydrodynamic parameters, the methods for their quantitative assessment remain underdeveloped. Already in the early works [Woodcock, 1965], CML self-ordering was noted, which manifests itself in the existence of convective cells. However, the question of the degree of such ordering and regularity of the CML spatial structure remains largely open. The issues related to estimates of turbulence parameters and transport coefficients in this layer are even less developed.

These problems are largely due to the complexity of setting up experiments in ice-covered lakes. Especially noteworthy are difficulties in identifying the horizontal structure of the CML: in most cases, the experiments are carried out in the moorings format, when vertical profiling of a certain parameter is carried out only in the place where the device is located. In recent years, significant progress in overcoming these problems has been achieved with the use of underwater gliders [Austin, 2019; Forrest et al., 2008]. This has enabled, in particular, esti-

mating the size of the cells and their aspect number.

The use of acoustic Doppler current profilers (ADCP) also has a significant potential in dealing with these problems. In particular, these instruments can be used to calculate averaged currents in the CML and some parameters of turbulence. The results of such calculations based on measured data from Petrozavodsk Bay of Lake Onega (March 2017) and from Lake Vendyurskoe (March and April of 2016 and 2020) are presented e.g. in [Bogdanov et al., 2019, 2021; Bouffard et al., 2019; Volkov et al., 2019].

This paper presents further elaboration of the methods developed in the above-mentioned publications and associated, first of all, with the analysis of progressive vector diagrams and the calculation of pulsation parameters from the data of two rigidly coupled ADCPs. The primary goal was to design an upgraded method for calculating all six com-

ponents of the turbulent stress matrix < 𝑢𝑖𝑢𝑗 >the main parameters that determine the large-scale

structure of turbulence.

Study Site and Experimental Setup. Materials and Methods

The CML studies were carried out in the period from March 27 to April 6, 2020 in a small and shallow Lake Vendyurskoe (6212′54′′ N, 3315′50′′E, surface area is 10.4 km2, average depth is 5.3 m). The instruments were located on ice, on the distance of 350 m from the northern shore, in an area with 7 m depth. The measuring complex included three «star-shaped» pyranometers, a thermistor chain with 13 temperature sensors and two ADCPs. The specificity of the experiment consisted in the use of two rigidly coupled (sideby-side) ADCPs operating in asynchronous regime with a lag of 30 s. The detailed experimental setup, the instruments’ regimes, and the measuring procedures were described by Bogdanov et al. [2021]. In the first stage of the experiment (SC setup — Single Crossing — from 17:00 on March 27 to 09:30 on March 30), the ADCPs emitters were located at a distance of 𝑙1 = 1.5 m from each other, and at a depth of ℎ = 𝑙1/(2 tg 𝛼0) 1.61 m two beams intersected. In the second stage (Double Crossing, DC, from 10:20 am on March 30 to 10:00 am

Figure 1. The DC experimental setup. Panel (a) presents the frame of reference and the numbering of beams for the second ADCP. Panel (b) illustrates the relative position of the devices by both side view and above view. The velocity measurements were carried out at points A, B, C, D. Dashed lines correspond to the beam-crossing depth.

on April 6), the distance between the emitters was half that, and two pairs of beams intersected at the same depth (1.61 m). The scanning range in depth was 2.875 m (115 cells with a size of 25 mm), while the points of intersection of the beams were at a depth of 1.61 m for both SC and DC setups. The main features of the DC setup are shown in Figure 1.

The data from two horizontally separated ADCPs were used to directly check the horizontal homogeneity of the currents — the main condition ensuring the adequacy of the vertical profiling of the hydrodynamic parameters. The main advantage of using two instruments consists, however, in the possibility of calculating turbulent stresses directly. A note to be given from the start is that when using single ADCPs (even the 4-beam and 5-beam versions, «Janus» and «Janus +«), it is not possible to extract stresses from the data on directly measured «beam» velocity components: the number of independent equations is less than six. Moreover, in the general case, it is possible to obtain only relations between turbulent stresses, but not explicit expressions for them. For example, when using a 3-beam device, it is not possible to obtain an explicit expression for any element of the stress tensor. In turn, a 4-beam device allows one to obtain

only two explicit expressions («variance method» [e.g., Howarth and Souza, 2005; Nystrom et al., 2007]) for the off-diagonal components of the matrix, for example:

< 𝑏′2 > − < 𝑏′2 >

advantages of both methods, is presented in detail below, but for the sake of clarity, only the DC setup is considered.

The beam velocities 𝑏𝑖 of both devices are directly expressed in terms of the orthogonal veloc-

< 𝑢′1𝑢′3 >=

1 3.

2 sin 2𝛼0

ity components at the same point, for example, for beam 4 (see Figure 1):

Here, 𝑏′1 and 𝑏′3 are the pulsation components of

«radial» velocities on opposite beams 1 and 3 of a

Janus-type device, 𝛼0 — is the angle of deflection of the beam from the vertical.

To overcome these difficulties two alternative approaches were proposed in [Vermeulen et al., 2011] and [Bogdanov et al., 2021]. Both are based on

𝑏4 = 𝑢𝑥(𝐴) sin 𝛼0 + 𝑢𝑧(𝐴) cos 𝛼0. (1)

This transformation does not involve a loss of information. Moreover, due to linearity, relation (1) and similar ones for the rest of the beam velocities retain their form also for averaged (< 𝑏𝑖 >) and pulsation (𝑏′ = 𝑏𝑖− < 𝑏𝑖 >) quantities.

{ }

the use of two coupled ADCPs. Vermeulen et al. 𝑖 ′2

[2011] used two 4-beam instruments, one of them («master») installed vertically, and the axis of the second one («slave») was tilted. The last condition is necessary when using the method: with the vertical orientation of both devices, the system of equations is singular: only five out of eight equations are linearly independent. This feature of the instrumental setup imposes significant restrictions on the application of the method: the scanning area in the horizontal plane increases significantly, and acceptable results are obtained only at rather large (more than 20) tilt angles for the slave ADCP, which is not recommended by the device manufacturer.

In the second approach [Bogdanov et al., 2021], these difficulties can be avoided, since the axes of

The intensities 𝑏𝑖 of pulsations of all six

beam velocities, as well as correlations of the form

< 𝑏′2𝑏′6 > and < 𝑏′3𝑏′5 >, can be directly calculated from experimental data. On the other hand, all these quantities can be expressed in terms of the

sought-for components of the stress matrix using six relations of the form (1). The resulting (overdetermined) system of eight equations for six unknowns can be represented in the following compact form:

𝐵𝑖 = 𝑠2M𝑖𝑗𝑅𝑗/4, 𝑖 = 1, . . . 8; 𝑗 = 1, 6. (2)

For convenience, we introduce the following «vectors»:

𝐵 = (< 𝑏′2 >, < 𝑏′2 >, < 𝑏′2 >, < 𝑏′2 >,

both instruments remain parallel. In this case, 𝑖

5 6 2 6 3 5

however, a rigorous calculation of all six compo-

1 2 3 4

nents of the stress tensor is possible only at a depth corresponding to the point of intersection of the

< 𝑏′2 >, < 𝑏′2 >, < 𝑏𝑏′ >, <𝑏𝑏



𝑅 = (< 𝑢′2 >, < 𝑢𝑢

>, < 𝑢𝑢′ >,

In addition, the system of equations for stresses

turns out to be overdetermined: in the DC setup

𝑖 𝑥

𝑥 𝑦

𝑥 𝑧

e.g., eight equations are obtained from the two ADCPs data (six intensities of beam velocity pulsations and two correlations of beam components at the intersection points of the beams) for six unknowns. Formally, the problem can be solved by

< 𝑢′2 >, < 𝑢𝑢′ >, < 𝑢′2 >).

𝑦 𝑦 𝑧 𝑧

The matrix M in (2) is determined only by the constructive angle 𝛼0:

removing any pair of «extra» equations from con-

4 √0 8𝑐 0 0√



sideration (as was done in [Bogdanov et al., 2021]),

but this procedure leads to the loss of useful infor-

mation. To minimize such a loss of information,

1 2 3 −4𝑐 3 −4𝑐 3 4𝑐
M =

√0 −8𝑐 0 0√




√ √

1 −2 3 −4𝑐 3 4𝑐 3 4𝑐2


one can, following [Vermeulen et al., 2011], seek an approximate solution of the overdetermined system of equations using the least squares method. Such

a procedure, which, in our opinion, combines the

1 2 3 4𝑐 3 4𝑐 3 4𝑐

√ √

1 −2 3 4𝑐 3 −4𝑐 3 4𝑐


⎜⎝ −1 0 0 3 −4𝑐√3 4𝑐2

−1 0 0 3 4𝑐

3 4𝑐2

In formulas (2) and (3) the notations 𝑠 = sin 𝛼0,

𝑐 = ctg 𝛼0 is used for brevity.

Using the standard procedure for finding an approximate solution to the overdetermined system

(2) by the least squares method, we obtain the following explicit expression for the «vector» of stresses 𝑅𝑖:

of velocity components on time calculated for each of the devices practically coincide. This fact can be taken as direct evidence of the horizontal flow homogeneity. It should also be noted that the dynamic curves show elements of quasi periodicity. The oscillations are not sinusoidal, but rather follow a relaxation, discontinuous pattern.


= 4V


𝐵𝑗/𝑠2, 𝑖 = 1, . . . 6, 𝑗 = 1, . . . 8. (4)

All the features noted above are clearly manifested in progressive vector diagrams (PVD). Ele-


Here, V = (M𝑇 M)−1M𝑇 is a matrix of coefficients with dimensions 6 8; the upper index 𝑇 is used to denote the transposed matrix.

The calculation of errors was carried out according to the method developed in [Stacey et al., 1999; Williams and Simpson, 2004], using the Gaussian approximation for the distribution of beam velocity pulsations, as well as the assumption that there are no correlations of these velocities for different beams.


During the experiment, the maximum values of the under-ice solar radiation flux varied from 100 W/m2 in the first four days of measurements to 50–60 W/m2 on April 3–6. During April 1 and 2, radiation practically did not penetrate under the ice, due to cloudy weather and snowfall. Accordingly, the CML depth varied over the measurement period non-monotonically, from 3 to 6 meters, while the scanning zone of the device was always within this layer.
Mean Velocity

The magnitude and direction of the mean velocity also varied significantly over time intervals of several hours, with the horizontal components reaching 6 mm/s, and the vertical 3 mm/s. At a first glance, the temporal variation of horizontal velocity components had a generally rather irregular character, in particular, no significant correlation with the radiation flux was found. Meanwhile, the root-mean-square values of the pulsation velocities were comparable with the values of the averaged ones.

At the same time, the analysis of the dynamics of velocity components reveals some features that indicate the quasi-ordered nature of the mean velocity field. Firstly, the curves of the dependence

ments of self-organization are revealed when analyzing the sets of these diagrams corresponding to different depths (Figure 2). The family of lines for different depths demonstrates a certain consistency of changes in velocity components with time, which may be due to the presence of quasi-deterministic structures of the CML.

Special attention should be given to one more feature of the PVD plotted in the horizontal plane (Figure 2c): at all the scanned depths, there was a systematic «drift» with a speed of  1 mm/s, and the direction close to NW. The presence of such a drift can be explained by the presence of large-scale geostrophic circulation [Palshin et al., 2017]. The calculated direction of such a geostrophic drift is approximately along the shoreline, and this direction, taking into account the location of the measuring complex, is consistent with the sign (counterclockwise) of such a circulation.

PVDs in two vertical planes (Figure 2a, Figure 2b) also illustrate a systematic drift, but its sign and magnitude change monotonically with depth. A drift of this kind can be caused by the asymmetry of convective cells and, accordingly, the difference in the area of the upand downwelling zones. Options for quantifying such drift and a discussion of its possible nature are presented below.

In general, the above features of the dynamics of the mean velocity components confirm the conclusions of [Bogdanov et al., 2019]: the CML can be considered as a system of quasi-deterministic convective cells or energy-containing eddies. Moreover, as shown in [Volkov et al., 2019], the geometric and energy parameters of these cells can be estimated based on a statistical analysis of jumps in the magnitude and direction of the mean velocity. Such an analysis is largely enabled by the presence of an average horizontal drift. It is due to such a drift that the initial Euler characteristics can be converted, at some approximation, to Lagrangian ones, and, accordingly, the PVD and ve-

Figure 2. Progressive vector diagrams (PVD) for data from the first device; a — on East — Downward plane, b — on North — Downward plane, c — on East-North horizontal plane. Panel d — PVD on the horizontal plane for both devices, depth = 1.6 m. Panels a, b, c present the depth range from 0.12 m to 2.87 m, with 0.25 m discreetness.

locity hodographs can be considered as images of quasi-deterministic cells obtained as a result of their «scanning» by the device in a reference frame related to the average drift velocity. In such a frame of reference, the water column as a whole is at rest, while the «velocity sensors» cross the cells and scan them, in a manner of speaking. In this case, the velocity jumps, in particular, serve as markers of cell boundaries.

Turbulent Stresses

All components of the turbulent stress matrix were calculated directly by formula (4). The averaging period (100 min) was chosen based on the

estimate of the turnover time 𝜏 ℎ/𝑤* of energycontaining structures in the CML, ℎ is the thickness of the CML, 𝑤* is the so-called convective velocity, calculated in a standard way from the buoyancy flux [Mironov et al., 2002].

The reliability of calculations, as noted in [Vermeulen et al., 2011], is largely governed by the value of the determinant 𝐷 of the coefficient matrix

𝑀𝑇 𝑀: for small values of 𝐷, the solutions, as a rule, are not physically consistent. Vermeulen et al. [2011] used a flexible configuration of two devices in which, due to the tilt of the second («slave») device, the value of 𝐷 could be varied. Physically acceptable results were obtained in this case only at large values of the tilt angle, which yielded

𝐷 values of the order 10−4. For the device con-

Figure 3. The calculated dynamics of turbulent stresses; a — pulsation intensities along axes 𝑌 and 𝑍, b — off-diagonal elements (shear stresses), c — pulsation intensities along axis 𝑋 and standard deviation 𝜎𝑥. Vertical lines correspond to noons. Time reading — from 00:00, March 30.

figuration considered in our work, such variation is impossible: the value of 𝐷 is determined only by the design angle 𝛼0. The calculated 𝐷 value, however, was 0.000687, which is almost an order of magnitude higher than the value presented in [Vermeulen et al., 2011]. As a consequence, physically adequate results were obtained for most of the considered time interval (see Figure 3a, and Figure 3c): positive definiteness of the diagonal components (intensities of pulsations along the orthogonal axes). As for the off-diagonal elements of the matrix («shear stresses»), the computation demonstrates their oscillating dynamics with a change in sign (Figure 3b).

For five days of the DC setup (March 30, 31 and April 3–5) with a high level of under-ice radiation, the intensities of horizontal velocity pulsations varied within 3–12 mm2/s2, which is in qualitative agreement with the calculations of Bogdanov et al. [2021]. The intensity of vertical pulsations was an order of magnitude lower: the average value of the anisotropy parameter 𝜁 =< 𝑢′2 > / < 𝑢′2 > over the entire measurement period with the DC setup was 0.1, with the ratio of root-mean-square pulsations along the 𝑍 and 𝑋 axes close to 0.28.


The errors in calculating the pulsation intensities along the 𝑋 and 𝑌 axes were 14% and 23%, respectively, and the accuracy turns out to be

Figure 4. Dynamics of anisotropy tensor eigenvalues. Dashed horizontal lines represent the physical limits for both eigenvalues. Vertical red and blue lines correspond to noons and midnights. Time reading — from 00:00, March 30.

somewhat higher than when using the calculation method based on an incomplete dataset [Bogdanov et al., 2021]. As for the intensity of vertical pulsations, the error was much greater in both methods, amounting to 30–40%. It should be noted, however, that such a large error is mainly due to the violation of the realizability conditions (positive definiteness of intensities and the fulfillment of the Cauchy-Schwarz inequalities < 𝑢𝑖𝑢𝑗 > 2 <

𝑢2 >< 𝑢2 >) at night, when the intensity of verti-

value. The calculated dynamics of these eigenvalues is shown in Figure 4.

Discussion and Conclusions

The paper presents a method (hereinafter: method II) which enables a fairly reliable estimation of all six components of the turbulent stress matrix.

𝑖 𝑗

cal pulsations was extremely low.

The calculated dependence of the stresses on time generally confirms other conclusions of [Bogdanov et al., 2021]: the daily correlation of the intensities of the pulsations with the radiation flux and the time lag ( 1.5 hours) of the system’s response to external forcing.

Yet another feature of the dynamics of pulsation intensities should be highlighted: its oscillating pattern with a period of several hours. These oscillations show no clear correlation with the radiation flux and are observed even at night. It is noteworthy that these oscillations are the most pronounced not for the intensities themselves, but for the eigenvalues of the stress matrix and its dimensionless deviatoric (trace-free) analog, the anisotropy tensor

𝑎𝑖𝑗 =< 𝑢𝑖𝑢𝑗 > / < 𝑢2 > −𝛿𝑖𝑗/3 (here, < 𝑢2 > is the full pulsation intensity). This tensor has only two independent eigenvalues, for which we have chosen those (𝜆1 and 𝜆2) that have a larger

Unlike a similar procedure (hereinafter: method I)

described in [Bogdanov et al., 2021], the use of method II does not entail a loss of useful information obtained from two ADCPs. In this case, however, the calculations are based on the analysis of an overdetermined system of equations and, accordingly, represent only its approximate solution. In general, both methods can be considered as complementary to each other. Each of them can be used as an additional tool for verification and refinement of the results. In particular, method II allows obtaining more accurate estimates for the intensity of pulsations along the 𝑌 axis. When using the alternative method I, just this component of the stress matrix that is the most problematic, due to the «unfavorable» direction of 𝑌 in the given instrument configuration. Besides, the method II is exceptionally easy to use: all stresses are calculated directly from experimental data using a single formula (4). As for vertical pulsations, the errors in calculating their intensities by the two methods

are very similar. This can be explained by the fact that the positive effect from the inclusion of «lost» data in method II is set back by the negative effect of an increase in error due to a greater number of terms in the final expression for pulsations intensity.

Also noteworthy is the question of the sign and magnitude of the ordered vertical drift detected during the construction of the PVD (panels a and b in Figure 2). The variability of such a drift over depths requires special consideration. As noted above, such a drift feature can be the basis for calculating the areas of the upand downwelling zones at different depths of the CML. Such estimates, in turn, are of great importance in the study of plankton dynamics. The ratio of the areas of these zones is an important parameter in the study of the depth distribution of plankton in general [Kelley, 1997] and its individual fractions (in size and mass) in particular.

The implementation of such estimations, however, is a rather difficult task. A direct indicator of the difference between upand downwelling zones is the skewness parameter associated with the third standardized moment of pulsations probability distribution. Experimental estimation of this parameter is fraught with great difficulties, and results of Large Eddy Simulation of this parameter are presented only in a few studies [Mironov et al., 2001]. To obtain a qualitative estimate, another approach is also possible, which, however, requires some refinement of the concepts of the CML structure and the concept of «measured average speed» of upand downwelling currents. The point is that the sign of the general vertical drift 𝑠 (here, we regard the direction «up» as positive) at some depth generally does not coincide with the sign of the average vertical velocity. In particular, according to the type of calculated PVD in vertical planes (Figure 2a and Figure 2b), it is still impossible to conclude that upwelling prevails in the upper layers (approximately to a depth of 1.5 m) or that downwelling prevails in the lower layers.

So, in the case of fixed deterministic cells at a given observation point, the sign of the vertical drift should be the same for all depths. In another limiting case — a chaotic structure of the convective layer — such a drift is absent. Both these extreme cases are thus ruled out according to the results of PVD calculations, and this is consistent with the

above-formulated concept of the CML as a system of quasi-deterministic cells, which in our case move relative to the stationary devices under the action of the geostrophic circulation.

In the absence of such circulation, vertical velocity measurements are actually made at a random set of points, and the number of points 𝑁+ and

𝑁− at which the vertical velocity is positive (𝑤+) or negative (𝑤−) is proportional to the area (𝐹+ and 𝐹−) of the corresponding zones in each cell. Having designated the sizes of the upand downwelling zones for (𝑙+) and (𝑙−), we can express the total vertical drift 𝑠 as: 𝑠 = (𝑤+𝑁+) − (𝑤𝑁−)


𝑤+𝑙2 − 𝑤𝑙2. As a result, taking into account the material balance 𝑤+𝐹+ = 𝑤𝐹− and the estimate

𝐹 𝑙2, we again obtain 𝑠 = 0.

Finally, we consider the situation where the cell system moves relative to the ADCPs, entrained by the circulating flow. In this case, the sign and magnitude of the vertical drift can be linked to the sizes and of the upand downwelling zones by using the following simple qualitative estimates.

In the presence of horizontal drift, the «residence time» of the velocity sensor in the upand downwelling zones is proportional to the sizes (not areas) of these zones, so that the vertical drift 𝑠 can be represented as: 𝑠 𝑤+𝑙+ − 𝑤𝑙− = 𝑤+𝑙+(1 −

𝑤𝑙−/𝑤+𝑙+). Again, given the ratio 𝐹 𝑙2, we


𝑠 = 𝑤+𝑙+(1 − 𝑙+/𝑙−).

It immediately follows from the last formula that

𝑠 > 0 if 𝑙+ < 𝑙− and vice versa. Thus, according to the sign of the vertical drift, it is possible to determine the depth ranges at which the areas of the downwelling zones are greater or less than the areas of upwelling. In our case, as follows from the analysis of the PVD, the prevalence of upwellings is typical for the deeper layers, where 𝑠 < 0. The opposite situation is typical of the upper layers.

Considering the detected differentiation of various CML strata according to the criterion «greater area of downwelling or upwelling», even where only the upper half of the CML was scanned, it cannot be ruled out that the generally accepted vision of the CML as a quasi-homogeneous layer requires a clarification. However, certain conclusions about the vertical structuring of the CML are possible only with the results of measurements in its entire thickness.
Acknowledgments. The work was carried out under state assignment to the Northern Water Problems Institute (Karelian Research Centre of the Russian Academy of Sciences). Field studies were carried out with financial support from Russian Foundation for Basic Research (Project No. 18-05-60261).



1. Austin, J. A. (2019), Observations of radiatively driven convection in a deep lake, Limnol. Oceanogr., 64, 2152-2160, Crossref

2. Bogdanov, S., et al. (2019), Structure and dynamics of convective mixing in Lake Onego under ice covered conditions, Inland Waters, 9, No. 2, 177-192, Crossref

3. Bogdanov, S., et al. (2021), Deriving of turbulent stresses in a convectively mixed layer in a shallow lake under ice by coupling two ADCPs, Fundamentalnaya i Prikladnaya Gidrofizica, 14,No. 2, 17-28, Crossref

4. Bouffard, D., A. W ̈uest (2018), Convection inlakes, Annu. Rev. Fluid Mech., 51, 189-215, Crossref

5. Bouffard, D., et al. (2019), Under ice convection dynamics in a boreal lake, Inland Waters, 9, No. 2, 142-161, Crossref

6. Bouillaut, V., S. Lepot, et al. (2019), Transition to the ultimate regime in a radiatively driven convection experiment, Journal of Fluid Mechanics, 861, R5, Crossref

7. Farmer, D. M. (1975), Penetrative convection in the absence of mean shear, Q. J. R. Meteorol.Soc., 101, 869-891, Crossref

8. Forrest, A., B. E. Laval, et al. (2008), Convectively driven transport in temperate lakes, Limnology and Oceanography, 53, 2321-2332, Crossref

9. Howarth, M. J., A. J. Souza (2005), Reynolds stress observations in continental shelf seas, Deep Sea Res. II, 52, 1075-1086, Crossref

10. Jewson, D., N. Granin, et al. (2009), Effect of snow depth on under ice irradiance and growth of Aulacoseira baicalensis in Lake Baikal, Aquatic Ecology, 43, 673-679, Crossref

11. Jonas, T., A. Terzhevik, et al. (2003), Radiatively driven convection in an ice covered lake investigated by using temperature microstructure technique, J. Geophys. Res., 108, No. C6,3183, Crossref

12. Kelley, D. (1997), Convection in ice covered lakes: effects on algal suspension, Journal of Plankton Research, 19, 1859-1880, Crossref

13. Lepot, S., S. Auma^ıtre, B. Gallet (2018), Radiative heating achieves the ultimate regime of thermal convection, Proc. Nat. Acad. Sci. USA,115, 8937-8941, Crossref

14. Maffioli, A., G. Brethouwer, E. Lindborg (2016),Mixing efficiency in stratified turbulence, Journal of Fluid Mechanics, 794, R3, Crossref

15. Mironov, D., S. Danilov, D. Olbers (2001), Large eddy simulation of radiatively driven convection in ice covered lakes, X. Casamitjana (ed.),Proc. Sixth Workshop on Physical Processes in Natural Waters, 27-29 June, 2001, Girona, Catalonia, Spain p. 71-75, Univ. of Girona, Spain.

16. Mironov, D., A. Terzhevik, et al. (2002),Radiatively driven convection in ice covered lakes: observations, scaling and a mixed layer model, J. Geophys. Res., 107, C4, Crossref

17. McNeill, A. R. (1990), Size, Speed and Buoyancy Adaptations in Aquatic Animals, American Zoologist, 30, No. 1, 189-196, Crossref

18. Nystrom, E. A., C. R. Rehmann, K. A. Oberg(2007), Evaluation of mean velocity and turbulence measurements with ADCPs, J. Hydraul. Eng., 133, 1310-1318, Crossref

19. Pernica, P., R. L. North, H. M. Baulch (2017),In the cold light of day: the potential importance of underice convective mixed layers to primary producers, Inland Waters, 7, No. 2, 138-150,Crossref

20. Palshin, N., et al. (2017), Geostrophic currents in the small ice covered lake, Advances in Current Natural Sciences, 11, 89-94. (URL, In Russian)

21. Palshin, N., et al.(2019), Effect of Under Ice Light Intensity and Convective Mixing on Chlorophylla Distribution in a Small Mesotrophic Lake, Water Resources, 46, 384-394, Crossref

22. Stacey, M. T., S. G. Monismith, J. R. Burau(1999), Measurements of Reynolds stress profiles in unstratified tidal flow, J. Geophys. Res.,104, 10,933-10,949, Crossref

23. Ulloa, H. N., A. W ̈uest, D. Bouffard (2018), Mechanical energy budget and mixing efficiency for a radiatively heated ice covered waterbody, Journal of Fluid Mechanics, 852, R1, Crossref

24. Vermeulen, B., A. J. F. Hoitink, M. G. Sassi(2011), Coupled ADCPs can yield complete Reynolds stress tensor profiles in geophysical surface flows, Geophys. Res. Lett., 38, L06406, Crossref

25. Volkov, S., et al. (2019), Large scale structure of convectively mixed layer in a shallow ice covered lake, Fundamentalnaya i Prikladnaya Gidrofizica, 12, No. 1, 30-39, (In Russian) Crossref

26. Williams, E., J. H. Simpson (2004), Uncertainties in Estimates of Reynolds Stress and TKE Production Rate Using the ADCP Variance Method, J. Atmos. Oceanic Technol., 21,347-357, Crossref

27. Woodcock, A. H. (1965), Melt patterns in ice over shallow waters, Limnology and Oceanography, 10, R290-R297, Crossref

28. Woodward, J. R., J. W. Pitchford, M. A. Bees(2019), Physical flow effects can dictate plankton population dynamics, J. R. Soc. Interface,20190247, Crossref

29. Yang, B., J. Young, et al. (2017), High Frequency Observations of Temperature and Dissolved Oxygen Reveal Under Ice Convection in a Large Lake, Geophys. Res. Lett., 44, No. 24,12,218-12,226,

Login or Create
* Forgot password?