Deriving Six Components of Reynolds Stress Tensor from Single-ADCP Data

Abstract

Acoustic Doppler current profilers (ADCP) are widely used in geophysical studies for mean velocity profiling and calculation of energy dissipation rate. On the other hand, the estimation of turbulent stresses from ADCP data still remains challenging. With the four-beam version of the device, only two shear stresses are derivable; and even for the five-beam version (Janus+), the calculation of the full Reynolds stress tensor is problematic currently. The known attempts to overcome the problem are based on the “coupled ADCP” experimental setup and include some hard restrictions, not to mention the essential complexity of performing experiments. In this paper, a new method is presented which allows to derive the stresses from single-ADCP data. Its essence is that interbeam correlations are taken into account as producing the missing equations for stresses. This method is applicable only for the depth range, for which the distance between the beams is comparable to the scales, where the turbulence is locally isotropic and homogeneous. The validation of this method was carried out for convectively-mixed layer in a boreal ice-covered lake. The results of computations turned out to be physically sustainable in the sense that realizability conditions were basically fulfilled. The additional verification was carried out by comparing the results, obtained by the new method and “coupled ADCPs” one.

1. Introduction

Acoustic Doppler current profilers (ADCP) are currently viewed as one of the most powerful tools for geophysical flows studies. The product family of these devices includes a lot of versions, which differ from each other by the number of beams, transducer head design, carrier acoustic frequency, cell sizes and time measurements settings. This variety of the device parameters provides flexibility in deployment and makes it possible to adjust the measurements to the broad range of research needs in meteorological (e.g., [1]), oceanological (e.g., [2,3,4]) and limnological (e.g., [5,6]) studies.

In particular, under the requirement of flow horizontal homogeneity, ADCPs are widely used for mean velocity profiling, as was originally designed. The later instrumental development (e.g., “burst” time settings, velocity measurements extending to ‘pulse coherent mode’) makes it possible to achieve higher resolution and better accuracy of the measurements, thus triggering the use of ADCP for enhanced studies of turbulence parameters. Within this new domain of ADCP applications, meaningful results have been obtained in fine-scale studies, including the estimations of dissipation rates ε [7,8,9].

At the same time some special methods have been developed for deriving the parameters of large-scale turbulence, with the special attention to the components $〈u_i^{\prime }{u}_j^{\prime }〉$ of Reynolds stress tensor ($u_i^{\prime } —$pulsation velocity components in orthogonal frame) (e.g., [3,4,10,11,12]). The equations for these components are derived from the directly available intensities $〈b_i^{\prime 2}〉$ of the “beam” velocities pulsations $b_i^{\prime }$. However, on the whole the problem remains challenging. Firstly, the number of ADCP beams usually varies from three to four (Janus configuration) and five (Janus+), and the system of equations is not complete. Secondly, in general case by calculating $〈b_i^{\prime 2}〉$ one can obtain only the relationships between the different required components $〈u_i^{\prime }{u}_j^{\prime }〉$, but not the explicit relations for each of them. For example, with three-beam ADCP, no explicit relations are available. As for the four-beam Janus configuration, after aligning the device axis with the mean velocity, by applying the so-called “variance method” one can derive two explicit relations for off-diagonal stresses (shear stresses), but that is all.

One of the ways to overcome the problem was presented in [13,14]. In both papers the main idea was based on ADCP coupling, when the experimental setup includes two rigidly connected ADCPs. In such a special setup, the design and implementation of the experiment become more complicated. Besides, in both cases the stresses derivation was conjugated with additional restrictions. The method suggested in [13] gave acceptable results only for the case when the axis of the second device was sufficiently (>20°) tilted to the vertical. However, this requirement is not recommended by the device manufacturers; besides, the tilting makes the horizontal homogeneity requirement tougher. In the method, presented in [14], the axes of both devices are vertical, and one pair of beams have the intersection point at some depth. However, in turn, this method also possesses restrictions: it is applicable only to a small range of depths, close to the depth of the intersection point.

