# Anomalous Expansion of the Copper-Apical Oxygen Distance in Superconducting
La2CuO4 $-$ La1.55Sr0.45CuO4 Bilayers

Hua Zhou1 Yizhak Yacoby2 Ron Pindak1<EMAIL_ADDRESS>Vladimir Butko1 Gennady
Logvenov1 Ivan Bozovic1 1Brookhaven National Laboratory, Upton, NY 11973 USA
2Racah Institute of Physics, Hebrew University, Jerusalem 91904, Israel

###### Abstract

We have introduced an improved X-ray phase-retrieval method with unprecedented
speed of convergence and precision, and used it to determine with sub-Ångstrom
resolution the complete atomic structure of an ultrathin superconducting
bilayer film, composed of La1.55Sr0.45CuO4 and La2CuO4 neither of which is
superconducting by itself. The results show that phase-retrieval diffraction
techniques enable accurate measurement of structural modifications in near-
surface layers, which may be critically important for elucidation of surface-
sensitive experiments. Specifically we find that close to the sample surface
the unit cell size remains constant while the copper-apical oxygen distance
shows a dramatic increase, by as much as 0.45 Å. The apical oxygen
displacement is known to have a profound effect on the superconducting
transition temperature.

###### pacs:

68.35.Ct, 61.05.Cm, 74.78.Fk, 74.72.Dn

††preprint: prepared for Arxiv Preprint

The exciting discovery of interface superconductivity in complex oxides Reyren
et al. (2007); Bozovic et al. (2002); Gozar et al. (2008); Ueno et al. (2008);
Yuli et al. (2008) has triggered intense debate about its origin and the
possibility to enhance Tc even further Kivelson (2002); Yunoki et al. (2007);
Okamoto and Maier (2008); Berg et al. (2008). The interfacial enhancement
Gozar et al. (2008) of the superconducting critical temperature (Tc) is
influenced by the crystal structure. The Z-axis lattice constant (c0) varies
significantly among La2CuO4 (LCO), La1.55Sr0.45CuO4 (LSCO), and bilayer
LSCO/LCO films, depending even on the deposition sequence, and it affects
superconductivity: Tc scales with c0 almost perfectly linearly Butko et al.
(2008). The reason for this is not understood at present, but notice that in
(La,Sr)2CuO4 the change in c0 goes together with the change in cA, the
distance between copper and the nearest apical oxygen, which some believe to
play a key role in the high temperature superconductivity (HTS) phenomenon
Fernandes et al. (1991); Ohta et al. (1991); Feiner et al. (1992); Lubrittoy
et al. (1996); Pavarini et al. (2001); Slezaka et al. (2008). In any case, it
is certain that (i) from one cuprate to another, cA varies more than any other
bond length, and (ii) at least in simple cuprates with a single CuO2 layer in
the unit cell it correlates with the maximal Tc \- the longer cA, the higher
Tc. At least, the first fact can be understood: cA is ’soft’ because apical
oxygen has no hard contact with the nearest copper ion; rather, it ”levitates”
on the electrostatic potential - a structural feature peculiar to certain
layered oxides with alternating ionic planes of opposite charge Radovic et al.
(2008). This makes apical oxygen prone to very large displacements - e.g., in
HgBa2CuO6 one finds cA $\approx$ 2.8 Å, longer by 0.9 Å than the in-plane Cu-O
bond; coincidentally, this compound has the highest Tc = 97 K among all
single-CuO2-layer cuprates Wagner et al. (1993).

It is thus important to find out what happens to the apical oxygen in LSCO-LCO
bilayers; however, standard X-ray diffraction (XRD) is not well suited for
this $-$ one needs to ”get inside the unit cell” and look for individual
atomic displacements. For this, the most suitable technique is the Coherent
Bragg Rod Analysis (COBRA) method Yacoby et al. (2002); Sowwan et al. (2002);
Willmott et al. (2007); Yacoby et al. (2008); Cionca et al. (2008); Fong et
al. (2005). However, COBRA is most effective for films that are just a few
unit cells (UCs) thick, and in the case of HTS compounds fabrication of
ultrathin films with bulk properties has been proven to be extremely
challenging. Fortunately, we had a technical solution at hand $-$ a unique
atomic layer-by-layer molecular beam epitaxy (ALL-MBE) system with proven
capability of fabricating ultrathin HTS layers Bozovic (2001); Bozovic et al.
(2003).

For this study we have synthesized by ALL-MBE a number of (n $\times$ LSCO + m
$\times$ LCO) bilayers, where (n,m) determine the thickness of the respective
layers expressed as the number of UCs. In this paper, we show the COBRA
results for two of these, (2.5,2.5) and (2,3). The films were deposited at T =
650∘C and p = 9 $\times$ 10-6 Torr of ozone and subsequently cooled down under
high vacuum to drive out all the interstitial oxygen. We used 10 $\times$ 10
$\times$ 1 mm 3 single-crystal LaSrAlO4 (LSAO) substrates polished with the
large surface perpendicular to the (001) direction. The substrate lattice
constants are a0 = b0 = 3.755 Å, c0 = 12.56 Å; the films are pseudomorphic
with LSAO and under compressive strain. The crystal structure of LCO is
illustrated in Fig. 1a. Atomic force microscopy scans over a large (10
$\times$ 10 $\mu$m2) area showed root-mean-square surface roughness of 0.25 nm
in the (2,3) and 0.11 nm in the (2.5,2.5) bilayer sample; this is
significantly less than the 1 UC step height which in LSCO is 1.32 nm.
Magnetic susceptibility was measured via two-coil mutual inductance technique
and revealed sharp superconducting transitions at Tc = 34 K in the (2,3), and
Tc = 36 K in the (2.5,2.5) bilayer, significantly higher than the values
reported for (n,m) bilayers in Ref. [3], which is remarkable given that these
films are only 5 UC thick. This was also confirmed by measuring the electric
resistance (see Fig. 1b) after the X-ray scattering experiments were completed
and gold pads were evaporated to enable four-point-contact measurements.

Figure 1: (Color online) A simplified structure model (one-half the
crystallographic unit cell) of La2CuO4 and the transport property for a
bilayer film. (a) At room temperature, the structure is tetragonal and the
space group is I4/mmm. Noted that the La(Sr)-O layers are strongly corrugated,
exaggerated in this sketch for clarity. (b) The electric resistance of the
(2.5,2.5) bilayer, measured by the four-point-contact technique, as a function
of temperature. Inset: a schematic of the bilayer on a LSAO substrate.

The atomic structure of the LSCO/LCO bilayer film was investigated at beamline
ID-33 of the Advanced Photon Source by measuring the diffraction intensities
along the substrate-defined Bragg rods. The sample and a PILATUS 100k photon-
counting pixel detector Schlepütz et al. (2005) were mounted on a six-circle
goniometer in Kappa geometry. The experimental set-up and procedures were
described in detail in previous works Sowwan et al. (2002); Yacoby et al.
(2008). Ten symmetry inequivalent Bragg rods were recorded with a maximum
value for the vertical reciprocal space coordinate of Lmax = 10.5 r.l.u.
(reciprocal lattice units) and a sampling density of 50 points per r.l.u.. The
X-ray flux after the Si (1,1,1) monochromator crystal was 3$\times 10^{12}$
photons/sec at a wavelength of $\lambda$ = 0.8266 Å. For all Bragg rod
measurements, except for the (0,0,L) rod, the angle of incidence had a fixed
value of 3.5∘. The X-ray beam was focused to 0.1 mm (V) $\times$ 0.2 mm (H),
resulting in a 2 mm long X-ray footprint. The background and diffuse X-ray
scattering contribution were removed efficiently and accurately using the
PILATUS detector images. The final results were then normalized by taking into
account the beam polarization and Lorentz factors. The results were
subsequently analyzed using the COBRA method Yacoby et al. (2002); Willmott et
al. (2007); Yacoby et al. (2008); Cionca et al. (2008); Fong et al. (2005). In
general, COBRA uses the measured diffraction intensities and the fact that the
complex structure factors (CSFs) vary continuously along the substrate-defined
Bragg rods to determine the diffraction phases and the CSFs. The CSFs are then
Fourier transformed into real space to obtain the 3-dimensional electron
density of the film and the substrate with sub-Ångstrom resolution.

Figure 2: (Color online) Representative Bragg rods of LSCO/LCO system (open
diamond) and calculated diffraction intensity obtained from COBRA-determined
electron density (solid line).

The experimental data of three representative Bragg rods are shown in Fig.
2(a). Notice that the diffraction intensity along the Bragg rods (excluding
the Bragg peaks) vary over more than 4 orders of magnitude with excellent
signal-to-noise ratio. The reference structures chosen as the starting point
for the COBRA analysis were the bulk LSAO structure and the tetragonal
LSCO/LCO bilayer with the nominal layer atomic positions. In our numerical
simulations, the topmost 4 UCs of the substrate were allowed to deform,
however the resulting deformations turned out to be very small. The COBRA
method uses the approximation that at two adjacent points along the Bragg rod
the change in CSFs contributed by the unknown part of electron density is
negligible compared to the change in CSFs contributed from the reference
structure Sowwan et al. (2002). The use of this approximation allows COBRA to
converge very quickly to approximately the right solution but not to the exact
one. To overcome this limitation we further refined the CSFs using the
Difference-Map algorithm introduced by Elser Elser (2003) and recently applied
to thin films Björck et al. (2008). Using the COBRA solution as the starting
point for the Difference-Map algorithm and using a proper filter program that
takes advantage of the fact that the CSFs vary continuously along the Bragg
rods, the Difference-Map algorithm converges after about 20 iterations; the
convergence accelerates by about two orders of magnitude. As seen in Fig. 2,
the final calculated and measured intensities are in very good agreement.
Similar agreement was found for all other Bragg rods and the overall X-ray
reliability factor $R=\frac{\Sigma||F_{0}|-|F_{c}||}{\Sigma|F_{0}|}=0.02$;
here, F0 and Fc are the observed and the calculated diffraction amplitudes,
respectively. To the best of our knowledge, so far there has been only one
attempt to determine the structure of a thin film using the Difference Map
method Björck et al. (2008). In that study, the atomicity constraint was
imposed and over 2,000 iterations were needed to achieve convergence. Our
analysis shows that the combined COBRA/Difference Map method combines the best
features of both methods and ensures rapid convergence to the correct solution
without the need to use the atomicity constraint.

