On application of empirical mode decomposition for turbulence analysis in open-channel flows

Large-scale coherent structures are key elements of open-channel flow turbulence, quantification of which remains elusive. In this work, we use empirical mode decomposition (EMD) to break down a velocity time series into different modes, denoted as “intrinsic mode functions” (IMFs). Analysis of velocity auto- and co-spectra indicates that large-scale (LSMs) and very large-scale (VLSMs) fluid motions are sufficiently represented by particular groups of IMFs. A correlation between LSMs and VLSMs, identified by the EMD analysis, was found to generate 7% of the Reynolds shear stresses. However, the EMD analysis of surrogate velocity signals with randomized spectral phases demonstrated that the revealed correlation is actually an artefact of the EMD approach and should not be interpreted physically.

The contributions of LSMs and VLSMs to turbulence quantities (e.g.turbulent kinetic energy and momentum fluxes) can be estimated using different methods.The most common and simplest approach involves an intuitively defined separation wavelength (or wavenumber) that divides LSMs and VLSMs in velocity spectra (e.g.Duan et al., 2020;Guala et al., 2006;Hutchins & Marusic, 2007), thereby allowing the detection of LSM and VLSM contributions to velocity signals and turbulence statistics.However, despite the attractiveness of this approach, a separation wavelength does not represent an actual border between turbulent structures, which are likely to overlap  in the spectral space over a range of wavenumbers.In an attempt to overcome this limitation, empirical mode decomposition (EMD, Huang et al., 1998) can be used to separate velocity signals into modes, known as "intrinsic mode functions" (IMFs).Indeed, unlike the sinusoidal Fourier components, each of these modes covers a range of wavenumbers, thereby allowing a coexistence of EMD-based modes in the spectral domain.Another velocity decomposition method is the proper orthogonal decomposition (POD, e.g.Holmes et al., 2012;Zampiron, Cameron, et al., 2020), which divides a velocity field into fully orthogonal (i.e.uncorrelated) modes.However, since the POD modes and the Fourier modes coincide for a 1D stationary velocity signal (which is a typical case of one-point measurements such as those used in this study), POD is not explored further.
EMD is gradually becoming a common tool for studying turbulent flows (e.g.Agostini & Leschziner, 2014, 2019;Franca & Lemmin, 2015;Peruzzi et al., 2021), and it is therefore beneficial to assess its advantages and limitations, particularly in relation to LSMs and VLSMs.The analysis presented in this work undertakes such an assessment using long-duration (8 h) acoustic-Doppler-velocimetry (ADV) records in rough-bed OCF.First, we will briefly introduce the experimental methodology (Section 2) and the EMD algorithm (Section 3).In Section 4, the results are outlined, including: examples of EMD modes (Section 4.1); the assessment of the contributions of LSMs and VLSMs to velocity auto-and co-spectra, and to turbulent stresses (Sections 4.2 and 4.3); and the analysis of surrogate velocity signals generated through phase randomization in the Fourier space, to evaluate the potential interrelations between LSMs and VLSMs identified by EMD (Section 4.4).Finally, the main outcomes of this work are summarized and discussed in Section 5.

Velocity measurements
ADV (Vectrino, Nortek) measurements were performed in the RS flume of the University of Aberdeen (0.4 m wide and 10.75 m long, e.g.Zampiron, Cameron, et al., 2020).The three velocity components were measured continuously for 8 h with a sampling frequency of 100 Hz at a relative elevation of z/H = 0.3 from the mean channel bed (z is wall-normal coordinate and H = 80.2 mm is mean flow depth).The measurements were located at 125H from the channel entrance, ensuring a good degree of flow development and the emergence of LSMs and VLSMs (Zampiron et al., 2023).The flow was steady and uniform for the total duration of the measurements.The bed of the channel was entirely covered by micro hooks (Pressogrip ® ) with a height of ≈ 1 mm, which ensured fully rough-bed flow conditions (Zampiron, Nikora, et al., 2020).The aspect ratio (i.e.width-to-depth) of 5 suggests quasi-two-dimensional flow conditions at the centre of the channel where measurements were made (e.g.Nezu & Nakagawa, 1993).Hydraulic parameters are listed in Table 1.The 8-hour velocity time series are provided in the online Supplemental Material.Further details on the measurements can be found in Zampiron et al. (2023).

