<figure><img src="image_1.png" /></figure>

HAL
open science

# A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves

Liad Paskin, Boris Conan, Yves Perignon, Sandrine Aubrun

► To cite this version:

Liad Paskin, Boris Conan, Yves Perignon, Sandrine Aubrun. A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves. Atmosphere, 2022, 13 (12), pp.2012. 10.3390/atmos13122012. hal-03946404

HAL Id: hal-03946404

https://hal.science/hal-03946404

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves

Liad Paskin, Boris Conan, Yves Perignon and Sandrine Aubrun

<figure><img src="image_1.png" /></figure>
<figure><img src="image_3.png" /></figure>

Article

# A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves

Liad Paskin, Boris Conan, Yves Perignon and Sandrine Aubrun

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>
Article

# A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves

Liad Paskin* and Boris Conan and Yves Perignon and Sandrine Aubrun

Nantes Université, École Centrale Nantes, CNRS, LHEEA, UMR 6598, F-44000 Nantes, France

* Correspondence: liad.paskin@ec-nantes.fr

**Abstract:** We performed large eddy simulations to study micro-scale wind-wave interactions under undisturbed freestream conditions. We identified that standard approaches lead to wave-related disturbances at the top boundary. Therefore, we developed a numerical strategy to maintain an undisturbed wind speed at the top, while considering arbitrary waves at the bottom. In a broader context, the method is capable of controlling the wind speed at any height in the domain, and may also be used to enhance atmospheric simulations over land. The method comprises an evolution equation that controls the dynamic evolution of the large-scale driving force, representing the geostrophic forcing from the meso- to the micro-scales. In flat-bottom applications, this guided the reference freestream velocities towards a certain target; convergence to a steady state regime was favored and self-similarity was ensured. In wavy bottom applications considering the prescription of a monochromatic wave, we were able to maintain a quasi-steady wind speed close to the target on the freestream. The wave-induced disturbances were then investigated as functions of varying wave age conditions. We performed a systematic wave age variation study by varying the reference wind speed, and evaluated wave-induced disturbances in the velocity, normal, and shear stress profiles.

**Keywords:** marine atmospheric boundary layer; large eddy simulation; air-sea fluxes; meso- to micro-scale forcing

## 1. Introduction

The marine atmospheric boundary layer (MABL) is often subject to significant wave-induced (WI) disturbances governing complex wind-wave interactions above the ocean. Ever since the pioneering theoretical work of [1–5], and up to modern experimental campaigns, such as those discussed in [6–10], the wave age (WA = c/U) has been used to classify wind-wave interactions in terms of the wave-phase (c) and mean wind (U) velocities. When a swell encounters light winds (large WA), the momentum transfers from the sea to the atmosphere in a situation referred to as old seas, as opposed to the wave-growing stage referred to as young seas (small WA).

A wave boundary layer (WBL) exists in the vicinity of the water surface, where the wave's influence directly alters the air-flow dynamics with the appearance of a WI flow [10,11], and notably a WI stress perturbation [12,13]. Neglecting viscous effects, the WI stress integrates into the form drag in the lower surface and is responsible for the vertical transfer of momentum between waves and wind [14,15]. If the waves travel sufficiently fast (large WA) with a significant wave slope, the upward transfer is such that the wave signature emerges as a low-level jet, characterized by a local maximum in the wind velocity profiles [16–18].

The definition of a WBL height is in debate [11,19,20], as the significance of the WI disturbance depends on the variables of interest [10]. Many authors considered the WBL height as being proportional to, and in the same order of magnitude as, the significant height [20,21]. However, it was observed that during swell conditions, the WBL can reach higher altitudes [10,19,22,23]. Moreover, WI disturbances were observed to decay with

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Paskin, L.; Conan, B.; Perignon, Y.; Aubrun, S. A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves. *Atmosphere* **2022**, *13*, 2012. https://doi.org/10.3390/atmos13122012

Academic Editor: Jimy Dudhia

Received: 31 October 2022

Accepted: 23 November 2022

Published: 30 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Article

# A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves

Liad Paskin* and Boris Conan and Yves Perignon and Sandrine Aubrun

Nantes Université, École Centrale Nantes, CNRS, LHEEA, UMR 6598, F-44000 Nantes, France

* Correspondence: liad.paskin@ec-nantes.fr

**Abstract:** We performed large eddy simulations to study micro-scale wind-wave interactions under undisturbed freestream conditions. We identified that standard approaches lead to wave-related disturbances at the top boundary. Therefore, we developed a numerical strategy to maintain an undisturbed wind speed at the top, while considering arbitrary waves at the bottom. In a broader context, the method is capable of controlling the wind speed at any height in the domain, and may also be used to enhance atmospheric simulations over land. The method comprises an evolution equation that controls the dynamic evolution of the large-scale driving force, representing the geostrophic forcing from the meso- to the micro-scales. In flat-bottom applications, this guided the reference freestream velocities towards a certain target; convergence to a steady state regime was favored and self-similarity was ensured. In wavy bottom applications considering the prescription of a monochromatic wave, we were able to maintain a quasi-steady wind speed close to the target on the freestream. The wave-induced disturbances were then investigated as functions of varying wave age conditions. We performed a systematic wave age variation study by varying the reference wind speed, and evaluated wave-induced disturbances in the velocity, normal, and shear stress profiles.

**Keywords:** marine atmospheric boundary layer; large eddy simulation; air-sea fluxes; meso- to micro-scale forcing

## 1. Introduction

The marine atmospheric boundary layer (MABL) is often subject to significant wave-induced (WI) disturbances governing complex wind-wave interactions above the ocean. Ever since the pioneering theoretical work of [1–5], and up to modern experimental campaigns, such as those discussed in [6–10], the wave age (WA = c/U) has been used to classify wind-wave interactions in terms of the wave-phase (c) and mean wind (U) velocities. When a swell encounters light winds (large WA), the momentum transfers from the sea to the atmosphere in a situation referred to as old seas, as opposed to the wave-growing stage referred to as young seas (small WA).

A wave boundary layer (WBL) exists in the vicinity of the water surface, where the wave's influence directly alters the air-flow dynamics with the appearance of a WI flow [10,11], and notably a WI stress perturbation [12,13]. Neglecting viscous effects, the WI stress integrates into the form drag in the lower surface and is responsible for the vertical transfer of momentum between waves and wind [14,15]. If the waves travel sufficiently fast (large WA) with a significant wave slope, the upward transfer is such that the wave signature emerges as a low-level jet, characterized by a local maximum in the wind velocity profiles [16–18].

The definition of a WBL height is in debate [11,19,20], as the significance of the WI disturbance depends on the variables of interest [10]. Many authors considered the WBL height as being proportional to, and in the same order of magnitude as, the significant height [20,21]. However, it was observed that during swell conditions, the WBL can reach higher altitudes [10,19,22,23]. Moreover, WI disturbances were observed to decay with

<figure><img src="image_6.png" /></figure>

check for
updates

**Citation:** Paskin, L.; Conan, B.; Perignon, Y.; Aubrun, S. A Dynamic Large-Scale Driving-Force to Control the Targeted Wind Speed in Large Eddy Simulations above Ocean Waves. *Atmosphere* **2022**, *13*, 2012. https://doi.org/10.3390/atmos13122012

Academic Editor: Jimy Dudhia

Received: 31 October 2022

Accepted: 23 November 2022

Published: 30 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_7.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
height at an exponential proportional to the wavelength [11,20]. Particularly in old seas, the waves can significantly disturb the wind profiles, but a clear parametrization of this disturbance from the sea-state conditions has not yet been achieved.

Atmospheric forecast systems [24,25], and predictive tools for design and opera-
tional purposes [26,27], rely on low-cost, semi-empirical numerical modeling, as discussed,
e.g., in [28]. The Monin–Obukhov similarity theory (MOST), [29,30], is expected to hold
sufficiently far above the wave’s influence [20,23,31,32], i.e., above the so-defined WBL.
Therefore, Charnock [33] developed a pioneer, widely employed parametrization of the WI
form drag in the upper part of the MABL. Charnock’s parametrization was followed and
improved by many [34–37]. Still, an insufficient understanding of wind–wave interaction
mechanisms prevents the accurate characterization of WI motions, or stresses, in these
semi-empirical models [11,13].

It is expected that high-fidelity computational fluid dynamics could drive the de-
velopment of more robust semi-empirical predictive tools. Fully resolved turbulence is
achievable through direct numerical simulations (DNS), within limited Reynolds numbers
representing laboratory scales [38–44]. Alternatively, large eddy simulations (LES) employ
a limited level of modeling that intervenes only at small scales, being able to resolve the
largest and most energetic scales of turbulence, at large Reynolds numbers representing
oceanic scales (c.f. [45]). Moreover, LES are often coupled to other complex dynamics
characteristic of engineering applications, such as, e.g., wind energy farming (c.f. [46]). A
pioneer LES over wavy surfaces was proposed in [47], and was followed by the develop-
ments of [48–59], with a recent review of the state-of-the-art methodology presented in [60].
The DNS model presented in Sullivan et al. [38] was adapted in [48–50], into the LES model
(from the NCAR, USA) employed herein. Further research is still required, however, before
numerous and systematically varying LES can improve the existing low-order closures
to the WBL reconstruction in semi-empirical models. In this context, we propose a new
methodology to expand LES capabilities in wind-wave interaction studies.

In the atmosphere, the flow is driven by a geostrophic wind (U_G), which results from
a balance between the large-scale pressure and Coriolis forces. U_G is assumed not to be
affected by the local terrain. To study large-scale dynamics, atmospheric LES approaches
use these physics to drive the flow in the numerical domain [48,60]. To avoid having to
model the entire ABL and above, a homogeneous longitudinal pressure gradient (∂P/∂x)
is widely used in micro-scale studies to drive the flow instead [47,49–57,59]. In this case,
unlike in reality, the flow travels along the pressure gradient direction as in a wind tunnel.
Combined with the case-dependent WI stresses, and with the free-slip condition (zero-
vertical-gradient of horizontal velocities) imposed in the upper boundary, both of these
strategies allow a varying wind speed profile across the vertical extension of the MABL. In
old seas and neutral atmospheric conditions, we reported the speed-up of the longitudinal
velocity up to the upper boundary of the numerical domain [52]. The WI disturbance of
the freestream is unlikely to occur in the physical MABL but is intrinsic to the numerical
formulations employed in the literature. This is observable, e.g., in [49,51,53,57]. Weakly
unstable stratification was considered in [49], and these references also report the slow-
down of the freestream velocity in young seas, or due to the presence of curved but
stationary bottoms. To keep a constant freestream velocity above a time-varying wave, we
develop herein a general strategy to control the wind speed in an atmospheric LES.

In previous studies of wind–wave interactions by LES, the large-scale driving force
(here ∂P /∂x) was usually prescribed as a constant [47–57,60]. As an exception, an ar-
bitrary time-varying force was recently employed to simulate a wind gust event [59].
Over land, however, a few applications considered transient large-scale forcing terms
in LES (c.f. [61–67]). Based on LiDAR measurements in the ABL, the overall accuracy was
improved through the prescription of a linear-time-varying, and spatially homogeneous,
pressure gradient [61]. The measurement of time-varying geostrophic winds drove the
ABL during one or multiple diurnal cycles [62,63]. Employing time-varying and height-
dependent forcing terms is in fact a common practice of data assimilation in meso- to
height at an exponential proportional to the wavelength [11,20]. Particularly in old seas, the waves can significantly disturb the wind profiles, but a clear parametrization of this disturbance from the sea-state conditions has not yet been achieved.

Atmospheric forecast systems [24,25], and predictive tools for design and opera-
tional purposes [26,27], rely on low-cost, semi-empirical numerical modeling, as discussed,
e.g., in [28]. The Monin–Obukhov similarity theory (MOST), [29,30], is expected to hold
sufficiently far above the wave’s influence [20,23,31,32], i.e., above the so-defined WBL.
Therefore, Charnock [33] developed a pioneer, widely employed parametrization of the WI
form drag in the upper part of the MABL. Charnock’s parametrization was followed and
improved by many [34–37]. Still, an insufficient understanding of wind–wave interaction
mechanisms prevents the accurate characterization of WI motions, or stresses, in these
semi-empirical models [11,13].

It is expected that high-fidelity computational fluid dynamics could drive the de-
velopment of more robust semi-empirical predictive tools. Fully resolved turbulence is
achievable through direct numerical simulations (DNS), within limited Reynolds numbers
representing laboratory scales [38–44]. Alternatively, large eddy simulations (LES) employ
a limited level of modeling that intervenes only at small scales, being able to resolve the
largest and most energetic scales of turbulence, at large Reynolds numbers representing
oceanic scales (c.f. [45]). Moreover, LES are often coupled to other complex dynamics
characteristic of engineering applications, such as, e.g., wind energy farming (c.f. [46]). A
pioneer LES over wavy surfaces was proposed in [47], and was followed by the develop-
ments of [48–59], with a recent review of the state-of-the-art methodology presented in [60].
The DNS model presented in Sullivan et al. [38] was adapted in [48–50], into the LES model
(from the NCAR, USA) employed herein. Further research is still required, however, before
numerous and systematically varying LES can improve the existing low-order closures
to the WBL reconstruction in semi-empirical models. In this context, we propose a new
methodology to expand LES capabilities in wind-wave interaction studies.

In the atmosphere, the flow is driven by a geostrophic wind ($U_G$), which results from a balance between the large-scale pressure and Coriolis forces. $U_G$ is assumed not to be affected by the local terrain. To study large-scale dynamics, atmospheric LES approaches use these physics to drive the flow in the numerical domain [48,60]. To avoid having to model the entire ABL and above, a homogeneous longitudinal pressure gradient ($\partial P/\partial x$) is widely used in micro-scale studies to drive the flow instead [47,49–57,59]. In this case, unlike in reality, the flow travels along the pressure gradient direction as in a wind tunnel. Combined with the case-dependent WI stresses, and with the free-slip condition (zero-vertical-gradient of horizontal velocities) imposed in the upper boundary, both of these strategies allow a varying wind speed profile across the vertical extension of the MABL. In old seas and neutral atmospheric conditions, we reported the speed-up of the longitudinal velocity up to the upper boundary of the numerical domain [52]. The WI disturbance of the freestream is unlikely to occur in the physical MABL but is intrinsic to the numerical formulations employed in the literature. This is observable, e.g., in [49,51,53,57]. Weakly unstable stratification was considered in [49], and these references also report the slow-down of the freestream velocity in young seas, or due to the presence of curved but stationary bottoms. To keep a constant freestream velocity above a time-varying wave, we develop herein a general strategy to control the wind speed in an atmospheric LES.

In previous studies of wind-wave interactions by LES, the large-scale driving force
(here ∂P/∂x) was usually prescribed as a constant [47–57,60]. As an exception, an ar-
bitrary time-varying force was recently employed to simulate a wind gust event [59].
Over land, however, a few applications considered transient large-scale forcing terms
in LES (c.f. [61–67]). Based on LiDAR measurements in the ABL, the overall accuracy was
improved through the prescription of a linear-time-varying, and spatially homogeneous,
pressure gradient [61]. The measurement of time-varying geostrophic winds drove the
ABL during one or multiple diurnal cycles [62,63]. Employing time-varying and height-
dependent forcing terms is in fact a common practice of data assimilation in meso- to
micro-scale coupling procedures [64,65]. A thorough discussion of the subject was pre-
sented in [66], who propose a dynamic modeler to control the wind speed at an arbitrary
height, in the context of wild-fire applications over canopies. To model wind direction
changes in LES, the governing equations were solved in a non-inertial rotating reference
frame [67].

To accommodate unknown, time-varying momentum stresses in the simulation, the present methodology proposes an original, dynamic $\partial P/\partial x$ modeler. An artificial evolution equation is built for $P_x = \partial P/\partial x$, so that a target velocity is reached at a certain height in the domain. $P_x$ remains spatially homogeneous but varies in time at periods much longer than the resolved motions. Our strategy differs from that of [61,62], because here, $P_x$ is modeled rather than prescribed. Considering here an evolution equation for $P_x$, it is an alternative to the data assimilation procedure presented in [64], or to the periodic adaptation employed in [66]. The evolution equation built is a Partial Differential Equation (PDE), inspired in the rigid-body (Newton) equation of motion, and acting as a dynamic filter to the fluid's (Navier-Stokes) momentum conservation equations.

The methodology employed is described in Section 2. The LES formulation is briefly
described in Section 2.1. The proposed $\partial P / \partial x$ modeler is described in Section 2.2. This
new methodology is compared with a reference case scenario in Section 3. Employing the
dynamic $\partial P / \partial x$ modeler, flat-bottom cases are built with varying wind speeds in Section 4,
ensuring self-similarity of the turbulent statistics, and building the initial conditions for
a varying WA comparison. Monochromatic waves are introduced in Section 5 to evalu-
ate the capability of the proposed methodology to provide consistent and efficient WA
parametrizations of WI disturbances in this canonical condition.

**2. Methods**

2.1. *Numerical Problem Formulation*

A volume domain is defined in the air (specific mass $\rho_\infty$) with sizes ($x_l, y_l, z_l$), bounded below at $S_{z-}$ by the air-water free surface, and above at $S_{z+}$ by the freestream, as exemplified in Figure 1. Wind and waves are aligned and propagate in the longitudinal direction (x axis). The free surface dynamics may be imposed by any single-valued function $\eta(x, y, t)$.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The wind tunnel-like simulation in Cartesian coordinates. Wind and waves are aligned to the x axis. The domain is bounded below by an arbitrary single-valued wave, and above by the freestream. The flow is driven by the large-scale pressure gradient $\partial P / \partial x$.</figcaption></figure>

In the basic configuration, a constant and homogeneous large-scale pressure gradient, $\partial P / \partial x$, drives the flow in the longitudinal (x axis) direction. In a flat-bottom, this gives the theoretical friction velocity $u_0^* = \sqrt{z_l} \partial P / \partial x / \rho_\infty$. We consider neutral stratification conditions, where the flow is defined with characteristic length $L_0$, characteristic velocity $u_0^*$, and the characteristic period $T_0 = L_0 / u_0^*$.

The canonical case of a single non-linear monochromatic wave is considered. Prevent-
ing the free surface evolution under the wind forcing, the wave is prescribed according
to a fifth-order Stokes solution [68]; with wavelength $\lambda = L_0$, wave number $k = 2\pi / \lambda$,
and non-dimensional amplitude $ka = 0.2$. The dispersion equation in deep water gives
a phase velocity $c = \sqrt{g/k}$, and a wave period $T_w = \lambda/c$, with g the gravitational accel-
eration. This strategy establishes the relation between characteristic length and velocity
micro-scale coupling procedures [64,65]. A thorough discussion of the subject was pre-
sented in [66], who propose a dynamic modeler to control the wind speed at an arbitrary
height, in the context of wild-fire applications over canopies. To model wind direction
changes in LES, the governing equations were solved in a non-inertial rotating reference
frame [67].

To accommodate unknown, time-varying momentum stresses in the simulation, the present methodology proposes an original, dynamic $\partial P/\partial x$ modeler. An artificial evolution equation is built for $P_x = \partial P/\partial x$, so that a target velocity is reached at a certain height in the domain. $P_x$ remains spatially homogeneous but varies in time at periods much longer than the resolved motions. Our strategy differs from that of [61,62], because here, $P_x$ is modeled rather than prescribed. Considering here an evolution equation for $P_x$, it is an alternative to the data assimilation procedure presented in [64], or to the periodic adaptation employed in [66]. The evolution equation built is a Partial Differential Equation (PDE), inspired in the rigid-body (Newton) equation of motion, and acting as a dynamic filter to the fluid's (Navier-Stokes) momentum conservation equations.

The methodology employed is described in Section 2. The LES formulation is briefly
described in Section 2.1. The proposed $\partial P / \partial x$ modeler is described in Section 2.2. This
new methodology is compared with a reference case scenario in Section 3. Employing the
dynamic $\partial P / \partial x$ modeler, flat-bottom cases are built with varying wind speeds in Section 4,
ensuring self-similarity of the turbulent statistics, and building the initial conditions for
a varying WA comparison. Monochromatic waves are introduced in Section 5 to evalu-
ate the capability of the proposed methodology to provide consistent and efficient WA
parametrizations of WI disturbances in this canonical condition.

**2. Methods**

2.1. *Numerical Problem Formulation*

