Impact of swell waves on atmospheric surface turbulence: wave–turbulence decomposition methods

To study turbulence properties, specifically vertical momentum fluxes under swell wave conditions, I investigate the impact of waves on the power spectrum and spectral coherence of turbulent wind across various spatial and temporal scales. I propose and apply a wave–turbulence decomposition method to split high-frequency surface wind data into distinct wind and wave components. Under the assumption of frozen turbulence, this method substitutes an empirically fitted spectrum for the observed or modelled wind spectrum within the wave-affected frequency range. I proceed to estimate time series of waves and turbulence through this decomposition technique. Using a few days of sonic anemometer wind measurements at 15 m height from 20 to 26 June 2015, the upward momentum transfer could be observed under high–steady (∼7 m s⁻¹) and decaying wind conditions. During the high and decaying winds, the atmospheric stability changes between unstable and stable conditions, blurring the wave signals due to the thermally and mechanically generated turbulence. The vertical wind spectra from selected episodes within the study period, acting as benchmarks, offer detailed insights into how waves affect energy elevation within the wave frequency band under low-wind, old-sea, and stable boundary layer conditions. These spectra also facilitate an effective performance assessment of the proposed decomposition method. Additionally, using a theoretical model derived from sonic anemometer measurements at heights of 15 and 20 m above the mean sea level, I parameterize the wave-contaminated coherence function, allowing for the synthetic generation of turbulent fluctuation spectra within the wave frequency band.

1 Introduction

Over the last several decades, a large number of laboratory studies and field experiments have shown modulation of turbulent momentum fluxes across a layer, the so-called wave boundary layer (WBL), on both sides of the air–sea interface (Chalikov, 1995; Rieder et al., 1994). The height of this sublayer, for example in the atmosphere, is approximately a few metres and on the order of the significant wave height. The atmospheric WBL is then limited from below by the air–sea interface and at the top by the atmospheric surface layer in which the Monin–Obukhov similarity theory (MOST) is applied. Stratification may have direct and indirect impacts on the wind–wave interaction within the WBL (Semedo et al., 2009). However, such effects depend on the height of the WBL with respect to the height of the dynamic sublayer, wave age, and the relative angle between wind and wave directions. Within the WBL, particularly under the influence of swell waves with low to moderate wind speeds, MOST or the logarithmic law should not be applied in order to estimate the drag coefficient or roughness length, and the wind profiles show a jet at the top of the WBL (Chalikov and Rainchik, 2011). This is because waves excite perturbations in this sublayer in addition to contributions from the buoyancy and shear productions. Although studies over the last few decades have significantly improved understanding of turbulent flows above the surface gravity waves, detailed knowledge of the WBL and its interaction with atmospheric turbulence under varying forcing conditions are key in better understanding turbulent processes and enhancing the accuracy of turbulent closure schemes used in oceanic, atmospheric, and climate models.

Fast-propagating swell waves generate wave-coherent structures in the WBL and challenge the widely accepted Monin–Obukhov scaling (Högström et al., 2013; Smedman et al., 2003; Semedo et al., 2009; Rieder and Smith, 1998). Wind velocity fluctuations are influenced by both turbulence and wave orbital motions, potentially leading to overestimated turbulence parameters due to significant wave-phase-dependent modulation of airflow. To mitigate this issue, various methods, such as phase averaging, linear transformation (Veron et al., 2008; Grare et al., 2013; Buckley and Veron, 2017), and orthogonal projection of wind onto the Hilbert space (Hristov et al., 2003; Wu et al., 2018), have been proposed and used for decomposing the turbulent velocity fields from instantaneous measurements of wind within the WBL. Many of these techniques rely on complex cross-spectra between horizontal and vertical air velocity fluctuations and on sea surface elevation to isolate the direct wave influence (Grare et al., 2013). It is worth mentioning that understanding the spatial coherence and the separation of wave-coherent structures from the wind measurements is crucial for generating turbulence data in structural load analysis models and wind turbine and farm control implementations and for simulating wind–wave interactions during swell waves (Bakhoday-Paskyabi et al., 2022).

In this study, I utilize a set of near-surface wind and wave measurements collected during the OBLEX-F1 campaign from the FINO1 offshore meteorological mast (met mast) in June 2015. Based on these data, I explore the turbulent structures within the WBL and the interactions of waves and wind stress under different atmospheric stability and sea-state conditions. I further calculate the swell-induced momentum fluxes from both available theories and high-frequency observational data. Specifically, the approach in this paper aims to unfold the following key aspects:

• identifying and assessing swell-related wind–wave imprints on atmospheric velocity and two-point coherence structures under specific swell-dominated conditions (proposing a theoretical model for representing wave-induced coherence),

• removing wave-induced peak effects from wind velocity spectra using a spectral technique and reconstructing turbulence and wave time series from the wave-affected sonic measurements, under mostly stable atmospheric conditions.

This paper is organized as follows. Section 2 will briefly explain the coherence spectrum within the wave-affected frequency band, the wave–turbulence decomposition, and the estimation of wave-induced momentum fluxes. Section 3 explains measurements of wind and waves at the FINO1 offshore met mast, and Sect. 4 describes the verification results of the suggested methods. Section 5 provides a brief discussion and summary of the work.

2 Methods

In the presence of swell waves, the total wind $\mathit{u}=(u,v,w)$ is linearly decomposed into the mean $\stackrel{\mathrm{‾}}{\mathit{u}}$, the turbulent u^′, and the wave-induced perturbation $\stackrel{\mathrm{̃}}{\mathit{u}}$ as

$$\begin{array}{}\text{(1)}& \mathit{u}=\stackrel{\mathrm{‾}}{\mathit{u}}+{\mathit{u}}^{\mathbf{\prime }}+\stackrel{\mathrm{̃}}{\mathit{u}},\end{array}$$

where u and v represent the horizontal velocity components and w represents the vertical velocity component. Based on this decomposition in this section, the WBL is studied through the wave–turbulence decomposition. Estimation of wave-induced motions $\stackrel{\mathrm{̃}}{\mathit{u}}=(\stackrel{\mathrm{̃}}{u},\stackrel{\mathrm{̃}}{v},\stackrel{\mathrm{̃}}{w})$ in Eq. (1) presents theoretical challenges, and it becomes even more intricate when dealing with measurements of wind over an undulating air–sea interface. To underscore the complexity in a realistic scenario, waves are often modelled as a continuous spectrum of monochromatic waves; e.g. refer to a finite expansion in Eq. (11). Each wave component may induce airflow perturbations travelling at its phase speed. This necessitates employing a space–time Fourier transform to differentiate components moving at the flow speed from those travelling at the wave speed (Ayet and Chapron, 2022). This multiscale wind–wave coupling, mediated by wave-coherent motions, is a responsible mechanism for variations in turbulent characteristics over the swell waves. For instance using Eq. (1), the total wind stress vector over the wavy surface is given as follows:

$$\begin{array}{}\text{(2)}& \mathit{\tau }={\mathit{\tau }}^{\mathbf{\prime }}+\stackrel{\mathrm{̃}}{\mathit{\tau }},\end{array}$$

where ${\mathit{\tau }}^{\mathbf{\prime }}=-{\mathit{\rho }}_{\mathrm{a}}(\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }},\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }})$ is the turbulent stress and $\stackrel{\mathrm{̃}}{\mathit{\tau }}=-{\mathit{\rho }}_{\mathrm{a}}(\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}},\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{v}\stackrel{\mathrm{̃}}{w}})$ denotes the wave-induced stress. Here, ρ_a indicates the air density. The total wind stress in Eq. (2) can be determined either through high-frequency measurements using the eddy covariance technique, to calculate the observed τ^′ and $\stackrel{\mathrm{̃}}{\mathit{\tau }}$, or by employing a bulk formula such as the one provided by COARE3.6 (Edson et al., 2013).

2.1 Wind–wave decomposition

To decompose wind and wave signals (or wave–turbulence decomposition), the energy spectrum of each velocity component in the inertial range at wavenumbers significantly above and below the wave band can be fitted using following the 1D Kaimal wavenumber spectrum:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{k{F}_{\mathit{\beta }\mathit{\beta }}\left(k\right)}{{\mathit{\sigma }}_{\mathit{\beta }}^{\mathrm{2}}}}={\displaystyle \frac{A(k/k_{\mathrm{0}\mathit{\beta }})}{\mathrm{1}+(k/k_{\mathrm{0}\mathit{\beta }}{)}^{\mathrm{5}/\mathrm{3}}}},\end{array}$$

where k denotes wavenumber; $\mathit{\beta }=u,v,w$; $A=\mathrm{5}\sin (\mathrm{3}\mathit{\pi }/\mathrm{5})/\left(\mathrm{6}\mathit{\pi }\right)$ is a constant; and k_0β and σ_β are two adjustable parameters describing the roll-off wavenumber (the length scales of turbulent eddies in the energy-containing subrange) and the standard deviation of β, respectively (Gerbi et al., 2009; Bakhoday-Paskyabi, 2019). Here, I perform a two-parameter least-squares fit of Eq. (3) to the observed data. This process allows us to estimate k_0β and σ_β, which characterize the variance and spatial scale, respectively, of the energy-containing eddies. The wavenumber spectrum is then converted to frequency scales by invoking Taylor's frozen turbulence hypothesis, $k=\mathit{\omega }/\stackrel{\mathrm{‾}}{u}$, where ω=2πf and $\stackrel{\mathrm{‾}}{u}$ is the mean (advection) wind speed:

$$\begin{array}{}\text{(4)}& \mathrm{d}f{E}_{\mathit{\beta }\mathit{\beta }}\left(f\right)=\mathrm{d}k{F}_{\mathit{\beta }\mathit{\beta }}\left(k\right).\end{array}$$

Here E_ββ(f) is the frequency spectrum of the wind β component, and the derivative $\mathrm{d}k/\mathrm{d}f$ is estimated by the use of the wave dispersion relation; see Appendix A. I then apply a two-parameter least-squares fitting of the model spectrum in Eq. (3) to the measured spectrum. For fitting (in log–log space to ensure equal weight is given to all parts of the model fit), the wave-affected band is first determined as $[\mathrm{0.6}{k}_{\mathrm{p}},k_{\mathrm{p}}+\mathrm{0.1}]$ (or $[\mathrm{0.6}{f}_{\mathrm{p}},f_{\mathrm{p}}+\mathrm{0.1}]$) at which k_p (f_p) denotes the peak wavenumber (frequency) measured from the recorded wave peak period, T_p. This interval is determined through a trial-and-error process using all datasets employed in this study and provides a reliable estimate of the frequency band for the entire campaign dataset too. The energy spectrum is then divided into three bands: below-wave-band (k<0.6k_p), wave-band, and above-wave-band ($k>k_{\mathrm{p}}+\mathrm{0.1}$) parts; see Fig. 1 and refer to Appendix 3 of Bakhoday-Paskyabi (2019). After discarding the wave band, the Kaimal spectrum of Eq. (3) is fitted over below-wave-band and above-wave-band wavenumbers and the wave-induced bump is replaced by the fitted curve. The wave-induced spectrum is then estimated as follows:

$$\begin{array}{}\text{(5)}& E_{\stackrel{\mathrm{̃}}{\mathit{\beta }}\stackrel{\mathrm{̃}}{\mathit{\beta }}}\left(f\right)=E_{\mathit{\beta }\mathit{\beta }}\left(f\right)-E_{{\mathit{\beta }}^{\prime }{\mathit{\beta }}^{\prime }}\left(f\right).\end{array}$$

Figure 1 Comparing the coherence of the vertical velocity component of wind turbulence using the Davenport empirical model (Term1) to the wind–wave coherence (depicted by the black line) as shown by combining Term1 and Term2 in Eq. (8) for all frequencies f≥0.6f_p. Here, we exclude the undulatory tail effects for $f\ge f_{\mathrm{p}}+\mathrm{0.1}$, as presented in Term3. Here $\stackrel{\mathrm{‾}}{u}=\mathrm{1}$ m s⁻¹, H_s=0.8 m, T_p=8 s, A₁=0.5, C₂=0.7, β₂=0.09, A₂=0.6, and A₃=0.1 for the above-the-wave frequencies (the grey-coloured area).

To estimate time series of turbulence and wave components in Eq. (5), I set $\mathit{u}\leftarrow \mathit{u}-\stackrel{\mathrm{‾}}{\mathit{u}}$ and transform Eq. (1) into the Fourier space in terms of Fourier coefficients of wind U_j, V_j, and W_j and waves ${\stackrel{\mathrm{̃}}{U}}_j$, ${\stackrel{\mathrm{̃}}{V}}_j$, and ${\stackrel{\mathrm{̃}}{W}}_j$; for example $E_{uu}=\left|U_j\right|$ and $E_{\stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{u}}=\left|{\stackrel{\mathrm{̃}}{U}}_j\right|$. These Fourier coefficients are expressed in phasor notation as follows (Bricker and Monismith, 2007):

$$\begin{array}{}\text{(6)}& {\displaystyle }{U}_j=\left|U_j\right|e^{i\mathrm{\angle }{U}_j},\text{(7)}& {\displaystyle }{\stackrel{\mathrm{̃}}{U}}_j=\left|{\stackrel{\mathrm{̃}}{U}}_j\right|e^{i\mathrm{\angle }{\stackrel{\mathrm{̃}}{U}}_j},\end{array}$$