Empirical mode decomposition (EMD) algorithm
EMD is used to decompose a generic time series θ (t) (which can be stationary or not, t is time) into a set of IMFs.The IMFs exhibit a dominant frequency that gradually decreases as the IMF order (m) increases from m = 1 to M, where M represents the total number of IMFs.An IMF must satisfy two conditions (Huang et al., 1998): (i) the count of extrema and zero-crossings must either be equal or differ by at most one; and (ii) the local mean values of the upper and lower envelopes must be zero or, in other words, the two envelopes must be symmetric about zero.It is important to note that IMFs may not necessarily be orthogonal and thus can exhibit some degree of correlation.Our implementation of EMD stems from Huang et al. (1998) and follows a sifting process to extract iteratively the mth IMF of a signal θ (IMF m θ ).The IMF order is initially set to m = 1.At the first iteration (n = 1, where n is iteration number), a provisional signal h is defined for the mth IMF as: where α is a summation index.A new provisional signal h m n+1 is computed for the next iteration n + 1 as: where e m n is the mean of the envelopes of h m n .Upper and lower envelopes are computed by fitting a cubic spline through the maxima and minima of h m n , respectively.However, this interpolation scheme is inefficient at the start and the end of the times series, resulting in unrealistic values that, through multiple iterations, can propagate throughout the signal.This issue was mitigated by the repetition of the extrema at the start and end of the signal (Huang et al., 1998).
Equation ( 2) is repeated iteratively, with n = n + 1 at each new iteration, until the following condition is satisfied: where ε is a threshold value assumed ε = 0.1 (Huang et al., 1998).Once the algorithm reaches convergence, the mth IMF of the signal θ is obtained as IMF m θ = h m n+1 and we may move to the (m + 1)th IMF.The process is repeated from Eq. ( 1) and continues until no further extrema can be detected in h m+1 n=1 , thus setting the total number of IMFs to M = m.The remaining monotonic signal h m+1 n=1 represents the residual r θ .In the case of stationary data, r θ is equal to the time-averaged value θ of the signal θ (overbar indicates time averaging).
Once all M IMFs have been obtained, EMD-filtered signals θ can be computed as: where m 1 and m 2 are chosen summation limits.The original signal θ is reconstructed as: A comprehensive review of the EMD approach and application examples are given in Huang and Shen (2014).

Results
In this work, we employ EMD to quantify the contributions of LSMs and VLSMs to velocity spectra and to turbulent stresses u u , w w and u w , where u and w are turbulent fluctuations of the streamwise (u) and vertical (w) velocity components, respectively, i.e. u = u − ū and w = w − w.Then, the obtained estimates are compared with those computed using the LSM-VLSM separation wavelength (Section 1) applied to the conventional velocity spectra.Considering that the length of both LSMs and VLSMs exceeds the flow depth H, we focus on streamwise wavelengths larger than H.The components in the streamwise (x), spanwise (y) and wall-normal (z) directions of the velocity vector are denoted as u, v and w, respectively.Note that the number of IMFs is not universal and may depend on the specifics of the EMD algorithm (e.g.convergence criteria).In Fig. 1, we show 20 s long samples for a subset of IMFs of the streamwise velocity u (IMF m u , with m = 4 to 6) together with their envelope curves.IMFs belonging to the other velocity components and/or of different order are qualitatively similar.As expected, the dominant frequency of the IMFs decreases with increase in m, while upper and lower envelopes exhibit a good level of symmetry (conditions of Huang et al., 1998).

