Wave and turbulence separation using dynamic mode decomposition

Wave-turbulence decomposition is a longstanding challenge that has been approached in a variety of ways. Simpler methods like moving average and band-pass filters are often used (Foster, 1997; Smyth et al., 2002; Williams et al., 2003; Zhu et al., 2016) due to their ease of use and time-resolved output; however, these filters do not distinguish the waves from the turbulence in the overlapping frequency range, resulting in the elimination of turbulent energy.

Many spectral decomposition techniques have been developed that can estimate the Reynolds stress without allowing for a time-series reconstruction; these methods often rely on the correlation between two separate measurements. For example, Benilov et al. (1974) proposed a linear filtration technique that leverages coherence between the pressure and velocity spectra. Similarly, Trowbridge (1998) and Shaw and Trowbridge (2001) proposed a method that leverages the coherence between two synchronized velocity measurements with a finite spacing. With only one velocity measurement, Bricker and Monismith (2007) proposed a method that linearly filters the energy spectrum by assuming an inertial subrange with a $-5/3$ slope. Whereas Gerbi et al. (2008) proposed fitting the measured spectrum to an empirical spectrum model to separate the waves and turbulence. While these methods are also frequently used due to their relative simplicity, they assume no wave-turbulence interactions and do not directly allow for a reconstructed time series (unless further approximations are made, e.g. see Cowherd et al. (2021)). Additionally, some of the methods require multiple synchronized measurements which further restricts their use.

Some researchers sought a more rigorous characterization of the interaction of wave and turbulence by using the nonlinear wave stream function (Dean, 1965), which imposed fewer assumptions and allowed for a time reconstruction of the filtered signal. Jiang et al. (1990) developed a method based on this work that consists of solving, in a least squares sense, the stream functions at the free surface. Similarly, the triple decomposition method from Thais and Magnaudet (1995, 1996) also used a stream function to decompose the instantaneous velocity into turbulent fluctuation, irrotational wave, and rotational wave components. Although these methods allow the construction of separate wave and turbulent velocity time series, they are most well-suited to laboratory conditions where dispersive effects can be neglected (Thais and Magnaudet, 1995).

More recently, data-driven dimensionality-reduction techniques have been proposed to separate waves and turbulence. Ensemble Empirical Mode Decomposition (EEMD) is the foundation of the method developed by Huang and Wu (2008). It is based on the empirical mode decomposition (EMD), which decomposes the signal into intrinsic mode functions (IMFs), which define a basis dictated by the data (Huang et al., 1998). While successful at separating a time series into different components, the method suffers because a single IMF can contain several oscillatory modes or a single mode can be split into several IMFs. Huang and Wu (2008) proposed the EEMD to overcome this issue by introducing noise in the original signal in order to properly organize the different scales in the time series. The limitations of this improvement are that a large number of iterations are needed to remove the effects of the noise and obtain an accurate decomposition. In addition, ensemble-averaging may affect the physical significance of the IMFs by artificially redistributing some energy among the IMFs. More recently, Bian et al. (2018) proposed the synchrosqueezed wavelet transform (SWT)-based method (Daubechies et al., 2011; Thakur et al., 2013), which decomposes the signal into a number of components separated in the time-frequency plane. This method forces the data to conform to a basis set determined a priori and tends to underestimate the turbulence in the wave frequency range. The EEMD and SWT are modal decomposition techniques that can be applied to a single time series and allow for a time-resolved signal reconstruction. For these reasons, we use them to compare the outputs of our DMD-based method in section 5.

To summarize, traditional decomposition methods applied in spectral space are limited to scenarios with no wave-turbulence interactions, often require more specific assumptions and do not directly give a time-series reconstruction. To address these limitations, methods based on modeling waves using stream functions were developed, but they are only effective under ideal conditions that are not commonly found in nature. Finally, newer modal decomposition methods have been successful in generating a time-resolved filtered signal, but they can have a significant impact on the turbulence in the wave frequency range.

In this section, we present a summary of relevant modal decomposition algorithms; for a more in-depth reference, see Schmid (2022). Modal decomposition refers to the identification and separation of modes, or pattern, in data from a dynamical system. This type of decomposition is favorable because the general motion of a system, even nonlinear and nonstationary systems, can be approximated by a linear superposition of its modes. In the context of fluid dynamics, modal decomposition refers to the identification and separation of features in flow data (e.g., generated from Particle Image Velocimetry). The singular value decomposition (SVD) is the basis for the implementation of several dimensionality-reduction techniques such as proper orthogonal decomposition (POD) (Lumley, 1967) and dynamic mode decomposition (DMD) (Schmid, 2010).