Figure 3: (Color online) The electron density variation, determined by COBRA,
along the [0, 0, Z] column of atoms as indicated by the dashed line in Figure
1a. Note that the topmost four unit cells of the substrate are included in the
structure refinement. The left and right dashed lines represent the nominal
LSAO/LSCO and LSCO/LCO interfaces, respectively.

The CSFs obtained have been Fourier transformed into real space yielding the
3D electron density (ED). As an example, we show in Fig. 3 the ED of a (2,3)
bilayer sample along the [0,0,Z] line that goes through the La(Sr), O, and
Cu(Al) atoms. Two points should be stressed. First, as seen the ED has almost
no negative parts. Together with the excellent agreement between the
calculated and measured diffraction intensities, this suggests that the ED is
very close to the correct one. Second, all the atoms including the oxygens can
be clearly identified and their positions determined with sub-Ångstrom
resolution. The small ED intensity fluctuation below -53 Å provides a measure
of the inaccuracy in the ED and as seen it is small even compared to the
oxygen ED.

The atomic positions in the Z direction were accurately determined by fitting
a Gaussian to each peak. We determined the size of the unit cell in the Z
direction by measuring the distance between consecutive pairs of La(Sr) and Cu
atoms. The results are shown in Fig. 4 inset. Each point corresponds to an
average of 4 La-La distances and 2 Cu-Cu distances. The measured lattice
constant of the bilayer film is 13.304 $\pm$ 0.016 Å and is larger by 0.148 Å
than that of the bulk LCO (c0 = 13.156 Å).

Figure 4: (Color online) Evolution of the inter-atomic distances. (a) The
measured Cu(Al)$-$ apical O distance, cA, varies as a function of the nominal
position of Cu(Al) atoms inside the refined structure. The data from two
representative bilayer samples and the average over the two are presented. The
lower and upper arrows represent the bulk values of cA for LSAO and for LCO,
respectively. (b) The comparison of cA, c1, and c2, averaged for each unit
cell, as a function of Z position. Inset: the lattice constant c0 as a
function of Z. The dotted line represents the bulk LSAO value. The horizontal
dashed line is the average value of c0 in bilayers extracted from the electron
density, as described in the text. In both (a), (b) and the inset, the
vertical dashed lines represent the nominal LSAO/LSCO interfaces,
respectively.

While the changes observed in the unit cell sizes are as expected from the
strain and the elastic parameters Ledbetter et al. (1992), the variations in
Cu-apical O and La-apical O distances are quite unexpected. The distances
cA,c1 and c2 are defined in Fig. 1a; the distances labeled
c${}_{A}^{{}^{\prime}}$,c${}_{1}^{{}^{\prime}}$ and c${}_{2}^{{}^{\prime}}$
would be their symmetry equivalents in bulk samples, but in thin films they
could differ in principle. For our samples, the measured values for the primed
and unprimed distances were in fact equal within the experimental error,
except at the LSAO/LSCO interface. The measured values averaged over cA and
c${}_{A}^{{}^{\prime}}$ are shown in Fig. 4a. The diamond and circular dots
represent the distances measured in (2,3) and (2.5,2.5) bilayers,
respectively. The triangular dots are averages over the two samples. Every
pair of triangular dots corresponds to one UC. The dashed vertical lines
represent the nominal LSAO/LSCO interfaces. The arrows on the right indicate
cA as measured in the bulk samples. The results show that, within the
experimental error, the values of cA in the substrate are equal to those in
the bulk but they are very different in the film. In both LSAO and LCO bulk
crystals Radovic et al. (2008), cA is equal to 2.41 Å. In the metal layer
closest to the substrate cA = 2.3 Å, and it then rises steadily all the way to
cA = 2.75 Å $-$ a change of 0.45 Å.

In Fig. 4b we display cA as well as the La-apical O distance, c2, and the La-
CuO2 plane distance, c1. Each point represents an average over the two bulk-
symmetry-equivalent distances and over the two measured samples. Notice that
c1 changes by less than 0.1 Å. On the other hand, cA increases by about 0.45
Å, while c2 decreases by about 0.25 Å. Close to the film surface, the apical
oxygen atoms are displaced away from the nearest Cu atoms. The La atoms are
displaced towards the closest CuO2 plane, but by a smaller amount, while the
separation between two adjacent CuO2 planes remains constant.

According to Ref. 15, cA = 2.7 Å should correspond to a Tc of 80 K at the
optimum doping. However, from Ref. 31 we know that the hole density drops
sharply on the I side of the interface and the screening length is equal to 6
$\pm$ 2 Å. This implies that on the I side and next to the M-I interface only
one or two CuO2 layers are doped via carrier accumulation while the others
remain insulating. Thus, unfortunately, we have a mismatch: in the optimally-
doped LCO layer, cA is close to its standard (bulk) value, while it is greatly
elongated only in insulating layers. It is tempting to speculate that one
could create LSCO-based samples with Tc much higher than 36 K, perhaps as high
as 80-90 K, if only one could achieve cA elongation and optimal doping in the
same LCO layer. An obvious avenue for further research is to try making I
layers even thinner, thus bringing the interface superconductivity closer to
the film surface. Another is to try engineering more sophisticated hetero-
structures and superlattices combining LCO with other metallic oxides
(nickelates, zincates, etc.).

In summary, we have used ALL-MBE to synthesize precise ultrathin bilayers
using metallic but non-superconducting LSCO and insulating LCO blocks, and
observed interface superconductivity with Tc = 34-36 K, significantly higher
than before. We have used synchrotron X-ray diffraction and the combined
COBRA/Difference-Map phase-retrieval method to determine accurately the atomic
structure and found the unit cell size to be constant despite dramatic atomic
displacements within the cell. In particular, near the surface the Cu$-$apical
O increases greatly, by as much as 0.45 Å, while it is known that variations
in the apical oxygen position strongly affect Tc. We conclude that in cuprates
the crystal structure can be modified in near-surface layers, and in such a
way that superconductivity properties can be dramatically altered. This result
amplifies the importance of high quality surface structure determination in
conjunction with surface sensitive probes of electronic states such as
scanning tunneling microscopy or angle-resolved photoemission spectroscopy.

This work was supported by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, Project MA-509-MACA, Contract No. DE-
AC02-98CH10886, and use of the Advanced Photon Source under Contract No. DE-
AC02-06CH11357. One of us (Y. Yacoby) would like to acknowledge with thanks
fruitful discussions with M. Björck.

## References

  * Reyren et al. (2007) N. Reyren et al., Science 317, 1196 (2007).
  * Bozovic et al. (2002) I. Bozovic et al., Phys. Rev. Lett. 89, 107001 (2002).
  * Gozar et al. (2008) A. Gozar et al., Nature 455, 782 (2008).
  * Ueno et al. (2008) K. Ueno et al., Nat. Mater. 7, 855 (2008).
  * Yuli et al. (2008) O. Yuli et al., Phys. Rev. Lett. 101, 057005 (2008).
  * Kivelson (2002) S. A. Kivelson, Physica B 318, 61 (2002).
  * Yunoki et al. (2007) S. Yunoki et al., Phys. Rev. B 76, 064532 (2007).
  * Okamoto and Maier (2008) S. Okamoto and T. A. Maier, Phys. Rev. Lett. 101, 156401 (2008).
  * Berg et al. (2008) E. Berg, D. Orgad, and S. A. Kivelson, Phys. Rev. B 78, 094509 (2008).
  * Butko et al. (2008) V. Butko et al., preprint (2008).
  * Fernandes et al. (1991) A. A. R. Fernandes et al., Phys. Rev. B 44, 7601 (1991).
  * Ohta et al. (1991) Y. Ohta, T. Tohyama, and S. Maekawa, Phys. Rev. B 43, 2968 (1991).
  * Feiner et al. (1992) L. F. Feiner, M. Grilli, and C. D. Castro, Phys. Rev. B 45, 10647 (1992).
  * Lubrittoy et al. (1996) C. Lubrittoy, K. Rosciszewski, and A. M. Olesz, J. Phys.: Condens. Matter 8, 11053 (1996).
  * Pavarini et al. (2001) E. Pavarini et al., Phys. Rev. Lett. 87, 047003 (2001).
  * Slezaka et al. (2008) J. A. Slezaka et al., Proc. Nat. Acad. Sci. 105, 3203 (2008).
  * Radovic et al. (2008) Z. Radovic, N. Bozovic, and I. Bozovic, Phys. Rev. B 77, 092508 (2008).
  * Wagner et al. (1993) J. L. Wagner et al., Physica C 210, 447 (1993).
  * Yacoby et al. (2002) Y. Yacoby et al., Nat. Mater. 1, 99 (2002).
  * Sowwan et al. (2002) M. Sowwan et al., Phys. Rev. B 66, 205311 (2002).
  * Willmott et al. (2007) P. R. Willmott et al., Phys. Rev. Lett. 99, 155502 (2007).
  * Yacoby et al. (2008) Y. Yacoby et al., Phys. Rev. B 77, 195426 (2008).
  * Cionca et al. (2008) C. N. Cionca et al., Appl. Phys. Lett. 92, 151914 (2008).
  * Fong et al. (2005) D. D. Fong et al., Phys. Rev. B 71, 144112 (2005).
  * Bozovic (2001) I. Bozovic, IEEE Trans. Appl. Supercond. 11, 2686 (2001).
  * Bozovic et al. (2003) I. Bozovic et al., Nature 422 422, 873 (2003).
  * Schlepütz et al. (2005) C. M. Schlepütz et al., Acta Crystallogr. Sect. A 61, 418 (2005).
  * Elser (2003) V. Elser, Acta Crystallogr. Sect. A 59, 201 (2003).
  * Björck et al. (2008) M. Björck et al., J. Phys.: Condens. Matter 20, 445006 (2008).
  * Ledbetter et al. (1992) H. Ledbetter, S. Kim, and A. Roshko, Z. Phys. B - Condensed Matter 89, 275 (1992).
  * Smadici et al. (2009) S. Smadici et al., Phys. Rev. Lett. 102, 107004 (2009).