A volume domain is defined in the air (specific mass $\rho_\infty$) with sizes ($x_l, y_l, z_l$), bounded below at $S_{z-}$ by the air-water free surface, and above at $S_{z+}$ by the freestream, as exemplified in Figure 1. Wind and waves are aligned and propagate in the longitudinal direction (x axis). The free surface dynamics may be imposed by any single-valued function $\eta(x, y, t)$.

<figure><img src="image_8.png" /><figcaption>**Figure 1.** The wind tunnel-like simulation in Cartesian coordinates. Wind and waves are aligned to the x axis. The domain is bounded below by an arbitrary single-valued wave, and above by the freestream. The flow is driven by the large-scale pressure gradient $\partial P / \partial x$.</figcaption></figure>

In the basic configuration, a constant and homogeneous large-scale pressure gradient, $\partial P / \partial x$, drives the flow in the longitudinal (x axis) direction. In a flat-bottom, this gives the theoretical friction velocity $u_0^* = \sqrt{z_l} \partial P / \partial x / \rho_\infty$. We consider neutral stratification conditions, where the flow is defined with characteristic length $L_0$, characteristic velocity $u_0^*$, and the characteristic period $T_0 = L_0 / u_0^*$.

The canonical case of a single non-linear monochromatic wave is considered. Prevent-
ing the free surface evolution under the wind forcing, the wave is prescribed according
to a fifth-order Stokes solution [68]; with wavelength $\lambda = L_0$, wave number $k = 2\pi / \lambda$,
and non-dimensional amplitude $ka = 0.2$. The dispersion equation in deep water gives
a phase velocity $c = \sqrt{g/k}$, and a wave period $T_w = \lambda/c$, with g the gravitational accel-
eration. This strategy establishes the relation between characteristic length and velocity
scales with nominal WA₀ = c(L₀)/u₀*. The turbulent characteristic time-scale becomes T₀ = λ/u₀* = WA₀ T_w. As an initialization procedure, the wave amplitude evolves linearly, from zero to ka, during 100T_w.

### 2.1.1. Governing Equations and Boundary Conditions

An incompressible and fully turbulent flow is considered, with the Boussinesq approximation for incompressibility (c.f. [69]) acting in the buoyancy terms of the momentum and turbulent kinetic energy (TKE) conservation equations. Coriolis and viscous forces are neglected for the sake of simplicity. The balance equations are filtered according to an LES approach. Let $\tilde{u}(x,t) = (\tilde{u}, \tilde{v}, \tilde{w})$, $\tilde{p}(x,t)$, and $\tilde{\theta}(x,t)$ be the spatially filtered velocity, pressure, and virtual temperature fields, respectively; $\hat{u}$, $\hat{p}$, and $\hat{\theta}$ are the corresponding subgrid-scale (SGS) fields.

The modified pressure is $p^*(x,t) = \tilde{p} - p_\infty + \rho_\infty g(z - z_\infty) + (2\hat{\theta}/3)$, accounting for: A reference value $p_\infty$; the static pressure $\rho_\infty g(z - z_\infty)$, dependent on the total average specific mass $\rho_\infty$, gravitational acceleration $g$, and height $(z - z_\infty)$, relative to a reference $z_\infty$; and the SGS TKE ($\hat{\theta}$) effect $2\hat{\theta}/3$. The (modified) specific mass is $\rho^*(x,t) = -\rho_\infty\beta_\theta(\hat{\theta} - \theta_\infty)$, dependent on the filtered virtual temperature $\tilde{\theta}$, the reference values [$\rho_\infty, \theta_\infty$], and the buoyancy parameter $\beta_\theta = g/\theta_\infty$. In neutral stratification, buoyancy forces are suppressed so $\rho^* = 0$. The filtered strain rate tensor is $\tilde{S} = [\nabla \tilde{u} + (\nabla \tilde{u})^T]/2$. The SGS-modeled shear stress tensor, $\tau^{SGS} = -2v_t \tilde{S}$, is defined within the scope of the eddy viscosity hypothesis, and is dependent on the turbulent kinematic viscosity $v_t$. The Kronecker delta is denoted by $\delta_{ij}$.

Mass and momentum balances are written as Equations (1) and (2). The SGS TKE is modeled by the Deardorff single Equation (3) [70].

$$ \frac{\partial \tilde{u}_j}{\partial x_j} = 0, \qquad (1) $$

$$ \rho_{\infty} \left[ \frac{\partial \tilde{u}_i}{\partial t} + \frac{\partial (\tilde{u}_j \tilde{u}_i)}{\partial x_j} \right] = - \rho^* \delta_{i3} - \frac{\partial p^*}{\partial x_i} - \rho_{\infty} \frac{\partial \tau_{ij}^{SGS}}{\partial x_j} - \frac{\partial \mathcal{P}}{\partial x_i} \delta_{i1}, \quad \text{and} \qquad (2) $$

$$ \frac{\partial \hat{e}}{\partial t} + \frac{\partial (\tilde{u}_i \hat{e})}{\partial x_i} = (2v_t \tilde{S}_{ij}) \tilde{S}_{ij} - \frac{g}{\theta_\infty} v_h \frac{\partial \tilde{\theta}}{\partial z} + \frac{\partial}{\partial x_j} \left( 2v_t \frac{\partial \hat{e}}{\partial x_j} \right) - \epsilon. \qquad (3) $$

The turbulent dissipation, $\epsilon = c_\epsilon \hat{e}^{3/2}/\Delta_f$, is determined according to the filter length scale $\Delta_f = [(3/2)^2 \Delta x \Delta y \Delta z]^{1/3}$. The turbulent kinematic viscosity and diffusivity are respectively $v_t = c_k l \hat{e}^{1/2}$ and $v_h = (1 + 2l/\Delta_f) v_t$, where l is here equal to $\Delta_f$.

The Cartesian physical space $(x, y, z)$ is mapped into the computational, terrain-following space, defined with $\zeta_1 = x$, $\zeta_2 = y$, and $\zeta_3 = \zeta_3(x, y, z, t)$. The governing Equations (1)–(3) are rewritten in the terrain-following space as demonstrated, e.g., in [14,41]. The fundamental unknowns $(\tilde{u}, \tilde{p}, \hat{e}, \hat{\theta})$ are located at the cell center. The full set of equations is presented in their discrete form and in the terrain-following space, and with the values of the constants given, e.g., in [48].

The geometrical conservation law introduced by [71], translates into Equation (4):

$$ \frac{\partial}{\partial t} \left( \frac{1}{J} \right) = \frac{\partial}{\partial \zeta_3} \left( \frac{\partial h}{\partial t} \right), \qquad (4) $$

relating the Jacobian of the grid transformation J, to the grid vertical position h.

A spectral discretization is adopted, so periodic Boundary Conditions (BCs) apply to the transversal and longitudinal boundaries.

### Upper Boundary Conditions

No-penetration and free-slip BCs are imposed as:
scales with nominal WA₀ = c(L₀)/u₀*. The turbulent characteristic time-scale becomes T₀ = λ/u₀* = WA₀ T_w. As an initialization procedure, the wave amplitude evolves linearly, from zero to ka, during 100T_w.

### 2.1.1. Governing Equations and Boundary Conditions

An incompressible and fully turbulent flow is considered, with the Boussinesq approximation for incompressibility (c.f. [69]) acting in the buoyancy terms of the momentum and turbulent kinetic energy (TKE) conservation equations. Coriolis and viscous forces are neglected for the sake of simplicity. The balance equations are filtered according to an LES approach. Let $\tilde{u}(x,t) = (\tilde{u}, \tilde{v}, \tilde{w})$, $\tilde{p}(x,t)$, and $\tilde{\theta}(x,t)$ be the spatially filtered velocity, pressure, and virtual temperature fields, respectively; $\hat{u}$, $\hat{p}$, and $\hat{\theta}$ are the corresponding subgrid-scale (SGS) fields.

The modified pressure is $p^*(x,t) = \tilde{p} - p_\infty + \rho_\infty g(z - z_\infty) + (2\hat{\epsilon}/3)$, accounting for: A reference value $p_\infty$; the static pressure $\rho_\infty g(z - z_\infty)$, dependent on the total average specific mass $\rho_\infty$, gravitational acceleration $g$, and height $(z - z_\infty)$, relative to a reference $z_\infty$; and the SGS TKE ($\hat{\epsilon}$) effect $2\hat{\epsilon}/3$. The (modified) specific mass is $\rho^*(x,t) = -\rho_\infty\beta_\theta(\hat{\theta} - \theta_\infty)$, dependent on the filtered virtual temperature $\tilde{\theta}$, the reference values [$\rho_\infty, \theta_\infty$], and the buoyancy parameter $\beta_\theta = g/\theta_\infty$. In neutral stratification, buoyancy forces are suppressed so $\rho^* = 0$. The filtered strain rate tensor is $\tilde{S} = [\nabla \tilde{u} + (\nabla \tilde{u})^T]/2$. The SGS-modeled shear stress tensor, $\tau^{SGS} = -2\nu_t \tilde{S}$, is defined within the scope of the eddy viscosity hypothesis, and is dependent on the turbulent kinematic viscosity $\nu_t$. The Kronecker delta is denoted by $\delta_{ij}$.

Mass and momentum balances are written as Equations (1) and (2). The SGS TKE is modeled by the Deardorff single Equation (3) [70].

$$ \frac{\partial \tilde{u}_j}{\partial x_j} = 0, \qquad (1) $$

$$ \rho_{\infty} \left[ \frac{\partial \tilde{u}_i}{\partial t} + \frac{\partial (\tilde{u}_j \tilde{u}_i)}{\partial x_j} \right] = -\rho^* \delta_{i3} - \frac{\partial p^*}{\partial x_i} - \rho_{\infty} \frac{\partial \tau_{ij}^{SGS}}{\partial x_j} - \frac{\partial \mathcal{P}}{\partial x_i} \delta_{i1}, \quad \text{and} \qquad (2) $$

$$ \frac{\partial \hat{\epsilon}}{\partial t} + \frac{\partial (\tilde{u}_i \hat{\epsilon})}{\partial x_i} = (2\nu_t \tilde{S}_{ij}) \tilde{S}_{ij} - \frac{g}{\theta_\infty} \nu_h \frac{\partial \tilde{\theta}}{\partial z} + \frac{\partial}{\partial x_j} \left( 2\nu_t \frac{\partial \hat{\epsilon}}{\partial x_j} \right) - \epsilon. \qquad (3) $$

The turbulent dissipation, $\epsilon = c_\epsilon \hat{\epsilon}^{3/2}/\Delta_f$, is determined according to the filter length scale $\Delta_f = [(3/2)^2 \Delta x \Delta y \Delta z]^{1/3}$. The turbulent kinematic viscosity and diffusivity are respectively $\nu_t = c_k l \hat{\epsilon}^{1/2}$ and $\nu_h = (1 + 2l/\Delta_f)\nu_t$, where l is here equal to $\Delta_f$.

The Cartesian physical space $(x, y, z)$ is mapped into the computational, terrain-following space, defined with $\zeta_1 = x$, $\zeta_2 = y$, and $\zeta_3 = \zeta_3(x, y, z, t)$. The governing Equations (1)–(3) are rewritten in the terrain-following space as demonstrated, e.g., in [14,41]. The fundamental unknowns $(\tilde{u}, \tilde{p}, \hat{\epsilon}, \hat{\theta})$ are located at the cell center. The full set of equations is presented in their discrete form and in the terrain-following space, and with the values of the constants given, e.g., in [48].

The geometrical conservation law introduced by [71], translates into Equation (4):

$$ \frac{\partial}{\partial t} \left( \frac{1}{J} \right) = \frac{\partial}{\partial \zeta_3} \left( \frac{\partial h}{\partial t} \right), \qquad (4) $$

relating the Jacobian of the grid transformation J, to the grid vertical position h.

A spectral discretization is adopted, so periodic Boundary Conditions (BCs) apply to the transversal and longitudinal boundaries.

### Upper Boundary Conditions

No-penetration and free-slip BCs are imposed as:
$$ \frac{\partial \tilde{\theta}}{\partial z} = 0, \quad \frac{\partial \tilde{\epsilon}}{\partial z} = 0, \quad \frac{\partial \tilde{u}_{\xi_1}}{\partial z} = 0, \quad \frac{\partial \tilde{u}_{\xi_2}}{\partial z} = 0, \quad \text{and} \quad \tilde{u}_{\xi_3} = 0, \quad \text{at} \quad S_{z+}. \qquad (5) $$

### Lower Boundary Conditions

The difference in the tangential velocities, $\Delta \tilde{u}_{[\xi_1, \xi_2]}$, is computed from the velocities resolved in the first cell, minus the ones prescribed at the lower surface. The wall fluxes $\tau_f^w$ at $[\xi_1, \xi_2]$ (at $S_{z-}$) are then obtained from the log law wall function:

$$ \Delta \tilde{u}_{[\xi_1, \xi_2]} = \frac{u^*_{[\xi_1, \xi_2]}}{\kappa} \ln\left(\frac{z}{z_0}\right), \quad \text{at } S_{z-}, \qquad (6) $$

specifying a Robin type of BC for the longitudinal momentum equations, with a roughness length $z_0 = 10^{-6} L_0$, and a von Kármán constant $\kappa = 0.4$.

The vertical velocity is imposed in the first cell as equal to the surface. The SGS TKE flux is null, and the turbulent kinematic viscosity is given by $v_t = |u^*|\kappa z$. The temperature can be specified either by its surface value or by the heat flux (here null). When the heat flux is other than null, Equation (6) is adapted by employing MOST [30].

$$ \Delta \tilde{u}_{\xi_3} = 0, \quad \frac{\partial \tilde{\theta}}{\partial z} = 0, \quad \text{and} \quad \frac{\partial \tilde{\epsilon}}{\partial z} = 0, \quad \text{at } S_{z-}. \qquad (7) $$

#### 2.1.2. Post-Processing

##### Surface Drag and Friction Velocity

The stresses are integrated and averaged at $S_{z-}$ to provide the bottom surface drag as a force per surface unit. The total drag is $F_t = F_p + F_f$, a combination of normal and tangential stresses. In a wavy surface, the projection of the modified pressure $p^*$ normal to the free surface $\eta(x, y, t)$ leads to a pressure stress vector of horizontal components $\tau_p^w = p^* \cdot (\partial\eta/\partial x_i)$, with $i = 1, 2$. The form drag is the magnitude ($| \cdot |$) of the space average of this vector $F_p = |\langle \tau_p^w \rangle|$. The friction drag is the magnitude of the averaged wall flux vector $F_f = |\langle \tau_f^w \rangle|$.

Equation (6) defines the point-wise wall flux and friction velocity vectors ($\tau_f^w = \rho_\infty u^{*2}$). The space-averaged friction velocity (scalar) is afterward obtained from the friction drag, as:

$$ u^* = \sqrt{F_f / \rho_\infty}, \qquad (8) $$

and $u^*$ is simply referred to as the (resolved) friction velocity hereafter.

##### Averages, Fluctuations, and Fluxes

Horizontal space averages $\langle \cdot \rangle$ are performed in the computational Cartesian grid, along $[\xi_1, \xi_2]$ planes of constant $\xi_3$. The space-averaged quantities are thus presented as a function of the averaged height $z = \langle h \rangle$, which remains constant for a given plane. Turbulent fluctuations of a general (filtered) field variable $\tilde{\phi}(x, t)$ are denoted and computed as $\phi' = \tilde{\phi} - \langle \tilde{\phi} \rangle$. The resolved TKE is $\tilde{e} = 0.5(u'u' + v'v' + w'w')$. Sliding time averages ($\cdot$) are afterward performed between $t - t_{avg}$ and $t$, with an average period of $t_{avg}$.

The wavy bottom surface invariably requires the problem to be formulated in surface- following coordinates, which leads to the pressure acting as a momentum flux mechanism in the Navier–Stokes equations [14]. Therefore, the total flux (stress) becomes a combination between pressure and friction (resolved and SGS) terms, written as $\tau_{13} = \tau_{13}^p + \tau_{13}^u + \tau_{13}^{SGS}$. The pressure term is $\tau_{13}^p = p^* \cdot (\partial h / \partial x_1)$, and equals the pressure drag ($F_p$) at $S_{z-}$, where $h = \eta(x, y, t)$. The resolved friction term is $\tau_{13}^u = u'w'$, and the SGS term is $\tau_{13}^{SGS} = -2v_t \tilde{S}_{13}$.
$$ \frac{\partial \tilde{\theta}}{\partial z} = 0, \quad \frac{\partial \tilde{\epsilon}}{\partial z} = 0, \quad \frac{\partial \tilde{u}_{\xi_1}}{\partial z} = 0, \quad \frac{\partial \tilde{u}_{\xi_2}}{\partial z} = 0, \quad \text{and} \quad \tilde{u}_{\xi_3} = 0, \quad \text{at } S_{z+}. \qquad (5) $$

### Lower Boundary Conditions

The difference in the tangential velocities, $\Delta \tilde{u}_{[\xi_1, \xi_2]}$, is computed from the velocities resolved in the first cell, minus the ones prescribed at the lower surface. The wall fluxes $\tau_f^w$ at $[\xi_1, \xi_2]$ (at $S_{z-}$) are then obtained from the log law wall function:

$$ \Delta \tilde{u}_{[\xi_1, \xi_2]} = \frac{u^*_{[\xi_1, \xi_2]}}{\kappa} \ln\left(\frac{z}{z_0}\right), \quad \text{at } S_{z-}, \qquad (6) $$

specifying a Robin type of BC for the longitudinal momentum equations, with a roughness length $z_0 = 10^{-6} L_0$, and a von Kármán constant $\kappa = 0.4$.

The vertical velocity is imposed in the first cell as equal to the surface. The SGS TKE flux is null, and the turbulent kinematic viscosity is given by $v_t = |u^*|\kappa z$. The temperature can be specified either by its surface value or by the heat flux (here null). When the heat flux is other than null, Equation (6) is adapted by employing MOST [30].

$$ \Delta \tilde{u}_{\xi_3} = 0, \quad \frac{\partial \tilde{\theta}}{\partial z} = 0, \quad \text{and} \quad \frac{\partial \tilde{\epsilon}}{\partial z} = 0, \quad \text{at } S_{z-}. \qquad (7) $$

### 2.1.2. Post-Processing

#### Surface Drag and Friction Velocity

The stresses are integrated and averaged at $S_{z-}$ to provide the bottom surface drag as a force per surface unit. The total drag is $F_t = F_p + F_f$, a combination of normal and tangential stresses. In a wavy surface, the projection of the modified pressure $p^*$ normal to the free surface $\eta(x, y, t)$ leads to a pressure stress vector of horizontal components $\tau_p^w = p^* \cdot (\partial\eta/\partial x_i)$, with $i = 1, 2$. The form drag is the magnitude ($| \cdot |$) of the space average of this vector $F_p = |\langle \tau_p^w \rangle|$. The friction drag is the magnitude of the averaged wall flux vector $F_f = |\langle \tau_f^w \rangle|$.

Equation (6) defines the point-wise wall flux and friction velocity vectors ($\tau_f^w = \rho_\infty u^{*2}$). The space-averaged friction velocity (scalar) is afterward obtained from the friction drag, as:

$$ u^* = \sqrt{F_f / \rho_\infty}, \qquad (8) $$

and $u^*$ is simply referred to as the (resolved) friction velocity hereafter.

#### Averages, Fluctuations, and Fluxes

Horizontal space averages $\langle \cdot \rangle$ are performed in the computational Cartesian grid, along $[\xi_1, \xi_2]$ planes of constant $\xi_3$. The space-averaged quantities are thus presented as a function of the averaged height $z = \langle h \rangle$, which remains constant for a given plane. Turbulent fluctuations of a general (filtered) field variable $\tilde{\phi}(x, t)$ are denoted and computed as $\phi' = \tilde{\phi} - \langle \tilde{\phi} \rangle$. The resolved TKE is $\tilde{e} = 0.5(u'u' + v'v' + w'w')$. Sliding time averages ($\cdot$) are afterward performed between $t - t_{avg}$ and $t$, with an average period of $t_{avg}$.

The wavy bottom surface invariably requires the problem to be formulated in surface- following coordinates, which leads to the pressure acting as a momentum flux mechanism in the Navier-Stokes equations [14]. Therefore, the total flux (stress) becomes a combination between pressure and friction (resolved and SGS) terms, written as $\tau_{13} = \tau_{13}^p + \tau_{13}^u + \tau_{13}^{SGS}$. The pressure term is $\tau_{13}^p = p^* \cdot (\partial h / \partial x_1)$, and equals the pressure drag ($F_p$) at $S_{z-}$, where $h = \eta(x, y, t)$. The resolved friction term is $\tau_{13}^u = u'w'$, and the SGS term is $\tau_{13}^{SGS} = -2v_t \tilde{S}_{13}$.
### 2.1.3. Numerical Strategy