where ∠ denotes the phase operator, $i=\sqrt{-\mathrm{1}}$ is the imaginary unit, and $|\cdot |$ represents the magnitude operator. By taking the inverse Fourier transformation of the above two-sided equations, the time series of turbulent and wave velocities are calculated (see Eqs. 10 and 11).

2.2 Wind–wave interaction: coherence and synthetic turbulence

The undulating surface of the ocean, as previously discussed, produces wave-coherent perturbations in the velocity (and pressure) fields. This has the potential to exert a dominant influence on turbulent properties within the WBL. In the case of low wind speeds and when vertical separation distances (for two-point turbulence problems) remain within typical turbulence length scales, turbulence can be considered frozen, allowing the application of Taylor's hypothesis using representations that may however deviate significantly from those in the existing literature. Moreover, the statistics of spatial structures at the microscale and the similarity of flow motions across different scales, as measured by the coherence between spatially separated data, may be influenced in the presence of a wavy surface.

Figure 2 (a) Times series of wind speed at 10 m height measured by FINO1's cup anemometer (black line) and significant wave height (red line) measured by a floating buoy operating in the close vicinity of the FINO1 platform; (b) wind (black markers) and wave (red markers) directions; and (c) the values of the stability parameter from which the Obukhov length, L, was calculated from sonic measurement at a height of 15 m above the mean sea level (a.m.s.l.) collected with sampling frequency of 25 Hz between 21 and 27 June 2015. The stability classes have been colour-coded in this figure. Vertical dashed grey lines highlight the study events of wind–wave interaction on 24 June at 13:00 UTC and 25 June at 00:00 UTC, respectively.

This study proposes a theoretical model for generating a wave-correlated wind field (or the wave-affected turbulence) in both time and frequency domains. The essence of the model lies in the representation of a coherence function for the fluctuating wind velocity by accounting for the impacts of surface wave processes. The model helps also to effectively isolate the wave contributions from the wind fluctuation signal, i.e. Eq. (9).

By relying on the Davenport representation of the coherence (Davenport, 1961), I suggest following a general representation for the vertical wave coherence formulation at the separation distance of Δz:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathit{\gamma }(\mathrm{\Delta }z,z,f)=\stackrel{\text{Term1}}{\overbrace{A_{\mathrm{1}}\exp \left(-C_{\mathrm{1}}{\displaystyle \frac{f\mathrm{\Delta }z}{\stackrel{\mathrm{‾}}{u}}}\right)}}\\ & +\stackrel{\text{Term2}}{\overbrace{A_{\mathrm{2}}\exp \left(-C_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{w}}{\displaystyle \frac{(f-f_{\mathrm{p}}{)}^{\mathrm{2}}z}{{\mathit{\beta }}_{\mathrm{2}}\left[\right(H_{\mathrm{s}}/T_{\mathrm{p}}{)}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{u}}^{\mathrm{2}}]}}\right)}}\\ & \stackrel{\text{Term3}}{\overbrace{\cos \left[{\displaystyle \frac{\mathrm{2}\mathit{\pi }fz}{{\mathit{\beta }}_{\mathrm{1}}{H}_{\mathrm{s}}/Tp}}+\mathit{\alpha }\mathit{\varphi }\right]}}+\stackrel{\text{Term4}}{\overbrace{A_{\mathrm{3}}}},\end{array}\end{array}$$

where C₁ is a dimensionless decay parameter and C₂, A₁, A₂, A₃, β₁, and β₂ are empirically determined coefficients, through fitting with observed coherence (i.e. Eq. 12), requiring a minimum of two high-frequency wind measurements at two heights with an appropriate separation distance. The variable z represents the average height over the sea surface for which coherence is calculated (for the vertical attenuation of wave effects). Oscillations of the wave component are controlled by setting a small value for the random phase shift of ϕ (i.e. α∼0.1) in Term3. Here, f is frequency; $\stackrel{\mathrm{‾}}{u}$ denotes the mean wind speed; and λ_w is the wavelength determined by utilizing the dispersion relation (Appendix A) and values of wave peak period, T_p, and the significant wave height, H_s. In Eq. (8), Term1 represents the original Davenport empirical model of the coherence without wave disturbances for the entire frequency range (the dashed red curve in Fig. 1). The separation distances (heights) are chosen to be within the maximum range where the Taylor hypothesis remains valid. Term2 together with Term1 is used in the wave-affected frequency band to represent the wind and wave coherence. The coherence in this band decreases with height z above the mean surface level and depends on H_s, T_p, wavelength, and $\stackrel{\mathrm{‾}}{u}$ (the yellow area in Fig. 1). To address the oscillatory behaviour in the tail of the coherence spectrum, observed in the data, I incorporate the sum of Term3 and Term4. It is important to note that this representation does not account for the effects of wind and wave misalignment and may require further adjustments for the vertical wave attenuation effects. Furthermore, this figure does not show the effects of Term3, and the undulating behaviour of Term3 can be seen in the shaded regions of Fig. 6a and b.

By estimating the squared coherence, γ², between vertical wind w and wave elevation η in the frequency domain, the spectral density of the wave-induced and turbulence fluctuations at the height of z is given according to Eq. (5) of Rieder et al. (1994):

$$\begin{array}{}\text{(9)}& E_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}\left(f\right)=(\mathrm{1}-{\mathit{\gamma }}^{\mathrm{2}})E_{\mathrm{ww}}\text{and}{E}_{\stackrel{\mathrm{̃}}{\mathrm{w}}\stackrel{\mathrm{̃}}{\mathrm{w}}}\left(f\right)={\mathit{\gamma }}^{\mathrm{2}}{E}_{\mathrm{ww}}.\end{array}$$

The random realizations of the wind and wave time series for vertical velocity are then estimated according to

$$\begin{array}{}\text{(10)}& {\displaystyle }{w}^{\prime }\left(t\right)={\mathrm{\Sigma }}_i\sqrt{\mathrm{2}{E}_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}\left(f\right)\mathrm{\Delta }f}\left[r_i^{\left(\mathrm{1}\right)}\cos \right(\mathit{\omega }t)+r_i^{\left(\mathrm{2}\right)}\sin (\mathit{\omega }t\left)\right],\text{(11)}& {\displaystyle }\stackrel{\mathrm{̃}}{w}\left(t\right)={\mathrm{\Sigma }}_i\sqrt{\mathrm{2}{E}_{\stackrel{\mathrm{̃}}{\mathrm{w}}\stackrel{\mathrm{̃}}{\mathrm{w}}}\left(f\right)\mathrm{\Delta }f}\left[q_i^{\left(\mathrm{1}\right)}\cos \right(\mathit{\omega }t)+q_i^{\left(\mathrm{2}\right)}\sin (\mathit{\omega }t\left)\right],\end{array}$$