LSM and VLSM contributions to velocity auto-and co-spectra
In this section, EMD is used to assess the contributions of LSMs and VLSMs to the pre-multiplied auto-spectra of the streamwise k x S uu (k x ) and vertical k x S ww (k x ) velocity components, and to the pre-multiplied co-spectrum k x C uw (k x ) of u and w, where k x = 2π / l x and l x are streamwise wavenumber and wavelength, respectively.In our estimates we assume the validity of Taylor's "frozen turbulence" hypothesis for OCFs, that gives l x = ū/f where f is frequency (e.g.Nikora & Goring, 2000).Spectra were computed for 2800 subsections 20.48 s long with a 50% overlap, which were then ensemble averaged.Pre-multiplied velocity spectra are particularly convenient for studying large-scale flow structures as they show energy distributions of the velocity fluctuations as a function of the spectral wavelength (unlike conventional power spectra that show energy density distributions).Auto-spectrum S vv exhibits low energy within the VLSM wavelength range and therefore is not considered in this paper.In Fig. 2, we show k x S uu (k x ) and k x S ww (k x ) of the original signals (dotted lines) together with the contributions from the IMFs.The spectrum of the original streamwise velocity signal (Fig. 2a) reveals two distinct "hills" at wavelengths of ≈ 3.5H and ≈ 17H, reflecting the presence of LSMs and VLSMs, respectively (e.g.Cameron et al., 2017;Kim & Adrian, 1999;Zampiron, Cameron, et al., 2020).As expected from Fig. 1, the dominant wavelength (i.e.corresponding to the spectral "hill") in the IMF spectra increases with m.The low-frequency IMFs 9 ≤ m ≤ 13 contain a small amount of energy at wavelengths exceeding the VLSM range (where our spectra are not resolved due to the limited size of the time subsections) and therefore are not visible.
In Fig. 3, the original auto-and co-spectra are compared to the contributions attributed to LSMs and VLSMs estimated using EMD.Considering that a zero-mean (i.e.fluctuating) velocity time series can be expressed as the sum of its IMFs (Eq.5), the auto-spectra of the streamwise (S uu ) and vertical (S ww ) velocity components can be expressed as (e.g.Bendat & Piersol, 2010, page 126): Similarly, for the velocity co-spectrum (C uw ) we have: The parts of the j th velocity component (u j ) representing LSMs ( u j LSM ) and VLSMs ( u j VLSM ) are reconstructed as the sums of different IMFs using Eq. ( 4).The remaining energy of S uu , S ww and C uw , not associated with LSMs, VLSMs, or their combination, is represented by the terms Φ uu , Φ ww and Φ uw , respectively (shown as grey areas in Fig. 3).
To define the LSM and VLSM signals, appropriate summation limits (m 1 and m 2 , defined in Eq. 4) are chosen based on seeking a good visual agreement between the shape and peaks in the original spectra, and the LSM and VLSM spectra computed from the constructed signals u j LSM and u j VLSM (Fig. 3), i.e. u LSM = 5 m=4 IMF m u , w LSM = 5 m=4 IMF m w , u VLSM = M u m=6 IMF m u and w VLSM = M w m=6 IMF m w .This approach is considered adequate due to the wide range of wavenumbers encompassed by each IMF (Fig. 2), which minimizes the subjectivity in selecting the IMFs.A similar approach was employed by Agostini andLeschziner (2014, 2019) when studying turbulence in closed-channel flows using two-dimensional EMD.
Similar to the streamwise auto-spectrum, velocity cospectrum also exhibits a strong signature of VLSMs, indicating their sizeable contribution to turbulent momentum fluxes (or shear stresses).Figure 3 suggests that a noticeable contribution to S uu and C uw arises from a combined effect of LSMs and VLSMs (i.e.LSM-VLSM correlation, Eqs 6 and 8).The LSM-VLSM contribution to S uu changes in sign with l x , while the contribution to -C uw remains positive throughout.Considering the potential physical significance of this combined contribution, its origin will be explored in Section 4.4.