Synthetic data allows us to test our wave-turbulence decomposition by comparing our decomposed velocities to the ground-truth signals directly. To create the dataset, we linearly add a wave $\tilde{u}$ and turbulence $u^{\prime}$ velocity time series. In this case, wave-turbulence interactions are exactly zero and $\tau/\rho={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}\langle u^{\prime}w^{\prime}\rangle$. The turbulence data come from a streamwise velocity time series collected with hot-wire anemometry in a wind tunnel sampled at 60 kHz (Castro and Vanderwel, 2021). The turbulence spectrum has a developed inertial range between 0.13 Hz to 0.6 Hz.

The wave amplitude spectrum $S(\omega)$ was generated using a Joint North Sea Wave Project (JONSWAP) model (Hasselmann et al., 1973)

where $\alpha=0.01$ is the intensity of the JONSWAP spectrum, $g$ is the acceleration of gravity, $\omega$ is the angular frequency, $\omega_{p}$ is the angular frequency of the peak of the spectrum, $\gamma=3.3$, and the factor $\beta=0.07$ for $\omega\leq\omega_{p}$ and $\beta=0.09$ for $\omega>\omega_{p}$. The peak period $T_{p}$ was chosen to be $5$ seconds, which results in $\omega_{p}=0.2$ Hz. At a point below the surface where $x=0$, and $z=-2$ m, linear wave theory was used to generate the velocity time series from the sum of $N_{w}=500$ wave components:

where $U_{i}=a_{i}\omega_{i}$ is the wave velocity amplitude vector, $a_{i}=\sqrt{2S(\omega_{i})d\omega}$ is the amplitude of the spectrum at frequency $\omega_{i}$, $k_{i}=\omega_{i}^{2}/g$ is the wave number that satisfies the deep-water linear dispersion relation, and $\phi_{i}$ is the phase which is randomly sampled from the interval $[0,\,2\pi)$.

The turbulence time series was decimated from 60 kHz to 10 Hz to match the sampling frequency of the JONSWAP waves time series. The total detrended horizontal component of velocity $u$ is obtained from adding the wave and turbulence components of horizontal velocity, defined in equation 1, as $u=u^{\prime}+\tilde{u}$. The power spectrum of $u$, shown in Fig. 2b, preserves the turbulence developed inertial range between 0.13 Hz to 0.6 Hz and shows that the wave energy is concentrated around the wave peak frequency at 0.2 Hz.

The field data consists of long waves measured in a turbulent boundary layer above the coastal shelf. The data comes from a published dataset in Reimers and Fogaren (2021). We are analyzing their acoustic Doppler velocimeter (ADV) data collected during a campaign to measure oxygen eddy covariance. The ADV captures a time series of the three perpendicular velocity components and pressure at a single point in space. We analyze only a time series of $u$ velocity components from their data. The measurements were collected from a lander deployment along the Oregon shelf above a sandy seafloor. The ADV probe’s sampling volume was 30 cm above the bottom of the sea floor in 80 m of water and was sampling at a frequency of $f_{s}=64$ Hz, which was then downsampled to $f_{s}=8$ Hz. Each deployment was broken down into 15-min (7200 observations) long bursts, and we consider only one burst in our study. The mean current of the analyzed burst time series is $0.11$ m/s, the peak wave period $T_{p}$ is $14.28$ seconds (see Fig. 8b) and for reference, the significant wave height was estimated to be 3.4 m. We chose this specific burst because it presented a mix of high and low amplitude waves, as shown in Fig. 8c. For further details about the data acquisition, please consult the original study (Reimers and Fogaren, 2021).

To contrast the bottom boundary layer field data, we collected laboratory data in a wind-driven surface boundary layer with short intermittently breaking waves. The data was collected during experiments conducted in the Washington Air-Sea Interaction Facility (WASIRF) which is a 12.2 m long, 0.91 m wide wind-wave tank. Fresh water was filled to a depth of 0.6 m leaving 0.6 m of headspace above for air circulation. To create wind, a suction fan propelled air through the test section’s headspace and recirculated it via an overhead duct. See Baker and DiBenedetto (2023) for more details on the facility.