# Wave and turbulence separation using dynamic mode decomposition

Julio Chávez-Dorado<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Washington, Seattle, WA 98195, USA  Isabel Scherl  Department of
Mechanical & Civil Engineering, California Institute of Technology, Pasadena,
CA 91125, USA  Michelle DiBenedetto  Department of Mechanical Engineering,
University of Washington, Seattle, WA 98195, USA

(February 29, 2024)

###### Abstract

Separating the effects of waves and turbulence in oceanographic time series is
an ongoing challenge because surface wave motion and turbulence fluctuations
can occur at overlapping frequencies. As a result, simple bandpass filters
cannot effectively separate their dynamics. While more advanced decomposition
techniques have been developed, they often entail restrictive assumptions
about the wave and turbulence interactions, require synchronized measurements,
and/or only decompose the signal spectrally without a time-series
reconstruction. To overcome these restrictions, we aim to apply modern data-
driven methods to this classical problem by using dynamic mode decomposition
(DMD). In this study, we present our new wave-turbulence decomposition
technique which uses DMD. The technique is designed to be signal-agnostic such
that it can be applied to any generic time series. Our approach requires
minimal tuning, where the main user input is the wave frequency range of
interest. To demonstrate the method, we apply it to synthetic, field, and
laboratory data, and compare the results to other published methods in the
literature. We conduct a sensitivity analysis on the synthetic data and find
that the most sensitive parameter to the accuracy is the rank truncation in
the DMD. Additionally, because of the accuracy of our decomposition technique,
we demonstrate how we are able to analyze the velocity autocorrelation of the
separated turbulence time series with minimal contamination by the waves.
Overall, we find that our decomposition method outperforms the other time-
series decomposition methods and provides for robust separation of the waves
and turbulence, demonstrating that it is a promising technique with wide
application to ocean signal processing.

††preprint: APS/123-QED

##  I Introduction

As observational advances increase the amount of ocean data available [1, 2],
effectively interpreting and using this information, both in post-processing
and real-time analysis, is imperative for gaining insights into the dynamics
of the ocean. One ongoing challenge is the interpretation of fluctuating data
that is influenced by both turbulence and surface gravity waves. This is
common in data obtained from the ocean surface and from the coastal ocean
where both turbulence and surface waves are expected to be strong. Separating
the turbulence and wave fluctuations in a signal is important for a variety of
analyses: e.g., isolating the turbulence fluctuations in any signal is
necessary for characterizing the turbulence Reynolds stress or for estimating
a scalar eddy covariance flux. And isolating the wave motion is necessary
e.g., for quantifying wave energy or for wave energy control strategies [3,
4]. We refer to this separation process as _wave-turbulence decomposition_.

The difficulties associated with wave-turbulence decomposition lie in the
overlapping frequency domains where waves and turbulence exist [5]—and that
both can manifest as broadband signals with correlation in space and time.
Even asserting that these signals can be effectively decomposed is an
assumption, given the nonlinear interactions that can occur between waves and
turbulence [6, 7, 8]. Over the years, numerous methods have been developed to
tackle this problem [9, 7, 10, 11, 12, 13, 14, 15, 16]. However, the various
methods are usually adapted to the specific data at hand [13, 17, 7], involve
restrictive assumptions (e.g., [9, 18]), and/or require multiple and often
complex synchronized measurements [11, 19, 17]. As a result, we still lack an
effective, universal technique for decomposing a generic time series.

Our goal is to develop a decomposition method that can be applied to a time
series of any type of data. However, for simplicity we will consider a
velocity time series as it is the most commonly decomposed signal in this
context. Consider the velocity decomposition

$u(x_{i},t)=\langle u(x_{i})\rangle+\tilde{u}(x_{i},t)+u^{\prime}(x_{i},t),$
(1)

where $u$ is the instantaneous horizontal velocity, $\langle u\rangle$ is the
time average, $\tilde{u}$ is the wave-orbital velocity and $u^{\prime}$ is the
turbulence velocity fluctuation. In two-dimensions, the result of Reynolds-
averaging the Navier-Stokes equations with this decomposition is the Reynolds
shear stress:

$-\frac{\tau}{\rho}=\langle\tilde{u}\tilde{w}\rangle+\langle\tilde{u}w^{\prime}\rangle+\langle
u^{\prime}\tilde{w}\rangle+\langle u^{\prime}w^{\prime}\rangle.$ (2)

The first term on the right hand side is identically zero for irrotational
waves described by linear wave theory but it can be non-zero if the coordinate
system is rotated, a common source of corrupted Reynolds stress measurements
[11]. The second and third terms are zero if there are no interactions between
the waves and turbulence. The fourth term is the turbulence Reynolds shear
stress, which is the only non-zero term as long as the aforementioned
assumptions regarding linear wave theory and wave-turbulence interaction hold.

Spectral filtering techniques [9, 18, 13, 15] are able to estimate bulk
turbulence quantities like the Reynolds stress with some accuracy, but they do
not allow for a time-series reconstruction of the decomposed velocity signals
because they are applied to the velocity power spectrum. Methods based on
nonlinear stream functions [20, 7, 10, 21] allow for a time-resolved
separation with fewer assumptions but only work for specialized scenarios.
More recent approaches involve dimensionality-reduction techniques [22, 23,
16] that make use of mode decomposition to capture the coherent, oscillating
nature of waves. However, it is common for these methods to remove turbulence
energy in excess from the wave frequency range. In addition, these methods
often require tuning parameters that lack a clear physical intuition.
Nevertheless, dimensionality reduction techniques such as proper orthogonal
decomposition (POD)[24], spectral proper orthogonal decomposition (SPOD) [25,
26], and dynamic mode decomposition (DMD) [27, 28] have become increasingly
popular in fluid dynamics research due to their adaptability and effectiveness
in extracting physically interpretable information from complex data. For
example, POD has been used in pattern recognition in chaotic flows [29], flow
control [30], and unsteady, high-Reynolds-number wake flows [31] and SPOD has
been employed in diverse flows such as boundary layers [32], mixing layers
[33], wakes [34], and jets [35], among others. Similarly, DMD is widely used
in the study of wakes [28], jets [28, 36, 37, 38], instabilities [39, 40], and
oscillations [36, 41, 29]. However, these methods have yet to be applied to
wave-turbulence decomposition, a problem for which they may be well-suited.
DMD shows promising filtering characteristics toward an accurate wave-
turbulence decomposition algorithm for single-point velocity measurements
since it captures the time dynamics of coherent structures present in the
signal.

Our approach is to develop a decomposition method that is as general as
possible. We seek a signal-agnostic framework, capable of handling one-
dimensional time-series data for a broad spectrum of flow parameters. Our
methodology does not rely on severe assumptions about wave nor turbulence
dynamics, making it adaptable to a variety of ocean and environmental
datasets. Thus, we will explore the effectiveness of DMD in a variety of
datasets: synthetic data generated with no wave-turbulence interactions, field
data collected in a swell-dominated bottom boundary layer, and laboratory data
collected in a surface boundary layer under wind-waves. The remainder of the
article has the following structure: section II gives relevant background on
other wave-turbulence decomposition methods and lays out the mathematical
basis for the mode decomposition technique used in our proposed method.
Section III outlines the procedure for our proposed decomposition method. In
section IV, we describe the datasets used for validation, and we present the
results of our decomposition in section V, including both discussion and a
sensitivity analysis. Finally, we summarize the study in section VI.

##  II Background

###  II.1 Wave-turbulence 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 [42, 43, 12, 44] 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._ [9] proposed a linear filtration technique that
leverages coherence between the pressure and velocity spectra. Similarly,
Trowbridge [11] and Shaw and Trowbridge [45] proposed a method that leverages
the coherence between two synchronized velocity measurements with a finite
spacing. With only one velocity measurement, Bricker and Monismith [18]
proposed a method that linearly filters the spectrum by assuming an inertial
range with a $-5/3$ slope. Whereas Gerbi  _et al._ [13] 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
[46]). 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 [20], which
imposed fewer assumptions and allowed for a time reconstruction of the
filtered signal. Jiang  _et al._ [7] 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 [10, 21] 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 wave-induced and
turbulent velocity time series, they are most well-suited to laboratory
conditions where dispersive effects can be neglected [10].

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
[14]. 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 [22]. 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 for time
series with time-varying dominant frequencies. Huang and Wu [14] 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. More recently, Bian  _et al._ [16] proposed the
synchrosqueezed wavelet transform (SWT)-based method [47, 48], 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 V.

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 have a significant impact on the turbulence in the wave frequency
range.