LSM and VLSM contributions to turbulent stresses
The portions of turbulent streamwise (u u ) and vertical (w w ) normal stresses, and of turbulent shear stresses (u w ) associated with LSMs and VLSMs are presented in Fig. 4. EMD estimates (Fig. 3) are compared to the separation wavelength method (introduced in Section 1).The separation wavelength is chosen at l x /H = 8, close to the intersection of the LSMs' and VLSMs' curves and the peak of the LSM-VLSM contribution to C uw (Fig. 3c).The lower bound wavelength of the LSMs is set equal to the flow depth (l x /H = 1), where the LSM curve approximately reaches zero.Contributions to turbulent stresses are computed as the integrals of the respective velocity spectra, i.e. u u Using EMD, we estimate that ∼ 24-26% of normal and shear stresses are associated with LSMs, while VLSMs account for 33% of u u , 12% of w w and 36% of u w .The inter-scale (correlation) contribution LSM-VLSM (Eqs 6-8) is essentially negligible for u u and w w (due to the change of sign, Fig. 3a  and b), whereas it accounts for 7% of the total shear stress u w .Contributions of comparable magnitude have been earlier reported by Agostini and Leschziner (2019).
Using the separation wavelength method, the contributions from LSMs to velocity variances u u and w w , and to velocity covariance u w are higher compared to the EMD estimates, primarily due to the minimum wavelength being set to l x /H = 1.The contributions of VLSMs to normal stresses closely match those obtained from EMD.However, a 42% VLSM contribution to u w exceeds the 36% from the EMD analysis.This discrepancy reflects the sharp transition between scales and the orthogonality of Fourier modes (i.e.LSM-VLSM cross-terms are expected to be zero due to the orthogonality).Cameron et al. (2020) and Duan et al. (2020) have reported contributions of VLSMs to u u of ∼ 40-50%, and to u w of ∼ 30-60%.
The above EMD results, together with those presented in Section 4.2, suggest the potential importance of the LSM-VLSM term as it may relate to interactions of LSMs and VLSMs and energy exchanges between them.The validity of such an interpretation is assessed in the next section.

The nature of LSM-VLSM correlation within the EMD approach
The nature of the correlated contribution from LSMs and VLSMs can be studied by randomizing the spectral phases of the velocity signals (Prichard & Theiler, 1994).Using phase randomization, we can generate surrogate signals (denoted by tilde) that preserve LSM and VLSM energy contributions to S uu , S ww and C uw (Eqs 6-8), while suppressing any cross-correlation between LSMs and VLSMs.Only the co-spectrum C uw will be considered, for brevity, but a qualitatively similar behaviour is also observed for the auto-spectra S uu and S ww .
The Fourier transform of a generic signal θ is Θ = A(f )e iφ(f ) , where A is the amplitude, i is the imaginary unit and φ is the phase angle.The phase of Θ can be randomized as Θ = A(f )e i[φ(f )+ϕ(f )] , where ϕ obeys a uniform random distribution within [0, 2π].Thus, the phase-randomized signal θ Figure 5 Co-spectra of (a) phase-randomized EMD signals u j LSM and u j VLSM ; and of (b) signal components u j LSM and u j VLSM obtained by applying EMD to the phase-randomized velocity signal u j (Section 4.4).Other lines, symbols and shadings are defined in Fig. 3 can be computed as the inverse Fourier transform of Θ.Following these steps, any velocity signal involved in Eqs (6-8) is phase-randomized.Using the same ϕ for different velocity components preserves the cross-correlation between them, which is otherwise destroyed (Prichard & Theiler, 1994).
In Fig. 5a, we show co-spectra of phase-randomized EMD velocity signals u j LSM and u j VLSM , with phase randomization applied after the EMD of the original velocity record.In this case, the cross-correlation between u and w is preserved, while the cross-correlation (if any) between LSM and VLSM signals is eliminated by phase randomization.This result is achieved by using a random phase angle ϕ that is the same for u and w, but different for their LSM and VLSM parts.The LSM and VLSM contributions in Fig. 5a are the same as in Fig. 3c, while the LSM-VLSM correlation term is suppressed by phase-randomization.As a result, the sum of the LSM and VLSM contributions (dashed line) deviates from the spectra of the original signal.
Then, EMD is applied a second time to a "surrogate" velocity signal u j reconstructed as the sum of all its phase-randomized constituent signals, and thus free of any correlation between LSMs and VLSMs.Using the same m 1 and m 2 as in Section 4.2, u j LSM and u j VLSM are obtained.Co-spectra of u j LSM and u j VLSM are shown in Fig. 5b.The co-spectrum computed as the sum of the different contributions (dashed line) is preserved by the EMD and therefore is the same as in Fig. 5a.However, we can observe the reappearance of the LSM-VLSM correlation term, which is balanced by a decrease in the LSM and VLSM co-spectra around l x /H ≈ 8 (Eq.8).
This result explicitly demonstrates that the LSM-VLSM correlation observed in the EMD outputs (Sections 4.2 and 4.3) is actually introduced by the EMD procedure and does not reflect physical mechanism(s).While a certain degree of correlation between IMFs is expected due to the lack of orthogonality inherent in EMD, the energy content of such correlation should be small compared to the other contributions (Huang et al., 1998).However, in our case, the LSM-VLSM correlated contribution to the uw cospectrum is quite noticeable, highlighting a certain flaw of EMD when used to assess the effects of multi-scale fluid motions on turbulent fluxes of momentum.Consequently, any attempts to separate the contributions of different scales to u w using EMD and their physical interpretation should be approached with caution.Instead, the use of alternative approaches like the wavelength separation method or other decompositions procedures such as the "dynamic mode decomposition" (Schmid, 2010) may be more appropriate.On the other hand, the contribution of the LSM-VLSM correlation to autospectra is much less significant, indicating that EMD remains a valuable tool in this case, capable of extracting LSM and VLSM contribution to velocity spectra with a realistic overlap across a range of wavenumbers.