In this paper, an alternative method for derivation of full Reynolds stresses from single-ADCP data is presented. As compared to the coupled-ADCP method, presented in [14], this new method is not restricted to stress computations for the special depth, and so gives the opportunity for stress profiling for a range of depths. The missing information is derived not from the additional beam data, but by taking into account the interbeam correlations of the velocity. It is worthy of note that usually these correlations are neglected, by suggesting the statistical independence of velocities at different beams. Meanwhile, in geophysical flows the integral scales of turbulence often are so large that the size of energy-containing eddies at some depths occurs commensurable with the distance between the beams. This is just the case, when beam velocities are correlated, and their covariance includes some “hidden” information, necessary for closing the equation system for $〈u_i^{\prime }{u}_j^{\prime }〉$.

The new method of turbulent stress derivation was applied to the case of the convectively mixed layer (CML) which develops in lakes during under-ice inhomogeneous solar heating of the water column.

2. Method Description

General Framework

In what follows below the simplest ADCP version with three azimuthally symmetric beams is considered (Figure 1). For the standard configuration the angle α₀ between the beam and vertical is 25°. As for the angle 2α between any pair of beams, its value may be determined by the following expression, derived from pure geometrical analysis:

$$\mathit{sin}\mathsf{\alpha }=\frac{\sqrt{3}}{2}\mathit{sin}{\mathsf{\alpha }}_0$$

In the orthogonal coordinate system XYZ, which is rigidly connected to the device, axis X is chosen in such a way that it lies on the same plane as beam 1 (Figure 1). In this frame of reference, the unit vectors ${\overrightarrow{n}}_1$, which identify the beams directions, have the following direction cosines:

$$\left\{\begin{array}{l}{\overrightarrow{n}}_1=\left(\mathit{sin}{\mathsf{\alpha }}_0,°0,°\mathit{cos}{\mathsf{\alpha }}_0\right),\\ {\overrightarrow{n}}_2=\left(-\frac{1}{2}\mathit{sin}{\mathsf{\alpha }}_0,°-\frac{\sqrt{3}}{2}\mathit{sin}{\mathsf{\alpha }}_0,°\mathit{cos}{\mathsf{\alpha }}_0\right),\\ {\overrightarrow{n}}_3=\left(-\frac{1}{2}\mathit{sin}{\mathsf{\alpha }}_0,°\frac{\sqrt{3}}{2}\mathit{sin}{\mathsf{\alpha }}_0,°\mathit{cos}{\mathsf{\alpha }}_0\right).\end{array}\right.$$

The trigonometric coefficients in the above expressions depend only on the basic angle α₀. Later on the correspondent matrix is denoted by N.

Each of the beam velocities {b_i}, measured directly at points A, B, C, is connected to the orthogonal projections u₁, u₂, u₃ of the velocity $\overrightarrow{u}$ at the same points by the linear relations, e.g., $b_1\left(A\right)=\left({\overrightarrow{n}}_1\stackrel{\rightharpoonup }{u}\left(A\right)\right)$, or, equivalently (note the summation over the repeated indexes):

$$b_1\left(A\right)=N_{1j}{u}_j\left(A\right)$$

(1)

For the flows, which are homogeneous in the horizontal plane, $〈u_i\left(A\right)〉=〈u_i\left(B\right)〉=〈u_i\left(C\right)〉 \equiv 〈u_i〉$ so the mean values three equations of the type (1) take the form: $〈b_i〉=N_{ij}〈u_j〉$. As the result, the mean velocity components at the given depth are obtained directly as the convolution of the so-called transformation matrix $\widehat{T}={\widehat{A}}^{-1}$ with the “vector” $\left(〈b_1\left(A\right)〉,〈b_2\left(B\right)〉 ,〈b_3\left(C\right)〉\right)$ of mean beam velocities.

The same homogeneity assumption makes it possible to express the beam velocity pulsations intensities $〈b_i^{\prime 2}〉$ through the turbulent stresses, e.g.,:

$$〈b_1^{\prime 2}〉=〈u_i^{\prime 2}〉{\mathit{sin}}^2{\mathsf{\alpha }}_0+〈u_3^{\prime 2}〉{\mathit{cos}}^2{\mathsf{\alpha }}_0+2〈u_1^{\prime }{u}_3^{\prime }〉\mathit{sin}{\mathsf{\alpha }}_{0}cos{\mathsf{\alpha }}_0$$

(2)

Expression (2) and two similar ones (for beams 2 and 3) represent three linear equations for six target components $〈u_i^{\prime }{u}_j^{\prime }〉$, so the system is incomplete. Moreover, the explicit relation for any stress through beam pulsation intensities is not available, as was mentioned in the introduction.

In the general case, the problem of yielding the turbulent stresses from three-beam ADCP data is highly problematic. For some flows, not only the external length scales, but also integral scales of turbulence are large enough as compared to the distance between beams. For such cases the correlations between beam velocities are not vanishing, and taking the values of $〈b_i{b}_j〉$ into account gives the opportunity to overcome the problem of missing equations.

To implement this opportunity, the structural function (SF) of the general type ${\tilde{D}}_{12}$ should be introduced into consideration. For beams 1 and 2, for example, this SF is defined as:

$${\tilde{D}}_{12}=〈\left(b_1\left(A\right)-b_1\left(B\right)\right)\left(b_2\left(A\right)-b_2\left(B\right)\right)〉$$

(3)

Under the assumption of local isotropy and homogeneity, the function ${\tilde{D}}_{12}$ is presented through the longitudinal SF $D_{LL}$ (see Appendix A): ${\tilde{D}}_{12}=\left(\frac{4}{3}co{s}^2\mathsf{\alpha }-si{n}^2\mathsf{\alpha }\right)D_{LL}$. So, after opening the brackets, and taking $〈b_1\left(A\right)b_2\left(A\right)〉=〈b_1\left(B\right)b_2\left(B\right)〉$ into account (horizontal homogeneity) Expression (3) takes the form:

$$\left(\frac{4}{3}co{s}^2\mathsf{\alpha }-si{n}^2\mathsf{\alpha }\right)D_{LL}=2〈b_1{b}_2〉-<b_1\left(A\right)b_2\left(B\right)>-<b_1\left(B\right)b_2\left(A\right)>$$

(4)

Here $〈b_1{b}_2〉=<b_1\left(A\right)b_2\left(A\right)>$.

Longitudinal SF $D_{LL}$ is derived directly by calculating along-beam velocity correlations. The second term in the r.h.s. of Equation (4) is also available directly from experimental data. As for the last term in the r.h.s., after taking into account both the horizontal homogeneity and refection invariance, the following expression is obtained in [15]: $<b_1\left(A\right)b_2\left(B\right)>=<b_1\left(B\right)b_2\left(A\right)>$. So, all except $〈b_1{b}_2〉$ terms in the Equation (4) are available from the experiment, and one may regard Equation (4) as the explicit expression for $〈b_1{b}_2〉$. With regard to the presentation $〈b_1{b}_2〉=〈b_1^{\prime }{b}_2^{\prime }〉+〈b_1〉〈b_2〉$ this expression may be also written as:

$$〈b_1^{\prime }{b}_2^{\prime }〉=\left(\frac{2}{3}{\mathit{cos}}^2\mathsf{\alpha }-\frac{1}{2}{\mathit{sin}}^2\mathsf{\alpha }\right)D_{LL}+〈b_1\left(A\right)b_2\left(B\right)〉-〈b_1〉〈b_2〉$$

(5)

The similar presentations are valid for two remaining pairs of beams (13 and 23).

On the other hand, with presentation (1) in mind, three “beam stresses” $〈b_i^{\prime }{b}_j^{\prime }〉\left(i\ne j\right)$ may be presented as the linear combination $N_{il}{N}_{jm}〈u_l^{\prime }{u}_m^{\prime }〉$ of the Reynolds tensor components $〈u_l^{\prime }{u}_m^{\prime }〉$ in the same way as Expression (2) for $〈b_i^{\prime 2}〉$ were derived. For example:

$$〈b_1^{\prime }{b}_2^{\prime }〉={\mathit{sin}}^2{\mathsf{\alpha }}_0(-〈u_1^{\prime 2}〉/2-\sqrt{3}〈u_1^{\prime }{u}_2^{\prime }〉/2+〈u_1^{\prime }{u}_3^{\prime }〉\left(\mathit{cot}{\mathsf{\alpha }}_0\right)/2-\sqrt{3}〈u_2^{\prime }{u}_3^{\prime }〉\left(\mathit{cot}{\mathsf{\alpha }}_0\right)/2+〈u_3^{\prime 2}〉{\mathit{cot}}^2{\mathsf{\alpha }}_0)$$

(6)

Expression (6) and two similar ones (for $〈b_1^{\prime }{b}_3^{\prime }〉$ and $〈b_2^{\prime }{b}_3^{\prime }〉$) one may regard as three missing equations for $〈u_i^{\prime }{u}_j^{\prime }〉$. Together with Equation (2) (and two similar expressions —for $〈b_2^{\prime 2}〉$ and $〈b_3^{\prime 2}〉$) they form the closed system of linear inhomogeneous equations.

To represent this system in the compact form, it is reasonable to introduce into consideration the following “vectors”:

$$B_i=\left(〈b_1^{\prime 2}〉,〈b_2^{\prime 2}〉,〈b_3^{\prime 2}〉,〈b_1^{\prime }{b}_2^{\prime }〉,〈b_1^{\prime }{b}_3^{\prime }〉,〈b_2^{\prime }{b}_3^{\prime }〉\right),$$

$$R_i=\left(〈u_1^{\prime 2}〉, 〈u_1^{\prime }{u}_2^{\prime }〉, 〈u_1^{\prime }{u}_3^{\prime }〉, 〈u_2^{\prime 2}〉, 〈u_2^{\prime }{u}_3^{\prime }〉, 〈u_3^{\prime 2}〉\right).$$

With these notations the system of equations becomes:

$$B_i=M_{ij}{R}_j, i,j=1\dots 6$$

(7)

The coefficient matrix M is derived directly from Equations (2) and (6) and similar ones:

$$M={\mathit{sin}}^2{\mathsf{\alpha }}_0\left(\begin{array}{cccccc}1& 0& 2\mathit{cot}{\mathsf{\alpha }}_0& 0& 0& {\mathit{cot}}^2{\mathsf{\alpha }}_0\\ 1/4& \sqrt{3}/2& -\mathit{cot}{\mathsf{\alpha }}_0& 3/4& -\sqrt{3}\mathit{cot}{\mathsf{\alpha }}_0& {\mathit{cot}}^2{\mathsf{\alpha }}_0\\ 1/4& -\sqrt{3}/2& -\mathit{cot}{\mathsf{\alpha }}_0& 3/4& \sqrt{3}\mathit{cot}{\mathsf{\alpha }}_0& {\mathit{cot}}^2{\mathsf{\alpha }}_0\\ -1/2& -\sqrt{3}/2& \mathit{cot}{\mathsf{\alpha }}_0/2& 0& -\sqrt{3}\mathit{cot}{\mathsf{\alpha }}_0/2& {\mathit{cot}}^2{\mathsf{\alpha }}_0\\ -1/2& \sqrt{3}/2& \mathit{cot}{\mathsf{\alpha }}_0/2& 0& \sqrt{3}\mathit{cot}{\mathsf{\alpha }}_0/2& {\mathit{cot}}^2{\mathsf{\alpha }}_0\\ 1/4& 0& -\mathit{cot}{\mathsf{\alpha }}_0& -3/4& 0& {\mathit{cot}}^2{\mathsf{\alpha }}_0\end{array}\right).$$

Summing up, it seems reasonable to stress some key points and the step by step procedure. The corresponding algorithm looks as follows:

1.

After proper choice of time averaging interval, the mean beam velocities $〈b_i〉$, pulsation intensities $〈b_i^{\prime 2}〉$ (Equation (2)) and correlations $〈b_i^{\prime }{b}_j^{\prime }〉$ (Equation (5), i ≠ j) are calculated directly from experimental data.

2.

For each beam the function $D_{LL}$ is calculated. After revealing the inertial interval, its extent is estimated, with the special attention to its upper scale limit l.

3.

The range of depths is chosen in such a way that the distance between beams does not exceed the scale l. The maximum depth h is derived from inequality AB < l (see Figure 1): $h<l/\left(\sqrt{3}\mathit{tan}{\mathsf{\alpha }}_0\right)$.

4.

For chosen depths, the turbulent stresses are calculated directly by solving the system (7):

$${\mathrm{R}}_i=M_{ij}^{-1}{}_{ij}{\mathrm{B}}_j, i,j=1\dots 6$$

(8)

Here the inverse matrix M⁻¹ looks like (here tan(α₀) is shortly denoted as t):

$$M^{-1}={\mathit{sin}}^{-2}{\mathsf{\alpha }}_0\left(\begin{array}{cccccc}4/9& 1/9& 1/9& -4/9& -4/9& 2/9\\ 0& 1/\left(3\sqrt{3}\right)& -1/\left(3\sqrt{3}\right)& -2/\left(3\sqrt{3}\right)& 2/\left(3\sqrt{3}\right)& 0\\ 2t/9& -t/9& -t/9& t/9& t/9& -2t/9\\ 0& 1/3& 1/3& 0& 0& -2/3\\ 0& -t/\left(3\sqrt{3}\right)& t/\left(3\sqrt{3}\right)& -t/\left(3\sqrt{3}\right)& t/\left(3\sqrt{3}\right)& 0\\ t^2/9& t^2/9& t^2/9& 2t^2/9& 2t^2/9& 2t^2/9\end{array}\right).$$

3. Experimental Setup and Results

Method validation was carried out with the velocity data, obtained from the special experiment on the shallow ice-covered lake Vendyurskoe (Karelia, Russia) between 27 March and 6 April 2020. The under-ice convection (most intense during 28–31 March and 4–6 April, when solar radiation was maximal) was clearly observed, with the convectively-mixed layer’s (CML) thickness varying from 3 to 6 m. The details of experimental setup are presented in [14].

The measurements were carried out near the northern shore of the lake, the location of the experimental complex is marked by a triangle on Figure 2a; the depth at this location was ~7 m. The instrumental complex included a thermistor chain with 13 temperature sensors (RBR Ltd., accuracy ± 0.002 °C, measurement interval 10 s). The vertical temperature profile is schematically presented in Figure 2b; it clearly demonstrates the splitting of the water body into three sublayers (thin underice gradient sublayer, CML, and the underlying stratified zone), which is typical for developed convection. The CML’s lower boundary was determined by thermistor chain data as a depth of the isotherm with a value exceeding the average temperature of the CML by 0.05 °C.

The most essential feature is connected with installing two rigidly connected down-looking ADCPs (2 MHz HR Aquadopp current velocity profiler, Nortek AS, Norway). Both devices were installed on the ice (Figure 2b) with emitters located 3 cm below the lower ice boundary (Figure 2c). The X-axes of both devices were aligned with the separation vector between two emitters, but were oppositely directed (Figure 1 and Figure 2c). Due to the choice of this specific configuration one (X-axes are oriented towards each other) or two (X-axes are oriented away from each other) pairs of beams have intersection points. Figure 2c illustrates the second variant, with the beams in question being 3, 5 and 2, 6. The first variant of the devices’ settings was realized from 17:00 on 27 March to 9:30 on 30 March, and the second between 10:00, 30 March and 10:00, 6 April 2020. In both cases the depth of the intersection points was the same (~1.6 m). By fitting the distance between the emitters. The presence of intersection points, as was shown in [14], is a key feature for deriving full stress tensor from coupled-ADCP data.

For both variants of coupled-ADCP setup, the signal discreteness was one minute (32 pulses with a frequency of 2 Hz) and the depth scanning range was 2.875 m (115 cells with a size of 25 mm). To exclude the mutual influence of the two ADCPs, the emitters were set in an asynchronous mode with a 30 s delay (Figure 2c). Then the radial velocities were averaged over 16 s active series; further processing was carried out using these averages, for which the same designations b_i were used. The root-mean-square error of b_i values varied in the range (0.1–0.5) mm/s.

Data processing was carried out in accordance to algorithm presented above. The details of the averaging procedure and the choice of averaging interval (100 min.) are presented in [14]. For specificity, only the results, which correspond to the active solar radiation period of 4–6 April and a depth of 1.6 m (corresponding to the position of beams intersection point) are presented below. For the rest of the data (another dates and depths) the results are similar.

The calculated dynamics of beam velocity intensities $〈b_i^{\prime 2}〉$ and interbeam correlations $〈b_i^{\prime }{b}_j^{\prime }〉$ are presented by Figure 3. The daytime maximums of the intensities reached the values (1–2) mm²/s², whereas interbeam correlations, remaining statistically significant, varied roughly from one third to one half of their limits.

The next step included the calculations of structure functions $D_{LL}$. During daylight time the structure function curves clearly demonstrate the presence of the inertial interval for all six beams. For illustration, the daytime sequence of the calculated SF for 4 April is presented by Figure 4 (the curves were averaged by all six beams). The upper bound of the inertial interval reached values up to 1 m, which is not much less as compared to the distance between beams at the depth under consideration. This fact gives grounds (see point 3 of the algorithm above) for involving Equation (6) as a missing equation. So Equation (8) may be used for stress estimations.

The calculations of the stresses were carried out by Equation (8) for both devices separately. The results of these independent computations demonstrate qualitative agreement, as Figure 5 illustrates.

Daily maximums of pulsation intensities along axes X, Y, Z achieved the values 7, 4 and 1 mm²/s² correspondently. The value of the anisotropy coefficient $〈u_3^{\prime 2}〉/〈u^{\prime 2}〉$ was subjected to irregular oscillations within the range (0.05–0.30).

The estimations of standard deviation for stresses were carried out at the same way as presented in [14]. The errors varied from 15% for $〈u_1^{\prime 2}〉$ to 25–30% for off-diagonal stresses. It is also worthy to note another criterion of physical sustainability of the results. This criterion—the so-called realizability condition—includes the positive definiteness of pulsations intensities and the Cauchy–Schwarz inequalities ${〈u_i^{\prime }{u}_j^{\prime }〉}^2 \le 〈u_i^{\prime 2}〉〈u_j^{\prime 2}〉$(here no summation over repeated indices). The fulfillment of both restrictions is the crucial point for low-energetic flow computations. In our case, the violations of this criterion were fixed during time intervals (presented by vertical red lines on the middle image of the top panel, Figure 5), which cover less than 5% of the whole observational period. Most of these intervals belong to nighttime, when the turbulence was sufficiently suppressed.

As was mentioned in the introduction, the same experimental data were used in [14], but the yielding of stresses was carried out by the coupled-ADCPs method. The comparison of these alternative computations may serve as additional verification of the new method presented in this paper. The results of this comparison are presented in

Table 1. Comparison of stress computations by two independent methods.

Stress Component	Correlation Coefficient, r	Coefficient of Determination, R2	Linear Regression Coefficient
$〈u_1^{\prime 2}〉$	0.96	0.92	1.31
$〈u_2^{\prime 2}〉$	0.98	0.96	1.03
$〈u_3^{\prime 2}〉$	0.92	0.85	0.61
$〈u_i^{\prime 2}〉$	0.99	0.98	1.14

and Figure 6. The computations by the new method are presented as the average values of the results, obtained for each device separately.

For all three pulsation intensities (diagonal components of the stresses matrix) the correlation coefficient r is higher than 0.9. At the same time, the new single-ADCP method gives the values 1.31 for $〈u_1^{\prime 2}〉$ and 0.61 for $〈u_3^{\prime 2}〉$, as compared to the coupled-ADCPs method. The values of the component $〈u_2^{\prime 2}〉$, calculated by both methods, are practically identical (with deviations within 3%). The best fitting (Figure 7) was observed for turbulent kinetic energy $〈u_i^{\prime 2}〉$(TKE) for which coefficient of determination R² achieved the value 0.98 (Table 1).

4. Discussion

The method for yielding turbulent stresses presented in this paper is not restricted to the case of three-beam ADCP. Moreover, with other conditions being equal, four- or five-beam devices, presumably, give some preferences for method implementation. First of all, with such devices, more equations of type (2) and (6) become available. The resulting system turns out to be overdetermined, but can be solved in a least-square sense, as was done in [13] for the case of coupled ADCPs. Though such a solution is only an approximation, one may expect that the procedure does not reduce the accuracy of stress estimation, due to the increase of information involved.

The additional advantage of Janus and Janus+ versions is that the angle 2α between beams is smaller, as compared to the three-beam device. As a result, the distance between beams for given depth becomes smaller too, so the depth range, where the correlations of the beams velocity increments satisfy the local isotropy and homogeneity requirements and Equation (6) are valid, becomes wider.