###  II.2 Modal decomposition of flow data

In the context of fluid dynamics, modal decomposition refers to the
identification and separation of _modes_ in high-dimensional flow data (e.g.,
generated from Particle Image Velocimetry). A _mode_ is extracted from a
dynamical system and characterizes a specific pattern. This type of
decomposition is favorable because the general motion of a system, even
nonlinear, nonstationary systems, can be approximated by a linear
superposition of its modes. The singular value decomposition (SVD) is the
basis for the implementation of several dimensionality-reduction techniques
such as proper orthogonal decomposition (POD) [24] and dynamic mode
decomposition (DMD) [28].

####  II.2.1 Singular value decomposition (SVD) and proper orthogonal
decomposition (POD)

The SVD decomposes a matrix, $\mathbf{M}\in\mathbb{R}^{m\times n}$ into the
product of three matrices

$\mathbf{M}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T},$

where $\mathbf{U}\in\mathbb{R}^{m\times m}$ and
$\mathbf{V}\in\mathbb{R}^{n\times n}$ are unitary matrices, and
$\mathbf{\Sigma}\in\mathbb{R}^{m\times n}$ is diagonal matrix with non-
negative real values. The singular values that make up the diagonal values
$\sigma_{k}\in\mathbf{\Sigma}$ are in descending order and weight the modes.
In other words, they are ordered based on how much of the total energy of the
original data matrix they possess. The columns of $\mathbf{U}$ (expressed as
$\mathbf{u}_{k}$) and the columns of $\mathbf{V}$ (expressed as
$\mathbf{v}_{k}$) are known as the left singular vectors and right singular
vectors of $\mathbf{M}$, respectively. A reduced-order approximation of the
matrix $\mathbf{M}$ can be achieved by truncation as

$\mathbf{M}_{r}\approx\sum_{k=1}^{r}\sigma_{k}\mathbf{u}_{k}\mathbf{v}^{T}_{k},$

where $\mathbf{M}_{r}$ is a rank $r$ approximation of $\mathbf{M}$ using the
$r\leq\min(m,n)$ most energetic singular modes.

The POD is simply an SVD applied to a specifically constructed $\mathbf{M}$
matrix. It was introduced in fluid mechanics by Lumley [24, 49], Sirovich [50,
51, 52], and Aubry [53]. The matrix $\mathbf{M}$ is constructed as follows:
given $n$ snapshots of flow data in a discrete spatial grid of $m$ points,
each consecutive snapshot is reshaped into a corresponding column vector of a
data matrix, $\mathbf{M}\in\mathbb{R}^{m\times n}$, where the $k$-th column of
$\mathbf{M}$, denoted as $\mathbf{x}_{k}$, contains all the spatial data of a
singular snapshot at time $k\Delta t$ or, similarly, each row contains a time
series of the flow data at a single point in space [54]. Thus we have,

$\mathbf{M}=\begin{bmatrix}|&|&|&&|\\\
\mathbf{x}_{1}&\mathbf{x}_{2}&\mathbf{x}_{3}&...&\mathbf{x}_{n}\\\
|&|&|&&|\end{bmatrix}.$

####  II.2.2 Dynamic mode decomposition (DMD)

Dynamic Mode Decomposition [28] is used to determine the temporal periodicity
of coherent structures present in flow data. DMD captures the system dynamics
in the eigenvectors (modes) and eigenvalues of a linear but infinite-
dimensional linear operator, known as the Koopman operator. DMD can separate
system dynamics into theses modes and can represent modal evolution over time.
Additionally, DMD has the advantage that it is tailored to recover oscillating
dynamics without any prior knowledge of the system dynamics.

To compute DMD, one begins with $n$ time steps of data across $m$ dimensions.
Each time step is shaped into a column vector of a present-state data matrix,
$\mathbf{X}\in\mathbb{R}^{m\times(n-1)}$ that contains time steps $1$ to
$n-1$, and a future-state data matrix,
$\mathbf{X}^{\prime}\in\mathbb{R}^{m\times(n-1)}$ that contains time steps $2$
to $n$. Thus we have,

$\mathbf{X}=\begin{bmatrix}|&|&|&&|\\\
\mathbf{x}_{1}&\mathbf{x}_{2}&\mathbf{x}_{3}&...&\mathbf{x}_{n-1}\\\
|&|&|&&|\end{bmatrix}\quad\text{and}\quad\mathbf{X}^{\prime}=\begin{bmatrix}|&|&|&&|\\\
\mathbf{x}_{2}&\mathbf{x}_{3}&\mathbf{x}_{4}&...&\mathbf{x}_{n}\\\
|&|&|&&|\end{bmatrix}.$

Similar to POD, the $k$-th column of either $\mathbf{X}$ or
$\mathbf{X}^{\prime}$, denoted as $\mathbf{x}_{k}$, contains all the data of a
singular time step at time $k\Delta t$ or, similarly, each row contains a time
series of the data for a single dimension. The DMD algorithm attempts to find
the best-fit linear operator
$\mathbf{A}:\mathbf{X}\mapsto\mathbf{X}^{\prime}\in\mathbb{R}^{m\times m}$
that advances the current state of the system, $\mathbf{X}$, one time step
into the future, $\mathbf{X}^{\prime}$, as

$\displaystyle\mathbf{X}^{\prime}=\mathbf{A}\mathbf{X}.$ (3)

The solution to equation 3 is

$\displaystyle\mathbf{A}=\mathbf{X}^{\prime}\mathbf{X}^{+},$ (4)

where $\mathbf{X}^{+}$ is the pseudo-inverse of the matrix $\mathbf{X}$.
Instead of directly finding the pseudo-inverse, the SVD of $\mathbf{X}$ yields
$\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}$ where $\mathbf{U}$,
$\mathbf{\Sigma}$, and $\mathbf{V}$ are unitary and, therefore, can be used to
express the pseudo-inverse as
$\mathbf{X}^{+}=\mathbf{V}\mathbf{\Sigma}^{-1}\mathbf{U}^{T}$, which yields

$\displaystyle\mathbf{A}\approx\mathbf{X}^{\prime}\mathbf{V}\mathbf{\Sigma}^{-1}\mathbf{U}^{T},$
(5)

where the eigenvalues and eigenvectors (DMD modes) of $\mathbf{A}$ capture the
dynamics of the system. The eigenvalues of $\mathbf{A}$ are complex numbers
that capture the growth/decay and oscillation of its corresponding DMD modes
in their real and imaginary component, respectively. The DMD modes of
$\mathbf{A}$ capture the spatial oscillatory patterns.

Calculating $\mathbf{A}$ directly would be computationally expensive. Instead,
a truncation of the matrices $\mathbf{U}$, $\mathbf{\Sigma}$, and $\mathbf{V}$
to retain only the first $r$ singular values results in the reduced-rank
matrices $\mathbf{U}_{r}$, $\mathbf{\Sigma}_{r}$, and $\mathbf{V}_{r}$ and
yields a best-fit linear operator $\tilde{\mathbf{A}}\in\mathbb{R}^{r\times
r}$ by projecting $\mathbf{A}$ onto the spatial modes $\mathbf{U}_{r}$

$\displaystyle\tilde{\mathbf{A}}=\mathbf{U}_{r}^{T}\mathbf{A}\mathbf{U}_{r}=\mathbf{U}_{r}^{T}\mathbf{X}^{\prime}\mathbf{V}_{r}\mathbf{\Sigma}_{r}^{-1}.$
(6)

The operator $\tilde{\mathbf{A}}$ has the same discrete-time eigenvalues,
$\mathbf{\Lambda}$, as $\mathbf{A}$. The eigenvectors, $\mathbf{W}$, of
$\tilde{\mathbf{A}}$ can be calculated as
$\tilde{\mathbf{A}}\mathbf{W}=\mathbf{\Lambda}\mathbf{W}$ and converted to the
high-dimensional eigenvalues, $\mathbf{\Phi}$, of $\mathbf{A}$, as follows
[55]

$\displaystyle\mathbf{\Phi}=\mathbf{X}^{\prime}\mathbf{V}_{r}\mathbf{\Sigma}_{r}^{-1}\mathbf{W}.$
(7)

From the leading eigenvalues, $\mathbf{\Lambda}$, and eigenvectors,
$\mathbf{\Phi}$, the solution of the system can be constructed as

$\displaystyle\mathbf{X}=\mathbf{\Phi}\exp^{\mathbf{\Omega}t}\mathbf{b}=\sum_{k=1}^{r}\phi_{k}\exp^{\omega_{k}t}b_{k},$
(8)

where $\mathbf{\Omega}=\log(\mathbf{\Lambda})/dt$ is the continuous-time
eigenvalue of $\mathbf{A}$, and
$\mathbf{b}=\mathbf{\Phi}\backslash\mathbf{X}(t_{0})$ is the initial condition
of the system.

##  III Methods