#### Discretization

Given the characteristic length $L_0$, the presented domain has sizes $(x_l, y_l, z_l) = (4, 2, 5)L_0$. The grid is structured. The domain is discretized with $(n_x, n_y, n_z) = (256, 128, 80)$, giving a total of 655,360 cells, with sizes $(\Delta x, \Delta y, \Delta z)$. The cells are uniformly distributed in $(x, y)$, but not in $z$. The first grid size in the $z$ direction is $z_{w1} = 0.5208 \times 10^{-2}L_0$, so that $z_{w1}/\Delta x = 1/3$, and it grows with a constant ratio of 1.05. Assuming for example that $\nu = 1.5 \times 10^{-5} \text{ m}^2/\text{s}$ (air at 15 °C), when $WA_0 = 60$ and $\lambda = 100 \text{ m}$, then $u^* = 0.210 \text{ m/s}$, and $z^+ = z_{w1} u^*/\nu = 7.3 \times 10^3$.

The spectral discretization applies to $(\xi_1, \xi_2)$ directions, and a second-order finite difference scheme applies to the $\xi_3$ direction. A third-order Runge-Kutta time-stepping scheme admits arbitrarily selected time steps $(\Delta t)$. An adaptive $\Delta t$ procedure ensures that $CFL = \Delta t \cdot \max[\tilde{u}/\Delta x, \tilde{v}/\Delta y, \tilde{w}/\Delta z] = 0.5$.

#### Initialization

We followed a standard procedure for the initialization of the neutral LES-generated MABL. First, the flow field was built from the mean theoretical solution for a turbulent boundary layer over a flat plate, superposed to artificial, randomly generated turbulent motions. In order to first generate resolved turbulence, we considered buoyant effects on the momentum conservation equations during the first thousand time steps. The test cases presented here were initialized from a converged case of constant $\partial \mathcal{P}/\partial x$, giving a nominal friction velocity $u_0^*$. The undisturbed, space-averaged (in $\xi_1$ and $\xi_2$) freestream velocity $\langle \tilde{u}_h \rangle$ was then resolved at $S_{z+}$, thus giving the reference freestream velocity $u_{ref0} = 52.11 u_0^*$.

### 2.2. A Large-Scale Dynamically Evolving Pressure Gradient Modeler

Isolating the local time derivative and the large-scale forcing in the first component (i = 1) of Equation (2), the other terms of the right-hand side are gathered in $\mathcal{R}^{[P]}$. Thus, the large-scale wind forcing is achieved in such an LES formulation through the longitudinal momentum balance equation:

$$ \rho_{\infty} \frac{\partial \tilde{u}}{\partial t} = \mathcal{R}^{[P]} - \frac{\partial \mathcal{P}}{\partial x}, \qquad (9) $$

which is evaluated at each time-step and Runge-Kutta iteration $n$.

The evaluation of Equation (9) leads to a non-zero, space-averaged local time derivative on the upper boundary $\partial \langle \tilde{u}_{h0} \rangle / \partial t$. As a motivation, and in a first attempt to keep a constant freestream velocity, we may consider the correction with an increment $\delta[\partial \mathcal{P}/\partial x]$:

$$ \rho_{\infty} \frac{\partial \langle \tilde{u}_{h1} \rangle}{\partial t} = \rho_{\infty} \frac{\partial \langle \tilde{u}_{h0} \rangle}{\partial t} - \delta \left[ \frac{\partial \mathcal{P}}{\partial x} \right], \qquad (10) $$

which is required to be null, leading to the iterative procedure:

$$ \overbrace{\left[ \frac{\partial \mathcal{P}}{\partial x} \right]_n}^{\mathcal{P}_{x,n}} = \overbrace{\left[ \frac{\partial \mathcal{P}}{\partial x} \right]_{n-1}}^{\mathcal{P}_{x,n-1}} + \delta \overbrace{\left[ \frac{\partial \mathcal{P}}{\partial x} \right]_n}^{\delta \mathcal{P}_{x,n}}, \quad \text{with} \quad \delta \left[ \frac{\partial \mathcal{P}}{\partial x} \right]_n = \rho_{\infty} \frac{\partial \langle \tilde{u}_{h0} \rangle}{\partial t} \bigg|_n. \qquad (11) $$

Through the adaptive procedure of Equations (10) and (11), $\langle \tilde{u}_h \rangle$ is constant at every instant, and $\partial \mathcal{P}/\partial x$ incorporates any turbulent disturbance in the upper boundary. This has to be avoided since by definition, $\partial \mathcal{P}/\partial x$ contains only large-scale dynamics. Instead, an evolution equation was built for $\mathcal{P}_x = \partial \mathcal{P}/\partial x$, in which it is dynamically filtered responding to the large-scale part of these equations.
### 2.1.3. Numerical Strategy

#### Discretization

Given the characteristic length $L_0$, the presented domain has sizes $(x_l, y_l, z_l) = (4, 2, 5)L_0$. The grid is structured. The domain is discretized with $(n_x, n_y, n_z) = (256, 128, 80)$, giving a total of 655,360 cells, with sizes $(\Delta x, \Delta y, \Delta z)$. The cells are uniformly distributed in $(x, y)$, but not in $z$. The first grid size in the $z$ direction is $z_{w1} = 0.5208 \times 10^{-2}L_0$, so that $z_{w1}/\Delta x = 1/3$, and it grows with a constant ratio of 1.05. Assuming for example that $\nu = 1.5 \times 10^{-5} \text{ m}^2/\text{s}$ (air at 15 °C), when $WA_0 = 60$ and $\lambda = 100 \text{ m}$, then $u^* = 0.210 \text{ m/s}$, and $z^+ = z_{w1} u^*/\nu = 7.3 \times 10^3$.

The spectral discretization applies to $(\xi_1, \xi_2)$ directions, and a second-order finite difference scheme applies to the $\xi_3$ direction. A third-order Runge-Kutta time-stepping scheme admits arbitrarily selected time steps $(\Delta t)$. An adaptive $\Delta t$ procedure ensures that $CFL = \Delta t \cdot \max[\tilde{u}/\Delta x, \tilde{v}/\Delta y, \tilde{w}/\Delta z] = 0.5$.

#### Initialization

We followed a standard procedure for the initialization of the neutral LES-generated MABL. First, the flow field was built from the mean theoretical solution for a turbulent boundary layer over a flat plate, superposed to artificial, randomly generated turbulent motions. In order to first generate resolved turbulence, we considered buoyant effects on the momentum conservation equations during the first thousand time steps. The test cases presented here were initialized from a converged case of constant $\partial \mathcal{P}/\partial x$, giving a nominal friction velocity $u_0^*$. The undisturbed, space-averaged (in $\xi_1$ and $\xi_2$) freestream velocity $\langle \tilde{u}_h \rangle$ was then resolved at $S_{z+}$, thus giving the reference freestream velocity $u_{ref0} = 52.11 u_0^*$.

### 2.2. A Large-Scale Dynamically Evolving Pressure Gradient Modeler

Isolating the local time derivative and the large-scale forcing in the first component (i = 1) of Equation (2), the other terms of the right-hand side are gathered in $\mathcal{R}^{[P]}$. Thus, the large-scale wind forcing is achieved in such an LES formulation through the longitudinal momentum balance equation:

$$ \rho_{\infty} \frac{\partial \tilde{u}}{\partial t} = \mathcal{R}^{[P]} - \frac{\partial \mathcal{P}}{\partial x}, \qquad (9) $$

which is evaluated at each time-step and Runge-Kutta iteration $n$.

The evaluation of Equation (9) leads to a non-zero, space-averaged local time derivative on the upper boundary $\partial \langle \tilde{u}_{h0} \rangle / \partial t$. As a motivation, and in a first attempt to keep a constant freestream velocity, we may consider the correction with an increment $\delta[\partial \mathcal{P}/\partial x]$:

$$ \rho_{\infty} \frac{\partial \langle \tilde{u}_{h1} \rangle}{\partial t} = \rho_{\infty} \frac{\partial \langle \tilde{u}_{h0} \rangle}{\partial t} - \delta \left[ \frac{\partial \mathcal{P}}{\partial x} \right], \qquad (10) $$

which is required to be null, leading to the iterative procedure:

$$ \overbrace{\left[ \frac{\partial \mathcal{P}}{\partial x} \right]}^{\mathcal{P}_{x,n}} = \overbrace{\left[ \frac{\partial \mathcal{P}}{\partial x} \right]}^{\mathcal{P}_{x,n-1}} + \delta \overbrace{\left[ \frac{\partial \mathcal{P}}{\partial x} \right]}^{\delta \mathcal{P}_{x,n}}, \quad \text{with} \quad \delta \left[ \frac{\partial \mathcal{P}}{\partial x} \right]_n = \rho_{\infty} \frac{\partial \langle \tilde{u}_{h0} \rangle}{\partial t} \bigg|_n. \qquad (11) $$

Through the adaptive procedure of Equations (10) and (11), $\langle \tilde{u}_h \rangle$ is constant at every instant, and $\partial \mathcal{P}/\partial x$ incorporates any turbulent disturbance in the upper boundary. This has to be avoided since by definition, $\partial \mathcal{P}/\partial x$ contains only large-scale dynamics. Instead, an evolution equation was built for $\mathcal{P}_x = \partial \mathcal{P}/\partial x$, in which it is dynamically filtered responding to the large-scale part of these equations.
The large-scale filter of Equations (10) and (11) is accomplished by a mass-spring-damper system represented in Figure 2. This determines the evolution of $P_x$ according to a second-order PDE, which is also known as the rigid body equation of motion [72]. Denoting the first and second-time derivatives of $P_x$ by $\dot{P}_x$ and $\ddot{P}_x$, respectively:

$$ \ddot{P}_x = -w_p^2 \delta P_x - 2\zeta w_p \dot{P}_x, \quad (12) $$

which advances with the same Runge–Kutta scheme applied to the other quantities in the LES code. The displacement from the equilibrium state is given by $\delta P_x$. The system is characterized by a natural angular-frequency $w_p$ (with natural period $T_p = 2\pi/w_p$), and a non-dimensional damping ratio $\zeta$.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The Mass-Spring-Damper system ruling the evolution of the longitudinal pressure gradient $P_x = \partial P / \partial x$, according to Equation (12).</figcaption></figure>

A property of Equation (12) is that $P_x$ is not excited by frequencies much higher than the natural, and another is that stability can be assured in the so-called critically damped regime [72,73]. Therefore, in order to ensure that the response is limited to large periods, we employ a large parametrization of $T_p$, and to preserve stability, we impose the critical damping ratio $\zeta$:

$$ T_p = 10^3 \frac{x_l}{u_{ref}}; \quad \text{and} \quad \zeta = 1. \quad (13) $$

Furthermore, we aim for a certain, non-null local time derivative with:

$$ \delta P_x = \rho_\infty \left[ \frac{\partial \langle \tilde{u}_{h0} \rangle}{\partial t} - \frac{\partial \langle \tilde{u}_{h1} \rangle}{\partial t} \right], \quad \text{and} \quad \frac{\partial \langle \tilde{u}_{h1} \rangle}{\partial t} = \frac{u_{ref} - \langle \tilde{u}_h \rangle}{T_p}, \quad (14) $$

so that seeking the equilibrium, $\langle \tilde{u}_h \rangle$ tends to the target $u_{ref}$ over the period $T_p$.

The pressure gradient thus defined is spatially homogeneous but varies instantaneously in time, so that the mean longitudinal velocity in the upper surface tends to the target. The methodology proposed can be easily generalized by associating the subscript $\cdot_h$ to an arbitrary height in the domain.

### 3. Pilot Application

The proposed methodology was built to mitigate the undesired variations of freestream velocities due to WI disturbances in LES (Section 1). Therefore to evaluate its performance, the dynamic $P_x$ modeler is here compared with the traditional, constant $P_x$ approach.

#### Test Cases

Three cases are presented, initializing with freestream velocity resolved as $u_{ref0} = 52.11 u_0^*$ (see initialization in Section 2.1.3). Case 01 is the reference, with flat-bottom BC and a constant $P_x = \rho_\infty u_0^2 / z_l$. Case 02 introduces the monochromatic wave (WA0 = 60), still keeping a constant $P_x$. Case 03 considers this same wave, but also the dynamic $P_x$ modeler here proposed, with target velocity ($u_{ref}$) equal to the initial $u_{ref0}$.
The large-scale filter of Equations (10) and (11) is accomplished by a mass-spring-damper system represented in Figure 2. This determines the evolution of $P_x$ according to a second-order PDE, which is also known as the rigid body equation of motion [72]. Denoting the first and second-time derivatives of $P_x$ by $\dot{P}_x$ and $\ddot{P}_x$, respectively:

$$ \ddot{P}_x = -w_p^2 \delta P_x - 2\zeta w_p \dot{P}_x, \qquad (12) $$

which advances with the same Runge-Kutta scheme applied to the other quantities in the LES code. The displacement from the equilibrium state is given by $\delta P_x$. The system is characterized by a natural angular-frequency $w_p$ (with natural period $T_p = 2\pi/w_p$), and a non-dimensional damping ratio $\zeta$.

<figure><img src="image_9.png" /><figcaption>Figure 2. The Mass-Spring-Damper system ruling the evolution of the longitudinal pressure gradient $P_x = \partial P / \partial x$, according to Equation (12).</figcaption></figure>

A property of Equation (12) is that $P_x$ is not excited by frequencies much higher than the natural, and another is that stability can be assured in the so-called critically damped regime [72,73]. Therefore, in order to ensure that the response is limited to large periods, we employ a large parametrization of $T_p$, and to preserve stability, we impose the critical damping ratio $\zeta$:

$$ T_p = 10^3 \frac{x_l}{u_{ref}}; \quad \text{and} \quad \zeta = 1. \qquad (13) $$

Furthermore, we aim for a certain, non-null local time derivative with:

$$ \delta P_x = \rho_\infty \left[ \frac{\partial \langle \tilde{u}_{h0} \rangle}{\partial t} - \frac{\partial \langle \tilde{u}_{h1} \rangle}{\partial t} \right], \quad \text{and} \quad \frac{\partial \langle \tilde{u}_{h1} \rangle}{\partial t} = \frac{u_{ref} - \langle \tilde{u}_h \rangle}{T_p}, \qquad (14) $$

so that seeking the equilibrium, $\langle \tilde{u}_h \rangle$ tends to the target $u_{ref}$ over the period $T_p$.

The pressure gradient thus defined is spatially homogeneous but varies instantaneously in time, so that the mean longitudinal velocity in the upper surface tends to the target. The methodology proposed can be easily generalized by associating the subscript $\cdot_h$ to an arbitrary height in the domain.

### 3. Pilot Application

The proposed methodology was built to mitigate the undesired variations of freestream velocities due to WI disturbances in LES (Section 1). Therefore to evaluate its performance, the dynamic $P_x$ modeler is here compared with the traditional, constant $P_x$ approach.

### Test Cases

Three cases are presented, initializing with freestream velocity resolved as $u_{ref0} = 52.11 u_0^*$ (see initialization in Section 2.1.3). Case 01 is the reference, with flat-bottom BC and a constant $P_x = \rho_\infty u_0^2 / z_l$. Case 02 introduces the monochromatic wave ($WA_0 = 60$), still keeping a constant $P_x$. Case 03 considers this same wave, but also the dynamic $P_x$ modeler here proposed, with target velocity ($u_{ref}$) equal to the initial $u_{ref0}$.
Results

Time histories are shown in Figure 3, by employing sliding-time statistics with long average periods $t_{avg} = 75T_0$. The normalized mean wind velocity in the last cell layer ($U_h/u_0^*$, with $U_h = |\langle \tilde{u}_h \rangle|$) is shown in Figure 3a. In the flat-bottom case (Case 01), the freestream velocity deviates ~2% from its initial value. The wave acts as an additional driving force dragging the wind above. Therefore, without a diminishing $P_x$ balancing the additional momentum transferred to the wind, the freestream velocity is greatly augmented by the introduction of the wave in Case 02, reaching values ~14% higher than $u_{ref}$. In Case 03, the dynamic $P_x$ algorithm is activated, and the target velocity $u_{ref}$ is shown as the horizontal line in the figure. With this approach, the freestream velocity deviates only ~1% from the target, when $t > 200T_0$ in Case 03.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time histories of non-dimensional: (a) wind speed at $S_{z+}$; (b) friction velocity; (c) large-scale pressure gradient; and (d) form drag. Averaged with $t_{avg} = 75T_0$. Case 01: Flat bottom with constant $P_x$, as the reference. Case 02: Wavy bottom with constant $P_x$. Case 03: Wavy bottom with the dynamic $P_x$ modeler.</figcaption></figure>

The resolved friction velocity ($\overline{u^*}$) is normalized by its nominal value in Figure 3b. Concerning Case 02, the wave first leads to a rapid drop in the friction velocity, which later returns towards $u_0^*$ over large periods. Employing the dynamic $P_x$ in Case 03, the friction velocity does not recover the nominal value but is kept at ~0.88 $u_0^*$.

Both the freestream and friction velocities of Figure 3a,b reveal a persistent transient regime for Case 02, which in fact has not completely converged after the 2 million time steps, and during the period of $250T_0$ shown in the figures. The transient regime is significantly mitigated in Case 03 with the dynamic $P_x$ method, which presents relatively slight velocity variations over the whole computation, reaching a statistically steady regime for $t > 200T_0$.

The homogeneous pressure gradient is shown in Figure 3c: Constant and equal to $1/z_l$ for Cases 01 and 02; dynamic and slowly evolving in Case 03, which converges towards $P_x \sim 0.6/z_l$ for $t > 200T_0$. The transient $P_x$ in Case 03 is responsible for keeping the quasi-steady wind and friction velocities.

Neglecting viscous effects, the WI stress integrates into the form drag in the lower surface, and is responsible for the vertical transfer of momentum between waves and wind [14,15]. If the waves travel sufficiently fast (large WA) with a significant wave slope, the upward transfer is such that the wave signature emerges as a low-level jet, characterized by a local maximum in the wind velocity profiles [16–18].

The momentum transfer between the wave and the wind is evaluated in the non-dimensional form drag (per surface unit) $\overline{F}_p/\overline{F}_t$, which also exhibits a transient behavior in Figure 3d. While Case 03 converges towards $\overline{F}_p/\overline{F}_t \sim -15\%$ for $t > 200T_0$, Case 02 continues to increase after its sudden drop and minimum occurring at $t \sim 10T_0$. Moreover in Case 03, $\overline{F}_p < 0$ and the momentum transfers exclusively upwards, but for Case 02, the form drag changes the sign, and upon becoming positive, it indicates a downward (wind-to-wave) momentum transfer for $t > 100T_0$.
Results

Time histories are shown in Figure 3, by employing sliding-time statistics with long average periods $t_{avg} = 75T_0$. The normalized mean wind velocity in the last cell layer ($U_h/u_0^*$, with $U_h = |\langle \tilde{u}_h \rangle|$) is shown in Figure 3a. In the flat-bottom case (Case 01), the freestream velocity deviates ~2% from its initial value. The wave acts as an additional driving force dragging the wind above. Therefore, without a diminishing $P_x$ balancing the additional momentum transferred to the wind, the freestream velocity is greatly augmented by the introduction of the wave in Case 02, reaching values ~14% higher than $u_{ref}$. In Case 03, the dynamic $P_x$ algorithm is activated, and the target velocity $u_{ref}$ is shown as the horizontal line in the figure. With this approach, the freestream velocity deviates only ~1% from the target, when $t > 200T_0$ in Case 03.

<figure><img src="image_10.png" /><figcaption>**Figure 3.** Time histories of non-dimensional: (a) wind speed at $S_{z+}$; (b) friction velocity; (c) large-scale pressure gradient; and (d) form drag. Averaged with $t_{avg} = 75T_0$. Case 01: Flat bottom with constant $P_x$, as the reference. Case 02: Wavy bottom with constant $P_x$. Case 03: Wavy bottom with the dynamic $P_x$ modeler.</figcaption></figure>

The resolved friction velocity ($\overline{u^*}$) is normalized by its nominal value in Figure 3b. Concerning Case 02, the wave first leads to a rapid drop in the friction velocity, which later returns towards $u_0^*$ over large periods. Employing the dynamic $P_x$ in Case 03, the friction velocity does not recover the nominal value but is kept at ~0.88 $u_0^*$.