where $r_i^{\left(\mathrm{1}\right)},r_i^{\left(\mathrm{2}\right)},q_i^{\left(\mathrm{1}\right)}$, and $q_i^{\left(\mathrm{2}\right)}$ are normal random numbers. In the synthetic time series, we determine the bulk wave parameters for the wind-sea condition using $H_{\mathrm{s}}=\mathrm{0.0248}|U_{\mathrm{10}}{|}^{\mathrm{2}}$ and $T_{\mathrm{p}}=\mathrm{0.729}\left|U_{\mathrm{10}}\right|$ (Carter, 1982) for simplicity and a clearer conceptual visualization; here U₁₀ indicates the wind speed at 10 m height. It is important to note that while Eq. (9) can also be employed for wave–turbulence decomposition using two-point coherence data (see the black curves in Fig. 6c and d), our focus in this study will be solely on the spectral technique presented in the next subsection because the proposed method relies exclusively on high-frequency sonic data at a single height.

Figure 3 (a) Coherence functions of the wind and wave for $A_{\mathrm{1}}/A_{\mathrm{2}}=\mathrm{7.5}$ in Eq. (8) (before normalizing the total coherence and only Term2) with peak wave frequencies of f_p=0.14 Hz and f_p=0.06 Hz corresponding to wind speeds of 10 and 20 m s⁻¹, respectively; (b) the energy spectra of the vertical velocity w estimated from the Kaimal spectra (Eq. 3) for two case studies, as shown by vertical dashed lines in Fig. 2; (c) the energy spectrum of the first case with U₁₀=10 m s⁻¹ and f_p=0.14 Hz (black line) and the corrected spectrum (blue line); and (d) the energy spectrum of the first case with U₁₀=20 m s⁻¹ and f_p=0.06 Hz (black line) and the corrected spectrum (blue line). The dashed red lines drawn in (c) and (d) show the spectral curves calculated from Eq. (3). Furthermore, the green-coloured areas in these panels represent the wave-affected frequency band with lower and upper frequencies of f_l=0.6f_p in hertz and $f_{\mathrm{u}}=f_{\mathrm{p}}+\mathrm{0.1}$ in hertz, respectively.

The observed coherence of vertical velocities is determined using the following relationship:

$$\begin{array}{}\text{(12)}& \mathit{\gamma }(z_{\mathrm{1}},z_{\mathrm{2}},f)={\displaystyle \frac{\left|{\text{Co}}_{z_{\mathrm{1}}{z}_{\mathrm{2}}}\right(f\left)\right|}{\sqrt{E_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}^{z_{\mathrm{1}}}\left(f\right)E_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}^{z_{\mathrm{2}}}\left(f\right)}}},\end{array}$$

where $E_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}^{z_{\mathrm{1}}}\left(f\right)$ and $E_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}^{z_{\mathrm{2}}}\left(f\right)$ are the power spectral density at heights z₁ and z₂, respectively. ${\text{Co}}_{z_{\mathrm{1}}{z}_{\mathrm{2}}}\left(f\right)$ denotes the two-point cross-power spectral density at heights z₁ and z₂.

2.3 Air–sea momentum flux

The wave boundary layer is a region where the non-static pressure distribution becomes apparent, with a height of impact corresponding to several significant wave heights. For medium waves (approximately from 2 to 4 m in height), the typical WBL height is a few metres, while for larger waves (more than approximately 4 m in height), it can extend up to say 20 m. The WBL interacts with the wavy air–sea interface below and merges with the Monin–Obukhov stratified boundary layer above. Within the WBL and according to Eq. (1), the momentum flux, to the leading order (compared to Eq. 2), can be alternatively written as follows:

$$\begin{array}{}\text{(13)}& \mathit{\tau }\left(z\right)={\mathit{\tau }}_{\mathit{\nu }}+{\mathit{\tau }}_f\left(z\right),\end{array}$$

where τ_ν and τ_f(z) are the viscous stress and form stress at the height of z, respectively (Donelan et al., 2012). Assuming that waves' impacts decay exponentially in the vertical, the wave form stress at z is defined by

$$\begin{array}{}\text{(14)}& {\mathit{\tau }}_f\left(z\right)={\mathit{\rho }}_{\mathrm{w}}\underset{k_{\min }}{\overset{k_{\max }}{\int }}\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}{e}^{-\mathrm{2}kz}{\mathit{\beta }}_{\mathrm{g}}(k,\mathit{\theta })\mathit{\omega }F(k,\mathit{\theta })\mathit{k}\mathrm{d}\mathit{\theta }\mathrm{d}k,\end{array}$$

where F(k,θ) is the 2D wave variance spectrum, θ denotes the wave direction, ρ_w is the water density, and k_min and k_max are the minimum and maximum wavenumbers. β_g is the wave growth rate as a function of wind speed at the height of $\mathit{\lambda }/\mathrm{2}$ (i.e. $u_{\mathit{\lambda }/\mathrm{2}}$ calculated by the logarithmic wind profile, where λ is the wavelength):

$$\begin{array}{}\text{(15)}& {\mathit{\beta }}_{\mathrm{g}}(k,\mathit{\theta })=A\mathit{\omega }{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{w}}}}{\displaystyle \frac{\left[u_{\mathit{\lambda }/\mathrm{2}}\cos \right(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{w}}\left)\right]\left|u_{\mathit{\lambda }/\mathrm{2}}\cos \right(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{w}}\left)\right|}{c^{\mathrm{2}}}},\end{array}$$

where θ_w indicates the wind direction, c is the wave phase speed, and the proportionality coefficient A is expressed by (Donelan et al., 2012)

$$A=\left\{\begin{array}{ll}\mathrm{0.11}& u_{\mathit{\lambda }/\mathrm{2}}\cos \left(\mathit{\theta }\right)>c\text{ (wind sea)}\\ \mathrm{0.01}& \mathrm{0}<u_{\mathit{\lambda }/\mathrm{2}}\cos \left(\mathit{\theta }\right)<c\text{ (fast running swell)}\\ \mathrm{0.1}& \cos \left(\mathit{\theta }\right)<\mathrm{0}\text{ (swell opposing the wind)}.\end{array}\right.$$

The form drag is calculated using the friction velocity ($u_{*}=\sqrt{{\mathit{\rho }}_{\mathrm{a}}^{-\mathrm{1}}\left|{\mathit{\tau }}^{\prime }\right|}$) and the wind speed at a reference height of z (i.e. u_z) as ${\mathrm{Cd}}_f=(u_{*}/u_z{)}^{\mathrm{2}}$. The viscous stress is expressed by

$$\begin{array}{}\text{(16)}& {\mathit{\tau }}_{\mathit{\nu }}={\mathit{\rho }}_{\mathrm{a}}{\mathrm{Cd}}_{\mathit{\nu }}^{\prime }\left|u_z\right|u_z,\end{array}$$

where the adjusted viscous drag by the form drag (sheltering effect) is given by

$$\begin{array}{}\text{(17)}& {\mathrm{Cd}}_{\mathit{\nu }}^{\prime }={\displaystyle \frac{{\mathrm{Cd}}_{\mathit{\nu }}}{\mathrm{3}}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathrm{Cd}}_{\mathit{\nu }}}{{\mathrm{Cd}}_{\mathit{\nu }}+{\mathrm{Cd}}_f}}\right).\end{array}$$