Figure 1 shows a diagram summarizing the algorithm workflow, which is
described below. For the method we propose, the input matrix contains a one-
dimensional time series of discrete, point-wise flow measurements such as one
velocity component ($u$, $v$, or $w$), or other variables like water elevation
$\eta$, scalar concentration $C$, etc. These measurements are collected at a
sampling frequency $f_{s}$ over a specified duration $T$, for $N=Tf_{s}$ total
observations in the given time series. Our approach is to use DMD to extract
the most coherent features of the signal, which we assume to be the waves, as
low-rank structures and subtract them from the raw signal to isolate the
turbulence. Although other order-reduction techniques exist, they have
limitations. For example, POD identifies correctly the rank of the system but
since this method is based entirely on energy (or contained variance), the
modes are not temporally orthogonal. In other words, POD can hierarchically
rank coherent structures in the signal based on their energy content but it
cannot separate the coherent temporal evolution in time. This is not ideal for
our problem because we expect waves to be represented as modes that have
coherent oscillations in time. In contrast, DMD can separate the dynamics in
the system and can identify their temporal oscillation frequency. This allows
us to identify a specific frequency per mode that we can associate with wave
motion [56].

In frequency space, waves can be narrow-banded when compared with turbulence,
but their energy can still span over multiple orders of magnitude in
frequency. The broader the wave frequency range, the larger the number of DMD
modes will be needed to represent the waves. However, the number of modes that
DMD can extract is limited by the rank of the input data matrix. Because we
want to extract a range of wave modes, we need to increase the rank of our
time series before we apply DMD.

Hence, we propose to address this problem and increase the rank of the input
data matrix by constructing a delay embedding of the time-series signals,
i.e.,

$\mathbf{H}=\left[\begin{array}[]{cccccc}x_{1}&x_{2}&x_{3}&\ldots&\ldots&x_{n}\\\
x_{2}&x_{3}&&&&\vdots\\\ x_{3}&&&&&\vdots\\\ \vdots&&&&&x_{N-2}\\\
\vdots&&&&x_{N-2}&x_{N-1}\\\
x_{m}&\ldots&\ldots&x_{N-2}&x_{N-1}&x_{N}\end{array}\right],$ (9)

where $m=N-n$ is the number of time lags, and $n$ is the number of columns.
Note that the matrix $\mathbf{H}$ can have multiple rows besides the time-
delay process if one has multiple data series (e.g., see [57, 58, 59]). The
introduction of the time delay serves a crucial role in increasing the rank of
the data by increasing the “information” content of the matrix [60].
Physically, this increase in rank can be analogous to Taylor’s frozen
hypothesis, since the time observations of the signals are converted to
spatial observations by construction of the time-delayed matrix $\mathbf{H}$.

In order to filter the signal using DMD, we first need to determine the wave
frequency range, the dimension of the time-delayed matrix $\mathbf{H}$, and
the rank truncation, $r$. Visual inspection of the signal power spectrum helps
determine the wave frequency range by identifying the points were the wave
peak meet the expected turbulence energy spectrum. For the construction of the
time-delayed matrix, the number of columns $n$ needs to be large enough to
capture the lowest wave frequencies of interest in each row. Based on our
sensitivity analysis described below, we find that a good rule of thumb is
within the range $n/N\approx[1/3,1/2]$; however, in our data the decomposition
is not as sensitive to this parameter as it is to the rank truncation $r$,
especially given a sufficiently long time series. In general, we suspect this
parameter might require some tuning in order to determine the appropriate
dimensions of $\mathbf{H}$. From this new configuration $\mathbf{H}$ of the
input data, the present and future-state matrices $\mathbf{X}$ and
$\mathbf{X}^{\prime}$ are constructed as indicated in section II.2.2 as input
for the DMD algorithm. To determine the rank truncation $r$ we perform a POD
decomposition of
$\mathbf{H}=\mathbf{U}_{H}\mathbf{\Sigma}_{H}\mathbf{V}_{H}^{T}$ and compute
the spectral density of each column $\mathbf{v}_{j}\in\mathbf{V}_{H}$ (known
as time expansion coefficients) as
$\mathbf{C}_{j}=|\mathcal{F}\\{\mathbf{v}_{j}\\}|^{2}/(f_{s}N)$. The vector
$\mathbf{C}_{j}$ contains the energy content at each frequency of the $j$-th
POD mode of $\mathbf{H}$ [25]. Plotting the energy content of each singular
value mode $\mathbf{C}$ as a function of frequency helps us visualize which
modes contain wave and/or turbulence energy. Next, we identify and discard the
modes that do not contain wave energy. The parameter $r$ consists of the
lowest modes that contain wave energy. This truncation allows us to remove the
higher-frequency turbulence from the data directly.

DMD is then applied $\mathbf{H}$, using the rank truncation $r$ that we just
determined. Following the procedure outlined in the previous section, we
obtain the eigenvalues $\mathbf{\Lambda}$, eigenvectors (or modes)
$\mathbf{\Phi}$, and amplitudes $\mathbf{b_{0}}$, which capture the coherent
structures in the matrix $\mathbf{H}$. The $j$-th DMD mode $\phi_{j}$ and its
corresponding eigenvalue $\lambda_{j}$ have a distinct oscillation frequency
$f_{j}$ that can be calculated as

$f_{j}=\frac{\mathrm{Im}[\log(\lambda_{j})]}{2\pi dt}.$ (10)

The modes with frequencies within the wave range are referred to as wave
modes, and are transformed to time domain as a wave time series using equation
8. The wave time series is then subtracted from the raw signal, leaving the
turbulence component of the signal.

Figure 1: Workflow of the proposed DMD wave-turbulence separation methodology.

##  IV Dataset descriptions

We test our method on three distinct datasets to show its broad applicability:
a synthetic dataset with known wave and turbulence components; a field dataset
collected on a swell-dominated coastal shelf in bottom boundary layer
published in Reimers and Fogaren [61], and laboratory data we collected in a
surface boundary layer in a wind-driven wave tank. Table 1 shows some flow and
DMD parameters used in our proposed decomposition method for each dataset.
Additionally, it shows the relative wave intensity,
$\tilde{u}_{\text{rms}}/u^{\prime}_{\text{rms}}$, estimated as the ratio of
the wave to the turbulence root-mean-square velocity obtained from the DMD-
based decomposition. Note that we cover a range of wave intensities and that
the field and laboratory waves are the most and least intense, respectively.

Data type | $N$ |  $f_{s}$ (Hz) |  $T_{p}$ (s) | Wave range (Hz) | $m$ | $n$ | $r$ | $r/n$ | $\psi_{\text{DMD}}$ | $\tilde{u}_{\text{rms}}/u^{\prime}_{\text{rms}}$
---|---|---|---|---|---|---|---|---|---|---
Synthetic | 10000 | 10 | 5 | $0.117-0.6~{}~{}~{}$ | 3000 | 5911 | 173 | 0.029 | 104 | 1.92
Field | 7200 | 8 | 14.29 | $0.031-0.175$ | 5700 | 1500 | 48 | 0.027 | 34 | 8.05
Laboratory, $u$ | 18000 | 30 | 0.5 | $1.02-3.16$ | 13000 | 5000 | 1000 | 0.200 | 506 | 0.30
Laboratory, $w$ | 18000 | 30 | 0.5 | $1.02-3.16$ | 12500 | 6500 | 1300 | 0.200 | 732 | 0.72
Table 1: Summary of dataset characteristics including some input and output
parameters from the wave-turbulence decomposition. Characteristics include
length of time series $N$, sampling frequency $f_{s}$, wave peak period
$T_{p}$, and wave frequency range. We also report the rows $m$ and columns $n$
of the time-delayed matrix $\mathbf{H}$ and the rank truncation $r$ used,
including the relative rank $r/n$ and the number of DMD modes
$\psi_{\text{DMD}}$ used to reconstruct the wave time series. The final column
reports the relative wave intensity
$\tilde{u}_{\text{rms}}/u^{\prime}_{\text{rms}}$ which is known for the
synthetic data, but calculated using the DMD decomposed velocities for the
field and laboratory data.

###  IV.1 Synthetic data

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=\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 [62]. The turbulence
spectrum has a developed inertial range in the 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 [63, 64]

$S(\omega)=\frac{\alpha
g^{2}}{\omega^{5}}\exp\left[-\frac{5}{4}\left(\frac{\omega_{p}}{\omega}\right)\right]\gamma^{\exp\left[-\frac{(\omega/\omega_{p}-1)^{2}}{2\beta^{2}}\right]},$
(11)

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:

$\displaystyle\tilde{u}(x,z,t)$
$\displaystyle=\sum_{i=1}^{N_{w}}U_{i}\cos(k_{i}x-\omega_{i}t+\phi_{i})\exp(k_{i}z),$
(12)

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 figure 2,
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.

###  IV.2 Field data

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 [61]. 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 figure 4) 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
figure 4. For further details about the data acquisition, please consult the
original study [61].

###  IV.3 Laboratory data

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. Tap 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 [65] 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 (figure 5 and 6). We post-processed the data using
the despiking method proposed by Goring and Nikora [66]. Even after despiking
the data, there is still some noise present in the data due to intermittent
wave-breaking which generated large bubbles and 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 complex nature
of the strongly wind-sheared flow and the intermittent wave-breaking, we see
some deviations from the expected -5/3 slope in the inertial range.

##  V Results and Discussion

In this section, we present the results obtained from analysing the synthetic,
laboratory, and field data with our proposed DMD method. We demonstrate our
DMD decomposition by applying it to all three datasets after first conducting
a sensitivity analysis on the synthetic data. Next, we compare our results to
the other two time-resolved mode-based decomposition methods: EEMD [67, 23]
and Synchrosqueezed Wavelet transform (SWT) [16]. After assessing the
decomposition in both spectral and time space, we analyze the temporal
autocorrelation of the decomposed turbulence signal, further demonstrating how
our DMD method outperforms the other methods.

###  V.1 Synthetic data decomposition

*[]*[]

*[]

*[]

*[]