Both the freestream and friction velocities of Figure 3a,b reveal a persistent transient regime for Case 02, which in fact has not completely converged after the 2 million time steps, and during the period of $250T_0$ shown in the figures. The transient regime is significantly mitigated in Case 03 with the dynamic $P_x$ method, which presents relatively slight velocity variations over the whole computation, reaching a statistically steady regime for $t > 200T_0$.

The homogeneous pressure gradient is shown in Figure 3c: Constant and equal to $1/z_l$ for Cases 01 and 02; dynamic and slowly evolving in Case 03, which converges towards $P_x \sim 0.6/z_l$ for $t > 200T_0$. The transient $P_x$ in Case 03 is responsible for keeping the quasi-steady wind and friction velocities.

Neglecting viscous effects, the WI stress integrates into the form drag in the lower surface, and is responsible for the vertical transfer of momentum between waves and wind [14,15]. If the waves travel sufficiently fast (large WA) with a significant wave slope, the upward transfer is such that the wave signature emerges as a low-level jet, characterized by a local maximum in the wind velocity profiles [16–18].

The momentum transfer between the wave and the wind is evaluated in the non-dimensional form drag (per surface unit) $\overline{F}_p/\overline{F}_t$, which also exhibits a transient behavior in Figure 3d. While Case 03 converges towards $\overline{F}_p/\overline{F}_t \sim -15\%$ for $t > 200T_0$, Case 02 continues to increase after its sudden drop and minimum occurring at $t \sim 10T_0$. Moreover in Case 03, $\overline{F}_p < 0$ and the momentum transfers exclusively upwards, but for Case 02, the form drag changes the sign, and upon becoming positive, it indicates a downward (wind-to-wave) momentum transfer for $t > 100T_0$.
Interestingly for Case 02, the freestream velocity converges towards the wave phase velocity in the last instants of Figure 3a, when $U_h/c = 99.5\%$. This means that Miles' critical height, where $U = c$ [3,74,75], is being brought inside the computational domain.

The mean wind and turbulence profiles in the last instants of the depicted evolution are presented in Figure 4. In the longitudinal velocities of Figure 4a, the speed-up of Case 02 occurs through the whole extension of the domain, while in Case 03 the speed-up is restricted to lower altitudes.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Height-dependent wind and turbulence profiles averaged with $t_{avg} = 75T_0$. Pilot application. (a) longitudinal velocity; (b) total shear stress; (c) total turbulent kinetic energy (TKE).</figcaption></figure>

Inferred from the velocity time histories in Figure 3a,b, the different wave ages between Cases 02 and 03 qualitatively justify the differences in the turbulent profiles of Figure 4b,c, where they are normalized by the resolved friction velocity $\overline{u}^*$. In the turbulent profiles, the disturbances previously reported in [52], for the TKE and the total shear stresses, are again observed in Figure 4, and increased in Case 03, where the WA is higher. A systematic WA comparison is presented in Section 5.

## 4. Flat-Bottom Application

In this section, we show that the proposed methodology preserves the self-similarity of turbulent statistics for flat-bottom cases while employing the dynamic $P_x$ modeler to build varying wind initial conditions prior to Section 5. These flat-bottom references provide guidance for the evaluation of WI disturbances later on.

### Test Cases

Five cases are presented, varying the freestream velocity from its initial value. The dynamic $P_x$ method is thus employed with varying target velocities: Cases [01, 02, 03, 04, 05] with target $u_{ref} = [0.75, 0.90, 1.00, 1.10, 1.25] u_{ref0}$, respectively, and with $u_{ref0} = 52.11 u_0$. The nominal friction velocity characterizing each case is thus $u_{ref}^* = u_{ref}/52.11$, and the corresponding characteristic time-scale is $T_{ref} = \lambda/u_{ref}^*$.

### 4.1. Evolution and Turbulent Profiles

Time histories are shown in Figure 5, by employing sliding-time statistics with short average periods $t_{avg} = 75 \times 10^{-3} T_0$, in order to assess the different turbulence levels and periods that characterize the evolution of different quantities. The mean wind velocity in the last cell layer ($U_h$) is shown in Figure 5a, showing that it tends to the specified values ($u_{ref}$) depicted by the horizontal lines in the figure. The low level of fluctuations revealed for $U_h$ is consistent with the low turbulence level imposed by the upper BC.

The dynamic $P_x$ is shown in Figure 5b. As expected, $P_x$ evolves in larger periods than $U_h$, with its natural period imposed as $T_P = 77T_0$. The friction velocity is normalized by $u_{ref}^*$ in Figure 5c. As the last cell velocity approaches the target, the non-dimensional solution converges between the cases attesting to their self-similarity.
Interestingly for Case 02, the freestream velocity converges towards the wave phase velocity in the last instants of Figure 3a, when $U_h/c = 99.5\%$. This means that Miles' critical height, where $U = c$ [3,74,75], is being brought inside the computational domain.

The mean wind and turbulence profiles in the last instants of the depicted evolution are presented in Figure 4. In the longitudinal velocities of Figure 4a, the speed-up of Case 02 occurs through the whole extension of the domain, while in Case 03 the speed-up is restricted to lower altitudes.

<figure><img src="image_11.png" /><figcaption>**Figure 4.** Height-dependent wind and turbulence profiles averaged with $t_{avg} = 75T_0$. Pilot application. (a) longitudinal velocity; (b) total shear stress; (c) total turbulent kinetic energy (TKE).</figcaption></figure>

Inferred from the velocity time histories in Figure 3a,b, the different wave ages between Cases 02 and 03 qualitatively justify the differences in the turbulent profiles of Figure 4b,c, where they are normalized by the resolved friction velocity $\overline{u}^*$. In the turbulent profiles, the disturbances previously reported in [52], for the TKE and the total shear stresses, are again observed in Figure 4, and increased in Case 03, where the WA is higher. A systematic WA comparison is presented in Section 5.

## 4. Flat-Bottom Application

In this section, we show that the proposed methodology preserves the self-similarity of turbulent statistics for flat-bottom cases while employing the dynamic $P_x$ modeler to build varying wind initial conditions prior to Section 5. These flat-bottom references provide guidance for the evaluation of WI disturbances later on.

### Test Cases

Five cases are presented, varying the freestream velocity from its initial value. The dynamic $P_x$ method is thus employed with varying target velocities: Cases [01, 02, 03, 04, 05] with target $u_{ref} = [0.75, 0.90, 1.00, 1.10, 1.25] u_{ref0}$, respectively, and with $u_{ref0} = 52.11 u_0$. The nominal friction velocity characterizing each case is thus $u_{ref}^* = u_{ref}/52.11$, and the corresponding characteristic time-scale is $T_{ref} = \lambda/u_{ref}^*$.

### 4.1. Evolution and Turbulent Profiles

Time histories are shown in Figure 5, by employing sliding-time statistics with short average periods $t_{avg} = 75 \times 10^{-3} T_0$, in order to assess the different turbulence levels and periods that characterize the evolution of different quantities. The mean wind velocity in the last cell layer ($U_h$) is shown in Figure 5a, showing that it tends to the specified values ($u_{ref}$) depicted by the horizontal lines in the figure. The low level of fluctuations revealed for $U_h$ is consistent with the low turbulence level imposed by the upper BC.

The dynamic $P_x$ is shown in Figure 5b. As expected, $P_x$ evolves in larger periods than $U_h$, with its natural period imposed as $T_p = 77T_0$. The friction velocity is normalized by $u_{ref}^*$ in Figure 5c. As the last cell velocity approaches the target, the non-dimensional solution converges between the cases attesting to their self-similarity.
<figure><img src="image_8.png" /><figcaption>Figure 5. Time histories of non-dimensional: (a) wind speed at $S_{z+}$, with target values denoted in the horizontal lines; (b) large-scale pressure gradient; (c) friction velocity. Averaged with $t_{avg} = 75 \times 10^{-3}T_0$. Flat-bottom application.</figcaption></figure>

The vertical wind and turbulent profiles are shown in Figure 6. They provide further evidence that self-similarity is preserved since the cases almost perfectly collapse in the non-dimensional, mean longitudinal velocity $\langle \bar{u} \rangle / \bar{u}^*$, turbulence intensity (TI) $\mathrm{TI}_{\phi\phi} = \sqrt{\langle \phi' \phi' \rangle / \langle \bar{u} \rangle}$, and total shear stresses $\langle \bar{\tau}_{13} \rangle / (\rho_{\infty} \bar{u}^{*2})$, shown respectively in Figure 6a–c.

<figure><img src="image_9.png" /><figcaption>Figure 6. Height-dependent wind and turbulence profiles averaged with $t_{avg} = 75T_0$. Flat-bottom application. (a) longitudinal velocity; (b) turbulence intensity (TI); (c) total shear stress. The normalizing friction velocity is resolved, as obtained from Equation (8).</figcaption></figure>

## 4.2. Wind Velocity Profiles and Fitting

We applied a logarithmic fitting to the wind profiles, by introducing the effective friction velocity ($\alpha_u \bar{u}^*$) and roughness length ($\alpha_z z_0$), in Equation (15):

$$ \bar{u}_{\mathrm{FIT}} = \frac{\alpha_u \bar{u}^*}{\kappa} \ln \left( \frac{z}{\alpha_z z_0} \right), \qquad (15) $$

thus recovering the log law of Equation (6) when $[\alpha_u = 1, \alpha_z = 1]$. With parameters $[\alpha_u, \alpha_z]$, the fitting minimizes the RMS difference $\mathrm{RMS}_{diff}$, defined in Equation (16) in an arbitrary region $z_{min} < z < z_{max}$.

$$ \mathrm{RMS}_{\mathrm{diff}} = \sqrt{\frac{1}{n} \sum_n \left[ \frac{(\bar{u}_{\mathrm{FIT}} - \langle \bar{u} \rangle)^2}{\langle \bar{u} \rangle^2} \right]}, \quad \forall \{n \mid z_{\mathrm{min}} \le z_n \le z_{\mathrm{max}}\}. \qquad (16) $$

This fitting is embedded in an adaptive procedure determining the log law region at $z_{min} \le z \le z_{max}$. Given the threshold $\mathrm{RMS}_c$, and starting at $z_{max}$, the lower limit ($z_{min}$) descends while $\mathrm{RMS}_{diff} < \mathrm{RMS}_c$. Thus, the log law region is defined as the most extensive portion of the ABL below $z_{max}$, in which the log fitting error is less than $\mathrm{RMS}_c = 5 \times 10^{-3}$. Here $z_{max} = 0.8L_0$ is fixed, in order to later guide the evaluation of WI disturbances in the lower part of the logarithmic region.

The logarithmic fitting is depicted in Figure 7, revealing different regions inside the domain. In the numerical inner surface layer, the modeled SGS stresses are significant. In
<figure><img src="image_12.png" /><figcaption>**Figure 5.** Time histories of non-dimensional: (a) wind speed at $S_{z+}$, with target values denoted in the horizontal lines; (b) large-scale pressure gradient; (c) friction velocity. Averaged with $t_{avg} = 75 \times 10^{-3}T_0$. Flat-bottom application.</figcaption></figure>

The vertical wind and turbulent profiles are shown in Figure 6. They provide further evidence that self-similarity is preserved since the cases almost perfectly collapse in the non-dimensional, mean longitudinal velocity $\langle \bar{u} \rangle / \bar{u}^*$, turbulence intensity (TI) $\mathrm{TI}_{\phi\phi} = \sqrt{\langle \phi' \phi' \rangle / \langle \bar{u} \rangle}$, and total shear stresses $\langle \bar{\tau}_{13} \rangle / (\rho_{\infty} \bar{u}^{*2})$, shown respectively in Figure 6a–c.

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Height-dependent wind and turbulence profiles averaged with $t_{avg} = 75T_0$. Flat-bottom application. (a) longitudinal velocity; (b) turbulence intensity (TI); (c) total shear stress. The normalizing friction velocity is resolved, as obtained from Equation (8).</figcaption></figure>

## 4.2. Wind Velocity Profiles and Fitting

We applied a logarithmic fitting to the wind profiles, by introducing the effective friction velocity ($\alpha_u \bar{u}^*$) and roughness length ($\alpha_z z_0$), in Equation (15):

$$ \bar{u}_{\mathrm{FIT}} = \frac{\alpha_u \bar{u}^*}{\kappa} \ln \left( \frac{z}{\alpha_z z_0} \right), \qquad (15) $$

thus recovering the log law of Equation (6) when $[\alpha_u = 1, \alpha_z = 1]$. With parameters $[\alpha_u, \alpha_z]$, the fitting minimizes the RMS difference $\mathrm{RMS}_{diff}$, defined in Equation (16) in an arbitrary region $z_{min} < z < z_{max}$.

$$ \mathrm{RMS}_{diff} = \sqrt{\frac{1}{n} \sum_n \left[ \frac{(\bar{u}_{\mathrm{FIT}} - \langle \bar{u} \rangle)^2}{\langle \bar{u} \rangle^2} \right]}, \quad \forall \{n \mid z_{min} \le z_n \le z_{max}\}. \qquad (16) $$

This fitting is embedded in an adaptive procedure determining the log law region at $z_{min} \le z \le z_{max}$. Given the threshold $\mathrm{RMS}_c$, and starting at $z_{max}$, the lower limit ($z_{min}$) descends while $\mathrm{RMS}_{diff} < \mathrm{RMS}_c$. Thus, the log law region is defined as the most extensive portion of the ABL below $z_{max}$, in which the log fitting error is less than $\mathrm{RMS}_c = 5 \times 10^{-3}$. Here $z_{max} = 0.8L_0$ is fixed, in order to later guide the evaluation of WI disturbances in the lower part of the logarithmic region.

The logarithmic fitting is depicted in Figure 7, revealing different regions inside the domain. In the numerical inner surface layer, the modeled SGS stresses are significant. In
the outer surface region, highly anisotropic turbulence dominates the shear flow, the SGS stresses become negligible with respect to the LES-resolved turbulence, and the longitudinal velocity is given by a log law profile; this so-called logarithmic region is observed above a numerical buffer region and below the outer layer. In the outer layer, the flow adapts to the freestream condition, turbulence vanishes, and the velocity profiles exhibit an almost linear behavior [52].

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Wind longitudinal velocity profiles, with different colors for Cases 01 to 05 of the flat-bottom application. The log law fitting to the outer surface layer (by Equation (15)) is shown with transparent lines.</figcaption></figure>

Considering the uncertainty as (±0.5) the maximum variation between the cases, the described fitting leads to $\alpha_u = 0.93 \pm 0.03$ and $\alpha_z = 0.43 \pm 0.16$, with $z_{min}/L_0 = 0.03$.

**5. Wave Age Variation**

The proposed methodology is employed here to perform a systematic WA variation study in old seas ($48 \le WA_{ref} \le 79$), by varying the wind speed in the range ±25%.

**Test Cases**

The wavy bottom surface replaces the flat-bottom cases previously discussed in Section 4. A single monochromatic wave is considered, with ka = 0.2 and λ = L₀. The varying wind conditions lead to nominal WA<sub>ref</sub> = c/u<sub>ref</sub> = [79, 66, 60, 54, 48], for Cases [01, 02, 03, 04, 05], respectively.

**5.1. Time Histories**

The normalized mean wind velocity in the upper cell layer ($U_h$) is shown in Figure 8a. The horizontal lines in the figure depict the target velocities ($u_{ref}$). The highest discrepancy observed between $U_h$ and $u_{ref}$ is 3.1% by the middle of the evolution depicted for Case 01 (highest WA). In the last steps of the evolution, the discrepancy in Case 01 is still the highest and drops to 1.9%.

Contrary to the flat-bottom cases of Section 4, here $\overline{u^*}/u^*_{ref}$ becomes a function of the WA as observed in Figure 8b. Except for Case 01, the friction velocities drop with respect to the reference values, and the discrepancy is larger with larger reference velocities, becoming the largest in Cases 03, 04, and 05, with 13% deviation in the last steps of Figure 8b. Cases 03, 04, and 05 converge to the same values of $\overline{u^*}/u^*_{ref}$, which suggests a saturation for this ratio at the lowest WA considered, but further investigations are needed with broader WA ranges to confirm this trend.

The dynamic pressure gradient decreases due to the introduction of the wave, evolving as reported in Figure 8c. The form drag evolution is depicted in Figure 8d, revealing an exponential trend between Cases 01 and 05 that is further evaluated in Section 5.3. The quantities presented in Table 1 were evaluated in the final time-step of the evolution depicted in Figure 8.
the outer surface region, highly anisotropic turbulence dominates the shear flow, the SGS stresses become negligible with respect to the LES-resolved turbulence, and the longitudinal velocity is given by a log law profile; this so-called logarithmic region is observed above a numerical buffer region and below the outer layer. In the outer layer, the flow adapts to the freestream condition, turbulence vanishes, and the velocity profiles exhibit an almost linear behavior [52].

<figure><img src="image_14.png" /><figcaption>**Figure 7.** Wind longitudinal velocity profiles, with different colors for Cases 01 to 05 of the flat-bottom application. The log law fitting to the outer surface layer (by Equation (15)) is shown with transparent lines.</figcaption></figure>

Considering the uncertainty as (±0.5) the maximum variation between the cases, the described fitting leads to $\alpha_u = 0.93 \pm 0.03$ and $\alpha_z = 0.43 \pm 0.16$, with $z_{min}/L_0 = 0.03$.

**5. Wave Age Variation**

The proposed methodology is employed here to perform a systematic WA variation study in old seas ($48 \le WA_{ref} \le 79$), by varying the wind speed in the range ±25%.

**Test Cases**

The wavy bottom surface replaces the flat-bottom cases previously discussed in Section 4. A single monochromatic wave is considered, with $ka = 0.2$ and $\lambda = L_0$. The varying wind conditions lead to nominal $WA_{ref} = c/u_{ref}^* = [79, 66, 60, 54, 48]$, for Cases [01, 02, 03, 04, 05], respectively.

### 5.1. *Time Histories*

The normalized mean wind velocity in the upper cell layer ($U_h$) is shown in Figure 8a. The horizontal lines in the figure depict the target velocities ($u_{ref}$). The highest discrepancy observed between $U_h$ and $u_{ref}$ is 3.1% by the middle of the evolution depicted for Case 01 (highest WA). In the last steps of the evolution, the discrepancy in Case 01 is still the highest and drops to 1.9%.

Contrary to the flat-bottom cases of Section 4, here $\overline{u^*}/u_{ref}^*$ becomes a function of the WA as observed in Figure 8b. Except for Case 01, the friction velocities drop with respect to the reference values, and the discrepancy is larger with larger reference velocities, becoming the largest in Cases 03, 04, and 05, with 13% deviation in the last steps of Figure 8b. Cases 03, 04, and 05 converge to the same values of $\overline{u^*}/u_{ref}^*$, which suggests a saturation for this ratio at the lowest WA considered, but further investigations are needed with broader WA ranges to confirm this trend.

The dynamic pressure gradient decreases due to the introduction of the wave, evolving as reported in Figure 8c. The form drag evolution is depicted in Figure 8d, revealing an exponential trend between Cases 01 and 05 that is further evaluated in Section 5.3. The quantities presented in Table 1 were evaluated in the final time-step of the evolution depicted in Figure 8.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Time histories of non-dimensional: (a) wind speed at $S_{z+}$, with target values denoted in the horizontal lines; (b) friction velocity; (c) large-scale pressure gradient; and (d) form drag. Averaged with $t_{avg} = 75T_0$. Wavy bottom application.</figcaption></figure>

<table><caption>**Table 1.** The parameters of Figure 8 probed in the final time-step, and the corresponding reference values.</caption><thead><tr><th>Case ID</th><th><i>u</i><sub>ref</sub>/<i>u</i><sub>ref0</sub></th><th><i>U</i><sub>h</sub> − <i>u</i><sub>ref</sub>/<i>u</i><sub>ref</sub></th><th><i>P</i><sub>x</sub> · <i>z</i><sub>l</sub>/<i>ρ</i><sub>∞</sub><i>u</i><sub>0</sub><sup>*2</sup></th><th><i>u</i><sup>*</sup>/<i>u</i><sub>ref</sub><sup>*</sup></th><th>WA<sub>ref</sub></th><th>WA</th></tr></thead><tbody><tr><td>1</td><td>0.75</td><td>1.85%</td><td>0.107</td><td>1.073</td><td>79</td><td>74</td></tr><tr><td>2</td><td>0.90</td><td>0.18%</td><td>0.361</td><td>0.926</td><td>66</td><td>72</td></tr><tr><td>3</td><td>1.00</td><td>1.24%</td><td>0.479</td><td>0.871</td><td>60</td><td>68</td></tr><tr><td>4</td><td>1.10</td><td>−0.04%</td><td>0.876</td><td>0.870</td><td>54</td><td>62</td></tr><tr><td>5</td><td>1.25</td><td>−0.49%</td><td>1.270</td><td>0.869</td><td>48</td><td>55</td></tr></tbody></table>