Here Cd_ν is the viscous drag coefficient.

Figure 4 (a) Time series of wave age $\mathit{\chi }=c_{\mathrm{p}}/U_{\mathrm{10}}$, in which U₁₀ denotes the wind speed at 10 m height and c_p is the wave phase speed. This figure contains the wind–wave misalignment time series θ−θ_wnd, where θ and θ_wnd denote the wave and wind directions, respectively. The dashed blue line (χ=1.2) represents the separation limit between the wind-sea condition and swell (i.e. χ≥1.2); (b) spectral energy evolution of surface wave elevation η measured from the AXYS buoy operating in very close vicinity of the FINO1 met mast; and (c) time evolution of energy spectra for the wind w component between 21 and 27 June 2015 calculated from the 15 m height sonic anemometer with a sampling frequency of 25 Hz. The dotted black markers are the wave peak frequencies calculated from the buoy-measured peak wave period T_p. Vertical dashed grey lines indicate the case studies.

3 Data

3.1 Dataset

The FINO1 measurement mast in the North Sea is located about 45 km north of the island of Borkum, Germany. Its geographical coordinates are 54°0^′53.5^′′ N, 6°35^′15.5^′′ E. The water depth at FINO1 is approximately 30 m, and the mast height is 100 m a.m.s.l. The site is exposed to an unlimited fetch area for northwesterly and northerly winds (Bakhoday-Paskyabi et al., 2018). The mast is equipped with different meteorological sensors such as cup anemometers to measure the velocity at 33, 40, 50, 60, 70, 80, 90, and 100 m and sonic anemometers with a sampling frequency of 10 Hz at 40, 60, and 80 m (Fig. 2b). During the OBLEX-F1 campaign between 2015 and 2016, two additional Gill R3-100 sonic anemometers were installed at 15 and 20 m a.m.s.l. with a sampling frequency of 25 Hz. The orientation of sonic anemometers was set at 135°, which means that the wind shadow zone extended approximately above 245°. While this paper's approach is primarily based on sonic anemometer wind data at 15 m height, I also make use of the 20 m sonic wind velocity data to improve the evaluation of the proposed theoretical coherence model in this study. This is crucial for a better understanding of the developed model's performance and applicability. Additionally, considering the significance of coherence function representation in relation to offshore wind turbine blade loads, this information is quite insightful for vertical separation distances that are not large. Furthermore, I use wave frequency spectra recorded by an AXYS wave buoy in the close vicinity of FINO1 platform during the study period.

Figure 5 (a, b) Scaled power spectra of vertical velocity fluctuations w^′ for two events drawn as vertical dashed grey lines in Fig. 2a and the corresponding corrected (decomposed) spectra (blue curves). The dashed red lines show the spectral curves calculated from Eq. (3). (c, d) Energy–time spectra of decomposed vertical velocity w^′ and the wave-induced vertical velocity $\stackrel{\mathrm{̃}}{w}$ fluctuations by the use of the suggested decomposition algorithm. The dotted black markers in (d) are the wave peak frequencies extracted from buoy measurements.

3.2 Data analysis

Figure 2 shows time series of wind, waves, and the stability parameter during the study period in June 2015. We do not exclude the effect of flow distortions by the FINO1 mast when the spectra show very clear wave-induced elevation (for wind directions between 245 and 360°); see one of the quality criteria in Appendix B. Sonic data at 15 m height scaled to 10 m height, i.e. U₁₀, using MOST show a range of wind speeds from moderate to high ($\mathrm{2}\le U_{\mathrm{10}}<\mathrm{11}$ m s⁻¹, Fig. 2a), and during the strong wind–wave interactions (after 24 June), the wind and wave directions are mostly misaligned (with direction differences larger than 100°). I study a few cases, specifically two cases representing the opposed wind and swell conditions (vertical dashed lines). To investigate coherence, I analyse two additional dates for which we have concurrent measurements of high-frequency wind at heights of 15 and 20 m, both displaying clear wave peaks in their spectra; see Fig. 6a and b.

Figure 6 (a, b) The coherence spectra of vertical velocity time series at two heights (15 and 20 m) are marked in red. In these plots, the thick black curves represent the empirical/theoretical coherence function, while the thin black and blue curves represent the fitted coherence derived from Eq. (8). The thin black curves illustrate the impact of wave-induced oscillations in the coherence tail (shaded regions). (c, d) Power spectra of the original vertical velocity fluctuations for two distinct dates under stable atmospheric conditions, depicted by the red curves. The black curves represent the spectra of turbulent components obtained directly through the coherence-based correction based on Eq. (9).

Figure 2c shows the atmospheric stability parameter $z/L$ at a height of z=15 m, where

$$L=-{\displaystyle \frac{u_{*}^{\mathrm{3}}{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{v}}}{g\mathit{\kappa }\left(\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}\right)}}$$

denotes the Obukhov length scale (in metres) and κ and g are the von Kármán constant and the gravitational acceleration, respectively. $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}$ is the flux of virtual potential temperature, and ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{v}}$ denotes the virtual potential temperature. The stability changes from stable (L>0) to unstable (L<0), and both coherence and decomposition study episodes represent, on average, stable conditions.

4 Results

In this section, I utilize three-dimensional wind measurement for calculating measured turbulence and wave-induced stresses. For clarity and brevity, the vertical wind component is used to study the performance of the wave–turbulence decomposition techniques. This further provides insights into vertical motion relevant for studying vertical momentum transport, vertical coherence, and turbulence variation with height. I initiate this section with idealized examples, based on the proposed coherence model, to establish a foundational understanding of the developed techniques under controlled parametric conditions.

Figure 7 (a) Time series of the wave-contaminated vertical velocity (red curve), the corrected w (black curve), and the wave component $\stackrel{\mathrm{̃}}{w}$ (blue curve) obtained by applying the developed decomposition method. (b) The contaminated and corrected coherence for a separation distance of 5 m using data from both 15 and 20 m sonic anemometers. The coherences are computed using Eq. (8). (c) The second-order structure function calculated from sonic anemometer measurements at 15 m height at the same time as illustrated in Fig. 6b.