Figure 2: DMD wave-turbulence decomposition of $u$-component velocity of
synthetic data. (a) POD spectrum of time-delayed matrix $\mathbf{H}$. (b)
Power spectra of the original data and DMD-separated wave and turbulence. A
-5/3 slope line is plotted for reference. In (c) we plot a portion of the
original time series, (d) the true and separated wave components using 104 DMD
modes, (e) and the true and separated turbulence components.

The first step of our decomposition method is constructing the time-delayed
matrix $\mathbf{H}$ and identifying the wave frequency range, as described in
section III and shown in table 1. Next, we determine the rank truncation $r$
by plotting the POD spectrum of $\mathbf{H}$ of the synthetic $u$ data (figure
2) 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 transfer from low to high frequencies, 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 low-frequency turbulence.

After determining the rank truncation, we are able to decompose the signal
using DMD. Figure 2 shows the spectra of the DMD-filtered wave and turbulence
time series which we compare with the true wave and 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 resolution only resolved 104 modes, we do not expect a perfect
reconstruction. In this method, the DMD is always discretizing a continuous
wave frequency range, so we will 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 less modes
compared to the other datasets.

(a) $\epsilon=0.1586$, at $r=173$, $n/N=0.5911$ ($N=10000$).

(b) $\epsilon=0.1894$, at $r=173$, $n/N=0.509$ ($N=7500$).

(c) $\epsilon=0.2109$, at $r=165$, $n/N=0.540$ ($N=5000$).

(d) $\epsilon=0.2737$, at $r=210$, $n/N=0.470$ ($N=2500$).

Figure 3: Sensitivity analysis to asses the separation performance of the DMD-
based algorithm for different time series length, $N$. A measurement of error,
$\epsilon$, was calculated for every combination of rank truncation, $r$, and
normalized number of columns, $n/N$, in the time-delayed matrix, $\mathbf{H}$,
and the optimal combination is indicated by a red cross.

Parameters | Length of raw time series, $N$
---|---
2500 |  5000 |  7500 |  10000
$r$ | 210 | 165 | 155 | 173
$n$ | 1175 | 2700 | 3900 | 5911
$n/N$ | 0.47 | 0.54 | 0.52 | 0.59
$\epsilon$ | 0.2737 | 0.2109 | 0.1894 | 0.1586

Method | $\quad\epsilon\quad$ |  $\quad\alpha$ (m/s) | $\quad\alpha/u^{\prime}_{\text{rms}}\quad$
---|---|---|---
DMD |  0.159 |  0.039 |  0.24
EEMD |  0.268 |  0.065 |  0.40
SWT |  0.325 |  0.079 |  0.49

Table 2: (Left) Summary of sensitivity analysis from figure 3 and (Right)
comparison of errors from DMD and other modal decomposition methods applied to
the longest time-series length of $N=10000$.

We conducted a sensitivity analysis on the decomposition based on the two main
tuning parameters: the rank truncation $r$ and the number of columns of the
time-delayed input matrix $n$ in order to identify the optimal $r$ and $n$ for
the decomposition. In addition, we also varied the length of the time series
$N$ and repeated the sensitivity analysis to assess how the optimal $r$ and
$n$ vary with $N$. We quantified the accuracy of the decomposition by
minimizing the normalized error $\epsilon$, which is defined as follows:

$\epsilon=\frac{\|u^{\prime}_{i}-u^{\prime}_{i,\text{DMD}}\|_{2}}{\|u^{\prime}_{i}\|_{2}},$
(13)

where $u^{\prime}_{i}$ and $u^{\prime}_{i,\text{DMD}}$ are the true and
decomposed turbulence time series (here the notation refers to DMD, but we
also calculate this error for the EEMD and SWT methods), respectively, at time
$t_{i}$, and $\|\cdot\|_{2}$ is the $L_{2}$ norm. For comparison, we also
measured the mean absolute error $\alpha$ defined as

$\alpha=\frac{1}{N}\sum_{i=0}^{N}|(u^{\prime}_{i}-u^{\prime}_{i,\text{DMD}})|.$
(14)

We conduct the sensitivity analysis across four different time-series lengths
to assess the importance of the length of the dataset. Figures 3(a-d), show
the results of the sensitivity analysis for $N=10000,7500,5000$, and $2500$,
respectively. We plot the norm-based error $\epsilon$ for different
combinations of $r$ and normalized number of columns $n/N$. The optimal
combination is indicated by a red cross. Overall, the optimal $n/N$ ratio is
not too sensitive and ranges between 0.5 and 0.6, while the optimal $r$
minimizes the error only in a narrow range of values. This indicates that, at
least for this simple case, the preferential dimension of the input time-
delayed matrix is an approximately square matrix. This might be due to the
fact that for a large dataset we have enough physical information to resolve
coherent features in the time domain, $n$, and the “spatial” domain, $m$, that
we constructed via time-delaying. However, when there is limited information
by the shortened time-series, as in the case of figure 3(d), the preferential
input matrix shape is a tall, skinny matrix. This suggests a preference for
“spatial” information over temporal information, and that meeting a critical
column length is more important than increased row length. We can observe that
the optimal $r$ shows a small decrease with $N$ (except for $N=2500$) but
overall is fairly stable. Selecting an appropriate rank truncation is
essential for accurately capturing the waves in the time series. If $r$ is too
low, some waves may be missed by the DMD algorithm, leading to an incomplete
reconstructed time series. Conversely, if $r$ is too high, it may lead to
overestimation of the wave energy.

Table 2 summarizes the different tuning parameters for all time-series length
$N$ and their associated errors. We see a decrease of the norm-based error
$\epsilon$ with $N$ which makes sense since the DMD is a data-based algorithm,
and more data generally improves performance. However, as the data increases
so does the computational cost of performing the decomposition. We compare the
DMD method to EEMD and SWT in table 2, where we see that given the same time-
series length, our DMD method outperforms the other modal decomposition
methods. While the EEMD outperforms the SWT, they both have higher errors than
the DMD method, even when DMD was applied to only half the data ($N=5000$).
Overall, our method seems to be perform the best when applied to this dataset.
At the end of this section, we compare the methods in more detail by applying
them to all of the datasets.

###  V.2 Field data decomposition

*[]*[]

*[]

Figure 4: DMD wave-turbulence decomposition of $u$-component velocity of field
data. (a) POD spectrum of time-delayed matrix $\mathbf{H}$. (b) Power spectra
of the original data and DMD-separated wave and turbulence. (c) A portion of
the time-series of the original and decomposed signals. The wave motion was
reconstructed using 34 DMD modes.

We plot the POD spectrum of the field data (figure 4), and find that similar
to the synthetic data in figure 2, there are two features: one associated with
the turbulence energy that is seen across the entire frequency range, and a
dip in the energetic modes associated with the the wave energy concentrated
between 0.04 and 0.16 Hz, marked with the vertical dotted lines. 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 figure 4 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 Figure 4, 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. We expect this to be true for all signals: that the higher the
wave intensity, the fewer modes relative to the total rank are needed.

###  V.3 Laboratory data decomposition

*[]*[]

*[]

Figure 5: DMD wave-turbulence decomposition of $u$-component velocity of
laboratory data. (a) POD spectrum of time-delayed matrix $\mathbf{H}$. (b)
Power spectra of the original data and DMD-separated wave and turbulence. (c)
A portion of the time-series of the original and decomposed signals. The wave
motion was reconstructed using 506 DMD modes.

*[]*[]

*[]

Figure 6: DMD wave-turbulence decomposition of $w$-component velocity of
laboratory data. (a) POD spectrum of time-delayed matrix $\mathbf{H}$. (b)
Power spectra of the original data and DMD-separated wave and turbulence. (c)
A portion of the time-series of the original and decomposed signals. The wave
motion was reconstructed using 732 DMD modes.

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$ (figures 5 and 6) 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 figures 5 and 6. 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 (figure 5), 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 figure 5 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 figure 6 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 figure 6
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 visualize the time
series of $\tilde{u}_{\textrm{DMD}}$ and $\tilde{w}_{\textrm{DMD}}$ together
in figure 7. 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 figure 7 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.

Figure 7: The $u$ and $w$ velocity components of DMD-separated wave time
series of the laboratory data.

###  V.4 Comparison with other decomposition methods

*[]*[]

*[]*[]

Figure 8: Comparison of modal decomposition methods, showing the spectra of
the original data and the decomposed turbulence data. (a) Synthetic data $u$
velocity. (b) Field data $u$ velocity. (c) and (d) Laboratory data $u$ and $w$
velocity, respectively.

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 Figure 8. 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 figures
8 and 8, respectively. We note though that the EEMD does not necessarily
conserve energy, and has actually added energy into the turbulence in the
field data, as seen in figure 8 between 0.2 to 0.4 Hz. Finally, the laboratory
data (figures 8 and 8) 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.

*[]*[]

*[]*[]

Figure 9: Normalized autocorrelation of raw (gray), DMD-filtered turbulence
(pink), EEMD- and SWT-filtered turbulence (green and purple, respectively)
signals. (a) Synthetic data $u$ velocity, the black line is the true
turbulence. (b) Field data $u$ velocity. (c) and (d) Laboratory data $u$ and
$w$ velocity, respectively.

We further compare the methods by computing the normalized autocorrelations
$R(\tau)$ of the filtered turbulence velocity across the different datasets,
as shown in figure 9. 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 (figure 9), we
know what the true turbulence signal is, and we see that while all the methods
follow the true turbulence autocorrelation closely, shown with a black line,
the SWT performs the worst, as shown by the small peak in the signal after the
initial decrease. In figure 9, 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 (figures 9 and 9), 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.