We collected a timeseries of 3-dimensional velocity at a height of 19 cm below the water surface with an ADV (Nortek Vectrino) under a wind speed of 16 m/s at a fetch of 7.5 m. The data was collected over 20 minutes at a sampling frequency $f_{s}=30$ Hz. The wave energy spans the frequency range of 1 to 3 Hz, with a peak at $2$ Hz (Figs. 9b and 10b). We post-processed the data using the despiking method proposed by Goring and Nikora (2002). Even after despiking the data, there is still some noise present in the data due to intermittent wave-breaking which generated bubbles which affected the data quality. This noisy data provides another challenge to the decomposition. We apply our decomposition method to both the horizontal and vertical components of velocity because they present contrasting conditions. The horizontal velocity data has high turbulence fluctuations when compared to the wave velocities, whereas the vertical velocity has wave and turbulence fluctuations with similar magnitude (see Table 1). Finally, we note that due the strong anisotropy in the surface boundary layer flow, we see deviation from the -5/3 slope in the inertial range (Fig. 9b).

The first step of our decomposition method is constructing the time-delayed matrix $\bm{\mathsf{H}}$ and identifying the wave frequency range, as described in section 3 and shown in Table 1. Next, we determine the rank truncation $r$ by plotting the POD spectrum of $\bm{\mathsf{H}}$ of the synthetic $u$ data (Fig. 2a) and selecting the modes that contain wave energy. We observe two strong features in this plot: (1) the higher modes contain energy in monotonically increasing frequencies across the entire domain due to turbulence, showing the energy decreasing from low to high frequencies consistent with the inertial range and (2) a coherent high energy signal in lower modes between the frequency range of 0.13 and 0.6 Hz indicated by the vertical dotted lines. This “dip” in the POD spectrum indicates an energetic range associated with the wave energy and accounts for a significant part of the variance in the time series. We use the location of the dip to inform the rank truncation. As indicated with the horizontal dashed line, we truncate at the $r=173$ mode; that is, we take only the first 173 modes for the next step in the DMD. We see that while most of the wave energy is contained in these first $r$ modes, there is also energy associated with low-frequency turbulence in these lowest modes. This clearly shows how POD alone fails at separating waves and turbulence since it can only isolate the high-frequency turbulence from the waves, not the lower-frequency turbulence.

After determining the rank truncation, we are able to decompose the signal using DMD. Figure 2b shows the spectra of the DMD-filtered wave and turbulence time series which we compare with the true wave and true turbulence spectra. Within the wave-frequency range, the DMD-filtered wave spectrum $\tilde{u}_{\textrm{DMD}}$ closely follows the wave peak. Outside the wave-frequency range, the $\tilde{u}_{\textrm{DMD}}$ spectrum falls sharply, capturing a similar shape to the true wave spectrum for multiple orders of magnitude before levelling off. The mismatch in the signals is clearer to see in the $u^{\prime}_{\textrm{DMD}}$ spectrum, which slightly underpredicts the turbulence in the wave-frequency range.

Figures 2(c-e) show the full $u$ time series and the decomposed $\tilde{u}$ and $u^{\prime}$ time series alongside their true counterparts. We see overall agreement between the decomposed and true time series, however some mismatch exists. While the turbulence time-series plot shows that the separated signal accurately captures the low frequency turbulence fluctuations, there is stronger disagreement in the fluctuations with frequencies in the wave-range. Given that the wave signal was generated with 500 discrete wave signals and that the DMD only resolved 104 modes, we do not expect a perfect reconstruction. In this method, the DMD is discretizing a wave frequency range that in reality is continuous, so we always expect some mismatch. Column $r/n$ in Table 1 shows the ratio of the number of modes used to truncate the input matrix normalized by the number of columns used to construct the input matrix. For the synthetic data, the low rank-to-column ratio indicates that the dominant wave patterns in the signal were reconstructed using fewer modes compared to the other datasets.