Figure 8 Time series of (a) the wave–turbulence strength ratio R; (b) momentum stresses at the surface calculated using the COARE3.6 algorithm (red line), the Donelan et al. (1999) parameterization of Eq. (13) (black line), and the wave-induced form stress estimated from Eq. (14) at the surface; (c) the momentum fluxes at 15 m height estimated using decomposed turbulent wind data based on the eddy covariance technique (black markers), and the estimated wave-induced form stress at measurement height z=15 m using Eq. (14), i.e. $\left|{\mathit{\tau }}_f\right(\mathrm{15}\left)\right|$ (blue markers) – here ${\stackrel{\mathrm{̃}}{u}}_{*}=\sqrt{{\mathit{\rho }}_{\mathrm{a}}^{-\mathrm{1}}\left|\stackrel{\mathrm{̃}}{\mathit{\tau }}\right|}$; and (d) the ratio of the wave-induced turbulence intensity and the corrected turbulence intensity, ${\mathrm{TI}}_{\stackrel{\mathrm{̃}}{\mathrm{w}}}/{\mathrm{TI}}_{\mathrm{w}}$.

4.1 Theoretical coherence

The theoretical (normalized) coherence spectra for two different wind speeds are obtained using Eq. (8) in an ideal setup. Figure 3a shows that the coherence functions have peak values at f_p with an exponential decay beyond the wave band, where f_p is the wave peak frequency in hertz. The decay coefficients for both wind and wave components in Eq. (8) are set to constant values for the sake of simplicity ($C_{\mathrm{1}}=C_{\mathrm{2}}=\mathrm{1}$) and do not increase by increasing the wind speed. From the procedure given in Sect. 2.2 and coherence information drawn in Fig. 3a, I obtain realizations of turbulent winds for two cases, with wind speeds of U₁₀=10 and 20 m s⁻¹, respectively. For simplicity, the wind-sea wave bulk parameters are treated as total bulk parameters (i.e. H_s and T_p). The spectral energy distributions of synthetically generated winds are shown in Fig. 3b, which identifies the impact of waves on the turbulence in the inertial subrange, particularly across the wave frequency band $\mathrm{0.6}{f}_{\mathrm{p}}\le f\le f_{\mathrm{p}}+\mathrm{0.1}$. Figure 3c and d depict the fitting of the model spectrum derived from Eq. (3) for vertical wind. This process is carried out to eliminate the wave peak and to synthetically generate the spectra of turbulent fluctuations within the wave frequency band (indicated by the blue curves, as discussed in Sect. 2.1). The parameters presented in this figure are idealized to enhance conceptual clarity and technical demonstration, rather than reflecting real-world values. They serve as a theoretical framework to facilitate detailed illustrations of the underlying concept and are not representative of practical measurements or actual conditions. It is important to note that the actual values and more detailed information for these parameters, obtained through fitting to the observed coherence data, can be found in Appendix C and Table C1, as well as Fig. 6a and b.

4.2 Measured wind–wave spectra

Figure 5a shows the time variation in the wave age, χ, during the study period covering both the mixed wind sea (i.e. χ<1.2) and the swell waves (i.e. χ≥1.2). Wind and waves are obviously aligned during the wind-sea conditions. For the conditions where the swell waves are dominant (see Fig. 4b), wind and waves are mainly misaligned with a difference of approximately more than 100°. This is particularly the case after 25 June, when the wind–wave misalignment shows an oscillation-like behaviour. The atmosphere experiences mainly stable and also to a lesser extent unstable conditions during this period (see Fig. 2c). The time evolution of the power spectral density of the vertical wind speed (i.e. the w component) is presented in Fig. 4c. This figure has been overlaid by the time series of the wave peak frequency, f_p. There is a good agreement between the buoy-measured f_p and the spectral peak of the measure vertical component wind at the wave frequency band, consistent with an increase in the values of wave age (i.e. values of χ greater than 1.2). The agreement is more pronounced under stable atmospheric conditions, with somewhat weaker agreement observed under unstable conditions, especially after 25 June.

To further investigate the ability of the suggested method in splitting the fluctuations in wave and turbulence, I use 30 min time series of sonic anemometer data at 15 m height for two study events, depicted in Fig. 5. These events correspond to strong swell–wind interaction, characterized by low wind, large values of χ, and spectrally enhanced energy within the wave-affected frequency band (i.e. around f_p). The spectra of the corrected vertical velocity fluctuations w^′ and the vertical wave orbital velocity $\stackrel{\mathrm{̃}}{w}$ are shown in Fig. 5c and d. It is observed how the decomposition detaches the wave-induced energy elevation from the inertial subrange of w. To assess the effectiveness and performance of the proposed decomposition technique, we compare the power spectra of the vertical velocity for three decomposition methods that depend exclusively on the time series of high-frequency wind at a single height. This comparison is presented in Fig. D1 (see Appendix D).

4.3 Spectral and coherence analyses

Figure 6 displays the coherence spectra of 15 min vertical wind data at heights of 15 and 20 m, with a separation distance of 5 m. The two cases depicted in this figure were measured under low-wind and stable atmospheric conditions. During the first period, the significant wave height was approximately 0.8 m, and the peak wave period was around 7 s (Fig. 6a and c). Consequently, this 15 min segment encompassed roughly 128 swell periods, providing a sufficiently robust estimate of coherence. In this scenario, we can clearly observe the presence of a wave-induced elevation. Figure 6c shows the power spectrum of vertical velocity fluctuations. Initially, the spectrum is primarily influenced by surface waves at frequencies within $\mathrm{0.6}{f}_{\mathrm{p}}-f_{\mathrm{p}}+\mathrm{0.1}$ Hz, corresponding to a wave peak period of 7 s, with a notable peak in this range. At higher frequencies, the spectrum closely followed the Kolmogorov $-\mathrm{5}/\mathrm{3}$ curve up to 5 Hz (not shown). Using the coherence-based method described in Eq. (9), it is possible to separate the total velocity fluctuations into wave and turbulent components, as depicted by black curves. In Fig. 6a and b, it is evident that scales at frequencies greater than f_p+0.1 Hz exhibit a sustained coherence value of approximately 0.1. This could be attributed to data characteristics and quality or to potential effects of residual flow distortion during low-wind conditions.

Figure 7a displays time series of original (red curve), corrected (black curve), and wave-induced (blue curve) vertical velocity during a segment on 26 June at 04:00 UTC (refer to Fig. 6b). The black and blue curves illustrate the corrected and estimated wave time series, respectively, extracted directly using the decomposition method developed in this paper (see Sect. 2.1). To examine the impact of corrections on turbulence statistics, the effects of corrections on the coherence spectrum and structure function are illustrated in Fig. 7b and c, respectively. Utilizing the decomposition method at heights of 15 and 20 m to compute the corrected coherence results in a significant reduction in coherence within the wave-affected frequency band, although some small wave-correlated points still remain in this frequency range (i.e. green area in Fig. 7b). Furthermore, we observe a more significant decline in coherence at frequencies above the wave band.

Starting from the time series of vertical wind velocities w(t), I define the second-order structure function as follows:

$$\begin{array}{}\text{(18)}& S^{\mathrm{2}}\left(t\right)=〈\left|w\right(t+{\mathit{\tau }}_{\mathrm{inc}})-w(t){|}^{\mathrm{2}}〉,\end{array}$$