##  VI Conclusion

While many wave-turbulence decomposition methods have been developed over the
years, they are typically limited in their use case. For example, they may be
restricted to the spectral domain precluding a time-series construction, they
may require multiple synchronized measurements, or they may have restrictive
assumptions such as no wave-turbulence interactions. In many scenarios, this
assumption can easily break, especially because waves and turbulence do
interact across different time scales [68]. Given these challenges, newer
modal decomposition methods have been developed, such as the EEMD and the SWT.
These methods still have their drawbacks. For example, the SWT uses a
predetermined basis to project the data whereas the EEMD introduces noise as a
filter to differentiate scales in the time series. These approaches are not
optimal for highly nonlinear data, resulting in either poor or excessive wave
energy separation. To address these limitations, we developed a new wave-
turbulence decomposition method using DMD, leveraging its data-driven adaptive
basis and minimal assumptions.

Our method applies DMD to a single time-series, assisted by a time-delay
embedding. It is able to isolate the wave-components without significantly
affecting the turbulent signal at similar frequencies. In the case of velocity
data, we show how the method is able to recover the inertial range scaling
below the wave peak. The main assumptions that underlie this method are that
the waves and turbulence can be separated and that waves are the most coherent
feature of the time-series, which are reasonable assumptions considering the
dynamics of the system. We applied our method to synthetic, field, and
laboratory data which covered a range of wave intensities. Our method clearly
outperformed the SWT and EEMD when applied to the synthetic and laboratory
data; whereas all methods performed similarly when applied to the field data.
This may be due to the fact that the field data had the largest wave
intensity, and therefore the wave motion was relatively easy to isolate from
the turbulence across all of the methods.

The proposed method does require some manual tuning: specifically the user
needs to input the shape of the time-delay matrix, the rank truncation, and
the wave frequency range of interest. However, we argue that the parameters
are physically related to the signal, especially compared to other methods
such as the SWT and EEMD which have somewhat arbitrary tuning parameters. The
rank truncation is discernible in the POD spectrum of the time-delayed input
matrix, and the wave-frequency range is clearly visible in the power spectrum
of the signal. While the number of time-delays is less clear, our sensitivity
analysis shows that it is not as important to the accuracy of the
decomposition, especially when the time-series length is long enough. We
recommend future work to examine this parameter more closely and relate it
directly to the wave spectrum. From the analyzed datasets, we found that the
ratio $r/n$ is inversely proportional to the wave intensity
$\tilde{u}_{\text{rms}}/u^{\prime}_{\text{rms}}$. The reported parameters in
table 1 can be a starting point for users of the method. Overall, the proposed
DMD-based decomposition method is successful at decomposing wave and
turbulence motion with minimal tuning and only one component of data, which is
valuable for time-series measurements of flow parameters.

We have demonstrated that DMD is powerful for this application, and some fine
tuning may need to occur as it is adopted for this problem. Future work can
focus on more robust DMD schemes and alternative ways of determining the
number of modes necessary to reconstruct the wave motion. A few examples
include the extended DMD [69] which attempts to circumvent certain constraints
of the (linear) DMD technique when decomposing data from nonlinear systems,
and sparsity-promoting DMD [70] that can reduce the complexity of the
calculations by identifying a smaller set of important DMD modes.

## Acknowledgements

We would like to thank Steve Brunton and Caroline Cardinale for helpful
comments and insights. We acknowledge support from the U.S. National Science
Foundation under grant No. CBET-2237550 and the Link Foundation Ocean
Engineering and Instrumentation Ph.D. Fellowship.