We plot the POD spectrum of the field data (Fig. 8a) and find that, similar to the synthetic data in Fig. 2a, there are two features: one associated with the turbulence energy that is seen across the frequency range and one that marks a dip in the energetic modes associated with the the wave energy concentrated between 0.03 and 0.2 Hz, demarcated with the vertical dotted lines. It is worth noting that in this case, there is not a well-defined turbulence signal in the wave frequency range. Wave and turbulence interactions present in real data, as compared with the synthetic data, may reduce the effectiveness of the POD modal separation in the overlapping frequency range. This further supports why POD alone cannot always isolate waves and turbulence. Based off the location of the dip, we rank-truncate to $r=48$ modes to preserve the energy in the wave frequencies and remove the high frequency turbulence. The final filtered velocity spectra are shown in Fig. 2b alongside the raw signal spectrum. We see that the DMD method is able to isolate the wave energy and flatten the turbulence spectrum under the wave peak. While the wave spectrum contains nonzero energy in the frequencies outside the wave range, the energy is small enough to not affect the turbulence spectrum.

In Fig. 2c, we present a portion of the separated time series over the raw signal. This portion of data was chosen to highlight a potentially difficult section to decompose where the wave amplitude sharply decreases over time. We expect the decomposition to perform best when the wave energy is large relative to that of the turbulence. While we see strong agreement with the wave signal during the first half of the plot, we see the mismatch increases once the wave amplitude decreases in the second half of the plot, however the decomposed wave velocity is still able to capture the transition moderately well and maintain the shape of the waves. When we later compare this decomposition to other methods, we find it has similar performance to the EEMD and SWT methods because this is the dataset with the largest relative wave energy, i.e., it is theoretically the easiest signal to decompose. Additionally, Table 1 indicates that this is the dataset with the highest relative wave intensity $\tilde{u}_{\text{rms}}/u^{\prime}_{\text{rms}}$ and also the lowest $r/n$ values. This indicates that the wave motion present in the signal can be captured using a smaller number of modes relative to the total rank.

The laboratory data represents a case with weak wave energy relative to the turbulence, and therefore we expect it to be the most challenging to decompose. The data is also more noisy than the field and synthetic data, but it is also the longest dataset of the three. We start by considering the POD spectra of both $u$ and $w$ (Figs. 9a and 10a) which show that the wave energy is contained within the frequency range of $1$ to $3$ Hz. What is interesting in these POD spectra, unlike in those of the synthetic and field data, is that the low number modes only contain low frequency turbulence, and that the wave energy is contained in intermediate modes. Based off the POD spectra, we truncate at $r=1000$ and $r=1300$ modes for $u$ and $w$ respectively.

The fully decomposed spectra are plotted in Figs. 9b and 10b. We see that the DMD method separated and removed most of the wave energy from the turbulence in both cases. First considering the $u$ spectrum (Fig. 9b), we see that the majority of the wave peak is removed from the turbulence spectrum, and that the wave spectrum has less energy than the turbulence spectrum at almost all frequencies, confirming that this is a turbulence-dominated flow. The time-series reconstruction in Fig. 9c again shows the small wave velocity amplitude compared to the turbulence fluctuations. We also see some noise spikes in the data which have been isolated to the turbulence time series.

Next, we consider the $w$ decomposition in Fig. 10b where we observe a cleanly separated turbulence and wave spectra. Note the difference between this and the $u$ spectra, specifically how they have similar wave energy but the $w$ spectra has much weaker turbulence. The time-series reconstruction in Fig. 10c show that the separated waves closely follow the trend of the raw data, and that both the high and low frequency turbulence fluctuations are isolated to the turbulence time series. In Table 1, we see that this is the dataset with the lowest wave intensity relative to the turbulence. This low wave intensity is a challenge for our decomposition, however by using a long enough time-series and a high $r/n$ truncation, the DMD method is able to give a better separation when compared with other methods.

In order to further assess the decomposition results, we compare the mean Reynolds stress $\langle u^{\prime}w^{\prime}\rangle$ from the raw signal with those from the decomposed turbulence and wave components. We expect the Reynolds stress associated with waves to be significantly lower than that generated by turbulence. Consistent with this expectation, our analysis shows that the wave Reynolds stress is over an order of magnitude smaller than the residual turbulence Reynolds stress. We find the Reynolds stress values to be $-5.1\times 10^{-4}$, $-4.5\times 10^{-4}$, and $-0.19\times 10^{-4}$ (m/s)² for the raw, decomposed turbulence, and decomposed wave signals respectively. This confirms that the wave momentum flux is inefficient compared to that of the turbulence, and that the DMD decomposition can successfully isolate wavelike motions from turbulence.

