1. Introduction
A weather forecasting system is composed of many subsystems: the atmospheric model, the wave model, the land-surface model and the ocean model. The coupling between these subsystems mainly relies on transfer coefficients that regulate the fluxes of momentum and heat between them (Charnock Reference Charnock1955). Of particular relevance is their behaviour at the ocean–atmosphere interface. The reason is given by the fact that the oceans cover 71 % of the Earth’s surface, making such interface phenomena very relevant for the atmospheric dynamics (Sullivan & McWilliams Reference Sullivan and McWilliams2010). The exchanges across the interface with oceans abruptly affect the atmosphere and the understanding of the driving mechanisms would be certainly fundamental to improving the reliability of the transfer coefficients and, hence, the weather and climate prediction capabilities. Several theoretical models were proposed to explain the wind–wave phenomena (Miles Reference Miles1957; Phillips Reference Phillips1957) and experimental data have provided valuable insights into the physical process (Csanady Reference Csanady2001; Sullivan & McWilliams Reference Sullivan and McWilliams2010; Buckley & Veron Reference Buckley and Veron2019). However, after decades of research efforts, our understanding of the wind–wave problem is still recognised as extremely elusive. The reason is the multiscale nature of the phenomena involved. Indeed, the large scales of the turbulent wind are affected by the smaller scales of water waves which in turn are influenced by the wind itself, thus forming a complex multiscale mechanism. Most of the previous attempts have focused on the position-dependent behaviour of the phenomena. However, the description in physical space alone is insufficient to capture the real dynamics. With this aim, a more general theoretical framework has been developed, the so-called generalised Kolmogorov equation (Hill Reference Hill2002). Such an equation provides a complete and formally precise description of both the scale- and position-dependent dynamics and will be used in the present work to unveil fundamental phenomena at the basis of the wind–wave interaction problem.
1.1. Approaches to the wind–wave problem
In weather prediction models, the momentum flux in the surface layer is represented as
\begin{equation} \tau =\rho _a u_{\tau _a}^2 ,\end{equation}
where
$\rho _a$
is the air density and
$u_{\tau _a}$
is the wind friction velocity. By using the law of the wall
\begin{equation} \frac {U(y)}{u_{\tau _a}}=\frac {1}{\kappa }\log \left ( \frac {y}{k_0} \right ) \; , \end{equation}
the value of
$u_{\tau _a}$
is derived from the neutral wind velocity (corrected for non-uniform temperature) at a reference height of 10 m,
$U_{10N} = U(y_{10})$
. Here,
$U(y)$
is the mean velocity profile,
$y$
is the height above the surface,
$\kappa$
is the von Kármán constant and
$k_0$
is the roughness length. The effect of stratification of the atmosphere is included by using the Monin–Obukhov similarity theory (Monin & Obukhov Reference Monin and Obukhov1954) which typically holds in a surface layer height of approximately 50 m. By using (1.2) we can rewrite (1.1) as
\begin{equation} \tau =\rho _a C_{DN} U_{10N}^2 ,\end{equation}
with
$C_{DN}=[\kappa /\ln (y_{10}/k_0)]^2$
being the neutral drag coefficient. The surface roughness
$k_0$
over land is considered constant or climatologically slowly varying with time following vegetative changes. Instead, over the ocean, the surface roughness is strongly coupled with the wind state conditions (Janssen Reference Janssen1989; Komen et al. Reference Komen, Cavaleri, Donelan, Hasselmann, Hasselmann and Janssen1994; Belcher & Hunt Reference Belcher and Hunt1998).
The characterisation of momentum transfer between wind and waves reveals three main factors that influence
$C_{DN}$
: the wind regime (low wind speed
$U_{10} \leqslant 4\,\mathrm{m\, s^{-1}}$
, high speed
$U_{10} \geqslant 20\,\mathrm{m\, s^{-1}}$
and moderate wind in between); the wind–wave coupling as expressed by the wave age parameter,
$\beta = c_p/u_{\tau _a}$
, with
$c_p$
the phase speed of waves at their peak frequency (
$\beta \approx 10$
for young waves,
$\beta \approx 30$
for old waves and fully developed sea conditions in between); and the alignment between wind and waves. Most of the uncertainties in the drag coefficient are found at low and high wind speeds. In low wind speed regimes a reduction of drag is often observed and this is often attributed to the presence of swell waves which pump momentum upward (wave-driven winds) (Hanley, Belcher & Sullivan Reference Hanley, Belcher and Sullivan2010; Edson et al. Reference Edson, Jampana, Weller, Bigorre, Plueddemann, Fairall, Miller, Mahrt, Vickers and Hersbach2013; Jiang et al. Reference Jiang, Sullivan, Wang, Doyle and Vincent2016). Also the misalignment between swell waves and the wind could explain the large scatter observed in swell-dominated conditions (Patton et al. Reference Patton, Sullivan, Kosović, Dudhia, Mahrt, Žagar and Marić2019; Cimarelli, Romoli & Stalio Reference Cimarelli, Romoli and Stalio2023). At high wind speed a saturation of the drag coefficient with increasing
$U_{10}$
is observed. This behaviour has been attributed to various mechanisms: sea-spray generation (Troitskaya et al. Reference Troitskaya, Kandaurov, Ermakova, Kozlov, Sergeev and Zilitinkevich2018), wave breaking, airflow separation (Donelan et al. Reference Donelan, Haus, Reul, Plant, Stiassnie, Graber, Brown and Saltzman2004; Scapin et al. Reference Scapin, Wu, Farrar, Chapron, Popinet and Deike2025) and local wind–wave misalignment (Zhou et al. Reference Zhou, Hara, Ginis, D’Asaro, Hsu and Reichl2022; Manzella, Hara & Sullivan Reference Manzella, Hara and Sullivan2024).
Overall, the complex interaction between wind and waves and the many flow properties involved in this problem (e.g. wave height, swell waves, wind fetch length, water depth) render field observations difficult to rationalise and a consensus among the different formulations proposed for the relevant scaling parameters in the wave boundary layer remains elusive (Edson et al. Reference Edson, Jampana, Weller, Bigorre, Plueddemann, Fairall, Miller, Mahrt, Vickers and Hersbach2013; Curcic & Haus Reference Curcic and Haus2020). In this context, numerical simulations have shown promising results, especially for low to moderate wind speed where the computational cost remains affordable (Li & Shen Reference Li and Shen2022; Wu, Popinet & Deike Reference Wu, Popinet and Deike2022; Cimarelli et al. Reference Cimarelli, Romoli and Stalio2023). The main difficulty is due to the numerical treatment of the high jump of fluid properties occurring at the air–water interface. To avoid such a problem the wind field is often decoupled from the wave evolution, which is prescribed by considering idealised (monochromatic) stationary or travelling waves (Sullivan, McWilliams & Moeng Reference Sullivan, McWilliams and Moeng2000; Kihara et al. Reference Kihara, Hanazaki, Mizuya and Ueda2007; Yang & Shen Reference Yang and Shen2010; Druzhinin, Troitskaya & Zilitinkevich Reference Druzhinin, Troitskaya and Zilitinkevich2012). On the contrary, a simple coupling between wind and waves can be added by considering the air pressure distribution over the water surface and solving the nonlinear waves using the Zakharov formulation (Zakharov Reference Zakharov1968). Such procedure has been used both in Direct Numerical Simulation (DNS) (Liu et al. Reference Liu, Yang, Guo and Shen2010; Wang et al. Reference Wang, Zhang, Hao, Huang, Shen, Xu and Zhang2020) and large-eddy equations (LES) (Yang, Meneveau & Shen Reference Yang, Meneveau and Shen2013; Hao & Shen Reference Hao and Shen2019). However, assuming the water motion as a potential flow requires the effects of viscosity, turbulence, surface tension and wave breaking to be negligibly small. Fully coupled numerical simulations with a direct solution of the air–water interface can overcome many problems of the above-mentioned approaches. However, only few attempts have been performed by Yang, Deng & Shen (Reference Yang, Deng and Shen2018) for studying breaking waves, by Wu et al. (Reference Wu, Popinet and Deike2022) for the early-stage transient growth of water waves and by Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023) for studying interaction phenomena at wind–wave equilibrium. In the following we focus on Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023), since this work represents a further deeper analysis of their simulation.
1.2. Direct numerical simulation of turbulent wind over water waves at equilibrium
To the best of our knowledge, the simulation performed in Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023) is the only example of a two-phase wind–wave interaction problem where the turbulent wind and the water waves reach an equilibrium. The flow case consists of a two-phase open channel where an air layer interacts with an underlying water column. The flow solution is obtained by solving the continuity and momentum equations
\begin{equation} \left \{ \begin{split} &\frac {\partial u_i^*}{\partial x_i} = 0 ,\\ &\frac {\partial \rho u_i^*}{\partial t} + \frac {\partial \rho u_i^* u_j^*}{\partial x_j} = -\frac {\partial p^*}{\partial x_i} + \frac {\partial \tau _{\textit{ij}}}{\partial x_j} + f_{\sigma _i} + \rho g_i ,\end{split} \right . \end{equation}
where
$u_i^*$
is the velocity field,
$p$
is the pressure field,
$\tau _{\textit{ij}} = 2 \mu S_{\textit{ij}}$
is the viscous stress tensor with
$S_{\textit{ij}}$
the strain rate tensor and
$\mu$
the dynamic viscosity,
$f_{\sigma _i}$
is the surface tension,
$g_i$
the acceleration due to gravity and
$\rho$
the density. In the following the indices
$i = 1, 2, 3$
correspond to the streamwise, vertical and spanwise
$(x,y,z)$
-directions and
$(u^*,v^*,w^*)$
-velocities and
$y=0$
is the average vertical position of the interface between the two fluids. The different physical properties of the two immiscible fluids are taken into account by solving a transport equation for the volume fraction (Hirt & Nichols Reference Hirt and Nichols1981). The properties of the two fluids are
$\rho _a = 1$
,
$\rho _w = 1 {\times}10^3$
,
$\mu _a = 1.48 {\times}10^{-5}$
and
$\mu _w = 1 {\times}10^{-3}$
where the subscripts
$w$
and
$a$
are used to denote quantities computed for water and air, respectively. Finally, a surface tension coefficient
$\sigma = 0.07$
is adopted. More details about the numerical setting, spatial resolution and domain dimensions can be found in Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023).
The evolution of the flow is studied starting from an initial condition where the two fluids are at rest and the air–water interface is flat. The flow is driven by a streamwise pressure gradient, thus allowing us to study the self-development of the wind–wave problem in a fully coupled first-principles framework. Hence, no wave pattern is imposed. After an initial transient of wave growth and turbulent wind evolution, the flow system reaches a statistically steady state where water waves are in equilibrium with the turbulent wind. Thanks to the use of periodic boundary conditions, such a statistical equilibrium state is reached for formally infinite fetch length. The obtained statistical equilibrium is characterised by a turbulent wind at friction Reynolds Number
${Re}_{\tau _a} \approx 320$
on top of an almost quiescent and laminar water layer at
${Re}_{\tau _w} \approx 100$
. A drag-reducing flow physics has been unveiled (Cimarelli et al. Reference Cimarelli, Romoli and Stalio2023). This evidence and the related flow statistics represent the starting point of the present work and are resumed in the next section for the sake of completeness. Hereafter, flow variables are made non-dimensional by using the wind friction length
$l_{\tau _a}=\nu _a/u_{\tau _a}$
and the wind friction velocity
$u_{\tau _a}$
. The resulting non-dimensional variables will be denoted with a superscript
$+$
. The value of
$u_{\tau _a}$
is evaluated at the location of the maximum of the mean velocity gradient, in analogy with wall turbulence. This location is found to be slightly above the water–wave pattern (
$y_0^+=0.37$
) with a corresponding mean velocity of
$U_0^+=0.75$
.
To better understand the effects the wind–wave interaction mechanisms have on the turbulent wind, we will compare the present wind–wave solution with that of a turbulent channel flow at a very similar friction Reynolds number,
${Re}_\tau = 300$
. The direct numerical simulation data of the turbulent channel have been already used in Cimarelli & De Angelis (Reference Cimarelli and De Angelis2011, Reference Cimarelli and De Angelis2012), to which the reader is referred to for the details about the numerical settings.
1.3. Wave-induced Stokes sublayer
The direct numerical simulation reveals that the flow reaches a statistical equilibrium consisting of a turbulent wind over an oblique pattern of water waves propagating upstream at an angle
$\gamma = 38^\circ$
with respect to the wind direction, see figure 1(a). The measured wavelength is
$\lambda ^+ =296$
that, due to the inclined wave pattern, leads to streamwise and spanwise wavelengths
$(\lambda _x^+, \lambda _z^+) = (475,380)$
. Despite the very small elevation of the water–wave pattern,
$\delta _w^+ = 0.3$
, the turbulent wind on top is found to be significantly affected, see figure 1(b). In particular, by comparing the mean velocity profile with that realised in wall turbulence (Nagaosa & Handler Reference Nagaosa and Handler2003), a significant drag reduction is observed in the form of an upward shift of the mean velocity profile
$\Delta B^+ = 1.6$
. Both the mean velocity and turbulent profiles agree with the presence of a near-interface weakening of turbulence.
Direct numerical simulation of the wind–wave interaction problem (Cimarelli et al. Reference Cimarelli, Romoli and Stalio2023). (a) Iso-contours of the wave elevation
$\eta ^+(x,z)$
. (b) Iso-surfaces of
$\lambda _2 = 2.5$
coloured with the streamwise velocity scaled with the top boundary velocity
$U_{top}$
.
Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023) argued that at the basis of the near-interface weakening of turbulence is the presence of a wave-induced Stokes sublayer. It consists of an oscillating spanwise forcing of the near-interface turbulent wind induced by the non-sheltering behaviour of the wind itself over the inclined water–wave pattern. Due to the low elevation of water waves, the wind is able to follow the water–wave shape, thus accelerating on the windward side and decelerating on the leeward side. This behaviour leads to a pressure field that is minimum above the wave crests and maximum within the trough region. Because the water–wave pattern is skewed with respect to the mean wind direction, these pressure variations give rise to periodically distributed pressure gradients in both the streamwise and spanwise directions. The latter is responsible for the generation of an oscillating spanwise forcing, thus inducing an alternating spanwise motion that takes the form of a streamwise travelling wave whose wavelength is
$\lambda _x^+ = 475$
and phase speed is
$c_x^+ = -10$
, see figure 2. This periodically distributed spanwise forcing emulates the flow behaviour of the so-called generalised Stokes layer that is widely recognised to reduce the momentum transfer and, hence, drag in wall-bounded turbulence (Quadrio & Ricco Reference Quadrio and Ricco2011). In particular, Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023) measures a penetration length of the wave-induced Stokes sublayer of the order of
$\ell _s^+ = 2$
with an induced spanwise motion of the order of
$|w|^+ \sim 10^{-2}$
. Let us point out that the measured reduction of streamwise momentum transport toward the interface is rather surprising. From the one hand, Ghebali, Chernyshenko & Leschziner (Reference Ghebali, Chernyshenko and Leschziner2017) tried to reproduce a similar Stokes sublayer by simulating the flow over skewed wavy walls but no clear evidence of drag reduction has been achieved. On the other hand, Nakanishi, Mamori & Fukagata (Reference Nakanishi, Mamori and Fukagata2012) found that drag reduction can be achieved with streamwise travelling wall deformations only when moving downstream, contrary to what we observe with upstream travelling wind waves. Despite the similarities with these works, the wind–wave problem exhibits substantially different flow properties that result in a net reduction of momentum transfer and, hence, that deserve a deep investigation.
Direct numerical simulation of the wind–wave interaction problem (Cimarelli et al. Reference Cimarelli, Romoli and Stalio2023). Iso-lines of the wave elevation
$\eta (x,z)$
(solid and dashed lines denote positive and negative values) superimposed on iso-contours of the spanwise velocity
$w^+ (x,z)$
.
1.4. Motivation and goals of the work
The discovered wave-induced Stokes sublayer unveils very basic flow phenomena that may explain some experimental evidence also in real wind–wave problems. Indeed, the condition for the development of the Stokes sublayer is a misalignment of the wind with respect to water waves. This condition often develops in field realisations because of the presence of swell waves moving in arbitrary directions with respect to the local wind. As a consequence, the change in the momentum flux at the air–water interface associated with the degree of wind–wave misalignment could explain the large scatter of drag coefficient data in field measurements. Because of the relevance of these arguments, the present work presents a further deeper analysis of these wind–wave interaction phenomena. In particular, given the intrinsic multiscale nature of the problem, we consider the theoretical framework provided by the generalised Kolmogorov equation (Hill Reference Hill2002). It consists of the exact evolution equation for the second-order structure function that allows us to undertake a detailed analysis of the turbulent processes among different scales and positions within the flow. The aim is twofold.
-
(i) To provide a deeper understanding of the drag-reducing phenomena related to the wave-induced Stokes layer. In particular: Why is the water interface less dissipative than a rough wall? Why are larger turbulent structures generated?
-
(ii) To establish a theoretical framework for large-eddy simulation of wind–wave flows.
2. The generalised Kolmogorov equation specialised to the wind–wave problem
Most previous works have concentrated their study of the wind–wave problem on classical single-point statistics. However, the understanding of the local coupling mechanisms between water waves and turbulent wind requires the use of more sophisticated tools, able to distinguish the different scales of the problem. Indeed, the small-scale interactions near the interface between the water and the turbulent wind are known to have a strong impact on the larger scales transporting momentum and heat in the outer part of the wind boundary layer. Hence, these multiscale cascade phenomena control the surface fluxes between the atmosphere and ocean and their understanding would have a strong impact on climate science. The generalised Kolmogorov equation (Hill Reference Hill2002) is known to address the multiscale issue in a unified and exact theoretical framework. It allows us to track the spatially evolving scale-by-scale cascade mechanisms propagating to and from the water–wave surface. Recently, the formalism has been extended to account for variable density flows (Lai, Charonko & Prestridge Reference Lai, Charonko and Prestridge2018) to study the geometrical properties of the liquid–gas interface in two-phase flows (Thiesset et al. Reference Thiesset, Duret, Ménard, Dumouchel, Reveillon and Demoulin2020; Gauding et al. Reference Gauding, Thiesset, Varea and Danaila2022). It is important to notice that the theoretical framework of the generalised Kolmogorov equation is strictly related to the coarse-grained approach to the Navier–Stokes equations with strong repercussions for closures, especially for large-eddy simulation (Cimarelli & De Angelis Reference Cimarelli and De Angelis2011, Reference Cimarelli and De Angelis2014; Togni, Cimarelli & De Angelis Reference Togni, Cimarelli and De Angelis2019). The theoretical framework is extended here for the first time to the wind–wave interaction problem.
In accordance with the two main goals of the present work, we are here interested in addressing the multiscale phenomena of the turbulent wind. Hence, the analysis can be conducted in a incompressible theoretical framework with constant fluid properties. By further considering that the water waves developed at the present low Reynolds number have a small elevation,
$\delta _w^+ = 0.3$
, we have that such an incompressible framework is valid down to the very-near interface region. Accordingly, in the present work we analyse the behaviour of the generalised Kolmogorov equation down to
$y^+ =1.5$
, i.e. sufficiently far from the variable properties induced by the wavy motion of the air–water interface but well within the viscous sublayer of the turbulent wind.
In the incompressible and constant properties framework, the generalised Kolmogorov equation represents the evolution equation of the second-order moment of the two-point velocity increment, the so-called second-order structure function
$\delta q^2 = \delta u_i \delta u_i$
, where
$\delta u_i = u_i (\boldsymbol{x'}, t)-u_i (\boldsymbol{x''}, t)$
and
$u_i = u_i^* - \langle u_i^* \rangle$
is the fluctuating velocity field with
$\langle \boldsymbol{\cdot }\rangle$
the ensemble average operator. Hereafter, we will refer to the second-order structure function as the scale energy even if such an interpretation is somewhat arguable (Hamba Reference Hamba2018). As shown in the appendix of Cimarelli et al. (Reference Cimarelli, Boga, Pavan, Costa and Stalio2024), in its more general form the equation reads
\begin{align} &\frac {\partial \langle \delta q^2 \rangle }{\partial t} + \frac {\partial \langle \delta q^2 \delta u_j \rangle }{\partial r_j} + \frac {\partial \langle \delta q^2 \rangle \delta U_j }{\partial r_j} - 2\nu \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial r_j \partial r_j} + \frac {\partial \langle \delta q^2 \widetilde {u}_j \rangle }{\partial x_{c_j}} + \frac {\partial \langle \delta q^2 \rangle \widetilde {U}_j }{\partial x_{c_j}}\nonumber\\&\quad + \frac {2}{\rho } \frac {\partial \langle \delta p \delta u_i \rangle }{\partial x_{c_i}} - \frac {\nu }{2} \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial x_{c_j} \partial x_{c_j}} = - 2 \langle \delta u_i \delta u_j \rangle \widetilde {\left ( \frac {\partial U_i}{\partial x_j} \right )} - 2 \langle \delta u_i \widetilde {u}_j \rangle \, \delta \! \left ( \frac {\partial U_i}{\partial x_j} \right ) - 4 \langle \widetilde {\epsilon } \rangle, \end{align}
where
$\delta \boldsymbol{\cdot}$
and
$\widetilde {\boldsymbol{\cdot }}$
denote the two-point difference and algebraic average operators and
$U_i = \langle u_i \rangle$
is the mean velocity field. The independent variables are time
$t$
, the two-point separation vector
$\boldsymbol{r} = \boldsymbol{x'}-\boldsymbol{x''}$
and the position vector of the mid-point
$\boldsymbol{x_c} = (\boldsymbol{x'}+\boldsymbol{x''})/2$
. Finally,
$p$
is the fluctuating pressure field and
$\epsilon = \nu (\partial u_i/\partial x_j)(\partial u_i/\partial x_j )$
is the turbulent pseudo-dissipation. By considering the statistical symmetries of the turbulent wind, the generalised Kolmogorov equation, (2.1), specialises to
\begin{align} &\frac {\partial \langle \delta q^2 \delta u_j \rangle }{\partial r_j} + \frac {\partial \langle \delta q^2 \rangle \delta U }{\partial r_x} - 2\nu \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial r_j \partial r_j} + \frac {\partial \langle \delta q^2 \widetilde {v} \rangle }{\partial y_c} + \frac {2}{\rho } \frac {\partial \langle \delta p \delta v \rangle }{\partial y_c} \nonumber\\&\quad - \frac {\nu }{2} \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial y_c \partial y_c} = - 2 \langle \delta u \delta v \rangle \widetilde {\left ( \frac {\partial U}{\partial y} \right )} - 2 \langle \delta u \widetilde {v} \rangle \, \delta \! \left ( \frac {\partial U}{\partial y} \right ) - 4 \langle \widetilde {\epsilon } \rangle, \end{align}
where the dependence on the mid-point positions
$x_c$
and
$z_c$
and on time
$t$
has been dropped because of statistical homogeneity. Equation (2.2) highlights that the multi-dimensional phenomena occurring in the turbulent wind can be described in a four-dimensional space, the three-dimensional space of scales
$\boldsymbol{r}$
plus the interface distance
$y_c$
.
This multi-dimensional nature of the problem is challenging for a clear understanding of the flow physics and reduction of the dimensionality is demanding. The most common solution is to consider the hyperplane
$(r_z, Y_c)$
for
$r_x=r_y=0$
because the space of spanwise scales is the one carrying most of the information about the space of scales in bounded flows (Cimarelli et al. Reference Cimarelli, Boga, Pavan, Costa and Stalio2024). In this hyperplane, the two-point difference operator applied to mean quantities gives a zero contribution while the two-point average operator does not alter mean quantities. Hence, (2.2) further simplifies and reads
\begin{align} \frac {\partial \langle \delta q^2 \delta u_j \!\rangle }{\partial r_j} -\! 2\nu \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial r_j \partial r_j} + \frac {\partial \langle \delta q^2 \widetilde {v} \rangle }{\partial y_c} + \frac {2}{\rho } \frac {\partial \langle \delta p \delta v \rangle }{\partial y_c} - \frac {\nu }{2} \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial y_c \partial y_c} = - 2 \langle \delta u \delta v \rangle \!\left ( \!\frac {\partial U}{\partial y} \!\right )\! - 4 \langle \epsilon \rangle .\end{align}
Equation (2.3) describes the turbulent wind as composed of inertial and diffusive transport mechanisms in the space of scales
\begin{equation} T_r = -\frac {\partial \langle \delta q^2 \delta u_j \rangle }{\partial r_j} ,\qquad D_r = 2\nu \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial r_j \partial r_j}, \end{equation}
and in the interface-normal direction
\begin{equation} T_c = - \frac {\partial \langle \delta q^2 \widetilde {v} \rangle }{\partial y_c} - \frac {2}{\rho } \frac {\partial \langle \delta p \delta v \rangle }{\partial y_c}, \qquad D_c = \frac {\nu }{2} \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial y_c \partial y_c} .\end{equation}
The source and sink of these transport mechanisms are provided by turbulence production by mean shear and by turbulent dissipation
\begin{equation} \varPi = - 2 \langle \delta u \delta v \rangle \left ( \frac {\partial U}{\partial y} \right ) ,\qquad E = 4 \langle \epsilon \rangle. \end{equation}
Hence, (2.3) is a budget that in a symbolic form reads
$T_r + D_r + T_c + D_c + \xi = 0$
, where
$\xi = \varPi - E$
is the net source/sink of scale energy in the flow.
In order to highlight the field of fluxes of scale energy occurring in the
$(r_x=r_y=0)$
-hyperplane, it is useful to recast (2.3) in a conservative form
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\phi } = \zeta ,\end{equation}
where the divergence operator and the scale-energy fluxes are defined as
\begin{equation} \boldsymbol{\nabla }= \left [ \begin{split} \frac {\partial }{\partial r_z}\\ \frac {\partial }{\partial y_c} \end{split} \right ] ,\qquad \boldsymbol{\phi } = \left [ \begin{split} \phi _{r_z} \\ \phi _c \end{split} \right ] = \left [ \begin{split} \partial \langle \delta q^2 \delta w \rangle - 2\nu \frac {\partial \langle \delta q^2 \rangle }{\partial r_z}\\ \langle \delta q^2 \widetilde {v} \rangle + \frac {2}{\rho } \langle \delta p \delta v \rangle - \frac {\nu }{2} \frac {\partial \langle \delta q^2 \rangle }{\partial y_c} \end{split} \right ], \end{equation}
while
$\zeta$
is an extended source term that take into account not only the net source/sink of scale energy
$\xi$
but also the exchanges of scale energy towards the
$r_x \neq 0$
and
$r_y \neq 0$
spaces
\begin{equation} \zeta = \xi - \frac {\partial \langle \delta q^2 \delta u \rangle }{\partial r_x} + 2\nu \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial r_x \partial r_x} - \frac {\partial \langle \delta q^2 \delta v \rangle }{\partial r_y} + 2\nu \frac {\partial ^2 \langle \delta q^2 \rangle }{\partial r_y \partial r_y} .\end{equation}
Hence, the conservative form (2.7) of the generalised Kolmogorov equation in the
$(r_x=r_y=0)$
-hyperplane highlights that a field of fluxes
$\boldsymbol{\phi }$
occurs in the two-dimensional space of spanwise scales and interface-normal positions
$(r_z, y_c)$
that is driven by an extended source term
$\zeta$
determined by the net source/sink of scale energy
$\xi$
plus the uncovered transports in the space of streamwise and interface-normal scales,
$r_x$
and
$r_y$
, respectively.
3. Large-scale asymptotic: turbulent kinetic energy budget
The generalised Kolmogorov equation restricted to the
$(r_x=r_y=0)$
-hyperplane exhibits a well-defined asymptotic behaviour at large separations
$r_z / \ell \gg 1$
, where
$\ell$
is the relevant correlation length of the flow. As an example, the second-order structure function can be rewritten as
$\langle \delta q^2 \rangle (r_z,y_c) = 4 \langle k \rangle (y_c) - 2 R_{ii} (r_z,y_c)$
, where
$k=u_i u_i/2$
is the turbulent kinetic energy and
$R_{\textit{ij}} (r_z, y_c) =\langle u_i(x,y_c,z+r_z/2,t) \, u_j(x,y_c,z-r_z/2,t) \rangle$
is the two-point correlation tensor. It is then clear that, for large separations
$r_z / \ell \gg 1$
, the trace of the two-point correlation tensor vanishes and the second-order structure function reduces to four times the turbulent kinetic energy,
$\langle \delta q^2 \rangle (r_z/\ell \gg 1, y_c) = 4 \langle k \rangle (y_c)$
. As shown in Marati, Casciola & Piva (Reference Marati, Casciola and Piva2004), the same reasoning applies to different terms of the generalised Kolmogorov (2.3), thus showing that it reduces within a factor of 4 to the turbulent kinetic energy equation
\begin{equation} \frac {{\rm d} \langle k v \rangle }{{\rm d} y} + \frac {1}{\rho } \frac {{\rm d} \langle p v \rangle }{{\rm d} y} - \nu \frac {d^2 \langle k \rangle }{{\rm d} y^2} = - \langle u v \rangle \frac {{\rm d}U}{{\rm d}y} - \langle \epsilon \rangle .\end{equation}
3.1. Sources and spatial fluxes
In figure 3(a), the profiles of turbulence production
$- \langle u v \rangle ({\rm d}U/{\rm d}y)$
and dissipation
$- \langle \epsilon \rangle$
are reported. Analogously to what happens in wall-bounded flows, turbulence production is mostly active in an intermediate layer, the so-called buffer layer for
$7 \lt y^+ \lt 50$
. In this region, the production rate exceeds the dissipation rate. Hence, the buffer layer is a net source region of turbulence. On the contrary, the very near-interface region, the so-called viscous sublayer for
$y^+ \lt 7$
, is a net sink region of turbulence, being the site of the highest levels of viscous dissipation and of a vanishingly small production. Another sink region for turbulence, although weaker, is the outer flow for
$y^+ \gt 0.6 {Re}_\tau$
. Here, the dissipation rate is small but prevails over the vanishing scaling of turbulence production, thus leading to a weak sink region. Intermediate between the buffer layer and the outer region, an overlap layer takes place for
$50 \lt y^+ \lt 0.6 {Re}_\tau$
, where production and dissipation almost balance each other. This equilibrium layer is also the site of the logarithmic region for the mean velocity profile.
(a) Profiles of turbulence production
$- \langle u v \rangle ^+ ({\rm d}U/{\rm d}y)^+$
(solid line) and dissipation
$- \langle \epsilon \rangle ^+$
(dashed line). (b) Profiles of spatial flux
$\langle \psi \rangle ^+$
. Wind–wave data are reported in black while channel data are reported in red.
This behaviour of the energy source/sink mechanisms qualitatively conforms with that commonly observed in wall-bounded flows, in agreement with laboratory experiments under low wind conditions (Yousefi, Veron & Buckley Reference Yousefi, Veron and Buckley2021). The main difference is indeed of a quantitative nature. In particular, we found that the peak of turbulence production
$- \langle u v \rangle ({\rm d}U/{\rm d}y)$
is slightly increased and shifted away, being located at
$y^+ = 13.8$
against the value of
$y^+ = 11.5$
occurring in the channel. Turbulence production in the wind–wave problem is found to be slightly higher than in channels in the outer part of the buffer layer and smaller in the inner part. On the contrary, the rate of turbulence dissipation is always slightly smaller in the wind–wave problem. Of particular interest is its behaviour in the viscous sublayer, as better shown in the inset of figure 3(a). There, the classical channel flow configuration exhibits a drastic increase of the rate of turbulent dissipation
$\langle \epsilon \rangle$
by moving towards the wall where its maximum value is reached. On the contrary, in the wind–wave problem such an increase is absent. In particular, an almost flat scaling of dissipation is observed for
$y^+ \lt 5$
where it remains almost constant around the value
$\langle \epsilon \rangle ^+ = 0.135$
. As a consequence, the maximum rate of dissipation is significantly reduced with respect to wall-bounded turbulence, where dissipation reaches a maximum value at the wall of the order of
$\langle \epsilon \rangle ^+ = 0.195$
. The physical origin of the near-interface scaling of turbulent dissipation is reported in § 3.2.
By rewriting the turbulent kinetic energy budget (3.1) in a conservative form
\begin{equation} \frac {{\rm d} \langle \psi \rangle }{{\rm d}y} = - \langle u v \rangle \frac {{\rm d}U}{{\rm d}y} - \langle \epsilon \rangle ,\end{equation}
it is clear that the observed wind–wave effects on the source term
$- \langle u v \rangle ({\rm d}U/{\rm d}y) - \langle \epsilon \rangle$
also affect the spatial flux of turbulent kinetic energy
\begin{equation} \langle \psi \rangle = \langle k v \rangle + \frac {1}{\rho } \langle p v \rangle - \nu \frac {{\rm d} \langle k \rangle }{{\rm d}y} .\end{equation}
As shown in figure 3(b), the wind–wave problem as much as wall turbulence is characterised by a spatial flux taking its origin from the peak of energy excess provided by the source term
$- \langle u v \rangle ({\rm d}U/{\rm d}y) - \langle \epsilon \rangle$
in the buffer layer. In particular, the spatial flux redistributes the turbulent energy excess of the buffer layer towards inner
$\langle \psi \rangle \lt 0$
and outer
$\langle \psi \rangle \gt 0$
locations, thus feeding the corresponding sink regions. Even if not shown for brevity reasons, the flux directed towards the interface is mostly driven by viscous diffusion,
$\langle \psi \rangle \sim - \nu {\rm d} \langle k \rangle /{\rm d}y$
for
$y^+ \lt 13.8$
while the flux directed towards the bulk of the flow is mostly given by inertial mechanisms,
$\langle \psi \rangle \sim \langle k v \rangle$
for
$y^+ \gt 13.8$
. The pressure flux is always directed towards the interface
$\langle p v \rangle /\rho \lt 0$
, but its magnitude is almost always negligible
$|\langle p v \rangle /\rho |/|\langle \psi \rangle | \ll 1$
.
The main differences in the spatial redistribution of turbulent kinetic energy of the wind–wave problem with respect to wall turbulence can be attributed to the main differences already observed for the source term. In particular, the lower dissipative nature of the interface region with respect to the wall region is associated with a reduced intensity of the spatial flux feeding it. On the contrary, the spatial flux feeding the outer region is enhanced. This outer enhancement of the flux can be related to the fact that the source
$- \langle u v \rangle ({\rm d}U/{\rm d}y) - \langle \epsilon \rangle$
is not reduced in the wind–wave problem with respect to wall turbulence. Hence, the energy excess provided by the buffer layer that has not been conveyed towards the interface region is instead redirected towards to outer flow regions.
3.2. Near-interface scaling of turbulent dissipation
The lower dissipative nature of the viscous sublayer in the wind–wave problem with respect to wall turbulence is inherently related to the low-drag feature of the former and, hence, to the wave-induced Stokes sublayer. In order to support this conjecture, it is relevant to explicitly write the near-wall scaling of dissipation in wall turbulence
\begin{equation} \langle \epsilon \rangle (y) = \langle \epsilon \rangle _w + \left (\frac {{\rm d} \langle \epsilon \rangle }{{\rm d} y} \right )_w \!\! y + \mathcal{O}(y^2) ,\end{equation}
where the subscript w denotes variables evaluated at the wall and
\begin{equation} \left (\frac {{\rm d} \langle \epsilon \rangle }{{\rm d} y} \right )_w = \frac {2}{\rho } \left ( \langle \frac {\partial p}{\partial x} \frac {\partial u}{\partial y} \rangle + \langle \frac {\partial p}{\partial z} \frac {\partial w}{\partial y} \rangle \right ) \lt 0 .\end{equation}
In wall turbulence, the high levels of the two covariances in (3.5) are shown in figure 4(a) and can be associated with a flow pattern sustained by the turbulent motion populating the bulk of the flow. Indeed, in the viscous sublayer, turbulence cannot self-sustain. Hence, the pattern taken by the fluctuating pressure and velocity fields is externally induced as a near-wall footprint of the most coherent structures of the buffer layer. The coherence of this flow pattern reflects to a significant degree of anti-correlation between horizontal pressure gradients and vertical gradients of horizontal velocities, think of the pattern of
$(\partial p/\partial x)$
and of
$(\partial u/\partial y)$
induced at the wall by streamwise velocity streaks and of
$(\partial p/\partial z)$
and of
$(\partial w/\partial y)$
induced at the wall by quasi-streamwise vortices. To measure this degree of anti-correlation, we compute the following correlation functions:
\begin{equation} \beta _x = \frac {\langle (\partial p/\partial x) (\partial u/\partial y) \rangle }{\sqrt {\langle (\partial p/\partial x)^2 \rangle \langle (\partial u/\partial y)^2 \rangle }}, \qquad \beta _z = \frac {\langle (\partial p/\partial z) (\partial w/\partial y) \rangle }{\sqrt {\langle (\partial p/\partial z)^2 \rangle \langle (\partial w/\partial y)^2 \rangle }}. \end{equation}
As shown in figure 4(b), the degree of anti-correlation is significant for
$\beta _z$
while is almost negligible for
$\beta _x$
, thus showing that the magnitude of
${\rm d} \langle \epsilon \rangle / {\rm d}y$
in wall turbulence is essentially the result of the near-wall velocity and pressure patterns induced by quasi-streamwise vortices in the buffer layer.
(a) Profiles of the covariances
$\langle (\partial p/\partial x) (\partial u/\partial y) \rangle$
(solid line) and
$\langle (\partial p/\partial z) (\partial w/\partial y) \rangle$
(dashed line) contributing to the scaling (3.5). (b) Profiles of the correlation functions
$\beta _x$
(solid line) and
$\beta _z$
(dashed line) from (3.6). The vertical grey line denotes the outer limit of the wave-induced Stokes sublayer. Wind–wave data are reported in black while channel data are reported in red.
In the wind–wave problem the scenario is modified since we measure
${\rm d} \langle \epsilon \rangle / {\rm d}y \approx 0$
. As shown in Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023), a very thin wave-induced Stokes sublayer is found to dominate the flow dynamics of the interface region for
$y^+ \lt 2$
. There, the velocity gradients
$\partial w / \partial y$
are found to be correlated with the alternating pressure gradients
$\partial p / \partial z$
induced by the non-sheltering behaviour of the wind over the skewed water–wave pattern. Hence, we might expect that, in the wave-induced Stokes sublayer, both the covariance
$\langle (\partial p/\partial z) (\partial w/\partial y) \rangle$
and the correlation
$\beta _z$
are not negligible, thus leading to a non-negligible
${\rm d} \langle \epsilon \rangle / {\rm d}y$
. On the contrary, outside the Stokes sublayer for
$y^+ \gt 2$
, the alternating pressure gradients induced by the wave pattern are still present but not the associated velocity field, as demonstrated in Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023). As a consequence, we might expect that the viscous sublayer for
$2\lt y^+\lt 7$
is dominated by pressure fluctuations induced by the wave pattern and velocity fluctuations induced from above by the coherent structures of the buffer layer. Since, the wave pattern and the structures of the buffer layer are not correlated, we might expect that, for
$2\lt y^+\lt 7$
, the correlations
$\beta _x$
and
$\beta _z$
and the covariances
$\langle (\partial p/\partial x) (\partial u/\partial y) \rangle$
and
$\langle (\partial p/\partial z) (\partial w/\partial y) \rangle$
are all negligible, thus explaining the observed
${\rm d} \langle \epsilon \rangle / {\rm d}y \approx 0$
. As shown in figure 4(a), within the wave-induced Stokes sublayer
$y^+ \lt 2$
,
${\rm d} \langle \epsilon \rangle / {\rm d}y$
is not negligible because of the covariance
$\langle (\partial p/\partial z) (\partial w/\partial y) \rangle$
induced by the Stokes sublayer itself, as demonstrated by the increasing value of anti-correlation
$\beta _z$
, see figure 4(b). On the other hand, outside the Stokes sublayer for
$2\lt y^+\lt 7$
, both the covariances
$\langle (\partial p/\partial x) (\partial u/\partial y) \rangle$
and
$\langle (\partial p/\partial z) (\partial w/\partial y) \rangle$
are negligible, thus leading to
${\rm d} \langle \epsilon \rangle / {\rm d}y \approx 0$
, see again figure 4(a). The uncorrelation between the pressure and velocity fields for
$2\lt y^+\lt 7$
is also demonstrated by the negligible values of
$\beta _x$
and
$\beta _z$
, see again figure 4(b).
In conclusion, we conjecture that the low dissipative nature of the wind–wave interface region is associated with the presence of a wave-induced Stokes sublayer whose penetration length,
$l_s^+=2$
, is within the viscous sublayer depth,
$y^+=7$
. In the wave-induced Stokes sublayer, velocity and pressure gradients, although weak, are correlated, both being induced by the same phenomenon of wind–wave interaction. In the outer part of the viscous sublayer,
$2 \lt y^+ \lt 7$
, wave-induced pressure fluctuations interact with bulk-induced velocity fluctuations in an incoherent manner. This lack of correlation of the velocity and pressure solutions in the outer viscous sublayer is at the basis of the low dissipative nature of the wind–wave interface region and, hence, of its low-drag feature. This finding agrees with laboratory experiments by Yousefi et al. (Reference Yousefi, Veron and Buckley2021). They found that wind turbulence and wave-induced motions interact when the Stokes sublayer height, evaluated as
$\delta _s^+= [ \lambda _x^+/(\pi c_x^+)]^{1/2}$
, exceeds the viscous sublayer. In their experiments, this condition is verified for sufficiently high wind speeds. By substituting the wave parameters of our DNS at low wind speed, we measure
$\delta _s^+ \approx 3.8$
which is well within the viscous sublayer depth. Hence, the observed lack of correlation between wave-induced motions and wind turbulence is also supported by experimental findings.
4. The multi-dimensional compound space of scales and positions
4.1. Scale energy
We start the analysis of the compound
$(r_z,yc)$
-space by considering the distribution of scale energy
$\langle \delta q^2 \rangle$
, as shown in figure 5(a). The buffer layer is recognised as the region of highest activity of turbulence. The peak activity is found to occur in a well-defined range of intermediate scales that clearly matches the length scales of the turbulent structures involved in the self-sustaining mechanisms of this region of the flow. From this region of the flow, the intensity of the turbulent fluctuations shows a rapid decrease by moving towards the interface while a more smooth decay is observed by moving towards the bulk of the flow. This behaviour qualitatively resembles the one observed in wall turbulence and reported in figure 5(b). The main differences are indeed of quantitative nature and consist of an increased intensity and a shift towards larger scales and interface distances of scale energy in the wind–wave problem with respect to wall turbulence. In particular, we measure that the peak intensity of scale energy is
$\langle \delta q^2 \rangle _{max}^+ = 29$
and is located at
$(r_z^+, y_c^+) = (85, 19)$
in the wind–wave problem with respect to
$\langle \delta q^2 \rangle _{max}^+ = 19$
occurring at
$(r_z^+, y_c^+) = (70, 16)$
in wall turbulence. This increase intensity and shift towards larger scales and positions conforms with the presence of a more coherent flow pattern, as classically observed in drag-reducing flows (White & Mungal Reference White and Mungal2008). Phenomenologically speaking, in drag-reducing flows turbulence is weakened near the boundary and the turbulence self-sustaining mechanisms of the buffer layer are moved outwards. As a consequence, the quasi-streamwise vortices and streaks involved in the self-sustaining cycle of turbulence (Jiménez & Pinelli Reference Jiménez and Pinelli1999), being located further away from the boundary, can be larger and more coherent, as demonstrated here by the higher intensity and shift of the iso-contours of scale energy.
Scale-energy
$\langle \delta q^2 \rangle ^+$
distribution in the compound space of scales and positions
$(r_z^+, y_c^+)$
for the wind–wave problem (a) and for the turbulent channel (b).
4.2. Scale-energy paths and sources
The description of the energetics of the flow given by the turbulent kinetic energy budget in § 3, is here extended by removing the asymptotic limit of the generalised Kolmogorov equation, i.e. by distinguishing the different scales involved in the processes. In particular, we start by considering the energetic behaviour of the flow described by the generalised Kolmogorov equation in its conservative form (2.7). In figure 6, the field of fluxes
$\boldsymbol{\phi }$
transporting scale energy among different spanwise scales and positions
$(r_z, y_c)$
is shown on top of iso-contours of the source term feeding (
$\xi \gt 0$
) or absorbing (
$\xi \lt 0$
) them. Again, both the wind–wave and the channel solutions are reported for comparison in figures 6(a) and 6(b), respectively. The buffer layer is again recognised as the scale-energy source region of the flow where the source term reaches its maximum. Accordingly, the field of fluxes takes its origin from a singularity point
$|\boldsymbol{\phi }| = 0$
located within the buffer layer and spreads, feeding turbulence fluctuations at different sizes of both the bulk and interface regions. The branch of fluxes flowing towards the bulk of the flow, while ascending
$\phi _c \gt 0$
, exhibits both forward
$\phi _{r_z} \lt 0$
and reverse
$\phi _{r_z} \gt 0$
cascades, see also figure 7. The picture consists of spatially ascending fluxes taking their origin from a well-defined spanwise scale of the buffer layer that feed wider turbulence fluctuations before reaching the
$y_c$
-distributed dissipative range at small scales; see Cimarelli et al. (Reference Cimarelli, De Angelis and Casciola2013, Reference Cimarelli, De Angelis, Jiménez and Casciola2016, Reference Cimarelli, Boga, Pavan, Costa and Stalio2024) for a detailed analysis of the physics of such an energetic pattern. On the contrary, the branch of fluxes flowing towards the interface region,
$\phi _c \lt 0$
, does not exhibit cascade phenomena in the space of scales,
$|\phi _{r_z}| / |\phi _c| \ll 1$
. Hence, the
$r_z$
-distributed scale-energy excess of the buffer layer is simply transported towards the interface without any significant modulation in scale space. As a result, the scale-energy sink of the viscous sublayer is performed by a viscous dissipation whose velocity gradients are the near-interface footprint of the variety of scales populating the outer flow regions, see Cimarelli et al. (Reference Cimarelli, Boga, Pavan, Costa and Stalio2024) for a deeper physical analysis.
Streamlines of the flux of scale energy
$(\phi _{r_z}^+, \phi _c^+)$
in the compound space of scales and positions
$(r_z^+, y_c^+)$
superimposed on iso-contours of the source term
$\xi ^+$
for the wind–wave problem (a) and for the turbulent channel (b).
Iso-contours of scale-energy flux
$\phi _{r_z}^+$
in the compound space of scales and positions
$(r_z^+, y_c^+)$
for the wind–wave problem (a) and for the turbulent channel (b).
The overall multi-dimensional pattern described so far qualitatively conforms with the present knowledge in wall-bounded flows, see e.g. Cimarelli et al. (Reference Cimarelli, De Angelis and Casciola2013, Reference Cimarelli, De Angelis, Jiménez and Casciola2016, Reference Cimarelli, Boga, Pavan, Costa and Stalio2024). In fact, comparing the two plots in figures 6(a) and 6(b), almost the same qualitative picture emerges. However, quantitatively speaking the wind–wave problem exhibits some distinctive features that are worth analysing. Overall, they consist of a larger-scale structure of the flow and enhanced reverse cascade phenomena, as shown in the following.
Starting from the buffer layer, it can be noted that the compound region of scales and positions involved in the source of scale energy is wider and thicker in the wind–wave problem. In particular, the peak of the source term
$\xi _{max}$
occurs at
$r_z^+ = 60$
for
$y_c^+ = 14.8$
in the wind–wave problem while in the turbulent channel it is at
$r_z^+ = 40$
for
$y_c^+ = 11.5$
. Note that this shift towards larger coherent scales and higher distances from the interface is a classical feature commonly observed in drag-reducing flows, as already pointed out in § 4.1
This displacement of the energy processes related to the production of turbulence towards larger scales and interface distances is also reflected in the field of fluxes
$\boldsymbol{\phi }$
. In particular, we measure that the singularity point at which the scale-energy paths in the wind–wave problem take their origin is located at
$r_z^+ = 82$
for
$y_c^+ = 16.2$
while in the turbulent channel it is at
$r_z^+ = 63$
for
$y_c^+ = 15.5$
. This shift towards larger scales and upward positions is a direct consequence of the similar shift observed for the scale-energy sourcing mechanisms
$\xi _{max}$
. The difference in the field of fluxes between the wind–wave and the channel problems is, however, not limited to this shift. Indeed, the branch of fluxes moving towards the bulk exhibits a more divergent pattern in the wind–wave problem. This is the result of the presence of more intense reverse cascade mechanisms, as clearly shown also in figure 7. Hence, it is possible to conjecture that the larger coherent structures observed in § 4.1 to populate the wind–wave problem that are also typical of drag-reducing flows are the results of enhanced reverse energy cascade processes and of sourcing mechanisms at wider scales. In particular, the shift of the sourcing mechanisms can be related to the generation of larger quasi-streamwise vortices and streaks locally to the buffer layer. On the other hand, the larger flow structures observed in the bulk of the flow can be related to a non-local effect of transport of the scale-energy excess of the buffer layer to the outer region where the reverse cascade is significantly more intense.
Let us finally notice that the source term, although weak, is positive also in the overlap layer and outer regions, see again figure 6. As shown in Hwang & Cossu (Reference Hwang and Cossu2010), the overlap layer is also the site of self-sustaining mechanisms of turbulence, although the amount of production there locally almost balances the rate of dissipation. As demonstrated in Cimarelli et al. (Reference Cimarelli, De Angelis, Schlatter, Brethouwer, Talamelli and Casciola2015), these phenomena are at the basis of the appearance of a second weaker peak of scale-energy source in the overlap layer whose scaling supports the hypothesis of an overlap layer dominated by attached eddies (Townsend Reference Townsend1976). As shown by the iso-contours of the source term in figure 6, this outer source region is also present in the wind–wave problem. The main difference is the lack of a sink region in between the outer and inner sources and the extension of the source to very high distances from the interface. The lack of a clear separation between the inner and outer sources can be attributed to the wave-induced shift of the inner source towards larger scales and positions. On the other hand, the increased extension of the outer source can be used to explain the intensification of the reverse cascade processes in the wind–wave problem. In particular, we argue that the deeper penetration of the outer source in the bulk of the turbulent wind has a more efficient repulsive effect for the fluxes, thus promoting their divergence and sustaining the reverse cascade.
4.3. Direct effects of wind–wave interactions
The wind–wave interactions have been analysed so far by considering their effect on the multi-dimensional compound space of scales and positions overall. We restrict ourselves here to the analysis to the very-near interface region where phenomena directly related to wind–wave interactions can be observed. As shown in Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023), the direct effect of the interaction of the wind with the wave pattern can be observed in the velocity field only in the very-near interface region for
$y^+ \le 2$
, where the wave-induced Stokes sublayer takes place. On the other hand, the footprint of the wind–wave interactions on the pressure field penetrates more deeply into the wind boundary layer with non-negligible effects up to
$y^+ \approx 60$
. The analysis of the scale-energy fluxes in the interface region confirms this observation. In fact, only the spatial flux carried out by the pressure field
$2\langle \delta p \delta v \rangle /\rho$
is found to be clearly modulated by the water–wave pattern. In contrast with the behaviour of the overall spatial flux
$\phi _c$
that is directing scale energy from the source peak in the buffer towards both the wall
$\phi _c \lt 0$
and the bulk of the flow
$\phi _c \gt 0$
, the contribution of the pressure field is a flux that is always directed from the bulk towards the wall, see figure 8(b). This scenario is changed in the near-interface region of the wind–wave problem. Indeed, as shown in figure 8(a), the spatial flux carried out by the pressure field exhibits a positive sign
$2\langle \delta p \delta v \rangle /\rho \gt 0$
for
$y^+ \lt 9$
at scales corresponding to a multiple of half the water wavelength,
$r_z = \lambda _z/2$
,
$3 \lambda _z /2$
and
$5 \lambda _z/2$
with
$\lambda _z^+ = 380$
. This is clearly a direct effect of the wind–wave interaction and represents a periodically distributed pumping of scale energy into the turbulent wind that otherwise would be simply transferred toward the interface and then dissipated; compare figures 8(a) and 8(b). We conjecture that this upward pumping of scale energy from the water–wave interface is the phenomenon responsible for the upward shift of the self-sustaining production mechanisms of turbulence in the buffer layer and hence for the generation of larger flow structures together with more intense reverse cascade phenomena. Hence, the Stokes sublayer pumping of scale energy can be understood as the energetic counterpart of the drag-reducing mechanism of the wave-induced Stokes sublayer.
Iso-contours of scale-energy flux
$(2\langle \delta p \delta v \rangle /\rho )^+$
in the compound space of scales and positions
$(r_z^+, y_c^+)$
for the wind–wave problem (a) and for the turbulent channel (b).
5. Turbulence decomposition and coarse-grained approaches
Almost every approach to turbulence is based on a decomposition of the turbulent field into different elements. The general aim is to provide a description of turbulence that is simpler than that given by the full Navier–Stokes equations. Famous examples are the Reynolds decomposition in the mean and fluctuating fields and the scale decomposition in large and small size motions obtained either via spectral operators or via physical operators such as the two-point difference used here. The study of the interactions between the different elements composing turbulence has certainly been fundamental for the development of the current state of theoretical and physical knowledge as we tried to do here for the wind–wave problem by using the generalised Kolmogorov equation. However, the turbulence decomposition has been fundamental also for the computation of numerical solutions of turbulence problems that otherwise would be prohibitive; think of the coarse-grained approaches given by the Reynolds-average equations and by the LES. In this context, the formalism of the generalised Kolmogorov equation is strictly related to the filtering approach of LES. As shown in Germano (Reference Germano2012), the two-point average and the two-point difference are the simplest operators for the large-scale filtering and the associated fluctuation, respectively. In particular, the second-order structure function
$\langle \delta q^2 \rangle$
can be directly related to the trace of the subgrid-stress tensor
$\langle \tau _{ii} \rangle /2$
that in turn represents the mean statistical value of turbulent energy at subgrid level (Germano Reference Germano2007a
,
Reference Germanob
). On the other hand, the terms of the generalised Kolmogorov equation, (2.1), can be understood as the energetic processes that can be captured by a coarse-grained approach to the Navier–Stokes equations at a spatial resolution
$\boldsymbol{r}$
(Cimarelli & De Angelis Reference Cimarelli and De Angelis2011, Reference Cimarelli and De Angelis2012).
5.1. Cross-over scales
In accordance to the previous arguments, the generalised Kolmogorov equation can be used to identify the characteristic scales at which different relevant processes dominate, thus allowing us to identify a priori the spatial resolution required in coarse-grained approaches to accurately resolve the wind–wave problem. In this context, it is well known that one of the most important feature of coarse-grained approaches is to accurately resolve the energy-containing motions where turbulence is produced. Furthermore, most of turbulence closures are based on the assumption of a filter length placed well within the inertial subrange of turbulence where the scale-energy flux is from large to small scales and isotropy is recovered. In order to identify these physically relevant ranges we introduce here three cross-over scales. The first one is the cross-over scale between forward and reverse cascade, defined as
\begin{equation} \phi _{r_z} (\ell _b, y_c) = 0 .\end{equation}
The cross-over scale
$\ell _b$
allows us to split the turbulence spectrum into scales that are dominated by reverse energy cascade processes for
$r_z \gt \ell _b$
from those that are involved in more classical forward cascade processes for
$r_z \lt \ell _b$
. The second relevant cross-over scale splits turbulence into a range dominated by production by mean shear and that dominated by inertial transport and is defined as
\begin{equation} T_r (\ell _s, y_c) = \varPi (\ell _s, y_c) .\end{equation}
The cross-over scale
$\ell _s$
is physically related to the more classical shear scale that from dimensional arguments is defined as
$L_s = \sqrt {\langle \epsilon \rangle / ({\rm d}U/{\rm d}y)^3}$
(Cimarelli et al. Reference Cimarelli, Boga, Pavan, Costa and Stalio2024). Physically speaking, the cross-over scale
$\ell _s$
allows us to distinguish the range of turbulence that is strongly anisotropic, being dominated by production by mean shear for
$r_z \gt \ell _s$
, from the range where an isotropy recovery is expected, being dominated by inertial transport phenomena for
$r_z \lt \ell _s$
. Finally, the third relevant cross-over scale is defined as
\begin{equation} D_r (\ell _\nu , y_c) = T_r (\ell _\nu , y_c) .\end{equation}
The cross-over scale
$\ell _\nu$
splits turbulence into a range of scales dominated by inertial transport for
$r_z \gt \ell _\nu$
and that where viscous diffusion prevails for
$r_z \lt \ell _\nu$
. Hence,
$\ell _\nu$
represents the viscous cutoff scale and is strictly related to the more classical Kolmogorov scale that from dimensional arguments is defined as
$\eta = (\nu ^3 /\langle \epsilon \rangle )^{1/4}$
. Clearly, the range of scales in between the viscous and shear scales quantitatively defines the inertial subrange of turbulence, being the dynamic of turbulence dominated by an inertial transport with isotropy recovery for
$\ell _\nu \lt r_z \lt \ell _s$
. By replacing the two cross-over scales
$\ell _\nu$
and
$\ell _s$
with their dimensional counterparts
$\eta$
and
$L_s$
, it is possible to derive asymptotic scaling for the inertial subrange. Indeed, by invoking an equilibrium assumption and a logarithmic law for the mean velocity profile in the overlap layer of wall turbulence, we have that
$\langle \epsilon \rangle = \langle u v \rangle ({\rm d}U/{\rm d}y) = u_\tau ^3/(\kappa y)$
, where
$\kappa$
is the von Kármán constant. By replacing this estimate of turbulent dissipation in the definition of the Kolmogorov and shear scales we have
\begin{equation} \eta ^+ \sim (\kappa y^+)^{1/4}, \qquad L_s \sim \kappa y^+ ,\end{equation}
thus showing that, asymptotically, the width of the inertial subrange scales as
$\kappa y^+ - (\kappa y^+)^{1/4}$
.
5.2. Scalings
The behaviour of the cross-over scale between forward and reverse cascades
$\ell _b$
is reported in figure 9(a). In the wind–wave problem, a minimum value
$\ell _b^+ = 81$
is found to occur at
$y_c^+ =14$
where the maximum value of the source
$\xi _{max}$
is attained. From this minimum, the scales involved in reverse cascade processes increase by moving away from the interface region, thus reaching a maximum value
$\ell _b^+ = 150$
for
$y_c^+ = 90$
. Indeed, in the overlap layer a slight decrease of
$\ell _b$
is observed, thus implying that the phenomena of reverse cascade penetrate more and more inside the turbulent spectrum in this region of the flow. This decrease is interrupted in the outer region of the flow where reverse cascade are finally absent and
$\ell _b$
formally goes to infinity. The main difference with respect to wall turbulence is given by the fact that the increase of the cross-over scale
$\ell _b$
in the inner region is faster, i.e. a smaller portion of the turbulent spectrum is dominated by reverse cascade processes. This observation, together with the already observed increase of the intensity of the reverse energy fluxes, further supports the idea that reverse cascade mechanisms have more impact in the flow dynamics of the wind–wave problem. These results suggest that the spanwise spatial resolution should be
$\unicode{x1D6E5} _z \lt \ell _b$
in coarse-grained approaches when using classical turbulence closures that strictly reproduce a forward cascade, such as eddy viscosity models. Otherwise, when using
$\unicode{x1D6E5} _z \gt \ell _b$
in order to reduce the computational costs, more sophisticated nonlinear closures able to capture a net backscatter from the small unresolved motion are demanded (Cimarelli & De Angelis Reference Cimarelli and De Angelis2014; Cimarelli, Abbà & Germano Reference Cimarelli, Abbà and Germano2019).
Cross-over scales for the wind–wave problem (black lines) and for the turbulent channel (red lines). (a) Scaling of the cross-over scale
$\ell _b^+$
between the forward and reverse cascades. (b) Scaling of the cross-over scale
$\ell _s^+$
between the production- and cascade-dominated range of scales (solid lines) and
$\ell _\nu ^+$
between the cascade- and viscous-dominated range of scales (dashed lines). The dotted lines denote the theoretical scalings
$\ell _s^+ = \kappa y_c^+$
and
$\ell _\nu ^+ = 7 (\kappa y_c^+)^{1/4}$
, where
$\kappa = 0.41$
.
The behaviour of the inertial subrange is reported in figure 9(b) where both cross-over scales
$\ell _s$
and
$\ell _\nu$
are shown. In the wind–wave problem, it is found that, for
$y_c^+ \lt 35$
, an inertial subrange cannot be realised because
$\ell _\nu \gt \ell _s$
. Hence, the anisotropic self-sustaining production processes of the buffer layer are superimposed on viscous diffusion mechanisms. For
$y_c^+ \gt 35$
, both cross-over scales increase but the higher rate of
$\ell _s$
with respect to
$\ell _\nu$
allows for the generation of an inertial subrange of increasing width. However, for
$y_c^+ \gt 100$
the behaviour of
$\ell _s$
flattens while
$\ell _\nu$
still increases, thus eroding the inertial subrange. For
$y_c^+ \lt 100$
, almost the same behaviour is observed in wall turbulence that indeed agrees with the expected asymptotic scaling (5.4). On the contrary, for
$y_c^+ \gt 100$
, the shear scale
$\ell _s$
in wall turbulence is found to maintain an increasing behaviour that, compared with a weaker increase of the viscous scale
$\ell _\nu$
, allows for the generation of an increasingly larger inertial subrange, contrary to the wind–wave problem. This different behaviour in the outer region can be associated with the different boundary conditions experienced by the wind–wave problem in the top boundary with respect to the turbulent channel that, as already mentioned, leads to enhanced sourcing mechanisms in the outer region of the wind–wave problem. In coarse-grained approaches, the present results suggest the use of spanwise spatial resolution
$\unicode{x1D6E5} _z \lt \ell _s$
in order to resolve the production-dominated energy-containing motions. By comparing figures 9(a) and 9(b), this condition is more stringent than that given by the scale
$\ell _b$
and in general suggests that the spanwise spatial resolution in the buffer layer, where most of the energy is produced, cannot exceed the value
$\unicode{x1D6E5} _z^+ \approx 20$
, analogously to classical recommendations in wall turbulence (Cimarelli & De Angelis Reference Cimarelli and De Angelis2012).
6. Conclusions
In Cimarelli et al. (Reference Cimarelli, Romoli and Stalio2023), the interaction of a turbulent wind with a water–wave surface has been solved by means of a fully coupled numerical framework. Starting from a quiescent initial condition characterised by a flat air–water interface, the flow solution has been evolved until a statistical equilibrium between the turbulent wind and the water–wave pattern is reached. The solution is characterised by a skewed pattern of water waves of low elevation and by a turbulent wind where streamwise momentum transfer is reduced by the presence of a wave-induced spanwise oscillating Stokes sublayer. The strongly inhomogeneous and multiscale features of these phenomena are addressed here by using, for the first time, the Kolmogorov equation generalised to the wind–wave problem.
We found that the turbulent wind is characterised by the presence of more coherent flow structures compared with wall turbulence. We recognise that, at the basis of the generation of these larger and more intense flow structures, is an upward shift of the self-sustaining production processes of turbulence and the presence of stronger reverse energy cascade phenomena. The upward shift of turbulence production is here recognised to be the result of a periodically distributed pumping of scale energy from the interface region induced by the pressure field associated with the wave-induced Stokes sublayer. Contrary to wall turbulence, where the spatial flux carried out by pressure is always pointing towards the wall, we found a reverse of sign of the spatial flux driven by the pressure field in the interface region for scales that match those of the water–wave pattern. This upward pressure flux of scale energy is here conjectured to be at the basis of the upward shift of the production mechanisms of turbulence. The generated flow structures such as streamwise velocity streaks and quasi-streamwise vortices, are then larger and more intense with respect to wall turbulence, being less constrained by the distance from the interface. The energy associated with these larger flow structures, while transported in the bulk of the turbulent wind, is conveyed to even larger flow structures by reverse energy cascade phenomena that are found to be more intense than those observed in wall turbulence, thus explaining the more coherent flow structures populating the entire turbulent wind.
These results have strong consequences for wind boundary layer theories but also for coarse-grained approaches to their solution. Indeed, these results unveil that turbulence closures must be able to reproduce a flow of energy from nowhere in the space of unresolved scales when the level of coarsening is such that the reverse energy cascade lies in the subgrid-scale motion. This scenario is very critical for the development of high-fidelity coarse-grained approaches deemed essential to catch the high Reynolds number regime typical of atmospheric realisations. In this context, we have shown that the generalised Kolmogorov equation provides a complete theoretical framework for assessing the spatial resolution and the associated best suited turbulence closure for large-eddy simulation. By introducing the cross-over scales
$\ell _s$
and
$\ell _\nu$
, we have quantitatively defined the scaling of the inertial subrange of turbulence within which the spatial resolution should lie in order to directly resolve the anisotropic energy-containing motions. As rigorously shown by the scaling of the split energy cascade scale
$\ell _b$
, this range of scales is dominated by a forward energy cascade, thus suggesting the use of classical eddy viscosity closures when the grid resolution lies within it.
Finally, a theoretical analysis of the low dissipative nature of the interface region is introduced and reveals that the scaling of dissipation is strictly related to the degree of correlation between pressure and velocity gradients. We found that pressure fluctuations generated by the wave-induced Stokes sublayer interact in an incoherent manner with the wind-induced velocity fluctuations in the external part of the interface region. The resulting lack of correlation acts as a layering of the flow dynamics in the interface region and is recognised to be at the basis of its low dissipative nature and, hence, of its low-drag feature. The lack of coherence between waves and turbulent wind is expected to exists as long as the Stokes sublayer depth is thinner than the viscous sublayer depth. Such physical explanation of the lower dissipative nature of wind waves at low Reynolds number, as well as the periodically distributed pumping of scale energy as found by the analysis of the generalised Kolmogorov equation, offer a possible interpretation of the observed drag reduction in low wind speed regimes as observed in field experiments.
Acknowledgements
We acknowledge the Italian supercomputing centre CINECA that, under the ISCRA B project WWW (turbulent Wind and Water Waves interaction: multiscale analysis and modelling) and the ISCRA C project WWIF (simulations of Wind–Wave Interfacial Flows), provided the high-performance computing resources for the present work.
Funding
This work was supported by the Italian Ministry of University and Research (MUR), through the PRIN 2022 PNRR Project P20229KW2R – SEAPLANE – ‘SimulatIon and modelling of interface fluxes in wind-wave flows for an improved climate science’, CUP E53D23017010001, funded by the National Recovery and Resilience Plan (PNRR), Italy, Mission 04 Component 2 Investment 1.1 – NextGenerationEU. Financial support has been provided also by the project ECOSYSTER – SPOKE 6 under the National Recovery and Resilience Plan, Mission 04 Component 2 Investment 1.5 – NextGenerationEU, Call for tender no. 3277 dated 30 December 2021.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available under request.
Open access