## 5.2. Wind Velocity Profiles

The normalized wind profiles vary substantially between the cases, as shown in Figure 9a. The log law functions were fitted to the outer surface layer, according to the procedure described in Section 4.2. The fitting is here depicted in Figure 9b, with the solution reported in Table 2. The flat cases presented in Section 4.2 appear in the tables with the case-averaged values, and the uncertainties are taken as (±0.5) the maximum variation between the flat-bottom cases.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Longitudinal velocity profiles averaged with $t_{avg} = 75T_0$. (a) Normalized by the resolved friction velocity. (b) Normalized by the reference friction velocity; with logarithmic profiles fitted to the outer region, depicted in black transparent lines.</figcaption></figure>

In Figure 9b, it can be seen that in the (numerical) inner surface layer, the effective roughness length increases with increasing WA, and the buffer region extends higher in the domain. As a consequence in Table 2, one observes a clear trend of increasing $z_{min}$ with increasing WA. The logarithmic layer is thus confined to higher altitudes, suggesting an increased WBL height scaling with the WA.
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Time histories of non-dimensional: (a) wind speed at $S_{z+}$, with target values denoted in the horizontal lines; (b) friction velocity; (c) large-scale pressure gradient; and (d) form drag. Averaged with $t_{avg} = 75T_0$. Wavy bottom application.</figcaption></figure>

<table><caption>**Table 1.** The parameters of Figure 8 probed in the final time-step, and the corresponding reference values.</caption><thead><tr><th>Case ID</th><th><i>u</i><sub>ref</sub>/<i>u</i><sub>ref0</sub></th><th><i>U</i><sub>h</sub> − <i>u</i><sub>ref</sub>/<i>u</i><sub>ref</sub></th><th><i>P</i><sub>x</sub> · <i>z</i><sub>l</sub>/<i>ρ</i><sub>∞</sub><i>u</i><sub>0</sub><sup>*2</sup></th><th><i>u</i><sup>*</sup>/<i>u</i><sub>ref</sub><sup>*</sup></th><th>WA<sub>ref</sub></th><th>WA</th></tr></thead><tbody><tr><td>1</td><td>0.75</td><td>1.85%</td><td>0.107</td><td>1.073</td><td>79</td><td>74</td></tr><tr><td>2</td><td>0.90</td><td>0.18%</td><td>0.361</td><td>0.926</td><td>66</td><td>72</td></tr><tr><td>3</td><td>1.00</td><td>1.24%</td><td>0.479</td><td>0.871</td><td>60</td><td>68</td></tr><tr><td>4</td><td>1.10</td><td>−0.04%</td><td>0.876</td><td>0.870</td><td>54</td><td>62</td></tr><tr><td>5</td><td>1.25</td><td>−0.49%</td><td>1.270</td><td>0.869</td><td>48</td><td>55</td></tr></tbody></table>

## 5.2. Wind Velocity Profiles

The normalized wind profiles vary substantially between the cases, as shown in Figure 9a. The log law functions were fitted to the outer surface layer, according to the procedure described in Section 4.2. The fitting is here depicted in Figure 9b, with the solution reported in Table 2. The flat cases presented in Section 4.2 appear in the tables with the case-averaged values, and the uncertainties are taken as (±0.5) the maximum variation between the flat-bottom cases.

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Longitudinal velocity profiles averaged with $t_{avg} = 75T_0$. (a) Normalized by the resolved friction velocity. (b) Normalized by the reference friction velocity; with logarithmic profiles fitted to the outer region, depicted in black transparent lines.</figcaption></figure>

In Figure 9b, it can be seen that in the (numerical) inner surface layer, the effective roughness length increases with increasing WA, and the buffer region extends higher in the domain. As a consequence in Table 2, one observes a clear trend of increasing $z_{min}$ with increasing WA. The logarithmic layer is thus confined to higher altitudes, suggesting an increased WBL height scaling with the WA.
<table><caption>**Table 2.** Log law-fit of Equation (15) to the outer surface layer: Optimal [$\alpha_u, \alpha_z$], for RMS\_diff < 3 × 10⁻³, in [$z_{min}, z_{max}$], with $z_{max} = 0.8\lambda$. [$\alpha_u, \alpha_z$] = [1, 1] recovers the BC of Equation (6).</caption><thead><tr><th>iCase</th><th>WA<sub>ref</sub></th><th>α<sub>u</sub></th><th>α<sub>z</sub></th><th>z<sub>min</sub>/λ</th></tr></thead><tbody><tr><td>1</td><td>79</td><td>0.94</td><td>0.03</td><td>0.18</td></tr><tr><td>2</td><td>66</td><td>1.40</td><td>4.16</td><td>0.14</td></tr><tr><td>3</td><td>60</td><td>1.48</td><td>9.41</td><td>0.11</td></tr><tr><td>4</td><td>54</td><td>1.57</td><td>20.09</td><td>0.08</td></tr><tr><td>5</td><td>48</td><td>1.58</td><td>29.01</td><td>0.06</td></tr><tr><td>Flat</td><td>-</td><td>0.93 ± 0.03</td><td>0.43 ± 0.16</td><td>0.03</td></tr></tbody></table>

Contrary to the inner region, but as expected from the literature [34], in the outer surface layer the effective roughness length decreases with increasing WA, as seen in Table 2. Except for Case 01 (highest WA), there are strong friction velocity disturbances ($\alpha_u$) due to the wave incidence; between Cases 2 and 5 these vary from 40% to 58%, compared to only −6% in Case 01 and the uncertainty of ±3%.

### 5.3. Form Drag Parametrization

In the literature [12,38,53,76], the form drag $F_p$ (evolution in Figure 8d), or the growth rate parameter $\beta = 2 F_p / [u^* (ka)^2]$, are usually assumed as dependent on the WA. Nevertheless, different definitions are then employed to define such a WA parameter [76]. We discussed different definitions of friction velocity. So in order to evaluate these relations, Figure 10a presents $\bar{F}_p / \bar{F}_t$ as either dependent on the resolved WA = $c/u^*$ (in red circles), or the nominal WA$_{ref}$ = $c/u_{ref}^*$ (in blue dots). Defining the nominal WA$_{ref}$, $u_{ref}^*$ scales linearly to the target velocity $u_{ref}$, and $u_{ref}$ approximates the resolved freestream velocity ($U_h$) with better than 3% accuracy.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Wave age dependency (Resolved WA, or Nominal WA$_{ref}$) and fitting of the normalized form drag. The LES results are shown in dots and circles, and the fittings are given by full and dot-dashed lines. (a) Form drag and the exponential fitting in linear scale. (b) The negative signs of the slopes in the log scale are evaluated between the intervals of (a). A linear fit applies to the log of the slopes in (b), and integrates into the exponential in (a). The fitting quality was evaluated by the $R^2$ criterion.</figcaption></figure>

Exponential functions are fitted to the data. The fitting of constants A and B from Equation (17) applies to the slopes of the LES results shown in Figure 10a, evaluated in the middle of the intervals by a first-order finite difference, as denoted by the dots and circles in Figure 10b. The constant part of the approximation function (C in Equation (17)) is set to match exactly the solution in Case 03.

The fitting is denoted by the curves, and the quality of the procedure is evaluated by the $R^2$ criterion in Figure 10a,b. It is thus seen that in the present formulation, the form drag scales better with the nominal WA$_{ref}$ ($R^2 = 0.99$), which characterizes the
<table><caption>**Table 2.** Log law-fit of Equation (15) to the outer surface layer: Optimal [$\alpha_u, \alpha_z$], for RMS<sub>diff</sub> < 3 × 10<sup>-3</sup>, in [$z_{min}, z_{max}$], with $z_{max} = 0.8\lambda$. [$\alpha_u, \alpha_z$] = [1, 1] recovers the BC of Equation (6).</caption><thead><tr><th>iCase</th><th>WA<sub>ref</sub></th><th>α<sub>u</sub></th><th>α<sub>z</sub></th><th>z<sub>min</sub>/λ</th></tr></thead><tbody><tr><td>1</td><td>79</td><td>0.94</td><td>0.03</td><td>0.18</td></tr><tr><td>2</td><td>66</td><td>1.40</td><td>4.16</td><td>0.14</td></tr><tr><td>3</td><td>60</td><td>1.48</td><td>9.41</td><td>0.11</td></tr><tr><td>4</td><td>54</td><td>1.57</td><td>20.09</td><td>0.08</td></tr><tr><td>5</td><td>48</td><td>1.58</td><td>29.01</td><td>0.06</td></tr><tr><td>Flat</td><td>-</td><td>0.93 ± 0.03</td><td>0.43 ± 0.16</td><td>0.03</td></tr></tbody></table>

Contrary to the inner region, but as expected from the literature [34], in the outer surface layer the effective roughness length decreases with increasing WA, as seen in Table 2. Except for Case 01 (highest WA), there are strong friction velocity disturbances ($\alpha_u$) due to the wave incidence; between Cases 2 and 5 these vary from 40% to 58%, compared to only −6% in Case 01 and the uncertainty of ±3%.

### 5.3. Form Drag Parametrization

In the literature [12,38,53,76], the form drag $F_p$ (evolution in Figure 8d), or the growth rate parameter $\beta = 2 F_p / [u^* (ka)^2]$, are usually assumed as dependent on the WA. Nevertheless, different definitions are then employed to define such a WA parameter [76]. We discussed different definitions of friction velocity. So in order to evaluate these relations, Figure 10a presents $\bar{F}_p / \bar{F}_t$ as either dependent on the resolved WA = $c/u^*$ (in red circles), or the nominal WA<sub>ref</sub> = $c/u_{ref}^*$ (in blue dots). Defining the nominal WA<sub>ref</sub>, $u_{ref}^*$ scales linearly to the target velocity $u_{ref}$, and $u_{ref}$ approximates the resolved freestream velocity ($U_h$) with better than 3% accuracy.

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Wave age dependency (Resolved WA, or Nominal WA<sub>ref</sub>) and fitting of the normalized form drag. The LES results are shown in dots and circles, and the fittings are given by full and dot-dashed lines. (a) Form drag and the exponential fitting in linear scale. (b) The negative signs of the slopes in the log scale are evaluated between the intervals of (a). A linear fit applies to the log of the slopes in (b), and integrates into the exponential in (a). The fitting quality was evaluated by the R<sup>2</sup> criterion.</figcaption></figure>

Exponential functions are fitted to the data. The fitting of constants A and B from Equation (17) applies to the slopes of the LES results shown in Figure 10a, evaluated in the middle of the intervals by a first-order finite difference, as denoted by the dots and circles in Figure 10b. The constant part of the approximation function (C in Equation (17)) is set to match exactly the solution in Case 03.

The fitting is denoted by the curves, and the quality of the procedure is evaluated by the R<sup>2</sup> criterion in Figure 10a,b. It is thus seen that in the present formulation, the form drag scales better with the nominal WA<sub>ref</sub> (R<sup>2</sup> = 0.99), which characterizes the
non-disturbed flow, and thanks to the $P_x$ modeler used here it scales with the freestream velocity. Therefore, the form drag dependency on the $WA_{ref}$ is given here by the fitting denoted as blue full lines in Figure 10a, obtained with parameters [$A = 1.07 \times 10^{-1}$, $B = 6.75 \times 10^{-6}$, $C = 2.04 \times 10^{-1}$].

$$ \frac{F_p}{F_t} = - \frac{B \exp[A \cdot WA]}{A} + C \qquad (17) $$

## 5.4. Turbulent Profiles

### Turbulent Normal Stresses

The TI for longitudinal and vertical velocities are shown, respectively, in Figure 11a,b. In these figures, the colored lines depict the wavy Cases 01 to 05 as described in the legend, and the light gray lines depict the reference, case averaged flat-bottom solution. As usual, the variances are strongly enhanced close to the surface due to the introduction of the wave. Farther away from the surface, on the contrary, TI diminishes. As the WA increases, the disturbances are amplified, the variances increase close to the surface, and diminish in the outer region.

A striking feature revealed in Figure 11 is that for each $TI_{uu}$ and $TI_{ww}$, the point where the cases intersect with each other and the flat case solution remains almost constant: $z/\lambda = 0.1$ for $TI_{uu}$ in Figure 11a, and $z/\lambda = 0.26$ for $TI_{ww}$ in Figure 11b. Clearly, the enhancement in TI close to the surface is due to the existence of WI motions, shown in Paskin et al. [52] to propagate higher in the WBL for the vertical velocities. TKE budget analyses reported in the literature (c.f. [19,77]) indicate that the decrease in TI on the outer region is due to an imbalance between the turbulent and pressure transport terms in the lower region, bringing turbulent motions from the outer layer into the WBL.

<figure><img src="image_14.png" /><figcaption>Figure 11. Turbulence intensity (TI) profiles averaged with $t_{avg} = 75T_0$. TI for longitudinal and vertical velocity fluctuations, in (a) and (b), respectively. The light gray lines represent the case averaged flat-bottom solution. The zoomed regions depict the intersections between the curves.</figcaption></figure>

### Turbulent Shear stresses

The total shear stresses ($\overline{\langle \tau_{13} \rangle}$) are shown in Figure 12, where they exhibit a mostly linear behavior with varying slopes between the cases. A linear fitting was applied to each profile and is shown with transparent lines in the figure. The theoretical flat-bottom solution is a straight line, denoted in black, varying from $\overline{\langle \tau_{13} \rangle} = \rho_{\infty}u^{*2}$ at the bottom ($S_{z-}$), to $\overline{\langle \tau_{13} \rangle} = 0$ at the top ($S_{z+}$). Due to the waves, and a non-negligible WI form drag, the total shear stress profiles depart from other values rather than $\rho_{\infty}u^{*2}$ in the bottom, and evolve until vanishing in the upper boundary. The fitted profiles evolve from the surface with the values given in Table 3, and are null when $z = z_l$. The stresses deviate from the linear behavior in the (numerical) buffer region, zoomed in the lower right corner of Figure 12.
non-disturbed flow, and thanks to the $P_x$ modeler used here it scales with the freestream velocity. Therefore, the form drag dependency on the WA$_{ref}$ is given here by the fitting denoted as blue full lines in Figure 10a, obtained with parameters [$A = 1.07 \times 10^{-1}$, $B = 6.75 \times 10^{-6}$, $C = 2.04 \times 10^{-1}$].

$$ \frac{F_p}{F_t} = - \frac{B \exp[A \cdot WA]}{A} + C \qquad (17) $$

## 5.4. Turbulent Profiles

### Turbulent Normal Stresses

The TI for longitudinal and vertical velocities are shown, respectively, in Figure 11a,b. In these figures, the colored lines depict the wavy Cases 01 to 05 as described in the legend, and the light gray lines depict the reference, case averaged flat-bottom solution. As usual, the variances are strongly enhanced close to the surface due to the introduction of the wave. Farther away from the surface, on the contrary, TI diminishes. As the WA increases, the disturbances are amplified, the variances increase close to the surface, and diminish in the outer region.

A striking feature revealed in Figure 11 is that for each TI$_{uu}$ and TI$_{ww}$, the point where the cases intersect with each other and the flat case solution remains almost constant: $z/\lambda = 0.1$ for TI$_{uu}$ in Figure 11a, and $z/\lambda = 0.26$ for TI$_{ww}$ in Figure 11b. Clearly, the enhancement in TI close to the surface is due to the existence of WI motions, shown in Paskin et al. [52] to propagate higher in the WBL for the vertical velocities. TKE budget analyses reported in the literature (c.f. [19,77]) indicate that the decrease in TI on the outer region is due to an imbalance between the turbulent and pressure transport terms in the lower region, bringing turbulent motions from the outer layer into the WBL.

<figure><img src="image_18.png" /><figcaption>Figure 11. Turbulence intensity (TI) profiles averaged with $t_{avg} = 75T_0$. TI for longitudinal and vertical velocity fluctuations, in (a) and (b), respectively. The light gray lines represent the case averaged flat-bottom solution. The zoomed regions depict the intersections between the curves.</figcaption></figure>

### Turbulent Shear stresses

The total shear stresses ($\overline{\langle \tau_{13} \rangle}$) are shown in Figure 12, where they exhibit a mostly linear behavior with varying slopes between the cases. A linear fitting was applied to each profile and is shown with transparent lines in the figure. The theoretical flat-bottom solution is a straight line, denoted in black, varying from $\langle \tau_{13} \rangle = \rho_{\infty}u^{*2}$ at the bottom ($S_{z-}$), to $\langle \tau_{13} \rangle = 0$ at the top ($S_{z+}$). Due to the waves, and a non-negligible WI form drag, the total shear stress profiles depart from other values rather than $\rho_{\infty}u^{*2}$ in the bottom, and evolve until vanishing in the upper boundary. The fitted profiles evolve from the surface with the values given in Table 3, and are null when $z = z_l$. The stresses deviate from the linear behavior in the (numerical) buffer region, zoomed in the lower right corner of Figure 12.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Total shear stress, averaged with $t_{avg} = 75T_0$. Linear fittings are shown as transparent lines, parameterized according to Table 3. The black line represents the theoretical flat-bottom solution. The fitting is less accurate in the near-wall region, zoomed in the lower right corner.</figcaption></figure>

The form drag ($\overline{F}_p$) was previously discussed in Figure 8d. In Cases 01 and 02, the drag is negative, and the waves are expected to grow. However, in Case 02, the value is very close to zero, and not surprisingly the shear stresses closely approach the flat-bottom solution in Figure 12. The pressure drag is clearly the determining factor that imposes the slope of the profiles in Figure 12. To clarify this relationship, Figure 13 compares the dependence of $\overline{F}_t$ on the WA$_{ref}$, as obtained from either: $\overline{F}_t = \overline{F}_p + \overline{F}_f$, with $\overline{F}_f = \rho_{\infty}u^{*2}$, resolved $u^*$ reported in Table 1, and $\overline{F}_p/\overline{F}_t$ depicted in Figure 10a; or from the fitted value of $\overline{F}_t = \langle\tau_{13}\rangle(z = 0)$, reported in Table 3. There is a very good agreement between the predictions, for which linear fittings appear as dotted lines in the figures. Nevertheless, there is a small but persistent bias between these estimations, probably induced by the linear fitting employed for the $\langle\tau_{13}\rangle$ profile (in Figure 12, leading to Table 3).

<table><caption>**Table 3.** Linear fitting for $\langle\tau_{13}\rangle(z)$, from Figure 12. $\langle\tau_{13}\rangle(z_l) = 0$. The table gives the values of $\langle\tau_{13}\rangle(0)$.</caption><thead><tr><th>Case ID</th><th>WA<sub>ref</sub></th><th>⟨τ<sub>13</sub>⟩(0)/u<sup>*2</sup></th></tr></thead><tbody><tr><td>1</td><td>79</td><td>-0.21</td></tr><tr><td>2</td><td>66</td><td>-0.62</td></tr><tr><td>3</td><td>60</td><td>-0.80</td></tr><tr><td>4</td><td>54</td><td>-0.98</td></tr><tr><td>5</td><td>48</td><td>-1.10</td></tr></tbody></table>

The fitting coefficients for $\overline{F}_t/u^{*2}$ and $\langle\tau_{13}\rangle(0)/u^{*2}$ in Figure 13 are given in Table 4. The offset between the curves in Figure 13 is observed in the 2% discrepancy in coefficient B, which represents either: an offset of 3% (1.84 of 60) in the WA$_{ref}$; or an offset of 5% (0.05 of 1) in $\overline{F}_t/u^{*2}$.

<table><caption>**Table 4.** The total drag ($\overline{F}_t$) as a function of the WA$_{ref}$ is approximated as $\overline{F}_t/u^{*2} = A \cdot WA_{ref} + B$, with coefficients given below.</caption><thead><tr><th>Case ID</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Resolved <span style="text-decoration: overline;">F</span><sub>t</sub></td><td>2.84×10<sup>-2</sup></td><td>-2.54</td></tr><tr><td>Fitted from ⟨τ<sub>13</sub>⟩(z)</td><td>2.84×10<sup>-2</sup></td><td>-2.48</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>**Figure 12.** Total shear stress, averaged with $t_{avg} = 75T_0$. Linear fittings are shown as transparent lines, parameterized according to Table 3. The black line represents the theoretical flat-bottom solution. The fitting is less accurate in the near-wall region, zoomed in the lower right corner.</figcaption></figure>