Conclusions
We used EMD to decompose velocity time-series into different modes (IMFs), which, unlike the Fourier components, can coexist over a range of spectral wavenumbers.Appropriately selected groups of IMFs were found to sufficiently represent the LSMs' and VLSMs' contributions to velocity times series, velocity auto-spectra and turbulent normal stresses (Section 4.2).The EMD analysis suggested that up to 7% of the turbulent shear stress may be due to the interactions (correlation) of LSMs and VLSMs (Section 4.3), with potentially important implications for scale-to-scale momentum exchanges and turbulence dynamics overall.However, non-zero values of the LSM-VLSM correlation terms in Eqs ( 5) and ( 6) could also reflect the lack of orthogonality of the IMFs.By randomizing the spectral phases of the velocity signals (Section 4.4), we demonstrated that the "correlated" contributions of LSMs and VLSMs to velocity spectra (Section 4.2) and turbulent stresses (Section 4.3) are an artefact of the EMD algorithm.In light of these considerations, any outcomes of the EMD analysis related to statistics across different velocity scales (e.g. in relation to Reynolds stresses and velocity co-spectra) should be treated with caution.On the other hand, our analysis suggests that EMD is an effective tool to isolate the signatures of LSMs and VLSMs on velocity signals and their contributions to auto-spectra and turbulent normal stresses for which the LSM-VLSM correlation term is negligible.

Figure 1
Figure 1 Selected IMFs (m = 4,5,6) of the streamwise velocity IMF m u (solid lines).Dashed and dash-dotted lines represent upper and lower envelopes, respectively

Figure 2
Figure 2 Pre-multiplied auto-spectra of the mth IMF of the (a) streamwise u and (b) vertical w velocity components as functions of the normalized streamwise wavelength l x /H.Dotted lines show the spectra of the original velocity signals

Figure 3
Figure3EMD contributions (Eqs 6-8) to pre-multiplied (a) streamwise velocity auto-spectrum, (b) vertical velocity auto-spectrum, and (c) uw co-spectrum.Vertical dashed lines show the separation wavelength l x / H = 8 used in Section 4.3.Grey areas show energy that is not associated with LSMs, VLSMs or their combination (Φ uu , Φ ww and Φ uw in Eqs 6-8, respectively)

Figure 4
Figure 4 Estimated relative contribution of LSMs and VLSMs to streamwise (u u ) and vertical (w w ) normal stresses, and to shear stress (u w ) using EMD (empty bars) and a separation wavelength of l x /H = 8 (filled bars)

Table 1
Flow conditions for the ADV measurements BH ) is bulk flow velocity, Q is flowrate, B is channel width; R e = UH /ν is bulk Reynolds number, ν is fluid kinematic viscosity; F = U/ √ gH is Froude number; B/H is flow aspect ratio; H /Δ is relative submergence, Δ is roughness height; H + = u * H /ν is friction Reynolds number, superscript + denotes normalization by the viscous length ν/u * ; and Δ + = u * Δ/ν is roughness Reynolds number.