where τ_inc indicates a time lag (or an incremental time step) and 〈⋅〉 denotes an average over all time lags. By having knowledge of S²(t), it is possible to transform this into a function in space S²(r) using the parameter $r={\mathit{\tau }}_{\mathrm{inc}}\stackrel{\mathrm{‾}}{u}$, where $\stackrel{\mathrm{‾}}{u}$ represents the mean wind speed averaged over the entire time series, based on Taylor's hypothesis. The values of r vary between the Taylor length scale and the integral length scale. In Fig. 7c, the non-corrected structure function exhibits wave-induced oscillations that gradually dampen at larger spatial scales. The shape of S²(r) indicates that wave orbital contamination influences the slope of the second-order structure function for ranges before the onset of the oscillating tail. The application of decomposition to generate corrected time series in computing the structure function eliminates not only these wave-induced oscillations but also the slope enhancement from the second-order structure function in this figure.

4.4 Momentum flux estimation

Regarding the calculation of momentum flux, splitting between wave and wind fluctuations is not reasonable if there is no obvious footprint of waves in the measured velocity spectra at 15 m height. As a result, the decomposition algorithm is applied if the ratio of the energy variances in the wave frequencies and model spectrum, R, is larger than 1:

$$\begin{array}{}\text{(19)}& R={\displaystyle \frac{{\int }_{f_{\mathrm{l}}}^{f_{\mathrm{u}}}{E}_{\mathrm{ww}}^{\mathrm{o}}\left(f\right)\mathrm{d}f}{{\int }_{f_{\mathrm{l}}}^{f_{\mathrm{u}}}{E}_{\mathrm{ww}}^{\mathrm{m}}\left(f\right)\mathrm{d}f}}>\mathrm{1},\end{array}$$

where $E_{\mathrm{ww}}^{\mathrm{o}}$ and $E_{\mathrm{ww}}^{\mathrm{m}}$ are the sonic-based energy spectrum and the model spectrum of the vertical wind speed given by Eq. (3), respectively. Figure 8a shows a measure to assess the strength of wavy structures in the observed velocity spectra. I apply the decomposed turbulence time series when $R>\stackrel{\mathrm{‾}}{R}$ (i.e. when there exists a well-pronounced energy elevation around f_p). Here, $\stackrel{\mathrm{‾}}{R}$ represents the average value of R. In Fig. 8b, I compare the total wind stress at the surface obtained from Eq. (13) (black curve), with the bulk estimation (red curve) derived from the COARE3.6 algorithm (see Eq. 2). The two stress estimates align consistently when $R<\stackrel{\mathrm{‾}}{R}$. Additionally, it is noteworthy that the wave-induced form stress at the surface, shown by blue markers, undergoes transitions from positive to negative for swells moving opposite to the wind direction (i.e. when $R>\stackrel{\mathrm{‾}}{R}$). Figure 8c shows that the estimated form stress at z=15 m according to Eq. (14) is approximately in acceptable agreement with the measured $\left|\stackrel{\mathrm{̃}}{\mathit{\tau }}\right|$ from the sonic data, using the eddy covariance technique according to Eq. (2). In the estimated form stress at 15 m height, the dimensionless function for the vertical decay, i.e. Eq. (14), plays a significant role in the vertical distribution of the wave-induced momentum flux. Moreover, the decomposition method outlined in Sect. 2.1 has been applied to all three components of measured wind velocities to estimate the observed wave-induced stress $\left|\stackrel{\mathrm{̃}}{\mathit{\tau }}\right|$ represented in Eq. (2). Figure 8d illustrates that the ratio of wave-induced (turbulence) intensity (the standard deviation of $\stackrel{\mathrm{̃}}{w}$ over the mean wind speed) to corrected turbulence intensity (the standard deviation of the corrected w over the mean wind speed) is most pronounced when wave elevations are clearly visible around the peak frequency f_p (in agreement with Fig. 8a).

5 Conclusions

I have suggested a wind–wave decomposition algorithm for the turbulent airflow over the ocean in the presence of swell waves based on high-frequency data recorded from a sonic anemometer at 15 m a.m.s.l. The wave–turbulence decomposition method proposed in this study possesses some key characteristics: (1) it relies solely on sonic wind velocity data, eliminating the need for simultaneous high-frequency wave measurements in the decomposition process. It assumes turbulence field stability during transformation into wavenumber space and disregards velocity fluctuations within the wave band. (2) The method uniquely adopts a statistical approach, employing a turbulence spectral model to effectively bridge the gap between high- and low-frequency sections in the observed spectra. This allows for the estimation of the variance attributed to turbulent velocity fluctuations within the wave frequency band, relying solely on the energy spectrum of the corresponding wind component. (3) Notably, this method provides wind-corrected and wave time series, crucial data for structural analysis that, to the best of my knowledge, are not available through many known methods.

Furthermore, I have introduced a theoretical formulation for coherence that considers both wind fluctuations and the influence of swell waves. This model not only has been utilized to create idealized wind time series under swell conditions but also can be adjusted using observational data. In this case, I additionally used sonic data at 20 m (along with sonic data at 15 m) to estimate observed coherence and determine the fitting coefficients presented in Eq. (8). Furthermore, I quantified the wave-induced stresses by assuming both a certain growth rate and a vertical decay function. A few real cases during the opposed wind–wave conditions were selected to demonstrate physical aspects of wind and waves regimes. Under the wind-sea conditions, it was found that the wind stresses were quite similar to those observed under atmospheric neutral conditions. This is because the height of the wave boundary layer was below the height of measurement, and waves could contribute to a very small fraction of the total wind stress.

Although I have introduced a theoretical coherence function with a small number of fitting parameters, there is potential for improving the model and its efficiency by reducing the number of adjustable parameters and the model mathematical representation. Further details on this, as well as the reconstruction of the turbulence box for assessing the structural loads of offshore wind turbines over swell waves, are discussed in a separate independent study.

Appendix A: Calculation of $\mathrm{d}f/\mathrm{d}k$ based on the linear dispersion relation

The wavenumber and frequency spectra are interrelated through the dispersion relation (I assume here the linear dispersion relation):

$${\mathit{\omega }}^{\mathrm{2}}=gk\tanh \left(kd\right),$$

where ω is the angular frequency and d denotes the water depth. The spectral variance, whether expressed in frequency or wavenumber spectra, can be determined accordingly:

$${\mathit{\sigma }}^{\mathrm{2}}=\int F\left(k\right)\mathrm{d}k=\int E\left(f\right)\mathrm{d}f,$$

where E(f) and F(k) are the frequency and wavenumber spectra, respectively. Assuming linear dispersion, we can estimate $\mathrm{d}k/\mathrm{d}f$, essential for the transformation between these two spectra, as