The form drag ($\overline{F}_p$) was previously discussed in Figure 8d. In Cases 01 and 02, the drag is negative, and the waves are expected to grow. However, in Case 02, the value is very close to zero, and not surprisingly the shear stresses closely approach the flat-bottom solution in Figure 12. The pressure drag is clearly the determining factor that imposes the slope of the profiles in Figure 12. To clarify this relationship, Figure 13 compares the dependence of $\overline{F}_t$ on the WA$_{ref}$, as obtained from either: $\overline{F}_t = \overline{F}_p + \overline{F}_f$, with $\overline{F}_f = \rho_{\infty}u^{*2}$, resolved $u^*$ reported in Table 1, and $\overline{F}_p/\overline{F}_t$ depicted in Figure 10a; or from the fitted value of $\overline{F}_t = \langle\tau_{13}\rangle(z = 0)$, reported in Table 3. There is a very good agreement between the predictions, for which linear fittings appear as dotted lines in the figures. Nevertheless, there is a small but persistent bias between these estimations, probably induced by the linear fitting employed for the $\langle\tau_{13}\rangle$ profile (in Figure 12, leading to Table 3).

<table><caption>**Table 3.** Linear fitting for $\langle\tau_{13}\rangle(z)$, from Figure 12. $\langle\tau_{13}\rangle(z_l) = 0$. The table gives the values of $\langle\tau_{13}\rangle(0)$.</caption><thead><tr><th>Case ID</th><th>WA<sub>ref</sub></th><th>⟨τ<sub>13</sub>⟩(0)/u<sup>*2</sup></th></tr></thead><tbody><tr><td>1</td><td>79</td><td>-0.21</td></tr><tr><td>2</td><td>66</td><td>-0.62</td></tr><tr><td>3</td><td>60</td><td>-0.80</td></tr><tr><td>4</td><td>54</td><td>-0.98</td></tr><tr><td>5</td><td>48</td><td>-1.10</td></tr></tbody></table>

The fitting coefficients for $\overline{F}_t/u^{*2}$ and $\langle\tau_{13}\rangle(0)/u^{*2}$ in Figure 13 are given in Table 4. The offset between the curves in Figure 13 is observed in the 2% discrepancy in coefficient B, which represents either: an offset of 3% (1.84 of 60) in the WA$_{ref}$; or an offset of 5% (0.05 of 1) in $\overline{F}_t/u^{*2}$.

<table><caption>**Table 4.** The total drag ($\overline{F}_t$) as a function of the WA$_{ref}$ is approximated as $\overline{F}_t/u^{*2} = A \cdot WA_{ref} + B$, with coefficients given below.</caption><thead><tr><th>Case ID</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Resolved <span style="text-decoration: overline;">F</span><sub>t</sub></td><td>2.84×10<sup>-2</sup></td><td>-2.54</td></tr><tr><td>Fitted from ⟨τ<sub>13</sub>⟩(z)</td><td>2.84×10<sup>-2</sup></td><td>-2.48</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Normalized total drag as a function of the reference $WA_{ref}$. Comparison between the resolved value, and the one extrapolated from the linear fitting depicted in Figure 12. Here, a linear fitting applies to $\overline{F_t}(WA_{ref})$, and is shown with dotted lines in the figure. The fitting coefficients are given in Table 4.</figcaption></figure>

# 6. Discussion

## 6.1. Reviewing the Mathematical Formulation of a Physical MABL in LES

We considered the mathematical formulation of a marine atmospheric boundary layer (MABL) in large eddy simulation(s) (LES), in a simplified framework summarized in Section 6.2. The MABL is bounded below by ocean waves traveling in the free surface, and above by the freestream. These physics translate in LES as a moving bottom with no-slip-no-penetration, and an upper surface with free-slip-no-penetration boundary conditions (BCs). Ocean waves carry momentum upwards or downwards in the MABL, but wave-induced (WI) disturbances are supposedly limited to the free surface vicinity [10–13]. Nevertheless, case-dependent WI stresses combined with a free-slip upper BC in LES allow a varying Wind Speed (WS) profile across the vertical extension of the domain. The WI disturbance of the freestream is unlikely to occur in the MABL but is intrinsic to the mathematical formulations previously employed in LES (c.f. [49,51,53,57]).

Atmospheric Reynolds averaged numerical simulations (RANS) employ a no-slip upper BC to drive the flow in the domain, which precisely imposes the mean freestream velocity in the upper boundary (c.f. [78]). Imposing no-slip upper BCs in DNS (c.f. [38,42]) or LES (c.f. [58]) formulations differ from an atmospheric setup, as a secondary boundary layer develops attached to the top boundary of the domain. Instead, in atmospheric LES the flow is driven by large-scale forces. These large-scale forces exist in the MABL coupling the meso- to the micro-scales, both as a large-scale pressure gradient and the Coriolis force (c.f. [66]). In the vast majority of previous LES studies of micro-scale wind-wave interactions the large-scale forces were considered as a constant [47–51,53–57,60]. However, there is no physical reason for this assumption (c.f. [59,62,65]), since the large-scale conditions could well vary in scales much larger than the computational domain.

## 6.2. A Simplified Framework

The simplified numerical framework exploited in this paper and detailed in Section 2.1 can be summarized as follows: Wind and waves propagate in the same direction, in a prismatic domain representing a wind tunnel; a single monochromatic wave is prescribed and not resolved, in order to prevent the free surface evolution under the wind forcing; old seas are considered, as the waves propagate rapidly compared to the wind; buoyancy forces are neglected, thus representing neutral stratification conditions; turbulence is fully developed in the air, so that the molecular viscosity can be neglected; Coriolis forces are neglected, as we focus on the airflow close to the sea surface; the test cases are exploited once the solution has converged to a statistically stationary regime.
<figure><img src="image_20.png" /><figcaption>**Figure 13.** Normalized total drag as a function of the reference $WA_{ref}$. Comparison between the resolved value, and the one extrapolated from the linear fitting depicted in Figure 12. Here, a linear fitting applies to $\overline{F_t}(WA_{ref})$, and is shown with dotted lines in the figure. The fitting coefficients are given in Table 4.</figcaption></figure>

# 6. Discussion

## 6.1. Reviewing the Mathematical Formulation of a Physical MABL in LES

We considered the mathematical formulation of a marine atmospheric boundary layer (MABL) in large eddy simulation(s) (LES), in a simplified framework summarized in Section 6.2. The MABL is bounded below by ocean waves traveling in the free surface, and above by the freestream. These physics translate in LES as a moving bottom with no-slip-no-penetration, and an upper surface with free-slip-no-penetration boundary conditions (BCs). Ocean waves carry momentum upwards or downwards in the MABL, but wave-induced (WI) disturbances are supposedly limited to the free surface vicinity [10–13]. Nevertheless, case-dependent WI stresses combined with a free-slip upper BC in LES allow a varying Wind Speed (WS) profile across the vertical extension of the domain. The WI disturbance of the freestream is unlikely to occur in the MABL but is intrinsic to the mathematical formulations previously employed in LES (c.f. [49,51,53,57]).

Atmospheric Reynolds averaged numerical simulations (RANS) employ a no-slip upper BC to drive the flow in the domain, which precisely imposes the mean freestream velocity in the upper boundary (c.f. [78]). Imposing no-slip upper BCs in DNS (c.f. [38,42]) or LES (c.f. [58]) formulations differ from an atmospheric setup, as a secondary boundary layer develops attached to the top boundary of the domain. Instead, in atmospheric LES the flow is driven by large-scale forces. These large-scale forces exist in the MABL coupling the meso- to the micro-scales, both as a large-scale pressure gradient and the Coriolis force (c.f. [66]). In the vast majority of previous LES studies of micro-scale wind-wave interactions the large-scale forces were considered as a constant [47–51,53–57,60]. However, there is no physical reason for this assumption (c.f. [59,62,65]), since the large-scale conditions could well vary in scales much larger than the computational domain.

## 6.2. A Simplified Framework

The simplified numerical framework exploited in this paper and detailed in Section 2.1 can be summarized as follows: Wind and waves propagate in the same direction, in a prismatic domain representing a wind tunnel; a single monochromatic wave is prescribed and not resolved, in order to prevent the free surface evolution under the wind forcing; old seas are considered, as the waves propagate rapidly compared to the wind; buoyancy forces are neglected, thus representing neutral stratification conditions; turbulence is fully developed in the air, so that the molecular viscosity can be neglected; Coriolis forces are neglected, as we focus on the airflow close to the sea surface; the test cases are exploited once the solution has converged to a statistically stationary regime.
### 6.3. Evaluating Standard LES Formulations of the MABL

Considering the simplified framework adopted, the results for a standard LES formulation are presented in Section 3 as Cases 01 (flat-bottom reference) and 02 (wavy bottom reference). Case 02 reveals a certain number of artificial features unlikely to occur in the physical MABL. Notably, the mean velocities increased to the upper boundary, and the freestream velocity is greatly augmented due to the incidence of the wave. The freestream velocity increases up to a moment where it drags back the propagating wave. The downward momentum transfer that occurs is characteristic of young seas, even though we attempted to reproduce old seas conditions in this simulation. The normalized total shear stress converges towards the flat-bottom case. These processes reveal a long transient regime, while the flow adapts all the way up to the top boundary. Most interestingly and hardly a coincidence, the freestream self-adapts towards a condition where its velocity matches the wave phase velocity ($U_h/c = 99.5\%$), thus bringing Miles' critical height (c.f. [3,74,75]) inside the numerical domain.

Considering the dynamic pressure gradient proposed herein and employed in Case 03 of Section 3, these artificial features disappear. The friction velocity self-adapts to the appearance of the wave, but the freestream velocity is kept constant. A steady state regime is rapidly reached. The momentum flux occurs exclusively upwards. The WI disturbance in the wind velocity profile is kept confined to lower altitudes.

### 6.4. Large-Scale Dynamically Evolving Pressure Gradient Modeler

In the physical MABL and as long as the boundary layer height is sufficiently distant from the ocean (in neutral or unstable stratification conditions), it is more likely that the large-scale pressure gradient (traveling instantaneously in incompressible flows) will adapt to maintain a freestream velocity, rather than that, the freestream velocity will vary to maintain a constant large-scale pressure gradient condition.

In order to maintain an undisturbed freestream, while resolving arbitrary WI stresses, we developed a general method to control the wind speed at a certain height in the numerical domain (Section 2.2). To achieve this, an artificial evolution equation was built for the spatially-homogeneous large-scale pressure gradient: $P_x = \partial P/\partial x$. As an innovation in the wind-wave interaction studies, $P_x$ is thus allowed to evolve in time, restricted to scales much larger than the resolved ones. Unlike other solutions developed for inland applications [61,62,64,66], this methodology does not require previous knowledge about the instantaneous characteristics of a wind speed and considers a continuous evolution of $P_x$. This method was employed here to control the wind speed on the upper boundary but can be easily generalized to arbitrary heights in the domain.

The evolution equation built is a second-order PDE, also known as the rigid-body equation of motion, employed to dynamically filter the fluid momentum conservation equations. The mass-spring-damper system is critically damped ($\zeta = 1$) to preserve stability. Large natural periods ($T_P$) ensure that $P_x$ responds to motions occurring at scales much larger than the numerical domain. In flat-bottom applications (Section 4), this guided the convergence of freestream velocities to a certain target $u_{ref}$, over periods comparable to $T_P$. The convergence to a steady state regime was favored, and the normalized turbulence statistics were shown to be invariant to different values of $u_{ref}$ (setup ranging by $\pm 25\%$).

### 6.5. WI Disturbances and WA Parametrization

In the wavy bottom applications exploited here in Section 5, the quasi-steady state was maintained in the freestream with better than 3% accuracy. This allowed a systematic WA variation study ($48 \le WA_{ref} \le 79$) by varying the wind speed, in old seas, and under the prescription of a single monochromatic wave with moderate to high steepness ($ka = 0.2$). This canonical scenario was chosen for demonstration purposes, as it provides a critical evaluation of WI momentum stresses in the atmosphere.

The WI disturbances are usually assumed as a function of the WA parameter. However, many different definitions have been employed in the literature to define the WA, either
### 6.3. Evaluating Standard LES Formulations of the MABL

Considering the simplified framework adopted, the results for a standard LES formulation are presented in Section 3 as Cases 01 (flat-bottom reference) and 02 (wavy bottom reference). Case 02 reveals a certain number of artificial features unlikely to occur in the physical MABL. Notably, the mean velocities increased to the upper boundary, and the freestream velocity is greatly augmented due to the incidence of the wave. The freestream velocity increases up to a moment where it drags back the propagating wave. The downward momentum transfer that occurs is characteristic of young seas, even though we attempted to reproduce old seas conditions in this simulation. The normalized total shear stress converges towards the flat-bottom case. These processes reveal a long transient regime, while the flow adapts all the way up to the top boundary. Most interestingly and hardly a coincidence, the freestream self-adapts towards a condition where its velocity matches the wave phase velocity ($U_h/c = 99.5\%$), thus bringing Miles' critical height (c.f. [3,74,75]) inside the numerical domain.

Considering the dynamic pressure gradient proposed herein and employed in Case 03 of Section 3, these artificial features disappear. The friction velocity self-adapts to the appearance of the wave, but the freestream velocity is kept constant. A steady state regime is rapidly reached. The momentum flux occurs exclusively upwards. The WI disturbance in the wind velocity profile is kept confined to lower altitudes.

### 6.4. Large-Scale Dynamically Evolving Pressure Gradient Modeler

In the physical MABL and as long as the boundary layer height is sufficiently distant from the ocean (in neutral or unstable stratification conditions), it is more likely that the large-scale pressure gradient (traveling instantaneously in incompressible flows) will adapt to maintain a freestream velocity, rather than that, the freestream velocity will vary to maintain a constant large-scale pressure gradient condition.

In order to maintain an undisturbed freestream, while resolving arbitrary WI stresses, we developed a general method to control the wind speed at a certain height in the numerical domain (Section 2.2). To achieve this, an artificial evolution equation was built for the spatially-homogeneous large-scale pressure gradient: $P_x = \partial P/\partial x$. As an innovation in the wind-wave interaction studies, $P_x$ is thus allowed to evolve in time, restricted to scales much larger than the resolved ones. Unlike other solutions developed for inland applications [61,62,64,66], this methodology does not require previous knowledge about the instantaneous characteristics of a wind speed and considers a continuous evolution of $P_x$. This method was employed here to control the wind speed on the upper boundary but can be easily generalized to arbitrary heights in the domain.

The evolution equation built is a second-order PDE, also known as the rigid-body equation of motion, employed to dynamically filter the fluid momentum conservation equations. The mass-spring-damper system is critically damped ($\zeta = 1$) to preserve stability. Large natural periods ($T_P$) ensure that $P_x$ responds to motions occurring at scales much larger than the numerical domain. In flat-bottom applications (Section 4), this guided the convergence of freestream velocities to a certain target $u_{ref}$, over periods comparable to $T_P$. The convergence to a steady state regime was favored, and the normalized turbulence statistics were shown to be invariant to different values of $u_{ref}$ (setup ranging by $\pm 25\%$).

### 6.5. WI Disturbances and WA Parametrization

In the wavy bottom applications exploited here in Section 5, the quasi-steady state was maintained in the freestream with better than 3% accuracy. This allowed a systematic WA variation study ($48 \le WA_{ref} \le 79$) by varying the wind speed, in old seas, and under the prescription of a single monochromatic wave with moderate to high steepness ($ka = 0.2$). This canonical scenario was chosen for demonstration purposes, as it provides a critical evaluation of WI momentum stresses in the atmosphere.

The WI disturbances are usually assumed as a function of the WA parameter. However, many different definitions have been employed in the literature to define the WA, either
in terms of the friction velocity or the wind speed measured at an arbitrary height [79]. We have shown that when $P_x$ is kept constant, the friction velocity converges towards undisturbed values in LES. Employing the dynamic $P_x$ modeler, on the contrary, the friction velocity exhibits a strong dependency on the WA. Therefore if the WA is defined as a function of the resolved friction velocity, one establishes a circular dependency between these parameters. This problem is likely to occur in the physical MABL, as there is no reason to assume that the friction velocity remains constant with a changing sea-state condition.

A definite parametrization of WI disturbances remains beyond the scope of the present study, but we performed these studies to attest to the self-consistency of the methodology presented. To evaluate the behavior of different WA definitions, we fitted the form drag as either a function of the reference or the resolved WA. The $WA_{ref}$ characterizing the non-disturbed or freestream flow was found to be more accurate to describe the observations. When evaluating the turbulent profiles, this parametrization also provided a reasonable linear approximation of the height-dependent shear stresses.

Other observations are in accordance with the literature and attest to the robustness of the methodology proposed. The region of validity of the logarithmic wind profiles displaces upwards with increasing WA [19,80]. In the outer surface layer the effective roughness length decreases with increasing WA [34,37]. Turbulence intensities (TI) are increased in the WBL due to the WI fluctuations and diminish above the WBL due to the downward transport of turbulent fluctuations [19,77].

## 7. Conclusions

We employed LES to study micro-scale wind-wave interactions, proposing a new method to control the wind speed above arbitrary stresses in the bottom. When the sea-state condition is modified in LES, the wind speed at the top of the domain ($U_h$) self-adapts to the changes in wave-induced (WI) stresses that carry momentum downwards (wind-to-wave) or upwards (wave-to-wind) through the domain. However, a constant $U_h$ is desired to correctly compare different test cases submitted to an undisturbed freestream condition.

In order to maintain an undisturbed freestream while resolving arbitrary WI stresses, we developed a general method to control the wind speed at a certain height in the numerical domain. To achieve this, an evolution equation was built for the spatially-homogeneous large-scale pressure gradient: $P_x$. In flat-bottom applications, this guided the convergence of freestream velocities to a certain target $u_{ref}$. The convergence to a steady state regime was favored, and the normalized turbulence statistics were shown to be invariant to different values of $u_{ref}$ (setup ranging by ±25%). In wavy bottom applications, the quasi-steady state was maintained in the freestream with better than 3% accuracy. This allowed a systematic WA variation study ($48 \le WA_{ref} \le 79$) by varying the wind speed, in old seas, and under the prescription of a single monochromatic wave. This canonical scenario provides a critical evaluation of WI momentum stresses in the atmosphere. This methodology provides self-consistent WA parametrizations from LES and enables valuable analyses of wind-wave interactions.

The proposed parametrizations present a qualitative description of physical processes and trends. They are not intended to be taken in a quantitative manner due to the limited resources in terms of grid resolution. Future work will address more realistic conditions and finer grid resolutions, in order to present more quantitative descriptions of WI disturbances from LES. Moreover, in future work, the dynamic large-scale forcing modeler will be employed to enhance the comparison of physical and numerical experiments, by employing the wind speed control at the measurement height, rather than at the upper boundary.

**Author Contributions:** Conceptualization, methodology, validation, investigation, L.P., B.C. and Y.P.; formal analysis, data curation, L.P., B.C.; resources, B.C. and Y.P.; writing—original draft preparation, L.P.; writing—review and editing, L.P., B.C., Y.P. and S.A.; project administration, funding acquisition, S.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was carried out within the framework of the WEAMEC, West Atlantic Marine Energy Community, and with funding from the Pays de la Loire Region. This work was developed
in terms of the friction velocity or the wind speed measured at an arbitrary height [79]. We have shown that when $P_x$ is kept constant, the friction velocity converges towards undisturbed values in LES. Employing the dynamic $P_x$ modeler, on the contrary, the friction velocity exhibits a strong dependency on the WA. Therefore if the WA is defined as a function of the resolved friction velocity, one establishes a circular dependency between these parameters. This problem is likely to occur in the physical MABL, as there is no reason to assume that the friction velocity remains constant with a changing sea-state condition.

A definite parametrization of WI disturbances remains beyond the scope of the present study, but we performed these studies to attest to the self-consistency of the methodology presented. To evaluate the behavior of different WA definitions, we fitted the form drag as either a function of the reference or the resolved WA. The $WA_{ref}$ characterizing the non-disturbed or freestream flow was found to be more accurate to describe the observations. When evaluating the turbulent profiles, this parametrization also provided a reasonable linear approximation of the height-dependent shear stresses.