## References

  * Smith  _et al._ [2019] S. R. Smith, G. Alory, A. Andersson, W. Asher, A. Baker, D. I. Berry, K. Drushka, D. Figurskey, E. Freeman, P. Holthus, T. Jickells, H. Kleta, E. C. Kent, N. Kolodziejczyk, M. Kramp, Z. Loh, P. Poli, U. Schuster, E. Steventon, S. Swart, O. Tarasova, L. P. de la Villéon, and N. Vinogradova-Shiffer, Ship-Based Contributions to Global Ocean, Weather, and Climate Observing Systems, Frontiers in Marine Science 6 (2019).
  * Rosa  _et al._ [2021] T. L. Rosa, A. M. Piecho-Santos, R. Vettor, and C. Guedes Soares, Review and Prospects for Autonomous Observing Systems in Vessels of Opportunity, Journal of Marine Science and Engineering 9, 366 (2021).
  * Li  _et al._ [2012] G. Li, G. Weiss, M. Mueller, S. Townley, and M. R. Belmont, Wave energy converter control by wave prediction and dynamic programming, Renewable Energy 48, 392 (2012).
  * Perez  _et al._ [2020] L. Perez, R. Cossu, C. Couzi, and I. Penesis, Wave-Turbulence Decomposition Methods Applied to Tidal Energy Site Assessment, Energies 13, 1245 (2020).
  * Lumley and Terray [1983] J. L. Lumley and E. A. Terray, Kinematics of Turbulence Convected by a Random Wave Field, Journal of Physical Oceanography 13, 2000 (1983).
  * Grant and Madsen [1979] W. D. Grant and O. S. Madsen, Combined wave and current interaction with a rough bottom, Journal of Geophysical Research: Oceans 84, 1797 (1979).
  * Jiang  _et al._ [1990] J.-Y. Jiang, R. L. Street, and S. P. Klotz, A study of wave-turbulence interaction by use of a nonlinear water wave decomposition technique, Journal of Geophysical Research: Oceans 95, 16037 (1990).
  * Magnaudet and Thais [1995] J. Magnaudet and L. Thais, Orbital rotational motion and turbulence below laboratory wind water waves, Journal of Geophysical Research: Oceans 100, 757 (1995).
  * Benilov  _et al._ [1974] A. Yu. Benilov, O. A. Kouznetsov, and G. N. Panin, On the analysis of wind wave-induced disturbances in the atmospheric turbulent surface layer, Boundary-Layer Meteorology 6, 269 (1974).
  * Thais and Magnaudet [1995] L. Thais and J. Magnaudet, A triple decomposition of the fluctuating motion below laboratory wind water waves, Journal of Geophysical Research: Oceans 100, 741 (1995).
  * Trowbridge [1998] J. H. Trowbridge, On a Technique for Measurement of Turbulent Shear Stress in the Presence of Surface Waves, Journal of Atmospheric and Oceanic Technology 15, 290 (1998).
  * Williams  _et al._ [2003] J. J. Williams, P. S. Bell, and P. D. Thorne, Field measurements of flow fields and sediment transport above mobile bed forms, Journal of Geophysical Research: Oceans 108, 10.1029/2002JC001336 (2003).
  * Gerbi  _et al._ [2008] G. P. Gerbi, J. H. Trowbridge, J. B. Edson, A. J. Plueddemann, E. A. Terray, and J. J. Fredericks, Measurements of Momentum and Heat Transfer across the Air–Sea Interface, Journal of Physical Oceanography 38, 1054 (2008).
  * Huang and Wu [2008] N. E. Huang and Z. Wu, A review on Hilbert-Huang transform: Method and its applications to geophysical studies, Reviews of Geophysics 46, 10.1029/2007RG000228 (2008).
  * Young and Webster [2018] D. L. Young and D. R. Webster, Assessment of single-instrument techniques for removing wave bias from Reynolds stress estimates, Limnology and Oceanography: Methods 16, 35 (2018).
  * Bian  _et al._ [2018] C. Bian, Z. Liu, Y. Huang, L. Zhao, and W. Jiang, On Estimating Turbulent Reynolds Stress in Wavy Aquatic Environment, Journal of Geophysical Research: Oceans 123, 3060 (2018).
  * Feddersen and Williams [2007] F. Feddersen and A. J. Williams, Direct Estimation of the Reynolds Stress Vertical Structure in the Nearshore, Journal of Atmospheric and Oceanic Technology 24, 102 (2007).
  * Bricker and Monismith [2007] J. D. Bricker and S. G. Monismith, Spectral Wave–Turbulence Decomposition, Journal of Atmospheric and Oceanic Technology 24, 1479 (2007).
  * Doron  _et al._ [2001] P. Doron, L. Bertuccioli, J. Katz, and T. R. Osborn, Turbulence Characteristics and Dissipation Estimates in the Coastal Ocean Bottom Boundary Layer from PIV Data, Journal of Physical Oceanography 31, 2108 (2001).
  * Dean [1965] R. G. Dean, Stream function representation of nonlinear ocean waves, Journal of Geophysical Research (1896-1977) 70, 4561 (1965).
  * Thais and Magnaudet [1996] L. Thais and J. Magnaudet, Turbulent structure beneath surface gravity waves sheared by the wind, Journal of Fluid Mechanics 328, 313 (1996).
  * Huang  _et al._ [1998] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis, Proceedings: Mathematical, Physical and Engineering Sciences 454, 903 (1998), 53161 .
  * Huang  _et al._ [2009] Y. Huang, F. G. Schmitt, Zhiming. Lu, and Y. Liu, Analysis of daily river flow fluctuations using empirical mode decomposition and arbitrary order Hilbert spectral analysis, Journal of Hydrology 373, 103 (2009).
  * Lumley [1967] J. L. Lumley, The structure of inhomogeneous turbulent flows, Atmospheric Turbulence and Radio Wave Propagation , 166 (1967).
  * Towne  _et al._ [2018] A. Towne, O. T. Schmidt, and T. Colonius, Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis, Journal of Fluid Mechanics 847, 821 (2018).
  * Schmidt and Colonius [2020] O. T. Schmidt and T. Colonius, Guide to Spectral Proper Orthogonal Decomposition, AIAA Journal 58, 1023 (2020).
  * Rowley  _et al._ [2009] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, Spectral analysis of nonlinear flows, Journal of Fluid Mechanics 641, 115 (2009).
  * Schmid [2010] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of Fluid Mechanics 656, 5 (2010).
  * Albidah  _et al._ [2021] A. B. Albidah, W. Brevis, V. Fedun, I. Ballai, D. B. Jess, M. Stangalini, J. Higham, and G. Verth, Proper orthogonal and dynamic mode decomposition of sunspot data, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 379, 20200181 (2021).
  * Gordeyev and Thomas [2013] S. V. Gordeyev and F. O. Thomas, A temporal proper decomposition (TPOD) for closed-loop flow control, Experiments in Fluids 54, 1477 (2013).
  * Durgesh and Naughton [2010] V. Durgesh and J. W. Naughton, Multi-time-delay LSE-POD complementary approach applied to unsteady high-Reynolds-number near wake flow, Experiments in Fluids 49, 571 (2010).
  * Tutkun and George [2017] M. Tutkun and W. K. George, Lumley decomposition of turbulent boundary layer at high Reynolds numbers, Physics of Fluids 29, 020707 (2017).
  * Braud  _et al._ [2004] C. Braud, D. Heitz, G. Arroyo, L. Perret, J. Delville, and J.-P. Bonnet, Low-dimensional analysis, using POD, for two mixing layer–wake interactions, International Journal of Heat and Fluid Flow Turbulence and Shear Flow Phenomena (TSFP-3), 25, 351 (2004).
  * Araya  _et al._ [2017] D. B. Araya, T. Colonius, and J. O. Dabiri, Transition to bluff-body dynamics in the wake of vertical-axis wind turbines, Journal of Fluid Mechanics 813, 346 (2017).
  * Heidt and Colonius [2023] L. Heidt and T. Colonius, Spectral proper orthogonal decomposition of harmonically forced turbulent flows (2023), arxiv:2305.05628 [physics] .
  * Seena and Sung [2011] A. Seena and H. J. Sung, Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations, International Journal of Heat and Fluid Flow 32, 1098 (2011).
  * Schmid [2011] P. J. Schmid, Application of the dynamic mode decomposition to experimental data, Experiments in Fluids 50, 1123 (2011).
  * Semeraro  _et al._ [2012] O. Semeraro, G. Bellani, and F. Lundell, Analysis of time-resolved PIV measurements of a confined turbulent jet using POD and Koopman modes, Experiments in Fluids 53, 1203 (2012).
  * Duke  _et al._ [2012] D. Duke, D. Honnery, and J. Soria, Experimental investigation of nonlinear instabilities in annular liquid sheets, Journal of Fluid Mechanics 691, 594 (2012).
  * Grilli  _et al._ [2012] M. Grilli, P. J. Schmid, S. Hickel, and N. A. Adams, Analysis of unsteady behaviour in shockwave turbulent boundary layer interaction, Journal of Fluid Mechanics 700, 16 (2012).
  * Massa  _et al._ [2012] L. Massa, R. Kumar, and P. Ravindran, Dynamic mode decomposition analysis of detonation waves, Physics of Fluids 24, 066101 (2012).
  * Foster [1997] D. L. Foster,  _Dynamics of the Nearshore Wave Bottom Boundary Layer_ , Ph.D. thesis, Oregon State University, United States – Oregon (1997).
  * Smyth  _et al._ [2002] C. Smyth, A. E. Hay, and L. Zedel, Coherent Doppler Profiler measurements of near-bed suspended sediment fluxes and the influence of bed forms, Journal of Geophysical Research: Oceans 107, 19 (2002).
  * Zhu  _et al._ [2016] Q. Zhu, B. C. van Prooijen, Z. B. Wang, Y. X. Ma, and S. L. Yang, Bed shear stress estimation on an open intertidal flat using in situ measurements, Estuarine, Coastal and Shelf Science 182, 190 (2016).
  * Shaw and Trowbridge [2001] W. J. Shaw and J. H. Trowbridge, The Direct Estimation of Near-Bottom Turbulent Fluxes in the Presence of Energetic Wave Motions, Journal of Atmospheric and Oceanic Technology 18, 1540 (2001).
  * Cowherd  _et al._ [2021] M. Cowherd, G. Egan, S. Monismith, and O. Fringer, Phase-Resolved Wave Boundary Layer Dynamics in a Shallow Estuary, Geophysical Research Letters 48, 10.1029/2020GL092251 (2021).
  * Daubechies  _et al._ [2011] I. Daubechies, J. Lu, and H.-T. Wu, Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool, Applied and Computational Harmonic Analysis 30, 243 (2011).
  * Thakur  _et al._ [2013] G. Thakur, E. Brevdo, N. S. Fučkar, and H.-T. Wu, The Synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications, Signal Processing 93, 1079 (2013).
  * Lumley [2007] J. L. Lumley,  _Stochastic Tools in Turbulence_ (Courier Corporation, 2007).
  * Sirovich [1987a] L. Sirovich, Turbulence and the Dynamics of Coherent Structures Part I: Coherent Structures, Quarterly of Applied Mathematics 45, 561 (1987a), 43637457 .
  * Sirovich [1987b] L. Sirovich, Turbulence and the Dynamics of Coherent Structures Part II: Symmetries and Transformations, Quarterly of Applied Mathematics 45, 573 (1987b), 43637458 .
  * Sirovich [1987c] L. Sirovich, Turbulence and the Dynamics of Coherent Structures Part III: Dynamics and Scaling, Quarterly of Applied Mathematics 45, 583 (1987c), 43637459 .
  * Aubry  _et al._ [1988] N. Aubry, P. Holmes, J. L. Lumley, and E. Stone, The dynamics of coherent structures in the wall region of a turbulent boundary layer, Journal of Fluid Mechanics 192, 115 (1988).
  * Scherl  _et al._ [2020] I. Scherl, B. Strom, J. K. Shang, O. Williams, B. L. Polagye, and S. L. Brunton, Robust principal component analysis for modal decomposition of corrupt fluid flows, Physical Review Fluids 5, 054401 (2020).
  * Tu  _et al._ [2014] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, On Dynamic Mode Decomposition: Theory and Applications, Journal of Computational Dynamics 1, 391 (2014), arxiv:1312.0041 [physics] .
  * Kutz  _et al._ [2016] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor,  _Dynamic Mode Decomposition_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016).
  * Filho and Lopes dos Santos [2019] E. V. Filho and P. Lopes dos Santos, A Dynamic Mode Decomposition Approach With Hankel Blocks to Forecast Multi-Channel Temporal Series, IEEE Control Systems Letters 3, 739 (2019).
  * Fujii  _et al._ [2019] K. Fujii, N. Takeishi, B. Kibushi, M. Kouzaki, and Y. Kawahara, Data-driven spectral analysis for coordinative structures in periodic human locomotion, Scientific Reports 9, 16755 (2019).
  * Lydon  _et al._ [2023] B. Lydon, B. Polagye, and S. Brunton, Data-driven modeling of an oscillating surge wave energy converter using dynamic mode decomposition (2023), arxiv:2306.03861 [physics] .
  * Takens [1981] F. Takens, Detecting strange attractors in turbulence, in  _Dynamical Systems and Turbulence, Warwick 1980_, Lecture Notes in Mathematics, edited by D. Rand and L.-S. Young (Springer, Berlin, Heidelberg, 1981) pp. 366–381.
  * Reimers and Fogaren [2021] C. E. Reimers and K. E. Fogaren, Bottom Boundary Layer Oxygen Fluxes During Winter on the Oregon Shelf, Journal of Geophysical Research: Oceans 126, e2020JC016828 (2021).
  * Castro and Vanderwel [2021] I. P. Castro and C. Vanderwel,  _Turbulent Flows: An Introduction_ (IOP Publishing, 2021).
  * Hasselmann  _et al._ [1973] K. Hasselmann, T. Barnett, E. Bouws, H. Carlson, D. Cartwright, K. Enke, J. Ewing, H. Gienapp, D. Hasselmann, P. Kruseman, A. Meerburg, P. Müller, D. Olbers, K. Richter, W. Sell, and H. Walden, Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP) (1973).
  * Farazmand and Sapsis [2019] M. Farazmand and T. Sapsis, Surface Waves Enhance Particle Dispersion, Fluids 4, 55 (2019).
  * Baker and DiBenedetto [2023] L. Baker and M. DiBenedetto, Large-scale particle shadow tracking and orientation measurement with collimated light, Experiments in Fluids 64, 52 (2023).
  * Goring and Nikora [2002] D. G. Goring and V. I. Nikora, Despiking Acoustic Doppler Velocimeter Data, Journal of Hydraulic Engineering 128, 117 (2002).
  * Qiao  _et al._ [2016] F. Qiao, Y. Yuan, J. Deng, D. Dai, and Z. Song, Wave–turbulence interaction-induced vertical mixing and its effects in ocean and climate models, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374, 20150201 (2016).
  * Guo and Shen [2013] X. Guo and L. Shen, Numerical study of the effect of surface waves on turbulence underneath. Part 1. Mean flow and turbulence vorticity, Journal of Fluid Mechanics 733, 558 (2013).
  * Williams  _et al._ [2015] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, A Data–Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition, Journal of Nonlinear Science 25, 1307 (2015).
  * Jovanović  _et al._ [2014] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, Sparsity-promoting dynamic mode decomposition, Physics of Fluids 26, 024103 (2014).