Additionally, we visualize the time series of $\tilde{u}_{\textrm{DMD}}$ and $\tilde{w}_{\textrm{DMD}}$ together in Fig. 11. We expect the horizontal and vertical wave velocities to have a 90° phase shift in time and for them to have similar amplitudes given that they were deep-water waves (short waves relative to the depth of the tank). The initial section of Fig. 11 shows a disagreement between both the expected wave magnitudes and the phase shift, which might be caused by the presence of strong turbulence; however, the rest of the time series shows good agreement with the expected wave characteristics from theory. Throughout the time series, the $\tilde{u}_{\textrm{DMD}}$ velocity amplitude is slightly underpredicted when compared with the $\tilde{w}_{\textrm{DMD}}$ signal; this is likely due to the $u$ decomposition underperforming due to the stronger turbulence. Overall, we find that the DMD method is able to decompose the signals fairly well, even given a noisy dataset with relatively weak, irregular waves under intermittent wave-breaking.

We compare our decomposition method to the EEMD and SWT, two other modal decomposition techniques that work with a single time series. To compare the methods, we plot the separated turbulence spectra using each method for each dataset in Fig. 12. Across all four plots, we see that the SWT tends to excessively remove energy from the turbulence in the wave range, causing significant dips in the turbulence spectra. In contrast, the EEMD generally performs as well as the DMD method for the synthetic and field data in Figs. 12a and 12b, respectively. We note that the EEMD has caused some turbulence energy redistribution, e.g., adding energy into the turbulence in the field data as seen in Fig. 12b between 0.2 to 0.4 Hz. Finally, the laboratory data (Figs. 12c and 12d) is the hardest to decompose due to its weak, irregular waves, and that is where we clearly see that the DMD outperforms the other two methods. The DMD filtered spectra is able to closely follow the expected slope of the turbulence spectra for both the $u$ and $w$ velocity spectra.

The methods show the most similar decomposition for the field data, which is the dataset with the highest wave intensity. This is unsurprising and suggests that the stronger the wave intensity, the easier the signal should be to decompose using any method. Even though the synthetic data has exactly no wave-turbulence interactions which should in theory mean it is relatively easy to decompose, its lower wave intensity causes issues for the SWT method. Finally, the laboratory data has the lowest wave intensity, and this is where the benefits of the DMD method are the clearest. One reason the SWT may underperform is because it determines its basis a priori; the method tries to fit the data to a wavelet basis, which might not always be a good descriptor of the dynamics of the system. Alternatively, the EEMD uses IMFs as an a posteriori basis, that is, they adapt to the data being analyzed. However, it seems that the IMFs might not be separating the waves from the turbulence completely, resulting in mixed modes which contain both turbulence and wave energy. The utility of this method is further hindered by the fact that the number of IMFs is determined by the EEMD algorithm which removes a tuning parameter.

We further compare the methods by computing the normalized autocorrelations $R(\tau)$ of the filtered turbulence velocity across the different datasets, as shown in Fig. 13. The raw autocorrelations, plotted in gray for reference, show a clear periodic wave signal across all datasets. Successfully isolating the turbulence from the waves would remove all time periodicity and preserve the exponential decay expected as the turbulence fluctuations decorrelate over time. We see that in most cases, the methods are all able to remove the majority of the wave signal. With respect to the synthetic data (Fig. 13a), we know what the true turbulence signal is, and we see that while all the methods follow the true turbulence autocorrelation (shown with a blue line) closely, the SWT performs the worst, as shown by the small peak in the signal after the initial decrease. In Fig. 13b, the field data autocorrelations all remove the wave, and it is unclear which performs best here from visual inspection. Recall that the EEMD method had added energy into the turbulence spectra, so it’s likely that its higher autocorrelation values are inaccurate. Finally, in the laboratory data (Figs. 13c and 13d), we clearly see the DMD method performs the best. It is able to remove all wave peaks while still capturing the long-time temporal decorrelation of the turbulence. This analysis clearly shows the strength of the DMD-method across a variety of datasets and demonstrates how it can be used for more than just spectral analysis, particularly under conditions that pose challenges to other decomposition methods.