Other observations are in accordance with the literature and attest to the robustness of the methodology proposed. The region of validity of the logarithmic wind profiles displaces upwards with increasing WA [19,80]. In the outer surface layer the effective roughness length decreases with increasing WA [34,37]. Turbulence intensities (TI) are increased in the WBL due to the WI fluctuations and diminish above the WBL due to the downward transport of turbulent fluctuations [19,77].

## 7. Conclusions

We employed LES to study micro-scale wind-wave interactions, proposing a new method to control the wind speed above arbitrary stresses in the bottom. When the sea-state condition is modified in LES, the wind speed at the top of the domain ($U_h$) self-adapts to the changes in wave-induced (WI) stresses that carry momentum downwards (wind-to-wave) or upwards (wave-to-wind) through the domain. However, a constant $U_h$ is desired to correctly compare different test cases submitted to an undisturbed freestream condition.

In order to maintain an undisturbed freestream while resolving arbitrary WI stresses, we developed a general method to control the wind speed at a certain height in the numerical domain. To achieve this, an evolution equation was built for the spatially-homogeneous large-scale pressure gradient: $P_x$. In flat-bottom applications, this guided the convergence of freestream velocities to a certain target $u_{ref}$. The convergence to a steady state regime was favored, and the normalized turbulence statistics were shown to be invariant to different values of $u_{ref}$ (setup ranging by ±25%). In wavy bottom applications, the quasi-steady state was maintained in the freestream with better than 3% accuracy. This allowed a systematic WA variation study ($48 \le WA_{ref} \le 79$) by varying the wind speed, in old seas, and under the prescription of a single monochromatic wave. This canonical scenario provides a critical evaluation of WI momentum stresses in the atmosphere. This methodology provides self-consistent WA parametrizations from LES and enables valuable analyses of wind-wave interactions.

The proposed parametrizations present a qualitative description of physical processes and trends. They are not intended to be taken in a quantitative manner due to the limited resources in terms of grid resolution. Future work will address more realistic conditions and finer grid resolutions, in order to present more quantitative descriptions of WI disturbances from LES. Moreover, in future work, the dynamic large-scale forcing modeler will be employed to enhance the comparison of physical and numerical experiments, by employing the wind speed control at the measurement height, rather than at the upper boundary.

**Author Contributions:** Conceptualization, methodology, validation, investigation, L.P., B.C. and Y.P.; formal analysis, data curation, L.P., B.C.; resources, B.C. and Y.P.; writing—original draft preparation, L.P.; writing—review and editing, L.P., B.C., Y.P. and S.A.; project administration, funding acquisition, S.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was carried out within the framework of the WEAMEC, West Atlantic Marine Energy Community, and with funding from the Pays de la Loire Region. This work was developed
in the context of a PhD thesis, funded by research grants from the French Ministry of Research and Higher Education and École Centrale de Nantes.

**Data Availability Statement:** The data presented in this study are available upon request from the
corresponding author.

**Acknowledgments:** This work employs the atmospheric, large eddy simulation solver developed by Peter Sullivan, from the NCAR (National Center for Atmospheric Research) in Colorado, USA. This work was granted access to the HPC resources of CINES under the location A0090112000 made by GENCI (Grand Equipement National de Calcul Intensif).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

1. Thomson, W. Hydrokinetic Solutions and Observations. Philos. Mag. **1871**, 42, 326–378. [CrossRef]

2. Jeffreys, H. On the formation of water waves by wind. Proc. R. Soc. Lond. **1925**, 107, 189–206.

3. Miles, J.W. On The Generation Of Surface Waves By Shear Flows, PartI. J. Fluid Mech. **1957**, 3, 185–204. [CrossRef]

4. Phillips, O.M. On The Generation Of Waves By Turbulent Wind. J. Fluid Mech. **1957**, 2, 417–445. [CrossRef]

5. Belcher, S.E.; Hunt, J.C.R. Turbulent Shear Flow Over Slowly Moving Waves. J. Fluid Mech. **1993**, 251, 109–148. [CrossRef]

6. Edson, J.B.; Zappa, C.J.; Ware, J.A.; McGillis, W.R.; Hare, J.E. Scalar flux profile relationships over the open ocean. J. Geophys. Res. Ocean. **2004**, 109, C8. [CrossRef]

7. Grare, L.; Lenain, L.; Melville, W.K. Wave-Coherent Airflow and Critical Layers over Ocean Waves. J. Phys. Oceanogr. **2013**, 43, 2156–2172. [CrossRef]

8. Bourras, D.; Branger, H.; Reverdin, G.; Marié, L.; Cambra, R.; Baggio, L.; Caudoux, C.; Caudal, G.; Morisset, S.; Geyskens, N.; et al. A New Platform for the Determination of Air–Sea Fluxes (OCARINA): Overview and First Results. J. Atmos. Ocean. Technol. **2014**, 31, 1043–1062. [CrossRef]

9. Tamura, H.; Drennan, W.M.; Collins, C.O.; Graber, H.C. Turbulent Airflow and Wave-Induced Stress Over the Ocean. Bound.-Layer Meteorol. **2018**, 169, 47–66. [CrossRef]

10. Paskin, L.; Conan, B.; Perignon, Y.; Aubrun, S. Evidence of Ocean Waves Signature in the Space–Time Turbulent Spectra of the Lower Marine Atmosphere Measured by a Scanning LiDAR. Remote Sens. **2022**, 14, 3007. [CrossRef]

11. Hristov, T. Mechanistic, empirical and numerical perspectives on wind-waves interaction. Procedia IUTAM **2018**, 26, 102–111. [CrossRef]

12. Harris, J.A.; Belcher, S.E.; Street, R.L. Linear dynamics of wind waves in coupled turbulent air-water flow. Part 2. Numerical model. J. Fluid Mech. **1996**, 308, 219–254. [CrossRef]

13. Ayet, A.; Chapron, B. The Dynamical Coupling of Wind-Waves and Atmospheric Turbulence: A Review of Theoretical and Phenomenological Models. Bound.-Layer Meteorol. **2022**, 183, 1–33. [CrossRef]

14. Yousefi, K.; Veron, F. Boundary layer formulations in orthogonal curvilinear coordinates for flow over wind-generated surface waves. J. Fluid Mech. **2020**, 888, A11. [CrossRef]

15. Ayet, A. Air-Sea Momentum Fluxes in the Vicinity of the Sea Surface: A Theoretical Study of the Coupling between Turbulence and Wind-Waves. Ph.D. Thesis, Université de Bretagne Occidentale, Brest, France, 2020.

16. Smedman, A.; Högström, U.; Bergström, H.; Rutgersson, A.; Kahma, K.K.; Pettersson, H. A case study of air-sea interaction during swell conditions. J. Geophys. Res. Ocean. **1999**, 104, 25833–25851. [CrossRef]

17. Harris, D.L. The Wave-Driven Wind. J. Atmos. Sci. **1966**, 23, 688–693. [CrossRef]

18. Hanley, K.E.; Belcher, S.E. Wave-Driven Wind Jets in the Marine Atmospheric Boundary Layer. J. Atmos. Sci. **2008**, 65, 2646–2660. [CrossRef]

19. Sjöblom, A.; Smedman, A.S. The turbulent kinetic energy budget in the marine atmospheric surface layer. J. Geophys. Res. Ocean. **2002**, 107, 6–18. [CrossRef]

20. Cifuentes-Lorenzen, A.; Edson, J.B.; Zappa, C.J. Air–Sea Interaction in the Southern Ocean: Exploring the Height of the Wave Boundary Layer at the Air–Sea Interface. Bound.-Layer Meteorol. **2018**, 169, 461–482. [CrossRef]

21. Chalikov, D.; Babanin, A.V. Parameterization of Wave Boundary Layer. Atmosphere **2019**, 10, 686. [CrossRef]

22. Semedo, A.; Saetra, Ø.; Rutgersson, A.; Kahma, K.K.; Pettersson, H. Wave-Induced Wind in the Marine Boundary Layer. J. Atmos. Sci. **2009**, 66, 2256–2271. [CrossRef]

23. Liu, C.; Li, X.; Song, J.; Zou, Z.; Huang, J.; Zhang, J.A.; Jie, G.; Wang, J. Characteristics of the Marine Atmospheric Boundary Layer under the Influence of Ocean Surface Waves. J. Phys. Oceanogr. **2022**, 52, 1261–1276. [CrossRef]

24. Janssen, P.A.E.M.; Breivik, Ø.; Mogensen, K.; Vitart, F.; Alonso-Balmaseda, M.; Bidlot, J.R.; Keeley, S.; Leutbecher, M.; Magnusson, L.; Molteni, F. Air-sea interaction and surface waves. ECMWF Tech. Memo. **2013**, 712, 1–34.

25. Li, H.; Claremar, B.; Wu, L.; Hallgren, C.; Körnich, H.; Ivanell, S.; Sahlée, E. A sensitivity study of the WRF model in offshore wind modeling over the Baltic Sea. Geosci. Front. **2021**, 12, 101229. [CrossRef]
in the context of a PhD thesis, funded by research grants from the French Ministry of Research and Higher Education and École Centrale de Nantes.

**Data Availability Statement:** The data presented in this study are available upon request from the corresponding author.

**Acknowledgments:** This work employs the atmospheric, large eddy simulation solver developed by Peter Sullivan, from the NCAR (National Center for Atmospheric Research) in Colorado, USA. This work was granted access to the HPC resources of CINES under the location A0090112000 made by GENCI (Grand Equipement National de Calcul Intensif).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## References

1. Thomson, W. Hydrokinetic Solutions and Observations. *Philos. Mag.* **1871**, *42*, 326–378. [CrossRef]

2. Jeffreys, H. On the formation of water waves by wind. *Proc. R. Soc. Lond.* **1925**, *107*, 189–206.

3. Miles, J.W. On The Generation Of Surface Waves By Shear Flows, PartI. *J. Fluid Mech.* **1957**, *3*, 185–204. [CrossRef]

4. Phillips, O.M. On The Generation Of Waves By Turbulent Wind. *J. Fluid Mech.* **1957**, *2*, 417–445. [CrossRef]

5. Belcher, S.E.; Hunt, J.C.R. Turbulent Shear Flow Over Slowly Moving Waves. *J. Fluid Mech.* **1993**, *251*, 109–148. [CrossRef]

6. Edson, J.B.; Zappa, C.J.; Ware, J.A.; McGillis, W.R.; Hare, J.E. Scalar flux profile relationships over the open ocean. *J. Geophys. Res. Ocean.* **2004**, *109*, C8. [CrossRef]

7. Grare, L.; Lenain, L.; Melville, W.K. Wave-Coherent Airflow and Critical Layers over Ocean Waves. *J. Phys. Oceanogr.* **2013**, *43*, 2156–2172. [CrossRef]

8. Bourras, D.; Branger, H.; Reverdin, G.; Marié, L.; Cambra, R.; Baggio, L.; Caudoux, C.; Caudal, G.; Morisset, S.; Geyskens, N.; et al. A New Platform for the Determination of Air-Sea Fluxes (OCARINA): Overview and First Results. *J. Atmos. Ocean. Technol.* **2014**, *31*, 1043–1062. [CrossRef]

9. Tamura, H.; Drennan, W.M.; Collins, C.O.; Graber, H.C. Turbulent Airflow and Wave-Induced Stress Over the Ocean. *Bound.-Layer Meteorol.* **2018**, *169*, 47–66. [CrossRef]

10. Paskin, L.; Conan, B.; Perignon, Y.; Aubrun, S. Evidence of Ocean Waves Signature in the Space-Time Turbulent Spectra of the Lower Marine Atmosphere Measured by a Scanning LiDAR. *Remote Sens.* **2022**, *14*, 3007. [CrossRef]

11. Hristov, T. Mechanistic, empirical and numerical perspectives on wind-waves interaction. *Procedia IUTAM* **2018**, *26*, 102–111. [CrossRef]

12. Harris, J.A.; Belcher, S.E.; Street, R.L. Linear dynamics of wind waves in coupled turbulent air-water flow. Part 2. Numerical model. *J. Fluid Mech.* **1996**, *308*, 219–254. [CrossRef]

13. Ayet, A.; Chapron, B. The Dynamical Coupling of Wind-Waves and Atmospheric Turbulence: A Review of Theoretical and Phenomenological Models. *Bound.-Layer Meteorol.* **2022**, *183*, 1–33. [CrossRef]

14. Yousefi, K.; Veron, F. Boundary layer formulations in orthogonal curvilinear coordinates for flow over wind-generated surface waves. *J. Fluid Mech.* **2020**, *888*, A11. [CrossRef]

15. Ayet, A. Air-Sea Momentum Fluxes in the Vicinity of the Sea Surface: A Theoretical Study of the Coupling between Turbulence and Wind-Waves. Ph.D. Thesis, Université de Bretagne Occidentale, Brest, France, 2020.

16. Smedman, A.; Högström, U.; Bergström, H.; Rutgersson, A.; Kahma, K.K.; Pettersson, H. A case study of air-sea interaction during swell conditions. *J. Geophys. Res. Ocean.* **1999**, *104*, 25833–25851. [CrossRef]

17. Harris, D.L. The Wave-Driven Wind. *J. Atmos. Sci.* **1966**, *23*, 688–693. [CrossRef]

18. Hanley, K.E.; Belcher, S.E. Wave-Driven Wind Jets in the Marine Atmospheric Boundary Layer. *J. Atmos. Sci.* **2008**, *65*, 2646–2660. [CrossRef]

19. Sjöblom, A.; Smedman, A.S. The turbulent kinetic energy budget in the marine atmospheric surface layer. *J. Geophys. Res. Ocean.* **2002**, *107*, 6–18. [CrossRef]

20. Cifuentes-Lorenzen, A.; Edson, J.B.; Zappa, C.J. Air-Sea Interaction in the Southern Ocean: Exploring the Height of the Wave Boundary Layer at the Air-Sea Interface. *Bound.-Layer Meteorol.* **2018**, *169*, 461–482. [CrossRef]

21. Chalikov, D.; Babanin, A.V. Parameterization of Wave Boundary Layer. *Atmosphere* **2019**, *10*, 686. [CrossRef]

22. Semedo, A.; Saetra, Ø.; Rutgersson, A.; Kahma, K.K.; Pettersson, H. Wave-Induced Wind in the Marine Boundary Layer. *J. Atmos. Sci.* **2009**, *66*, 2256–2271. [CrossRef]

23. Liu, C.; Li, X.; Song, J.; Zou, Z.; Huang, J.; Zhang, J.A.; Jie, G.; Wang, J. Characteristics of the Marine Atmospheric Boundary Layer under the Influence of Ocean Surface Waves. *J. Phys. Oceanogr.* **2022**, *52*, 1261–1276. [CrossRef]

24. Janssen, P.A.E.M.; Breivik, Ø.; Mogensen, K.; Vitart, F.; Alonso-Balmaseda, M.; Bidlot, J.R.; Keeley, S.; Leutbecher, M.; Magnusson, L.; Molteni, F. Air-sea interaction and surface waves. *ECMWF Tech. Memo.* **2013**, *712*, 1–34.

25. Li, H.; Claremar, B.; Wu, L.; Hallgren, C.; Körnich, H.; Ivanell, S.; Sahlée, E. A sensitivity study of the WRF model in offshore wind modeling over the Baltic Sea. *Geosci. Front.* **2021**, *12*, 101229. [CrossRef]
26. Rezaeiha, A.; Pereira, R.; Kotsonis, M. Fluctuations of angle of attack and lift coefficient and the resultant fatigue loads for a large Horizontal Axis Wind turbine. *Renew. Energy* **2017**, *114*, 904–916. [CrossRef]

27. Leroy, V.; Gilloteaux, J.C.; Lynch, M.; Babarit, A.; Ferrant, P. Impact of aerodynamic modeling on seakeeping performance of a floating horizontal axis wind turbine. *Wind Energy* **2019**, *22*, 1019–1033. [CrossRef]

28. Kalvig, S.; Gudmestad, O.T.; Winther, N. Exploring the gap between 'best knowledge' and 'best practice' in boundary layer meteorology for offshore wind energy. *Wind Energy* **2014**, *17*, 161–171. [CrossRef]

29. Obukhov, A.M. Turbulence in an Atmosphere with a Non-Uniform Temperature. *Bound.-Layer Meteorol.* **1946**, *2*, 7–29. [CrossRef]

30. Monin, A.S.; Obukhov, A.M. Basic Laws of Turbulent Mixing in the Surface Layer of the Atmosphere. *Contrib. Geophys. Inst. Acad. Sci. USSR* **1954**, *24*, 163–187.

31. Edson, J.B.; Fairall, C.W. Similarity Relationships in the Marine Atmospheric Surface Layer for Terms in the TKE and Scalar Variance Budgets. *J. Atmos. Sci.* **1998**, *55*, 2311–2328. [CrossRef]

32. Hristov, T.; Ruiz-Plancarte, J. Dynamic Balances in a Wavy Boundary Layer. *J. Phys. Oceanogr.* **2014**, *44*, 3185–3194. [CrossRef]

33. Charnock, H. Wind stress on a water surface. *Q. J. R. Meteorol. Soc.* **1955**, *81*, 639–640. [CrossRef]

34. Donelan, M.A.; Dobson, F.W.; Smith, S.D.; Anderson, R.J. On the Dependence of Sea Surface Roughness on Wave Development. *J. Phys. Oceanogr.* **1993**, *23*, 2143–2149. [CrossRef]

35. Edson, J.B.; Jampana, V.; Weller, R.A.; Bigorre, S.P.; Plueddemann, A.J.; Fairall, C.W.; Miller, S.D.; Mahrt, L.; Vickers, D.; Hersbach, H. On the Exchange of Momentum over the Open Ocean. *J. Phys. Oceanogr.* **2013**, *43*, 1589–1610. [CrossRef]

36. Patton, E.G.; Sullivan, P.P.; Kosović, B.; Dudhia, J.; Mahrt, L.; Žagar, M.; Marić, T. On the Influence of Swell Propagation Angle on Surface Drag. *J. Appl. Meteorol. Climatol.* **2019**, *58*, 1039–1059. [CrossRef]

37. Porchetta, S.; Temel, O.; Muñoz Esparza, D.; Reuder, J.; Monbaliu, J.; van Beeck, J.; van Lipzig, N. A new roughness length parameterization accounting for wind-wave (mis)alignment. *Atmos. Chem. Phys.* **2019**, *19*, 6681–6700. [CrossRef]

38. Sullivan, P.P.; McWilliams, J.C.; Moeng, C.H. Simulation of turbulent flow over idealized water waves. *J. Fluid Mech.* **2000**, *404*, 47–85. [CrossRef]

39. Yang, D.; Shen, L. Characteristics of coherent vortical structures in turbulent flows over progressive surface waves. *Phys. Fluids* **2009**, *21*, 125106. [CrossRef]

40. Yang, D.; Shen, L. Direct-simulation-based study of turbulent flow over various waving boundaries. *J. Fluid Mech.* **2010**, *650*, 131–180. [CrossRef]

41. Yang, D.; Shen, L. Simulation of viscous flows with undulatory boundaries. Part I: Basic solver. *J. Comput. Phys.* **2011**, *230*, 5488–5509. [CrossRef]

42. Druzhinin, O.A.; Troitskaya, Y.I.; Zilitinkevich, S.S. Direct numerical simulation of a turbulent wind over a wavy water surface. *J. Geophys. Res. Ocean.* **2012**, *117*, C11. [CrossRef]

43. Yang, Z.; Deng, B.Q.; Shen, L. Direct numerical simulation of wind turbulence over breaking waves. *J. Fluid Mech.* **2018**, *850*, 120–155. [CrossRef]

44. Deskos, G.; Ananthan, S.; Sprague, M.A. Direct numerical simulations of turbulent flow over misaligned traveling waves. *Int. J. Heat Fluid Flow* **2022**, *97*, 109029. [CrossRef]

45. Stoll, R.; Gibbs, J.A.; Salesky, S.T.; Anderson, W.; Calaf, M. Large-Eddy Simulation of the Atmospheric Boundary Layer. *Bound.-Layer Meteorol.* **2020**, *177*, 541–581.

46. Chanprasert, W.; Sharma, R.; Cater, J.; Norris, S. Large Eddy Simulation of wind turbine fatigue loading and yaw dynamics induced by wake turbulence. *Renew. Energy* **2022**, *190*, 208–222. [CrossRef]

47. Henn, D.S.; Sykes, R.I. Large-eddy simulation of flow over wavy surfaces. *J. Fluid Mech.* **1999**, *383*, 75–112. [CrossRef]