$${\displaystyle \frac{\mathrm{d}f}{\mathrm{d}k}}={\displaystyle \frac{g}{\mathrm{4}\mathit{\pi }\mathit{\omega }}}\left[\tanh \left(kd\right)+{\text{sech}}^{\mathrm{2}}\left(kd\right)kd\right].$$

Appendix B: Isotropic ratio

Here, I examine one of the quality criteria, which is the presence of the isotropic portion of the spectrum in the inertial subrange (outside the wave-contaminated frequency range). The cross-coherences, at zero separation, are not computed to assess isotropy, as not all of these cross-coherences may be exactly zero. To be more specific, I define the isotropy criteria as follows:

$${\displaystyle \frac{E_{u^{\prime }{u}^{\prime }}}{E_{{\mathrm{w}}^{\prime }{\mathrm{w}}^{\prime }}}}\sim {\displaystyle \frac{E_{u^{\prime }{u}^{\prime }}}{E_{v^{\prime }{v}^{\prime }}}}\sim {\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}.$$

It is worth noting that these ratios over the isotropic bandwidth converge to $\mathrm{3}/\mathrm{4}$ (Fig. B1, red-coloured areas). In Fig. B1, the power spectra and spectral ratio of the non-corrected time series are presented to assess the statistical isotropy of the three velocity components. It is evident that the non-corrected ratios in Fig. B1b approach a value of $\mathrm{3}/\mathrm{4}$ for frequencies larger than the wave-affected band. In some cases with wave peaks in the velocity spectrum during low-wind conditions, I notice increased anisotropic features (not shown), likely due to higher wave contamination beyond the wave frequency band (Bakhoday-Paskyabi, 2019) and partial effects of flow distortion.

Figure B1 (a) Three-dimensional velocity spectra f⋅E_ii(f) for $i=u,v,$ and w showing the spectral regions of the wave-affected peaks and the inertial subrange; (b) velocity spectra ratios with horizontal dashed yellow and black lines indicating the value of $\mathrm{4}/\mathrm{3}$ and $\mathrm{3}/\mathrm{4}$, respectively.

Appendix C: Theoretical coherence function in Eq. (8)

To obtain estimations for the fitting parameters in the theoretical coherence model outlined in Eq. (8), I provide the results of the fitting process in Table C1 for the scenario illustrated in Fig. 6b. The fitting model comprises only Term1 and Term2. To generate multiple ensembles using the derived parameters, I perturbed four critical parameters, including H_s, T_p, A₁, and A₂ (Fig. C1).

Figure C1 Theoretical coherences are generated based on the fundamental parameters listed in Table C1. To create 20 ensembles, I perturb H_s within the range of [0.5,1.5], T_p within [5,12], A₁ within [0.1,.5], and A₂ within [0.02,0.1]. It is important to note that in the theoretical model defined in Eq. (8), we only use Term1 and Term2 to produce this figure.

Table C1 Theoretical model parameters are determined by applying the least-squares estimation method to the observed coherence illustrated in Fig. 6b. The total wave bulk parameters used to generate this figure are $H_{\mathrm{s}}=H_{\mathrm{s}}^{\mathrm{sw}}+H_{\mathrm{s}}^{\mathrm{ws}}$ and $T_{\mathrm{p}}=T_{\mathrm{p}}^{\mathrm{sw}}+T_{\mathrm{p}}^{\mathrm{ws}}$, where $H_{\mathrm{s}}^{\mathrm{sw}}$ [m] and $T_{\mathrm{p}}^{\mathrm{sw}}$ are the significant wave height and peak period of swell waves and $H_{\mathrm{s}}^{\mathrm{ws}}$ and $T_{\mathrm{p}}^{\mathrm{ws}}$ are the wind-sea bulk parameters estimated using Carter's formula (Sect. 2.2).

Appendix D: Wind–wave decomposition methods

Figure D1 (a) Power spectra of the vertical velocity fluctuation w^′ (illustrated in Fig. 6d) and its decomposed intrinsic mode functions (IMFs; for n=10); (b) comparison of different wave–turbulence decomposition methods: stop-band technique (blue curve), spectral linear transformation (red line), the empirical mode decomposition (EMD) method after removing IMF4 (green curve), and the spectral decomposition method proposed in this paper (yellow curve). The coloured areas in this figure show the selected wave frequency band.

The choice of decomposition method in this paper is based on specific considerations related to the research objectives and the nature of the data. The approach solely utilizes sonic wind velocity data at a single height, omitting the need for concurrent high-frequency wave measurements in the decomposition process. In this appendix, I compare three existing wave–turbulence decomposition methods that rely only on the high-frequency measurement of wind at a single height. These methods are the stop-band (SB) method, spectral linear transformation, and an approach based on empirical mode decomposition (EMD). The first two are established filtering techniques, while the latter is a novel signal-processing method.

The stop-band filter eliminates frequency bands in which waves dominate in the sonic velocity time series (blue curve in Fig. D1b). In this method, we employ a second-order Butterworth filter. The stop-band frequency thresholds, which consist of lower and upper cutoff frequencies, are specifically designed to attenuate frequencies within the wave-affected band. The spectrum produced by the stop-band filter method exhibits a significant drop in energy at wave frequencies, resulting in an underestimation of turbulent energy. The second method, adapted from Rieder and Smith (1998), involves identifying the wave frequency band in the velocity spectrum. Subsequently, we remove the wave-correlated portion of the spectrum and replace the removed frequencies with a line, as illustrated in Fig. D1b (red line). While the effectiveness of these two methods relies heavily on the accurate choice of the wave frequency band, the line-fit method successfully eliminates orbital velocities from the analysed signal.

The third approach, known as the EMD method, decomposes the observed signal in the time domain into multiple intrinsic mode functions (IMFs). Each IMF, as a stationary stochastic process, characterizes a narrowband frequency–amplitude modulation typically associated with a specific physical process. I have adapted this technique from Qiao et al. (2016), who utilized the EMD method to extract wave signals from field velocity observations. In this method, the vertical velocity fluctuations, which include wave orbital velocities, are decomposed into n IMFs and a residual fluctuation $w_{\text{residual}}^{\prime }$ as follows:

$$\begin{array}{}\text{(D1)}& w^{\prime }=w_{\text{ IMF1}}^{\prime }+w_{\text{ IMF2}}^{\prime }+\mathrm{\cdots }+w_{\text{ IMFn}}^{\prime }+w_{\text{ residual}}^{\prime }.\end{array}$$

The IMFs can be categorized into two groups: wave-correlated components and non-wave-correlated components. In the vertical velocity time series, the peak frequencies of IMF4 and IMF5 in Fig. D1a fall within the wave frequency range. Thus, we specifically identify IMF4 as the wave component that is excluded from the sum in Eq. (D1) to remove significant wave contamination. IMF6–IMF10, on the other hand, represent high-frequency signals like turbulence.