48. Sullivan, P.P.; Edson, J.B.; Hristov, T.; McWilliams, J.C. Large-Eddy Simulations and Observations of Atmospheric Marine Boundary Layers above Nonequilibrium Surface Waves. *J. Atmos. Sci.* **2008**, *65*, 1225–1245. [CrossRef]

49. Sullivan, P.P.; McWilliams, J.C.; Patton, E.G. Large-Eddy Simulation of Marine Atmospheric Boundary Layers above a Spectrum of Moving Waves. *J. Atmos. Sci.* **2014**, *71*, 4001–4027. [CrossRef]

50. Hara, T.; Sullivan, P.P. Wave Boundary Layer Turbulence over Surface Waves in a Strongly Forced Condition. *J. Phys. Oceanogr.* **2015**, *45*, 868–883. [CrossRef]

51. Cathelain, M. Development of a Deterministic Numerical Model for the Study of the Coupling between an Atmospheric Flow and a Sea State. Ph.D. Thesis, Ecole Centrale de Nantes (ECN), Nantes, France, 2017.

52. Paskin, L.; Perignon, Y.; Conan, B.; Aubrun, S. Numerical study on the Wave Boundary Layer, its interaction with turbulence and consequences on the wind energy resource in the offshore environment. *J. Phys. Conf. Ser.* **2020**, *1618*, 062046. [CrossRef]

53. Yang, D.; Meneveau, C.; Shen, L. Dynamic modelling of sea-surface roughness for large-eddy simulation of wind over ocean wavefield. *J. Fluid Mech.* **2013**, *726*, 62–99. [CrossRef]

54. Yang, D.; Meneveau, C.; Shen, L. Effect of downwind swells on offshore wind energy harvesting—A large-eddy simulation study. *Renew. Energy* **2014**, *70*, 11–23. [CrossRef]

55. Yang, D.; Meneveau, C.; Shen, L. Large-eddy simulation of offshore wind farm. *Phys. Fluids* **2014**, *26*, 025101. [CrossRef]

56. Hao, X.; Shen, L. Wind-wave coupling study using LES of wind and phase-resolved simulation of nonlinear waves. *J. Fluid Mech.* **2019**, *874*, 391–425. [CrossRef]

57. Hao, X.; Cao, T.; Shen, L. Mechanistic study of shoaling effect on momentum transfer between turbulent flow and traveling wave using large-eddy simulation. *Phys. Rev. Fluids* **2021**, *6*, 054608. [CrossRef]
26. Rezaeiha, A.; Pereira, R.; Kotsonis, M. Fluctuations of angle of attack and lift coefficient and the resultant fatigue loads for a large Horizontal Axis Wind turbine. *Renew. Energy* **2017**, *114*, 904–916. [CrossRef]

27. Leroy, V.; Gilloteaux, J.C.; Lynch, M.; Babarit, A.; Ferrant, P. Impact of aerodynamic modeling on seakeeping performance of a floating horizontal axis wind turbine. *Wind Energy* **2019**, *22*, 1019–1033. [CrossRef]

28. Kalvig, S.; Gudmestad, O.T.; Winther, N. Exploring the gap between 'best knowledge' and 'best practice' in boundary layer meteorology for offshore wind energy. *Wind Energy* **2014**, *17*, 161–171. [CrossRef]

29. Obukhov, A.M. Turbulence in an Atmosphere with a Non-Uniform Temperature. *Bound.-Layer Meteorol.* **1946**, *2*, 7–29. [CrossRef]

30. Monin, A.S.; Obukhov, A.M. Basic Laws of Turbulent Mixing in the Surface Layer of the Atmosphere. *Contrib. Geophys. Inst. Acad. Sci. USSR* **1954**, *24*, 163–187.

31. Edson, J.B.; Fairall, C.W. Similarity Relationships in the Marine Atmospheric Surface Layer for Terms in the TKE and Scalar Variance Budgets. *J. Atmos. Sci.* **1998**, *55*, 2311–2328. [CrossRef]

32. Hristov, T.; Ruiz-Plancarte, J. Dynamic Balances in a Wavy Boundary Layer. *J. Phys. Oceanogr.* **2014**, *44*, 3185–3194. [CrossRef]

33. Charnock, H. Wind stress on a water surface. *Q. J. R. Meteorol. Soc.* **1955**, *81*, 639–640. [CrossRef]

34. Donelan, M.A.; Dobson, F.W.; Smith, S.D.; Anderson, R.J. On the Dependence of Sea Surface Roughness on Wave Development. *J. Phys. Oceanogr.* **1993**, *23*, 2143–2149. [CrossRef]

35. Edson, J.B.; Jampana, V.; Weller, R.A.; Bigorre, S.P.; Plueddemann, A.J.; Fairall, C.W.; Miller, S.D.; Mahrt, L.; Vickers, D.; Hersbach, H. On the Exchange of Momentum over the Open Ocean. *J. Phys. Oceanogr.* **2013**, *43*, 1589–1610. [CrossRef]

36. Patton, E.G.; Sullivan, P.P.; Kosović, B.; Dudhia, J.; Mahrt, L.; Žagar, M.; Marić, T. On the Influence of Swell Propagation Angle on Surface Drag. *J. Appl. Meteorol. Climatol.* **2019**, *58*, 1039–1059. [CrossRef]

37. Porchetta, S.; Temel, O.; Muñoz Esparza, D.; Reuder, J.; Monbaliu, J.; van Beeck, J.; van Lipzig, N. A new roughness length parameterization accounting for wind-wave (mis)alignment. *Atmos. Chem. Phys.* **2019**, *19*, 6681–6700. [CrossRef]

38. Sullivan, P.P.; McWilliams, J.C.; Moeng, C.H. Simulation of turbulent flow over idealized water waves. *J. Fluid Mech.* **2000**, *404*, 47–85. [CrossRef]

39. Yang, D.; Shen, L. Characteristics of coherent vortical structures in turbulent flows over progressive surface waves. *Phys. Fluids* **2009**, *21*, 125106. [CrossRef]

40. Yang, D.; Shen, L. Direct-simulation-based study of turbulent flow over various waving boundaries. *J. Fluid Mech.* **2010**, *650*, 131–180. [CrossRef]

41. Yang, D.; Shen, L. Simulation of viscous flows with undulatory boundaries. Part I: Basic solver. *J. Comput. Phys.* **2011**, *230*, 5488–5509. [CrossRef]

42. Druzhinin, O.A.; Troitskaya, Y.I.; Zilitinkevich, S.S. Direct numerical simulation of a turbulent wind over a wavy water surface. *J. Geophys. Res. Ocean.* **2012**, *117*, C11. [CrossRef]

43. Yang, Z.; Deng, B.Q.; Shen, L. Direct numerical simulation of wind turbulence over breaking waves. *J. Fluid Mech.* **2018**, *850*, 120–155. [CrossRef]

44. Deskos, G.; Ananthan, S.; Sprague, M.A. Direct numerical simulations of turbulent flow over misaligned traveling waves. *Int. J. Heat Fluid Flow* **2022**, *97*, 109029. [CrossRef]

45. Stoll, R.; Gibbs, J.A.; Salesky, S.T.; Anderson, W.; Calaf, M. Large-Eddy Simulation of the Atmospheric Boundary Layer. *Bound.-Layer Meteorol.* **2020**, *177*, 541–581.

46. Chanprasert, W.; Sharma, R.; Cater, J.; Norris, S. Large Eddy Simulation of wind turbine fatigue loading and yaw dynamics induced by wake turbulence. *Renew. Energy* **2022**, *190*, 208–222. [CrossRef]

47. Henn, D.S.; Sykes, R.I. Large-eddy simulation of flow over wavy surfaces. *J. Fluid Mech.* **1999**, *383*, 75–112. [CrossRef]

48. Sullivan, P.P.; Edson, J.B.; Hristov, T.; McWilliams, J.C. Large-Eddy Simulations and Observations of Atmospheric Marine Boundary Layers above Nonequilibrium Surface Waves. *J. Atmos. Sci.* **2008**, *65*, 1225–1245. [CrossRef]

49. Sullivan, P.P.; McWilliams, J.C.; Patton, E.G. Large-Eddy Simulation of Marine Atmospheric Boundary Layers above a Spectrum of Moving Waves. *J. Atmos. Sci.* **2014**, *71*, 4001–4027. [CrossRef]

50. Hara, T.; Sullivan, P.P. Wave Boundary Layer Turbulence over Surface Waves in a Strongly Forced Condition. *J. Phys. Oceanogr.* **2015**, *45*, 868–883. [CrossRef]

51. Cathelain, M. Development of a Deterministic Numerical Model for the Study of the Coupling between an Atmospheric Flow and a Sea State. Ph.D. Thesis, Ecole Centrale de Nantes (ECN), Nantes, France, 2017.

52. Paskin, L.; Perignon, Y.; Conan, B.; Aubrun, S. Numerical study on the Wave Boundary Layer, its interaction with turbulence and consequences on the wind energy resource in the offshore environment. *J. Phys. Conf. Ser.* **2020**, *1618*, 062046. [CrossRef]

53. Yang, D.; Meneveau, C.; Shen, L. Dynamic modelling of sea-surface roughness for large-eddy simulation of wind over ocean wavefield. *J. Fluid Mech.* **2013**, *726*, 62–99. [CrossRef]

54. Yang, D.; Meneveau, C.; Shen, L. Effect of downwind swells on offshore wind energy harvesting—A large-eddy simulation study. *Renew. Energy* **2014**, *70*, 11–23. [CrossRef]

55. Yang, D.; Meneveau, C.; Shen, L. Large-eddy simulation of offshore wind farm. *Phys. Fluids* **2014**, *26*, 025101. [CrossRef]

56. Hao, X.; Shen, L. Wind-wave coupling study using LES of wind and phase-resolved simulation of nonlinear waves. *J. Fluid Mech.* **2019**, *874*, 391–425. [CrossRef]

57. Hao, X.; Cao, T.; Shen, L. Mechanistic study of shoaling effect on momentum transfer between turbulent flow and traveling wave using large-eddy simulation. *Phys. Rev. Fluids* **2021**, *6*, 054608. [CrossRef]
58. Cao, T.; Shen, L. A numerical and theoretical study of wind over fast-propagating water waves. *J. Fluid Mech.* **2021**, *919*, A38.
[CrossRef]

59. Hao, X.; Shen, L. Large-eddy simulation of gusty wind turbulence over a travelling wave. *J. Fluid Mech.* **2022**, *946*, A8. [CrossRef]

60. Deskos, G.; Lee, J.C.Y.; Draxl, C.; Sprague, M.A. Review of Wind-Wave Coupling Models for Large-Eddy Simulation of the Marine Atmospheric Boundary Layer. *J. Atmos. Sci.* **2021**, *78*, 3025–3045. [CrossRef]

61. Pedersen, J.G.; Kelly, M.; Gryning, S.E.; Brümmer, B. The effect of unsteady and baroclinic forcing on predicted wind profiles in Large Eddy Simulations: Two case studies of the daytime atmospheric boundary layer. *Meteorologische Zeitschrift* **2013**, *22*, 661–674.
[CrossRef]

62. Basu, S.; Vinuesa, J.; Swift, A. Dynamic LES Modeling of a Diurnal Cycle. *J. Appl. Meteorol. Climatol.* **2008**, *47*, 1156–1174.
[CrossRef]

63. Nielson, J.; Bhaganagar, K. Capturing Day-to-day Diurnal Variations in Stability in the Convective Atmospheric Boundary Layer Using Large Eddy Simulation. *Open Atmos. Sci. J.* **2018**, *12*, 107–131. [CrossRef]

64. Allaerts, D.; Quon, E.; Draxl, C.; Churchfield, M. Development of a Time-Height Profile Assimilation Technique for Large-Eddy Simulation. *Bound.-Layer Meteorol.* **2020**, *176*, 329–348. [CrossRef]

65. Draxl, C.; Allaerts, D.; Quon, E.; Churchfield, M. Coupling Mesoscale Budget Components to Large-Eddy Simulations for Wind-Energy Applications. *Bound.-Layer Meteorol.* **2021**, *179*, 73–98. [CrossRef]

66. Pimont, F.; Dupuy, J.L.; Linn, R.R.; Sauer, J.A.; Muñoz-Esparza, D. Pressure-Gradient Forcing Methods for Large-Eddy Simulations of Flows in the Lower Atmospheric Boundary Layer. *Atmosphere* **2020**, *11*, 1343. [CrossRef]

67. Stieren, A.; Gadde, S.N.; Stevens, R.J. Modeling dynamic wind direction changes in large eddy simulations of wind farms. *Renew. Energy* **2021**, *170*, 1342–1352. [CrossRef]

68. Fenton, J.D. A fifth-order Stokes theory for steady waves. *J. Waterw. Port Coastal Ocean. Eng.* **1985**, *111*, 216–234. [CrossRef]

69. Spiegel, E.A.; Veronis, G. On The Boussinesq Approximation For A Compressible Fluid. *Astrophys. J.* **1960**, *131*, 442. [CrossRef]

70. Moeng, C.H. A Large-Eddy-Simulation Model for the Study of Planetary Boundary-Layer Turbulence. *J. Atmos. Sci.* **1984**, *41*, 2052–2062. [CrossRef]

71. Thomas, P.D.; Lombard, C.K. Geometric Conservation Law and Its Application to Flow Computations on Moving Grids. *AIAA J.* **1979**, *17*, 1030–1037. [CrossRef]

72. Kleppner, D.; Kolenkow, R. *An Introduction to Mechanics*, 1st ed.; Cambridge University Press: Cambridge, UK, 2010; pp. 410–440.

73. Georgi, H. *The Physics of Waves*, 1st ed.; Prentice Hall: Englewood Cliffs, NJ, USA, 1993; pp. 1–52.

74. Lighthill, M.J. Physical interpretation of the mathematical theory of wave generation by wind. *J. Fluid Mech.* **1962**, *14*, 385–398.
[CrossRef]

75. Hristov, T.S.; Miller, S.D.; Friehe, C.A. Dynamical coupling of wind and ocean waves through wave-induced air flow. *Nature* **2003**, *422*, 55–58. [CrossRef] [PubMed]

76. Donelan, M.; Babanin, A.; Young, I.; Banner, M. Wave-Follower Field Measurements of the Wind-Input Spectral Function. Part II: Parameterization of the Wind Input. *J. Phys. Oceanogr.* **2006**, *36*, 1672–1689. [CrossRef]

77. Högström, U.; Smedman, A.; Sahleé, E.; Drennan, W.M.; Kahma, K.K.; Pettersson, H.; Zhang, F. The Atmospheric Boundary Layer during Swell: A Field Study and Interpretation of the Turbulent Kinetic Energy Budget for High Wave Ages. *J. Atmos. Sci.* **2009**, *66*, 2764–2779. [CrossRef]

78. O'Sullivan, J.; Archer, R.; Flay, R. Consistent boundary conditions for flows within the atmospheric boundary layer. *J. Wind. Eng. Ind. Aerodyn.* **2011**, *99*, 65–77. [CrossRef]

79. Donelan, M.A.; Drennan, W.M.; Magnusson, A.K. Nonstationary Analysis of the Directional Properties of Propagating Waves. *J. Phys. Oceanogr.* **1996**, *26*, 1901–1914. [CrossRef]

80. Mastenbroek, C. Wind-Wave Interaction. Ph.D. Thesis, TU Delft, Delft, The Netherlands, 1996.
58. Cao, T.; Shen, L. A numerical and theoretical study of wind over fast-propagating water waves. *J. Fluid Mech.* **2021**, *919*, A38.
[CrossRef]

59. Hao, X.; Shen, L. Large-eddy simulation of gusty wind turbulence over a travelling wave. *J. Fluid Mech.* **2022**, *946*, A8. [CrossRef]

60. Deskos, G.; Lee, J.C.Y.; Draxl, C.; Sprague, M.A. Review of Wind-Wave Coupling Models for Large-Eddy Simulation of the Marine Atmospheric Boundary Layer. *J. Atmos. Sci.* **2021**, *78*, 3025–3045. [CrossRef]

61. Pedersen, J.G.; Kelly, M.; Gryning, S.E.; Brümmer, B. The effect of unsteady and baroclinic forcing on predicted wind profiles in Large Eddy Simulations: Two case studies of the daytime atmospheric boundary layer. *Meteorologische Zeitschrift* **2013**, *22*, 661–674.
[CrossRef]

62. Basu, S.; Vinuesa, J.; Swift, A. Dynamic LES Modeling of a Diurnal Cycle. *J. Appl. Meteorol. Climatol.* **2008**, *47*, 1156–1174.
[CrossRef]

63. Nielson, J.; Bhaganagar, K. Capturing Day-to-day Diurnal Variations in Stability in the Convective Atmospheric Boundary Layer Using Large Eddy Simulation. *Open Atmos. Sci. J.* **2018**, *12*, 107–131. [CrossRef]

64. Allaerts, D.; Quon, E.; Draxl, C.; Churchfield, M. Development of a Time-Height Profile Assimilation Technique for Large-Eddy Simulation. *Bound.-Layer Meteorol.* **2020**, *176*, 329–348. [CrossRef]

65. Draxl, C.; Allaerts, D.; Quon, E.; Churchfield, M. Coupling Mesoscale Budget Components to Large-Eddy Simulations for Wind-Energy Applications. *Bound.-Layer Meteorol.* **2021**, *179*, 73–98. [CrossRef]

66. Pimont, F.; Dupuy, J.L.; Linn, R.R.; Sauer, J.A.; Muñoz-Esparza, D. Pressure-Gradient Forcing Methods for Large-Eddy Simulations of Flows in the Lower Atmospheric Boundary Layer. *Atmosphere* **2020**, *11*, 1343. [CrossRef]

67. Stieren, A.; Gadde, S.N.; Stevens, R.J. Modeling dynamic wind direction changes in large eddy simulations of wind farms. *Renew. Energy* **2021**, *170*, 1342–1352. [CrossRef]

68. Fenton, J.D. A fifth-order Stokes theory for steady waves. *J. Waterw. Port Coastal Ocean. Eng.* **1985**, *111*, 216–234. [CrossRef]

69. Spiegel, E.A.; Veronis, G. On The Boussinesq Approximation For A Compressible Fluid. *Astrophys. J.* **1960**, *131*, 442. [CrossRef]

70. Moeng, C.H. A Large-Eddy-Simulation Model for the Study of Planetary Boundary-Layer Turbulence. *J. Atmos. Sci.* **1984**, *41*, 2052–2062. [CrossRef]

71. Thomas, P.D.; Lombard, C.K. Geometric Conservation Law and Its Application to Flow Computations on Moving Grids. *AIAA J.* **1979**, *17*, 1030–1037. [CrossRef]

72. Kleppner, D.; Kolenkow, R. *An Introduction to Mechanics*, 1st ed.; Cambridge University Press: Cambridge, UK, 2010; pp. 410–440.

73. Georgi, H. *The Physics of Waves*, 1st ed.; Prentice Hall: Englewood Cliffs, NJ, USA, 1993; pp. 1–52.

74. Lighthill, M.J. Physical interpretation of the mathematical theory of wave generation by wind. *J. Fluid Mech.* **1962**, *14*, 385–398.
[CrossRef]

75. Hristov, T.S.; Miller, S.D.; Friehe, C.A. Dynamical coupling of wind and ocean waves through wave-induced air flow. *Nature* **2003**, *422*, 55–58. [CrossRef] [PubMed]

76. Donelan, M.; Babanin, A.; Young, I.; Banner, M. Wave-Follower Field Measurements of the Wind-Input Spectral Function. Part II: Parameterization of the Wind Input. *J. Phys. Oceanogr.* **2006**, *36*, 1672–1689. [CrossRef]

77. Högström, U.; Smedman, A.; Sahleé, E.; Drennan, W.M.; Kahma, K.K.; Pettersson, H.; Zhang, F. The Atmospheric Boundary Layer during Swell: A Field Study and Interpretation of the Turbulent Kinetic Energy Budget for High Wave Ages. *J. Atmos. Sci.* **2009**, *66*, 2764–2779. [CrossRef]

78. O'Sullivan, J.; Archer, R.; Flay, R. Consistent boundary conditions for flows within the atmospheric boundary layer. *J. Wind. Eng. Ind. Aerodyn.* **2011**, *99*, 65–77. [CrossRef]

79. Donelan, M.A.; Drennan, W.M.; Magnusson, A.K. Nonstationary Analysis of the Directional Properties of Propagating Waves. *J. Phys. Oceanogr.* **1996**, *26*, 1901–1914. [CrossRef]

80. Mastenbroek, C. Wind-Wave Interaction. Ph.D. Thesis, TU Delft, Delft, The Netherlands, 1996.