††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Site and bond percolation on four-dimensional simple hypercubic lattice with
extended neighborhoods

Pengyu Zhao  Jinhong Yan  Zhipeng Xun<EMAIL_ADDRESS>Dapeng Hao  School of
Material Sciences and Physics, China University of Mining and Technology,
Xuzhou 221116, China  Robert M. Ziff<EMAIL_ADDRESS>Center for the Study of
Complex System and Department of Chemical Engineering, University of Michigan,
Ann Arbor, Michigan 48109-2800, USA

###### Abstract

The asymptotic behavior of the percolation threshold $p_{c}$ and its
dependence upon coordination number $z$ is investigated for both site and bond
percolation on four-dimensional lattices with compact neighborhoods. Simple
hypercubic lattices with extended neighborhoods, up to 9th nearest neighbors,
are studied to high precision by means of Monte-Carlo simulations based upon a
single-cluster growth algorithm. For site percolation, an asymptotic analysis
confirms the predicted behavior $zp_{c}\sim 16\eta_{c}=2.086$ for large $z$,
and finite-size corrections are accounted for by forms $p_{c}\sim
16\eta_{c}/(z+b)$ and $p_{c}\sim 1-\exp(-16\eta_{c}/z)$. For bond percolation,
the finite-$z$ correction is accounted for by $zp_{c}-1\simeq a_{1}z^{-x}$,
and data fitting gives $x\approx 3/4$. This result indicates that the relation
$x=(d-1)/d$, which has been shown to hold in two and three dimensions as a
surface-area effect, holds in four dimensions also.

###### pacs:

64.60.ah, 89.75.Fb, 05.70.Fh

##  I Introduction

Percolation on lattices with extended neighborhoods, which goes back to the
“equivalent neighbor model” of Dalton, Domb and Sykes in 1964
DaltonDombSykes64 ; DombDalton1966 ; Domb72 , has been a question of long-
standing interest in the percolation field, both because of its theoretical
significance and practical applications Kleinberg2000 ;
SanderWarrenSokolov2003 ; Ziff2021 ; KozaKondratSuszcaynski2014 ; KozaPola2016
. Following the work of Dalton, Domb and Sykes, many researchers have
performed extensive investigations, such as long-range site percolation on
compact regions in a diamond shape on a square lattice GoukerFamily83 , site
and bond percolation on body-centered cubic lattices with nearest and next-
nearest neighbors JerauldScrivenDavis1984 , site percolation on face-centered
cubic lattices up to fourth nearest neighbors GawronCieplak91 , site
percolation on all eleven of the Archimedian lattices with long-range
connections dIribarneRasigniRasigni95 ; dIribarneRasigniRasigni99 ;
dIribarneRasigniRasigni99b . The idea of “complex neighborhoods” where various
combinations of neighborhoods, not necessarily compact, was introduced by
Malarz and Galam MalarzGalam05 , and has been followed up by many subsequent
works in two, three, and four dimensions for both site percolation
MalarzGalam05 ; MajewskiMalarz2007 ; KurzawskiMalarz2012 ; Malarz2015 ;
KotwicaGronekMalarz19 ; Malarz2020 ; XunHaoZiff2021  and bond percolation
OuyangDengBlote2018 ; DengOuyangBlote2019 ; XunZiff2020 ; XunZiff2020b ;
XuWangHuDeng2021 .

Since the early days of the study of percolation, researchers have focused on
exploring the correlations between percolation thresholds $p_{c}$ and
properties of the lattice, especially the coordination number $z$ which gives
the number of nearest neighbors to a vertex, since $z$ (along with the
dimensionality) seems to be one of the main determinants of the threshold. In
the context of extended-range percolation, earlier investigations
GalamMauger1996 ; vanderMarck1998 ; KurzawskiMalarz2012 ; XunZiff2020 ;
XunZiff2020b  demonstrated that the percolation thresholds for lattices with
extended neighborhoods can be fitted well by a simple power given by
$p_{c}\sim z^{-a}$ or $p_{c}\sim(z-1)^{-a}$. It has been argued Domb72 ;
dIribarneRasigniRasigni99b ; KozaKondratSuszcaynski2014 ; KozaPola2016  that
for longer range site percolation, the asymptotic behavior for large $z$ could
be related to the continuum percolation threshold $\eta_{c}$ for objects of
the same shape as the neighborhood. In general $d$, these arguments imply
XunHaoZiff2021

$zp_{c}\sim 2^{d}\eta_{c},$ (1)

where $d$ is the dimension of the system. Equation (1) holds for systems with
compact neighborhoods, and is independent of the type of lattice. For finite
$z$, one expects asymptotic corrections of the form

$p_{c}\sim\frac{c}{z+b},$ (2)

where $c=2^{d}\eta_{c}$ and $b$ an empirical constant. Another proposed form
(without any additional constants) is XunHaoZiff2021

$p_{c}\sim 1-\exp\left(-\frac{2^{d}\eta_{c}}{z}\right).$ (3)

For bond percolation, the asymptotic behavior of $p_{c}$ tends to the Bethe-
lattice behavior

$p_{c}\sim\frac{1}{z-1},$ (4)

for $z\to\infty$, because for large $z$ and the corresponding low $p$, the
chance of hitting the same site twice is vanishingly small and the system
behaves basically like a tree. The finite-$z$ correction was theoretically
predicted FreiPerkins2016 ; XuWangHuDeng2021  to be

$zp_{c}-1\sim a_{1}z^{-x},$ (5)

where $x=(d-1)/d$ for $d=2,3$ (regarded as a surface area effect), while for
$d\geq 4$, it was suggested that $zp_{c}-1\sim a_{1}/z$, with logarithm
corrections in $d=4$ FreiPerkins2016 ; HofstadSakai2005 , implying that in
this case the asymptotic behavior of $zp_{c}$ is a bulk property. In two and
three dimensions, both the asymptotic behavior and the finite-$z$ correction
for site and bond percolation were confirmed by extensive simulation work
dIribarneRasigniRasigni99b ; OuyangDengBlote2018 ; DengOuyangBlote2019 ;
XunZiff2020b ; XunHaoZiff2021 ; XunHaoZiff2021b . However, more work is needed
to confirm the above asymptotic behavior for $d\geq 4$, and to understand the
difference of the correlation exponents in different dimensions.

In this paper, we focus on four-dimensional simple hypercubic (sc(4)) lattices
with extended neighborhoods. Both site and bond percolation on sc(4) lattices
with up to 9th nearest neighbors are investigated by employing Monte-Carlo
simulation, and we find site and bond percolation thresholds of these lattices
with high precision. For site percolation, both $z$ versus $1/p_{c}$ and $z$
versus $-1/\ln(1-p_{c})$ lead to the predicted asymptotic value of $zp_{c}\sim
16\eta_{c}=2.086$. For bond percolation, data fitting shows that the
thresholds satisfy $zp_{c}-1\sim a_{1}z^{-x}$, with $x\sim 3/4$, which is,
however, inconsistent with a recent theoretical suggestion FreiPerkins2016 ;
HofstadSakai2005 . In contrast, the value of exponent $x$ here indicates that
the relation of $x=(d-1)/d$, valid in two and three dimensions, can be
generalized to four dimensions.

The remainder of the paper is organized as follows. The simulation results and
discussions are given in Sec. II, including some simulation details (II.1),
results for site percolation (II.2) and bond percolation (II.3). In Sec. III,
we present our conclusions.

##  II Simulation results and discussions

###  II.1 Simulation details and basic theory

We use the notation sc(4)-$a,b,...$ to indicate a four-dimensional simple
hypercubic lattice with the $a$th nearest neighbors, $b$th nearest neighbors,
etc. In the Monte-Carlo simulations, by using a single-cluster growth
algorithm (see Refs. LorenzZiff1998 ; XunZiff2020 ; XunZiff2020b  for more
details), a site on the lattice of size $L\times L\times L\times L$ with
$L=128$ under periodic boundary conditions is chosen as the seed, and an
individual cluster begins to grow at that seeded site. We grow many samples of
individual clusters for each lattice, which are $5\times 10^{8}$ for sc(4)-3,
sc(4)-1,3, sc(4)-2,3, sc(4)-1,2, sc(4)-1,2,3, and sc(4)-1,2,3,4 lattices, and
$10^{8}$ for other five lattices. Clusters with different sizes are
distributed in bins of range of $(2^{n},2^{n+1}-1)$ for $n=0,1,2,\ldots$. An
upper size of the cluster, which is called the upper size cutoff, needs to be
set to halt the growth of clusters larger than that size, to avoid wrapping
around the boundaries. For clusters still growing when they reach an upper
size cutoff, we count them in the last bin. Here for site percolation, we set
the upper size cutoff to be $2^{15}$ occupied sites for sc(4)-1,2 lattice,
sc(4)-1,2,3 lattice; $2^{13}$ for the sc(4)-1,…,9 lattice; and $2^{14}$ for
other lattices. For bond percolation, we set the upper size cutoff to be
$2^{16}$ occupied (wetted) sites for all lattices.

If one defines $n_{s}(p)$ as the number of clusters (per site) containing $s$
occupied sites, as a function of the site or bond occupation probability $p$,
then in the scaling limit, in which $s$ is large and $(p-p_{c})$ is small such
that $(p-p_{c})s^{\sigma}$ is constant, $n_{s}(p)$ behaves as

$n_{s}(p)\sim A_{0}s^{-\tau}f[B_{0}(p-p_{c})s^{\sigma}],$ (6)

where $\tau$, $\sigma$, and $f(x)$ are universal, while $A_{0}$ and $B_{0}$
are lattice-dependent metric factors. When $p=p_{c}$ and for finite $s$, there
are corrections to Eq. (6)

$n_{s}(p_{c})\sim A_{0}s^{-\tau}(1+C_{0}s^{-\Omega}+\dots),$ (7)

where $\Omega$ is another universal exponent. Then the probability that a
point belongs to a cluster of size greater than or equal to $s$ is given by
$P_{\geq s}=\sum_{s^{\prime}=s}^{\infty}s^{\prime}n_{s^{\prime}}$, and it
follows XunZiff2020 ; XunZiff2020b

$s^{\tau-2}P_{\geq s}\sim
A_{1}[1+B_{1}(p-p_{c})s^{\sigma}+C_{1}s^{-\Omega}+\dots],$ (8)

where $A_{1}$, $B_{1}$ and $C_{1}$ are non-universal constants. Equation (8)
provides us two methods to determine the percolation threshold $p_{c}$, and we
will show them in detail in the next subsection, combined with our simulation
results.

Due to the universality of $\tau$, $\Omega$ and $\sigma$, these exponents
depend only on the system dimension, not the type of the lattice. We choose
$\tau=2.3135(5)$XunZiff2020 , $\Omega=0.40(3)$XunZiff2020 , and
$\sigma=0.4742$Gracey2015 , which are relatively accurate and acceptable, in
the numerical simulations. Then the number of clusters greater than or equal
to size $s$ could be found based on the data from our simulation, and the
quantity $s^{\tau-2}P_{\geq s}$ could be easily calculated.

###  II.2 Results of site percolation

From the behavior of Eq. (8), we can determine if we are above, near, or below
the percolation threshold. For large $s$ where the finite-size effect term
$s^{-\Omega}$ can be ignored, Eq. (8) becomes

$s^{\tau-2}P_{\geq s}\sim A_{1}[1+B_{1}(p-p_{c})].$ (9)

This implies that $s^{\tau-2}P_{\geq s}$ will convergence to a constant value
at $p_{c}$, while it deviates linearly from that constant value when $p$ is
away from $p_{c}$. Figure 1 shows the relation of $s^{\tau-2}P_{\geq s}$
versus $s^{\sigma}$ for the site percolation of the sc(4)-1,…,9 lattice under
probabilities $p=0.004828$, $0.004829$, $0.004830$, $0.004831$, $0.004832$,
and $0.004833$. For small clusters, $s^{\tau-2}P_{\geq s}$ shows a steep rise
due to the finite-size effect, while for large clusters, $s^{\tau-2}P_{\geq
s}$ shows a linear region. As $p$ tends to $p_{c}$, the linear part of
$s^{\tau-2}P_{\geq s}$ become more nearly horizontal. Based upon these
properties of the linear portions in Fig. 1, the central value of $p_{c}$ can
be determined by plotting

$\frac{\mathrm{d}(s^{\tau-2}P_{\geq s})}{\mathrm{d}(s^{\sigma})}\sim
B_{1}(p-p_{c}),$ (10)

as shown in the inset of Fig. 1, $p_{c}=0.0048301$ can be deduced from the $p$
intercept of the plot of the above derivative versus $p$.

If $p$ is very close to $p_{c}$, Eq. (8) can be reduced to

$s^{\tau-2}P_{\geq s}\sim A_{1}[1+C_{1}s^{-\Omega}],$ (11)

implying a linear relationship between $s^{\tau-2}P_{\geq s}$ and
$s^{-\Omega}$ for large $s$. In Fig. 2, we show the plot of $s^{\tau-2}P_{\geq
s}$ versus $s^{-\Omega}$ for the site percolation of the sc(4)-1,…,9 lattice
under probabilities $p=0.004828$, $0.004829$, $0.004830$, $0.004831$,
$0.004832$, and $0.004833$. It can be seen clearly that when $p$ is away from
$p_{c}$, the curves show an obvious deviation from linearity for large $s$,
while better linear behavior emerges if $p$ is very close to $p_{c}$. Then we
can conclude the range of the site percolation threshold of
$0.004830<p_{c}<0.004831$ here.

Figure 1: Plot of $s^{\tau-2}P_{\geq s}$ vs $s^{\sigma}$ with $\tau=2.3135$
and $\sigma=0.4742$ for the site percolation of the sc(4)-1,…,9 lattice under
different values of $p$. The inset indicates the slope of the linear portions
of the curves shown in the main figure as a function of $p$, and the central
value of $p_{c}=0.0048301$ can be calculated from the $p$ intercept. Figure 2:
Plot of $s^{\tau-2}P_{\geq s}$ vs $s^{-\Omega}$ with $\tau=2.3135$ and
$\Omega=0.40$ for the site percolation of the sc(4)-1,…,9 lattice under
different values of $p$.

Based upon the two methods indicated in Figs. 1 and 2, as well as the errors
for the values of $\tau=2.3135(5)$ and $\Omega=0.40(3)$, finally, the site
percolation threshold of the sc(4)-1,…,9 lattice can be deduced to be
$p_{c}=0.0048301(9)$, where the number in parentheses represents the estimated
error in the last digit. The figures for the other 7 lattices we simulated are
shown in the Supplementary Material XunZiff2020supplementary  in Figs. 1-14,
and the corresponding site percolation thresholds are summarized in Table 1.

In Table 1, we also show previous site percolation thresholds for sc(4)-1,2
lattice and sc(4)-1,2,3 lattice provided by Ref.KotwicaGronekMalarz19 , as can
be seen in the last column. For these two lattices, we get more precise site
percolation thresholds than previous results. We also get new values of
thresholds for other 6 lattices.

Table 1: Site percolation thresholds on four-dimensional simple hypercubic lattice with extended neighborhoods up to the 9th nearest neighbors. lattice | $z$ |  $p_{c}$(present) |  $p_{c}$(previous)
---|---|---|---
sc(4)-1,2 | 32 | 0.0617731(19) | 0.06190(23)KotwicaGronekMalarz19
sc(4)-1,2,3 | 64 | 0.0319407(13) | 0.03190(23)KotwicaGronekMalarz19
sc(4)-1,2,3,4 | 88 | 0.0231538(12) | —
sc(4)-1,…,5 | 136 | 0.0147918(12) | —
sc(4)-1,…,6 | 232 | 0.0088400(10) | —
sc(4)-1,…,7 | 296 | 0.0070006(6) | —
sc(4)-1,…,8 | 320 | 0.0064681(9) | —
sc(4)-1,…,9 | 424 | 0.0048301(9) | —

With regard to the asymptotic behavior between percolation thresholds $p_{c}$
and coordination number $z$ for the site percolation on lattices with compact
nearest neighborhoods, in four dimensions where $\eta_{c}$ for hyperspheres
equals $0.1304(5)$ TorquatoJiao2012 , one should expect from Eq. (1) that

$p_{c}=\frac{2.086(8)}{z}.$ (12)

The extended-range nearest-neighbor regions on a lattice are not exactly
spherical in shape, of course. A limit of a more non-spherical surface is a
system aligned hypercubes, where $\eta_{c}=0.1201(6)$ TorquatoJiao2012 ,
implying $c=zp_{c}=16\eta_{c}=1.922(10)$. Indeed, for some small values of
$z$, the neighborhood tends to a hypercube. However, for large $z$ where the
neighborhood becomes more spherical, one would expect that $zp_{c}$ should
approach the value 2.086.

By plotting the relation of $z$ versus $1/p_{c}$ and $z$ versus
$-1/\ln(1-p_{c})$, as shown in Figs. 3 and 4, respectively, we can see that
data fitting of both plots lead to the asymptotic value of $2.065$, which
agrees quite well with the theoretically predicted value for hyperspheres
above.

In addition, in two and three dimensions, there was a contradiction between
the choice of $z$ versus $1/p_{c}$ and $z$ versus $-1/\ln(1-p_{c})$, for we
found XunHaoZiff2021b  that, in two dimensions with largest $z=84$ and three
dimensions with largest $z=146$, $z$ versus $1/p_{c}$ is better than $z$
versus $-1/ln(1-p_{c})$. Meanwhile, we also found that the gap between the two
methods is smaller as $z$ increases. Here in four dimensions with largest
$z=424$, both $z$ versus $1/p_{c}$ and $z$ versus $-1/\ln(1-p_{c})$ nearly
give the same data fitting values, as one would expects since the two formulas
are asymptotically equal. These comparisons show that $z$ versus $1/p_{c}$
holds for almost $z$, while $z$ versus $-1/ln(1-p_{c})$ holds when $z$ is
large enough. One should pay more attention when use the two formulas in
further.

Figure 3: Plot of $z$ vs $1/p_{c}$ for the four-dimensional simple hypercubic
lattice with compact nearest neighborhoods shown in Table 1. Data fitting
gives the slope of $c=2.065$, compared with the prediction of $zp_{c}\sim
16\eta_{c}=2.086$, and the intercept of $-b=-1.290$. Figure 4: Plot of $z$
versus $-1/\ln(1-p_{c})$ for the four-dimensional simple hypercubic lattice
with compact nearest neighborhoods shown in Table 1. Data fitting gives the
slope of $2.065$, compared with the prediction of $zp_{c}\sim
16\eta_{c}=2.086$, and the intercept of $-0.250$.

###  II.3 Results of bond percolation

For the bond percolation simulated in this paper, Figs. 5 and 6 show the plots
of $s^{\tau-2}P_{\geq s}$ versus $s^{\sigma}$ and $s^{-\Omega}$, respectively,
for the sc(4)-1,…,9 lattice under probabilities $p=0.0024105$, $0.0024110$,
$0.0024115$, $0.0024120$, $0.0024125$, and $0.0024130$. Similar to the
discussion above, we conclude the bond percolation threshold of the
sc(4)-1,…,9 lattice here to be $p_{c}=0.0024117(7)$. The figures for the other
9 lattices we simulated are shown in the Supplementary Material
XunZiff2020supplementary  in Figs. 15-32, and the corresponding bond
percolation thresholds are summarized in Table 2.

Figure 5: Plot of $s^{\tau-2}P_{\geq s}$ versus $s^{\sigma}$ with $\tau=2.3135$ and $\sigma=0.4742$ for the bond percolation of the sc(4)-1,…,9 lattice under different values of $p$. The inset indicates the slope of the linear portions of the curves shown in the main figure as a function of $p$, and the central value of $p_{c}=0.0024117$ can be calculated from the $p$ intercept. Figure 6: Plot of $s^{\tau-2}P_{\geq s}$ versus $s^{-\Omega}$ with $\tau=2.3135$ and $\Omega=0.40$ for the bond percolation of the sc(4)-1,…,9 lattice under different values of $p$. Table 2: Bond percolation thresholds for the four-dimensional simple hypercubic lattice with extended neighborhoods up to the 9th nearest neighbors. lattice | $z$ | $p_{c}$ | $zp_{c}$
---|---|---|---
sc(4)-3 | 32 | 0.0338047(27) | 1.08175
sc(4)-1,3 | 40 | 0.0271892(22) | 1.08757
sc(4)-2,3 | 56 | 0.0194075(15) | 1.08682
sc(4)-1,2,3 | 64 | 0.0171036(11) | 1.09463
sc(4)-1,2,3,4 | 88 | 0.0122088(8) | 1.07437
sc(4)-1,…,5 | 136 | 0.0077389(9) | 1.05249
sc(4)-1,…,6 | 232 | 0.0044656(11) | 1.03602
sc(4)-1,…,7 | 296 | 0.0034812(7) | 1.03044
sc(4)-1,…,8 | 320 | 0.0032143(8) | 1.02858
sc(4)-1,…,9 | 424 | 0.0024117(7) | 1.02256

Table 2 also shows the values of $zp_{c}$ for each lattice. With the increase
of $z$, the value of $zp_{c}$ decreases and tends to the asymptotic value of
1. In Fig. 7, we show the relation of $zp_{c}$ versus $z^{-x}$, and further
study indicates that better linear fitting is found when $x\approx 3/4$. The
value of the exponent $x$ here is inconsistent with the theoretical suggestion
of $zp_{c}-1\simeq a_{1}/z$ (with logarithm corrections in $d=4$) in Refs.
FreiPerkins2016 ; HofstadSakai2005 . Here we believe our value of $x=3/4$
seems to make more sense. Firstly, In Ref. FreiPerkins2016 , Frei and Perkins
did a mapping to a random walk problem. The critical dimension for random
walks is four, so they get logarithmic corrections. However, in percolation,
the critical dimension is six. Secondly, the relation of $x=(d-1)/d$ (a
surface area effect) in two and three dimensions has already been proved by
extensive numerical simulations OuyangDengBlote2018 ; DengOuyangBlote2019 ;
XuWangHuDeng2021 ; XunHaoZiff2021b , so it is natural to expect that this
relation holds in higher but below the critical dimensions.

Figure 7: Plot of $zp_{c}$ vs $z^{-3/4}$ for the four-dimensional simple
hypercubic lattice with compact nearest neighborhoods shown in Table 2. Data
fitting gives the slope of $2.134$ and the intercept of $1.000$.

##  III Conclusion

To summarize, in this paper, in order to further explore the correlation
between percolation thresholds $p_{c}$ and coordination number $z$ in higher
dimensions, we have carried out extensive Monte Carlo simulations for the site
and bond percolation on four-dimensional simple hypercubic lattice with
extended neighborhoods up to 9th nearest neighbors. By employing an effective
single-cluster growth method, we found precise estimates of the percolation
thresholds for 18 systems.

For site percolation, the asymptotic behavior was investigated by plotting
both $z$ versus $1/p_{c}$ and $z$ versus $-1/\ln(1-p_{c})$, and both figures
tend to the theoretical value of $zp_{c}\sim 16\eta_{c}\approx 2.086$.

For bond percolation, the asymptotic value of $zp_{c}$ tends to the Bethe-
lattice behavior of unity with the increase of the coordination number $z$. We
accounted the finite-$z$ corrections by analyzing $zp_{c}-1\simeq
a_{1}z^{-x}$, and data fitting shows that $x=3/4$ gives the best linearity.
The value of the exponent $x$ indicated that the relation of $x=(d-1)/d$,
which holds in two and three dimensions, is still satisfied in four
dimensions. As the critical dimension for percolation is six, it will be
interesting to further check the validity of this relation in five dimensions.

##  IV Acknowledgments

The authors are grateful to the Advanced Analysis and Computation Center of
CUMT for the award of CPU hours to accomplish this work. This work is
supported by “the Fundamental Research Funds for the Central Universities”
under Grant No. 2020ZDPYMS31.

## References

  * [1] N. W. Dalton, C. Domb, and M. F. Sykes.  Dependence of critical concentration of a dilute ferromagnet on the range of interaction.  Proc. Phys. Soc., 83(3):496–498, 1964.
  * [2] C. Domb and N. W. Dalton.  Crystal statistics with long-range forces: I. The equivalent neighbour model.  Proc. Phys. Soc., 89(4):859–871, 1966.
  * [3] C. Domb.  A note on the series expansion method for clustering problems.  Biometrika, 59(1):209–211, 1972.
  * [4] Jon M. Kleinberg.  Navigation in a small world.  Nature, 406:845, 2000.
  * [5] L. M. Sander, C. P. Warren, and I. M. Sokolov.  Epidemics, disorder, and percolation.  Physica A, 325(1):1 – 8, 2003.
  * [6] Robert M. Ziff.  Percolation and the pandemic.  Physica A, 568:125723, 2021.
  * [7] Zbigniew Koza, Grzegorz Kondrat, and Karol Suszczyński.  Percolation of overlapping squares or cubes on a lattice.  J. Stat. Mech.: Th. Exp., 2014:P11005, 2014.
  * [8] Zbigniew Koza and Jakub Poła.  From discrete to continuous percolation in dimensions 3 to 7.  J. Stat. Mech.: Th. Exp., 10:103206, 2016.
  * [9] Mark Gouker and Fereydoon Family.  Evidence for classical critical behavior in long-range site percolation.  Phys. Rev. B, 28:1449–1452, 1983.
  * [10] G. R. Jerauld, L. E. Scriven, and H. T. Davis.  Percolation and conduction on the 3d Voronoi and regular networks: a second case study in topological disorder.  J. Phys. C: Solid State, 17(19):3429–3439, 1984.
  * [11] T. R. Gawron and Marek Cieplak.  Site percolation thresholds of fcc lattice.  Acta Phys. Pol. A, 80(3):461–464, 1991.
  * [12] C. d’Iribarne, G. Rasigni, and M. Rasigni.  Determination of site percolation transitions for 2d mosaics by means of the minimal spanning tree approach.  Phys. Lett. A, 209(1):95–98, 1995.
  * [13] C. d’Iribarne, M. Rasigni, and G. Rasigni.  Minimal spanning tree and percolation on mosaics: graph theory and percolation.  J. Phys. A: Math. Gen., 32(14):2611–2622, 1999.
  * [14] C. d’Iribarne, M. Rasigni, and G. Rasigni.  From lattice long-range percolation to the continuum one.  Phys. Lett. A, 263(1):65–69, 1999.
  * [15] Krzysztof Malarz and Serge Galam.  Square-lattice site percolation at increasing ranges of neighbor bonds.  Phys. Rev. E, 71:016125, 2005.
  * [16] M. Majewski and K. Malarz.  Square lattice site percolation thresholds for complex neighbourhoods.  Acta Phys. Pol. B, 38:2191, 2007.
  * [17] Łukasz Kurzawski and Krzysztof Malarz.  Simple cubic random-site percolation thresholds for complex neighbourhoods.  Rep. Math. Phys., 70(2):163 – 169, 2012.
  * [18] Krzysztof Malarz.  Simple cubic random-site percolation thresholds for neighborhoods containing fourth-nearest neighbors.  Phys. Rev. E, 91:043301, 2015.
  * [19] M. Kotwica, P. Gronek, and K. Malarz.  Efficient space virtualization for the Hoshen-Kopelman algorithm.  Int. J. Mod. Phys. C, 30(8):1950055, 2019.
  * [20] Krzysztof Malarz.  Site percolation thresholds on triangular lattice with complex neighborhoods.  Chaos, 30:123123, 2020.
  * [21] Zhipeng Xun, Dapeng Hao, and Robert M. Ziff.  Site percolation on square and simple cubic lattices with extended neighborhoods and their continuum limit.  Phys. Rev. E, 103:022126, 2021.
  * [22] Yunqing Ouyang, Youjin Deng, and Henk W. J. Blöte.  Equivalent-neighbor percolation models in two dimensions: Crossover between mean-field and short-range behavior.  Phys. Rev. E, 98:062101, 2018.
  * [23] Youjin Deng, Yunqing Ouyang, and Henk W. J. Blöte.  Medium-range percolation in two dimensions.  J. Phys.: Conf. Ser., 1163:012001, 2019.
  * [24] Zhipeng Xun and Robert M. Ziff.  Precise bond percolation thresholds on several four-dimensional lattices.  Phys. Rev. Research, 2:013067, 2020.
  * [25] Zhipeng Xun and Robert M. Ziff.  Bond percolation on simple cubic lattices with extended neighborhoods.  Phys. Rev. E, 102:012102, 2020.
  * [26] Wenhui Xu, Junfeng Wang, Hao Hu, and Youjin Deng.  Critical polynomials in the nonplanar and continuum percolation models.  Phys. Rev. E, 103:022127, 2021.
  * [27] Serge Galam and Alain Mauger.  Universal formulas for percolation thresholds.  Phys. Rev. E, 53:2177–2181, 1996.
  * [28] Steven C. van der Marck.  Calculation of percolation thresholds in high dimensions for fcc, bcc and diamond lattices.  Int. J. Mod. Phys. C, 9(4):529–540, 1998.
  * [29] Spencer Frei and Edwin Perkins.  A lower bound for $p_{c}$ in range-$r$ bond percolation in two and three dimensions.  Electron. J. Probab., 21:56, 2016.
  * [30] R. Hofstad and A. Sakai.  Critical points for spread-out self-avoiding walk, percolation and the contact process above the upper critical dimensions.  Probab. Theory Relat. Fields, 132:438–470, 2005.
  * [31] Zhipeng Xun, Dapeng Hao, and Robert M. Ziff.  Universal asymptotic behavior of the percolation on regular lattices with compact neighborhoods in 2 and 3 dimensions.  preprint, 2021.
  * [32] Christian D. Lorenz and Robert M. Ziff.  Precise determination of the bond percolation thresholds and finite-size scaling corrections for the sc, fcc, and bcc lattices.  Phys. Rev. E, 57:230–236, 1998.
  * [33] J. A. Gracey.  Four loop renormalization of ${\phi}^{3}$ theory in six dimensions.  Phys. Rev. D, 92:025012, 2015.
  * [34] See Supplemental Material at http://xxx for the simulation results.
  * [35] S. Torquato and Y. Jiao.  Effect of dimensionality on the continuum percolation of overlapping hyperspheres and hypercubes. II. Simulation results and analyses.  J. Chem. Phys., 137:074106, 2012.


# A sea surface-based drag model for Large Eddy Simulation of wind-wave
interaction

###### Abstract

Monin-Obukhov similarity theory (MOST) is a well tested approach for
specifying the fluxes when the roughness surfaces are homogeneous. For flow
over waves (inhomogeneous surfaces), phase-averaged roughness length scales
are often prescribed through models based on the wave characteristics and the
wind speed. However, such approaches lack generalizability over different wave
ages and steepnesses due to the reliance on model coefficients tuned to
specific datasets. In this paper, a sea surface-based hydrodynamic drag model
applicable to moving surfaces is developed to model the pressure-based surface
drag felt by the wind due to the waves. The model is based on the surface
gradient approach of Anderson and Meneveau (2010) applicable to stationary
obstacles and extended here to the wind-wave problem. The wave drag model
proposed specifies the hydrodynamic force based on the incoming momentum flux,
wave phase speed, and the surface frontal area. The drag coefficient
associated with the wind-wave momentum exchange is determined based on the
wave steepness. The wave drag model is used to simulate turbulent airflow
above a monochromatic wave train with different wave ages and wave
steepnesses. The mean velocity profiles and model form stresses are validated
with available lab-scale experimental data and show good agreement across a
wide range of wave steepnesses and wave ages. The drag force is correlated
with the wave surface gradient and out-of-phase with the wave height
distribution by a factor $\pi/2$ for the sinusoidal wave-train considered.
These results demonstrate that the current approach is sufficiently general
over a wide parameter space compared to wave phase-averaged models with a
minimal increase in computational cost.

## Significance Statement

Understanding the physics of wind-waves plays an important role in the context
of numerous geophysical and engineering applications. A drag-based model is
developed that characterizes the effect of the sea surface waves on the wind
above. The model is validated with existing experimental datasets and is shown
to be effective in predicting the average wind velocity and stress over waves
with varied steepnesses and phase speeds. The ease of implementation and low
computational cost of the model make it useful for studying turbulent
atmospheric-scale flows over the sea surface important in offshore wind energy
research as well as for modeling air-sea fluxes of momentum, heat, and mass.

##  1 Introduction

The physics of wind-wave interactions involve an interplay between atmospheric
turbulence and ocean waves and play an important role in the context of
numerous geophysical and engineering applications. Wind-wave interactions have
been a continuous topic of research spanning numerous experimental (Shemdin
and Hsu, 1967; Grare et al., 2013; Buckley and Veron, 2017; Buckley et al.,
2020), theoretical (Miles, 1993; Janssen, 1991; Belcher, 1999; Kudryavtsev and
Makin, 2004), and computational (Sullivan et al., 2000, 2008, 2014; Yang and
Shen, 2010; Yang et al., 2013; Cao et al., 2020; Cao and Shen, 2021) studies
devoted to quantifying the momentum exchange between the waves and the airflow
above. Many physical processes rely on accurately quantifying the heat, mass,
and momentum exchanges that take place at the air-sea interface and are
important for numerous practical applications such as weather and climate
modeling, safe steering of naval vessels, and offshore wind energy harvesting.

Measurements of air pressure in the laboratory and field have been carried out
to quantify the momentum exchange at the air-sea interface (Shemdin and Hsu,
1967; Snyder et al., 1981; Banner, 1990; Donelan et al., 2006; Peirson and
Garcia, 2008; Grare et al., 2013). The results and findings from these
experiments have been summarized by Grare et al. (2013) where they analyzed
the growth rate and associated form stress due to energy transfer between the
wind and the waves for different experimental data. They found that the
normalized wind  energy input rate values showed good collapse as a function
of wave steepness rather than the wave age that is typically used in the
literature for the datasets analyzed. Detailed experimental measurements of
the airflow structure and wind stress over the wave surface within the
airflow’s viscous sublayer were carried out by Buckley and Veron (2017) and
Buckley et al. (2020). They found that the form stress dominated the total
air-water momentum flux for high wave slopes while the viscous drag was
important for low wave slopes. The high resolution and quality of the wind and
wave profile measurements provide a useful database for comparison with models
and simulations.

The component of the total stress $\bm{\tau}$ imparted by the wind to the wave
field can be written as a sum of the pressure and tangential viscous stresses.
Considering a sloped surface ${\eta}(x,t)$ varying as a function of horizontal
coordinate $x$ and time $t$ representing the wave, the horizontal component of
the stress can be written as

$\tau_{xz}=\left\langle
p_{s}\eta_{x}\right\rangle+\left\langle\frac{\tau_{visc,xz}}{\sqrt{1+\eta_{x}^{2}}}\right\rangle,$
(1)

where $z$ is the vertical direction, $p_{s}$ is the surface pressure,
$\tau_{visc,xz}=\mu\partial u/\partial z$ is the local viscous surface
tangential shear stress, $\eta_{x}=\partial{\eta}/\partial x$ is the local
interface slope and the angular brackets denote averaging over one wavelength.
A more detailed discussion of the various stress components is provided in the
Appendix of Grare et al. (2013). Using a linear approximation $\eta_{x}<<1$,
the total  energy input rate can be related to the wave-coherent component of
the horizontal stress:

$S_{in}=dE/dt=\left\langle p_{s}\eta_{x}\right\rangle
c+\left\langle\tau_{visc,xz}u_{s}\right\rangle,$ (2)

where $u_{s}$ is the surface orbital velocity, $E=1/2\rho
g\langle\eta^{2}\rangle$ is the energy per unit area of sea surface averaged
over a wavelength, and $c$ is the wave phase speed.  When normalized by the
wave phase speed, the first term on the RHS of Equation 2 is defined as the
form stress $\tau_{form}$ and the second term is the wave-coherent tangential
stress $\tau_{visc,wc}$ (Longuet-Higgins, 1969). The total energy input is
then related to the total wave-coherent momentum flux as
$S_{in}=(\tau_{form}+\tau_{visc,wc})c$. Knowledge of the dependence of the
energy input rate on wave characteristics is important in modeling the wind-
sea momentum flux. Two common parameterizations for the  energy input rate
have been developed by Jeffreys (1925) and Miles (1957). The former takes the
form:

$S_{in}=\frac{1}{2}\rho_{a}s_{\lambda/2}(ak)^{2}c^{3}\left(\frac{U_{\lambda/2}}{c}-1\right)\left|\frac{U_{\lambda/2}}{c}-1\right|,$
(3)

where $\rho_{a}$ is the air density, $s_{\lambda/2}$ is the sheltering
coefficient evaluated at a height $z=\lambda/2$ above the water surface with
$\lambda$ being the wavelength, $k=2\pi/\lambda$ is the wavenumber, $a$ is the
amplitude of the wave, and $U_{\lambda/2}$ is the mean streamwise velocity at
$z=\lambda/2$. The latter parameterization by Miles (1957) is given by

$S_{in}=\beta\frac{\rho_{a}}{\rho_{w}}\left(\frac{u_{\ast}}{c}\right)^{2}\omega
E,$ (4)

where $E$ is the wave energy density, $\omega$ is the angular frequency,
$\rho_{w}$ is the water density, $u_{\ast}$ is the surface friction velocity,
and $\beta$ is the normalized wave growth coefficient. The angular frequency
and the wavenumber are related by the dispersion relation $\omega=\sqrt{gk}$
where $g$ is the gravitational acceleration.

Computational simulations provide an alternative approach to these theoretical
parameterizations. Direct Numerical Simulation (DNS) provides the highest
fidelity description of the problem. However, these simulations have been
restricted to low Reynolds number for air-flow over a wavy surface due to high
computational costs associated with resolving the viscous layer for high
Reynolds number (Sullivan et al., 2000; Lin et al., 2008; Yang and Shen, 2010;
Druzhinin et al., 2012). Large Eddy Simulation (LES) has been shown to yield
high-fidelity results for flow over complex terrain resolving the large- and
intermediate-scale turbulent motions and only requiring modeling of the
unresolved turbulence effects. The effect of waves in LES of the marine
atmospheric boundary layer (MABL) has been incorporated either using wave
phase-resolved simulations (Sullivan et al., 2008, 2014; Yang et al., 2014;
Sullivan et al., 2018; Hao et al., 2018) or wave phase-averaged models
(Charnock, 1955; Donelan, 1990; Toba et al., 1990; Drennan et al., 2005). The
former approach relies on constructing a terrain-following grid and having
high enough resolution to resolve the near wave dynamics. The approach has
been used successfully to simulate boundary layer flow with the bottom
boundary specified using either a sinusoidal wave train (Sullivan et al.,
2008) or the full 3D wave spectrum (Sullivan et al., 2014; Yang et al., 2014).
Yang et al. (2014) developed a dynamic coupling procedure to include the
effect of air flow pressure field on the wave motion as a kinematic boundary
condition. This approach is quite computationally expensive for high Reynolds
number atmospheric flows due to the large separation of scales and the high
grid resolution required near the surface to resolve the wave motion in the
vertical direction.

In high Reynolds number flow simulations, a wall model can instead be used to
specify the boundary fluxes (Piomelli and Balaras, 2002). The phase-averaged
approach models the total surface stress $\tau$ by specifying a surface
roughness for the waves and calculating the stress based on MOST (Moeng,
1984):

$\tau=\rho_{a}\left[\frac{\kappa}{\log{\left(z_{r}/z_{0}\right)}-\psi_{m}(z/L)}\right]^{2}u_{r}^{2},$
(5)

where $\kappa=0.4$ is the Von-Karman constant, $z_{0}$ is the surface
roughness, $u_{r}$ is the reference velocity at height $z_{r}$, and $\psi_{m}$
is an empirical function of the stability parameter, $\zeta=-z/L$, where $L$
is the Obukhov length, which accounts for stratification effects. In the
Monin-Obukhov similarity framework, the unresolved stresses, including the
viscous stresses, are absorbed into a single roughness length $z_{0}$. The sea
surface roughness elements are all in motion with waves of different
wavelengths propagating with different speeds according to the wave dispersion
relation. To account for this, numerous models have been proposed to
parameterize an equivalent effective roughness scale $z_{0}$. A widely used
model is one proposed by Charnock (1955) where the surface roughness $z_{0}$
is calculated as

$z_{0}=\frac{\alpha_{ch}u_{\ast}^{2}}{g},$ (6)

where $\alpha_{ch}=0.01-0.03$ is an empirical constant (Garratt, 1977). The
wide range of values for $\alpha_{ch}$ have been attributed to factors such as
the wind speed, wave age, wave steepness, and water depth (Donelan, 1990;
Smith et al., 1992; Taylor and Yelland, 2001; Fairall et al., 2003; Edson et
al., 2013; Jiménez and Dudhia, 2018). Alternative formulations have been
proposed based on the peak wave age $c_{p}/u_{\ast}$ (Donelan, 1990; Toba et
al., 1990; Smith et al., 1992):

$z_{0}=A(u_{\ast}/c_{p})^{B},$ (7)

where $A$ and $B$ are empirical parameters that vary across datasets (Deskos
et al., 2021). Taylor and Yelland (2001) proposed a parameterization based on
the characteristic wave steepness $H_{s}/\lambda_{p}$:

$z_{0}=1200H_{s}(H_{s}/\lambda_{p})^{3.4},$ (8)

where $H_{s}$ is the characteristic wave height, $\lambda_{p}$ is the peak
wavelength, and the factors $1200$ and exponent $3.4$ are determined using
data consisting of waves with larger fetch (“old waves”). Due to the modeling
uncertainties for the effective roughness in phase-averaged formulations,
these formulations may not be sufficiently general. There is a need to develop
a model that would be applicable to a range of wave steepnesses and wave ages
yet still retaining the computational efficiency of these wave phase-averaged
parameterizations. It is important to note that such a model would not resolve
the detailed near-surface flow around the waves but rather reproduce the
forces imposed by and flow above the waves in the atmospheric boundary layer.

In the atmospheric boundary layer (ABL), the effect of forest canopies and
vegetation has been represented using drag-based models termed canopy stress
models (Shaw and Schumann, 1992; Su et al., 1998; Arthur et al., 2019). The
hydrodynamic drag force due to obstacles in the flow is modeled as a function
of the obstacle frontal area and the incoming momentum flux. Brown et al.
(2001) applied the canopy stress models to flow over stationary sinusoidal
ridges. Anderson and Meneveau (2010) developed a surface-gradient drag (SGD)
model where the obstacle roughness was vertically unresolved but resolved
horizontally. Arthur et al. (2019) implemented the canopy stress model
framework in the Weather Research and Forecasting (WRF) model and applied it
to both resolved and unresolved roughness. The canopy stress approach provides
a simple framework in which effects of obstacles and topographies could be
included without additional computational complexity and without requiring
empirical correlations for some effective roughness.

In the current work, a wave drag formulation to represent the wave stress is
developed for the marine atmospheric boundary layer applicable to the wind-
wave problem. The ease of implementation of this model along with the lower
associated computational cost makes it an attractive alternative to more
expensive wave phase-resolved terrain following approaches. The rest of the
paper is organized as follows. The governing equations used for the LES are
described in §2. The wave drag model formulation is described in §3. The model
validation and results are presented in §4. Conclusions are drawn in §5.

##  2 Computational Framework

LES calculations were performed using NGA, which is a structured, finite
difference, low Mach number flow solver (Desjardins et al., 2008; MacArt and
Mueller, 2016). The wind velocity field is described using the filtered
Navier-Stokes equations in the incompressible limit:

$\nabla\cdot\mbox{\boldmath$\widetilde{u}$}{}=0,$ (9)
$\displaystyle\frac{\partial\mbox{\boldmath$\widetilde{u}$}{}}{\partial
t}+\mbox{\boldmath$\widetilde{u}$}{}\cdot\nabla\mbox{\boldmath$\widetilde{u}$}{}=$
$\displaystyle-\frac{1}{\rho}\nabla\mbox{$\widetilde{p}$}{}-\nabla\cdot{\bm{\widetilde{\tau}}}+{\bm{f}}_{d}.$
(10)

The equations are discretized on a Cartesian grid $(x,y,z)$, where $x$ and $y$
are the streamwise and spanwise coordinates, respectively, and $z$ is the
vertical coordinate. In Equations (9) and (10),
$\mbox{\boldmath$\widetilde{u}$}{}=(\mbox{$\widetilde{u}$}{},\mbox{$\widetilde{v}$}{},\mbox{$\widetilde{w}$}{})$
is the velocity vector with the tilde denoting variables filtered on the LES
grid;  ${\bm{f}}_{d}$ is the force per unit mass applied in the streamwise and
spanwise directions, representing the effects of the waves; and
$\widetilde{\tau}_{ij}=-2\nu\widetilde{S}_{ij}+\mbox{
$\widetilde{\tau}$}{}^{d}_{ij}$ is the total deviatoric stress (divided by
density), where $\nu$ is the molecular viscosity, $\widetilde{S}_{ij}$ is the
resolved strainrate tensor, and $\mbox{ $\widetilde{\tau}$}{}^{d}_{ij}$ is the
subfilter stress (SFS) tensor. The SFS tensor is modeled using a Lilly-
Smagorinsky type subfilter viscosity model $\mbox{
$\widetilde{\tau}$}{}_{ij}^{d}=-2\nu_{T}\widetilde{S}_{ij}$, where the
subfilter viscosity is computed using the Anisotropic Minimum Dissipation
(AMD) model (Rozema et al., 2015; Abkar and Moin, 2017). The subfilter
viscosity is given by

$\nu_{T}=\frac{-(\hat{\partial}_{k}\widetilde{u}_{i})(\hat{\partial}_{k}\widetilde{u}_{j})\widetilde{S}_{ij}}{(\partial_{l}\widetilde{u}_{m})(\partial_{l}\widetilde{u}_{m})},$
(11)

where $\hat{\partial}_{\gamma}=\sqrt{C}\Delta_{\gamma}\partial_{\gamma}$ (no
summation implied, $\gamma=1,2,3$) is the scaled gradient operator  (so
dimensionless), $\Delta_{\gamma}$ is the LES resolution in the respective
coordinate direction, and $C=1/3$ is the modified Poincare constant (Abkar and
Moin, 2017). The AMD subfilter model has been tested for large-scale
simulations of turbulent boundary layers and shown to be accurate and have a
lower computational cost compared to Lagrangian models (Gadde et al., 2021).

In LES of high Reynolds number flows, the viscous boundary layer at the bottom
boundary at the waves is unresolved. A wall-layer model is used to impose the
correct surface stress for the bottom boundary. This surface stress could be
formulated in terms of the total roughness scale $z_{0}$, which is generally
written as a sum of the smooth surface roughness $z_{0,s}$ and a wave-based
roughness scale $z_{0,w}$, giving a total roughness $z_{0}=z_{0,s}+z_{0,w}$
(Yang et al., 2013). The wave-based roughness scale is generally parameterized
using phase-averaged roughness models, as discussed in the Introduction.
However, in this work, the wave drag model will replace the wave-based
roughness scale $z_{0,w}$ to provide better generalizability and accuracy and
is discussed in the next section. For low wind speeds over smooth ocean
surfaces, the smooth roughness scale Reynolds number, defined based on the
friction velocity $u_{\ast}$ and the roughness scale $z_{0,s}$, approaches a
constant value $Re_{z_{0,s}}=0.11$ (Fairall et al., 1996). This microscale
smooth surface roughness length scale $z_{0,s}$ can then be calculated as

$z_{0,s}=0.11\frac{\nu_{a}}{u_{\ast}}.$ (12)

To apply the resulting wall stress locally, the stress is enforced using an
enhanced subfilter viscosity at the wall along with a Dirichlet boundary
condition for the velocity components. The subfilter viscosity is set to
enforce the correct wall stress given by (Moeng, 1984)

$\mbox{
$\widetilde{\tau}$}{}_{wall}=\left[\frac{\kappa\mbox{$\widehat{\widetilde{U}}$}{}_{avg}}{\log\left(\frac{\Delta_{z}/2-\widetilde{\eta}}{z_{0,s}}\right)}\right]^{2},$
(13)

where
$\mbox{$\widehat{\widetilde{U}}$}{}_{avg}=\sqrt{\mbox{$\widehat{\widetilde{u}}$}{}_{r}^{2}(x,y,t)+\mbox{$\widehat{\widetilde{v}}$}{}_{r}^{2}(x,y,t)}$
is the magnitude of the tangential wind velocity relative to the wave surface,
$\Delta_{z}/2$ is the first off-wall grid point, and $\widetilde{\eta}(x,y,t)$
is the sea surface elevation filtered at the LES resolution $\Delta$. The
relative surface velocities
$\mbox{$\widehat{\widetilde{u}}$}{}_{r},\mbox{$\widehat{\widetilde{v}}$}{}_{r}$
are calculated using the velocity at the first grid point and the wave surface
orbital velocity (Yang et al., 2013; Sullivan et al., 2014; Yang et al.,
2014):

$\mbox{$\widehat{\widetilde{u}}$}{}_{r}=\mbox{$\widehat{\widetilde{u}}$}{}(x,y,z_{1},t)-\mbox{$\widehat{\widetilde{u}}$}{}_{s}(x,y,t),\qquad\mbox{$\widehat{\widetilde{v}}$}{}_{r}=\mbox{$\widehat{\widetilde{v}}$}{}(x,y,z_{1},t)-\mbox{$\widehat{\widetilde{v}}$}{}_{s}(x,y,t),$
(14)

where
$\mbox{$\widehat{\widetilde{u}}$}{}_{s},\mbox{$\widehat{\widetilde{v}}$}{}_{s}$
are the wave orbital velocities. In the wall stress calculation, the
velocities are test-filtered using a box-filter with a filter-width equal to
twice the grid scale (LES resolution $\Delta$) and are denoted by
$(\widehat{\widetilde{...}})$. This additional filtering reduces velocity
fluctuations significantly and improves the local applicability of the
logarithmic law (Bou-Zeid et al., 2005). The subfilter viscosity is then set
such that

$\nu_{T}={\mbox{
$\widetilde{\tau}$}{}_{wall}}\frac{\Delta_{z}/2}{\mbox{$\widehat{\widetilde{U}}$}{}_{avg}}-\nu.$
(15)

A similar enhanced subfilter viscosity approach was shown to improve
predictions of statistics in the near-wall region of turbulent channel flows
(Bae and Lozano-Durán, 2021).

##  3 Wave drag model

The proposed drag model for the waves is based on a canopy stress-type
approach used for simulating flows over vegetation or forest canopies. The
model is based on applying a hydrodynamic drag force proportional to the
incoming momentum flux onto a frontal surface area (Brown et al., 2001). This
offers the advantage of incorporating the wave characteristics (surface
topology and wave speed) into the formulation without relying on empirical
roughness parameterizations. Typically, for LES where the wave-air interface
is unresolved, the effect of the waves is solely modeled by the subfilter
stress via a roughness parameter. In this section, a model is developed to
represent the pressure drag due to the waves and applied locally as a force in
the LES domain.

Figure 1:  Sketch of the simulation setup and the wave drag model. The waves
are vertically unresolved and lie beneath the first grid point.

The surface gradient drag model of Anderson and Meneveau (2010) is extended to
account for moving boundaries in order to model flow over waves. Figure 1
shows a sketch of the model setup. The waves are vertically unresolved and lie
below the first LES grid cell similar to Anderson and Meneveau (2010, 2011).
Consider the control volume encompassing the wave surface and the first grid
point, moving with a characteristic phase velocity $\bm{c}$. The incoming
momentum flux associated with the flow over a surface with an area $A$ is
given by
$\int\rho\mbox{\boldmath$\widetilde{u}$}{}(\mbox{\boldmath$\widetilde{u}$}{}_{r,c}\cdot{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{d}\bm{A}}})$,
where
$\mbox{\boldmath$\widetilde{u}$}{}_{r,c}=\mbox{\boldmath$\widetilde{u}$}{}-{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{e_{x}}c}}$
is the relative velocity between the incoming flow and the wave phase
velocity, where  $\bm{e_{x}}$ is the unit vector in the $x$-direction and
$\bm{\widetilde{u}}$ are the velocities evaluated at $z=\Delta_{z}/2$. It is
important to note that the phase velocity has been used in the current
formulation to preserve Galilean invariance. The model can be generalized to
waves represented by a wave spectrum. Each wave component imposes a drag
proportional to the relative velocity $\widetilde{\bm{u}}-\bm{e_{x}}{c}$ with
the total drag calculated as a summation of individual effects. This
formulation ensures that the drag for waves with different wavelengths
propagating with different wave speeds is incorporated in the formulation,
with the force contribution from each wavelength calculated separately.

The flow frontal area normal to the incoming flow can be calculated based on
the surface gradient as $\bm{\hat{n}_{u}}\cdot\nabla\eta\Delta_{x}\Delta_{y}$,
where $\bm{\hat{n}_{u}}=\mbox{\boldmath$\widetilde{u}$}{}_{r,c}/U^{\Delta}$
and $U^{\Delta}=\sqrt{(\widetilde{u}-c_{x})^{2}+(\widetilde{v}-c_{y})^{2}}$ is
the tangential velocity relative to the wave propagation speed. The momentum
flux can then be written as
$M_{wave}=\rho\mbox{\boldmath$\widetilde{u}$}{}U^{\Delta}\bm{\widehat{n}_{u}}\cdot\nabla\widetilde{\eta}\Delta_{x}\Delta_{y}$.
In the wind-wave problem, only a fraction of the incoming momentum flux is
imparted as drag and quantified in the model through a drag coefficient
$C_{D}$. The wave drag forcing term is included in the momentum equation by
dividing the momentum flux by the local grid volume. Additionally, the forcing
is only applied if $\bm{\widehat{n}\cdot\nabla\widetilde{\eta}}>0$, that is,
only on the ”front” of the wave with respect to the flow. The wave drag model
is then expressed as

$F_{d,i}=-{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}C_{D}}}\frac{M_{wave}}{\Delta_{x}\Delta_{y}\Delta_{z}}=-C_{D}\frac{\rho}{\Delta_{z}}\widetilde{u}_{i}U^{\Delta}\left(\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right){H\left\\{\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right\\}}\qquad i=x,y,$ (16)

where $H[x]$ is the Heaviside function that ensures that the drag force is
only applied when the flow is incident on the wave frontal area. The wave
surface height is not resolved on the LES grid, and its effects are completely
represented by Equation (16). However, the waves must be well resolved in the
horizontal direction to ensure adequate representation of the wave topology,
that is, surface gradients of $\widetilde{\eta}$.

The drag coefficient can be determined by considering the  energy input rate
from the the wind to the waves as described in Equations (3) and (4). The form
derived by Jeffreys (1925) gives a quadratic dependence of the  energy input
rate for both the wave steepness $ak$ and the inverse wave age
$c/U_{\lambda/2}$ (or equivalently $c/u_{\ast}$). The  energy input rate from
the model for a sinusoidal wave train can be calculated by substituting
$\widetilde{\eta}=a\cos{(k(x-ct))}$ in Equation (16):

$S_{in,LES}=F_{d}\Delta_{z}c=C_{D}\rho(ak)c^{3}\frac{\mbox{$\widetilde{u}$}{}_{\Delta_{z}/2}}{c}\left(\frac{\mbox{$\widetilde{u}$}{}_{\Delta_{z}/2}}{c}-1\right).$
(17)

Note that the  energy input rate is evaluated by integrating the Heaviside
function over one wavelength and assuming that the velocity at the first grid
point is approximately constant and any resulting prefactors are absorbed into
$C_{D}$. The model  energy input rate varies linearly with the wave steepness
and quadratically with wave age. Comparing Equations (3) and (17), the drag
coefficient can be set proportional to the wave steepness $C_{D}=C_{\rm
LES}\cdot ak$ to recover a scaling similar to the Jeffreys (1925)
parameterization for the wave  energy input rate. $C_{\rm LES}$ should be an
order unity coefficient assuming the sheltering coefficient to be a constant.
Peirson and Garcia (2008) analyzed the dependence of the sheltering
coefficient on the wave steepness and found that the values were constant for
$ak>0.1$. For lower wave steepness $ak<0.1$, the sheltering coefficients were
higher. In order to mimic the higher sheltering coefficients for lower wave
steepness, the constant $C_{\rm LES}$ can be adjusted to weakly depend on the
wave steepness. This allows for modulation of the drag coefficient for the low
and high $ak$ limits. For instance, Belcher (1999) developed a theoretical
framework that incorporated a wave coherent tangential stress contribution as
well as a modification to the turbulent stress adjacent to the interface. This
yielded an expression for the form stress that was proportional to $(ak)^{2}$
but modulated by a prefactor
$C_{Belcher}=(\beta_{f}+\beta_{t})/(2+\beta_{f}(ak)^{2})$, where $\beta_{f}$
and $\beta_{t}$ are contributions due to the form and tangential stress
respectively. For values of $\beta_{f}=12$ and $\beta_{t}=20$ (Peirson and
Garcia, 2008), the prefactor varies from $C_{Belcher}=16$ (equal to that
derived by Plant (1982)) at $ak=0$ to $C_{Belcher}=11$ at $ak=0.27$.
Incorporating a similar behavior into $C_{\rm LES}$, the drag coefficient
$C_{D}$ can be written as

$C_{D}=\frac{P}{1+Q(ak)^{2}}ak,$ (18)

where $P=1.2$ and $Q=6$ to be consistent with Belcher (1999). The drag force
per volume applied in LES at the first grid point can now be written as

${F}_{d,i}=-\frac{P}{1+Q(ak)^{2}}(ak)\frac{\rho}{\Delta_{z}}\widetilde{u}_{i}U^{\Delta}\left(\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right){H\left\\{\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right\\}}\qquad i=x,y.$ (19)

The force acts only in the horizontal directions, to remain consistent with
the MOST implementation (Equation (13)). This assumption neglects the effects
of a surface-normal velocity that could become important for swell conditions.
It is important to note that Equation (19) is written for moderate wave phase
speeds or “young-waves”. In the case of swell, where the wave age is much
higher, momentum is transported from the waves to the wind. The efficiency of
momentum transfer from the waves to wind would modify the drag coefficient
used in the model and is left for future studies.

The wave drag model applies a body force calculated based on the wave
characteristics at the first-grid point of the simulation. For vertically
unresolved surface waves, this force can be recast as a stress and applied as
a boundary condition. This formulation is equivalent to the body force
formulation. Casting the wave drag as a body force allows intuitive extension
of the current model to other problems of interest including wind-wave
misalignment and floating wind turbines where the hydrodynamic wave force
plays an important role. An equivalent surface stress and total roughness
formulation are presented here for the sake of completeness.

  1. 1.

Surface stress approach

The total stress at the wall can be written as a sum of the contribution from
the Monin-Obukhov similarity and the wave drag model:

$\widetilde{\tau}_{iz}=\widetilde{\tau}_{iz}|_{wall}(x,y)+\frac{1}{\rho}{F}_{d,i}(x,y,\Delta_{z}/2)\Delta_{z}\qquad
i=x,y.$ (20)

The first term on the right hand side represents the portion of the stress
modeled due to the smooth roughness, and the second term represents the stress
due to the resolved wave field calculated using the wave drag model.

  2. 2.

Equivalent roughness approach

These boundary conditions can also be collectively applied using the standard
logarithmic formulation using a total roughness including the waves, with a
spatially and temporally evolving wave roughness length. Using equations (13)
and (19) in (20) and rewriting for the total stress, the total stress is given
by

$\displaystyle\tau_{iz}^{T}$
$\displaystyle=\left(\left[\frac{\kappa}{\log\left(\frac{\Delta_{z}/2-\eta}{z_{0,s}}\right)}\right]^{2}\mbox{$\widehat{\widetilde{u}}$}{}_{r,i}\mbox{$\widehat{\widetilde{U}}$}{}_{avg}+C_{D}\left(\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right){H\left\\{\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right\\}}{\widetilde{u}_{i}U^{\Delta}}\right).$ (21)

The total roughness then can be calculated as:

$z_{0}^{WDM}=\left(\frac{\Delta_{z}}{2}-\eta\right)\exp{\left[-\kappa\left(\left(\frac{\kappa}{\log\left(\frac{\Delta_{z}/2-\eta}{z_{0,s}}\right)}\right)^{2}+C_{D}\left(\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right){H\left\\{\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right\\}}\right)^{-\frac{1}{2}}\frac{\widetilde{u}_{i}U^{\Delta}}{\mbox{$\widehat{\widetilde{u}}$}{}_{r,i}\mbox{$\widehat{\widetilde{U}}$}{}_{avg}}\right]}.$
(22)

This roughness can then be applied in the simulation to parameterize the
effects of the waves. The wave component of the total roughness depends on the
relative velocity of the airflow to the wave surface and implicitly depends on
the velocity at the first grid point in the LES domain. In the limit of very
slow waves
$\left({\widetilde{u}_{i}U^{\Delta}}\right)/\left(\mbox{$\widehat{\widetilde{u}}$}{}_{r,i}\mbox{$\widehat{\widetilde{U}}$}{}_{avg}\right)\approx
1$, a phase-averaged roughness length can be calculated using Equation 22
assuming the direction of wind propagation is known.

##  4 LES of wind over a sinusoidal wave-train

###  4.1 Simulation setup

The wave surface is specified by
$\widetilde{\eta}(x,t)=a\cos\left({k\left(x-ct\right)}\right)$, and the
orbital velocities are defined using the deep water wave solutions:

$u_{s}=a\omega\cos(k(x-ct));\qquad v_{s}=0;\qquad
w_{s}=a\omega\sin{(k(x-ct))}.$ (23)

The computational domain size in the horizontal directions is
$L_{x}=L_{y}=5\lambda$. The horizontal grid is discretized using sufficient
points in the streamwise and spanwise directions, ensuring that the wave is
horizontally well-resolved on the LES grid with  the grid resolution for each
case provided in Table 1. The vertical extent is of the order of a wavelength
in order to facilitate comparison with the experiments. The discretization
$N_{z}$ is chosen to ensure that the wave lies below the cell center of the
first grid point, $a=0.95\Delta_{z}/2$, and  that the grid aspect ratio is
ideal for wall-modeled LES, i.e $\Delta_{x}/\Delta_{z}\geq 1$ (Piomelli and
Balaras, 2002; Kawai and Larsson, 2012). The vertical resolution
$\Delta_{z}=0.02-0.06H$, where $H$ is the height of the boundary layer, is a
typical resolution for wall-modeled LES (Nicoud et al., 2001; Piomelli and
Balaras, 2002; Kawai and Larsson, 2012; Lee et al., 2013; Bae and Lozano-
Durán, 2021). Therefore, the wave is unresolved in the vertical direction, and
its effects are solely specified by the wave drag model described in the
previous section. Details of the grid and wave properties for the different
configurations are provided in Table 1. Periodic boundary conditions are
applied in the $x$ and $y$ directions, and the flow is driven with a constant
external pressure gradient. The pressure gradient results in a friction
velocity (or surface stress) $u_{\ast}^{2}=\rho^{-1}(\partial P/\partial x)H$.
A free-slip boundary condition is used for the top of the domain. The
simulations are run for a total 200 eddy turnover times $T_{eddy}=H/u_{\ast}$
with averaging carried out over the last 50 turnover times. The simulation
setup is depicted in Figure 1. The height $H$ can be calculated from the
values provided in Table 1 as the wavelength $\lambda$ can be determined using
$\lambda=2\pi/k$ where $k$ is the wavenumber, calculated using the relation
$k=g/c^{2}$.

Case | $ak$ | $c/u_{\ast}$ | $u_{\ast}$ | $N_{x}$ | $N_{z}$ | $H/\lambda$
---|---|---|---|---|---|---
Buckley et al. (2020) | $0.06$ | $6.57$ | $0.073$ | $128$ | 48 | 1
$0.12$ | $3.91$ | $0.167$ | $64$ | $28$ | 1
$0.2$ | $1.8$ | $0.538$ | $64$ | $32$ | 2
Buckley and Veron (2016, 2017) | $0.27$ | $1.4$ | $0.672$ | $48$ | $24$ | $2$
Different Wave Steepnesses | 0.1 | $4.2$ | $0.2$ | $64$ | $16$ | 1
$0.15$ | $4.2$ | $0.2$ | $64$ | $16$ | 1
$0.2$ | $4.2$ | $0.2$ | $64$ | $16$ | 1
Different Wave Ages | $0.1$ | $7.84$ | $0.111$ | $64$ | $16$ | 1
$0.1$ | $11.6$ | $0.076$ | $64$ | $16$ | 1
Table 1: Wave train simulation parameters.

###  4.2 Description of datasets

The model is validated using laboratory-scale experiments from Buckley and
Veron (2016, 2017, 2019) and Buckley et al. (2020). The validation datasets
consist of a wide parameter space with wave steepnesses ranging from $ak=0.06$
to $ak=0.27$ and wave ages from $c/u_{\ast}=1.4$ to $c/u_{\ast}=7$.

Buckley et al. (2020) carried out two-dimensional high resolution measurements
of the airflow above wind-generated waves. The experiment consisted of a tank
specifically designed for the study of air-sea interactions where wind-waves
were generated by a computer-controlled recirculating wind tunnel. The
apparatus and data processing techniques are described in detail in Buckley
and Veron (2017). A total of four cases are presented combining data from
Buckley and Veron (2017) and Buckley et al. (2020). The wave characteristics
of the cases considered are provided in Table 1. The wave characteristics are
chosen based on the peak frequencies from the experiments. The experimental
waves were generated using a wave paddle forced at a particular frequency
resulting in a narrow-banded spectrum.

###  4.3 Results and discussions

####  4.3.1 Velocity Statistics and Shear Stress

In this section, the emphasis is on the instantaneous and mean velocity
profiles obtained using the wave drag model for different wave steepnesses and
wave ages.

Figure 2: Instantaneous velocity profiles from LES-WDM for a) $ak=0.06$, b)
$ak=0.12$, c) $ak=0.21$, and d) $ak=0.27$. Additionally shown are the wave
height distribution in the black line. The effect of the wave phase is
captured in the instantaneous velocity profiles.

The $x-z$ contours at the spanwise mid-plane of the instantaneous streamwise
velocity are shown in Figure 2. A reduced value of streamwise velocity is
observed at the crests in each case. The decrease in streamwise velocity
extends to a larger vertical extent for the cases with higher steepnesses.
Figure 2a depicts the case with the lowest steepness and is closest to a rough
wall turbulent boundary layer with mild phase effects. For the steepest wave
case, the phase effects extend to about 20% of the boundary layer height. The
“phase-aware” nature of the model is well captured in the instantaneous
velocity contours.

To highlight the wave-coherent streamwise motions, a triple decomposition of
the velocity field is used (Hara and Sullivan (2015)). The field is first
split into an ensemble (Reynolds) average and a turbulent fluctuation,
$\phi=\overline{\phi}+\phi^{\prime}$, where the prime denotes the turbulent
fluctuation. The ensemble average $\overline{\phi}$ is further split into the
horizontal mean $\langle\phi\rangle$ and a wave fluctuation ${\phi^{w}}$. A
triple decomposition of the velocity field can then be written as

$\widetilde{u}=\overline{\tilde{u}}+\tilde{u}^{\prime}=\langle\tilde{u}\rangle+\tilde{u}^{w}+\tilde{u}^{\prime}.$
(24)

Figure 3 shows the normalized streamwise wave fluctuation velocity contours
near the wave surface for $ak=0.06$. The streamwise velocity fluctuation from
the model shows good qualitative agreement with the experiments of Buckley and
Veron (2019, see Figure 4a) above $kz=0.2$. At the wave surface there are
discrepancies since the current model does not include fluctuations in the
velocity due to the vertical wave motions.

Figure 3:  Contours of the normalized streamwise wave fluctuation velocity
$\tilde{u}^{w}$ near the wave surface for $ak=0.06$ and $U_{10}=2.19$ m/s. At
the first gridpoint above the wave surface the velocity decelerates as the
flow encounters the wave crest (due to the applied drag force) and accelerates
downwind. The phase aware nature of the model manifests in the form of the
alternate positive and negative wave fluctuation velocities.

In Figure 4, the mean velocity profiles are shown from the LES with the wave
drag model, hereafter referred to as LES-WDM, for the different wave
steepnesses and compared to the experimental data. As the wave steepness
increases and the wave age decreases, the mean wind velocity measured by the
experiments is reduced due to increasing resistance from the waves. The LES-
WDM accurately captures this effect with slower velocities predicted for
higher wave steepness and low wave age.

Figure 4a shows the comparison for $ak=0.06$. For small wave steepnesses the
deviation from the log-law is small, with $\kappa=0.4$. The LES-WDM
generalizes across the different experimental conditions and shows good
agreement for all cases considered. The effective roughness lengths from the
experiment and the models are compared in Table 2. The roughness length for
the experimental data is calculated by fitting a logarithmic profile to the
data. The roughness length from the simulation can be calculated similarly or
by averaging equation Equation (22) and using the velocity profile from the
LES. Additionally, for each case, the roughness lengths obtained by using the
Charnock model and the steepness-dependent roughness model, calculated using
Equations (6) and (8), respectively, are shown in Table 2. Figure 5a compares
the modeled experimental results of Table 2 in the form of a parity plot. The
total roughness, in general, can depend on the wind speed (or friction
velocity $u_{*}$), the wave age $c/u_{*}$, and/or the wave steepness $ak$. In
the experimental dataset considered, all three parameters are changing for the
different cases. An effective Charnock parameter can be calculated by
normalizing the total surface roughness by $\alpha_{ch,eff}=u_{*}^{2}/g$ and
is shown as a function of wave steepness in Figure 5b. The normalized
roughness for the Charnock model is a constant $\alpha_{ch}=0.11$ as the
length scale predicted by the Charnock model only depends on the friction
velocity and has no explicit dependence on the wave steepness or wave age. In
the Taylor model, the normalized roughness increases as a function of wave
steepness. It is evident from the experimental data that the effect of the
waves are more complicated as the normalized roughness depends on both the
wave age and the steepness. This trend is also captured in the current LES-
WDM. Modeling the form drag using the LES-WDM provides remarkable improvement
over the standard phase-averaged approaches (Eq. 6 & 8) as it incorporates the
effect of the wind speed, the wave age, and the wave steepness. The effective
roughness length scale predicted by the phase-averaged model of Charnock
(1955) is smaller than the corresponding experimental roughness, while the
roughness length calculated using the model proposed by Taylor and Yelland
(2001) is larger. The steepness dependent model only accounts for the wave
steepness effect, neglecting the wave propagation speed. The results from the
Charnock model could be improved by using a higher value of the parameter
$\alpha_{ch}$; however, this constant cannot be determined _a priori_ and
would vary with the other wave characteristics.

Figure 4: Mean velocity profiles from LES-WDM and experiments of Buckley et
al. (2020) for a) $ak=0.06$, b) $ak=0.12$, c) $ak=0.21$, and d) $ak=0.27$. The
wave phase-resolved LES from Husain et al. (2019) is shown in d).

Note that the LES-WDM formulation is developed for coarse grids and the focus
is on predicting the velocity profiles above the peak wave amplitude.
Therefore, for the higher wave steepnesses only a few points lie in the range
of the experimental data as the near-field of the wave is unresolved.

$ak$ | $c/u_{\ast}$ | $z_{0,Expt.}$ | $z_{0,LES}$ | $z_{0,Charnock}$ | $z_{0,Taylor}$
---|---|---|---|---|---
$0.06$ | $6.57$ | $5.5\times 10^{-5}$ | $6\times 10^{-5}$ | $8.15\times 10^{-6}$ | $6.69\times 10^{-6}$
$0.12$ | $3.91$ | $6.0\times 10^{-5}$ | $7\times 10^{-5}$ | $4.26\times 10^{-5}$ | $1.60\times 10^{-4}$
$0.2$ | $1.8$ | $2.3\times 10^{-4}$ | $4\times 10^{-4}$ | $4.42\times 10^{-4}$ | $4.43\times 10^{-3}$
$0.27$ | $1.4$ | $1.0\times 10^{-3}$ | $1.4\times 10^{-3}$ | $6.91\times 10^{-4}$ | $1.32\times 10^{-2}$
Table 2: Total roughness length calculated from the experimental data, the
phase-averaged roughness models, and the LES-WDM for the different cases.
Figure 5: Total roughness length calculated from the experimental data and
roughness models (Eqs. (6), (8) and (22)): a) parity plot of experimental and
modeled phase-averaged roughness length and b) variation of normalized total
roughness with wave steepness. The solid line in the left plane is where
$z_{0,expt}=z_{0,mod}$.

To highlight and isolate the effects of wave age and steepness on the airflow,
two additional sets of simulations are considered. First, the wave steepness
is fixed at $ak=0.1$, and the wave age varies as $c/u_{*}=4.2$,
$c/u_{*}=7.84$, and $c/u_{*}=11.5$. The second set of simulations has fixed
wave age $c/u_{*}=4.2$ and wave steepnesses $ak=0.1$, $ak=0.15$, and $ak=0.2$.
Figure 6a shows that the normalized mean velocity is faster for higher wave
ages. Based on Equation 19, increasing the wave phase speed reduces the wave
form drag applied by the model. The wave offers smaller resistance to the
wind, and the blocking effect on the airflow is diminished. In Figure 6b the
effect of only changing the wave steepness is highlighted. Larger steepness of
the wave results in higher blockage thereby applying more drag on the airflow.
Additionally, the reduction in the overall mean velocity profile is not
uniform between $ak=0.1$, $ak=0.15$ and $ak=0.2$ as the effect of the
steepness on the wave drag is not linear.

Figure 6: Normalized mean velocity profiles for a) different wave ages
$c/u_{\ast}=4.2,\ 7.84$, and $11.5$ with wave steepness fixed at $ak=0.1$ and
b) different wave steepnesses $ak=0.1$, $ak=0.15$, and $ak=0.2$ and wave age
fixed at $c/u_{*}=4.2$. The mean velocity is faster for higher wave ages and
lower wave steepnesses.

The averaged total shear stress consists of the averaged resolved Reynolds
stress
$\langle\mbox{$\widetilde{u}$}{}^{\prime}\mbox{$\widetilde{w}$}{}^{\prime}\rangle$,
where
$\mbox{$\widetilde{u}$}{}^{\prime}=\mbox{$\widetilde{u}$}{}-\overline{\mbox{$\widetilde{u}$}{}}$
is the fluctuating part of the resolved LES velocity obtained by subtracting
the time averaged part of the velocity, and the mean subfilter shear stress
$\langle{\mbox{ $\widetilde{\tau}$}{}}_{xz}\rangle$. The normalized resolved
and subfilter shear stresses from the four validation cases are shown in
Figure 7. For each case, the normalized resolved stress shown in Figure 7a
increases linearly from the top boundary towards the surface and is the
dominant contribution to the total stress. The mean normalized subfilter
stress depicted in Figure 7b is negligible far from the surface and starts to
increase rapidly below $z/H=0.1$. The total normalized stress shown in Figure
7c is similar to a pressure gradient driven turbulent boundary layer. The
total normalized shear stress does not reach $-1$ at the surface  as it does
not include the model form stress. At the surface, the total stress
$\widetilde{\tau}|_{wall}=\tau_{xz}+\rho^{-1}{F}_{d}\Delta z$ is the sum of
the form stress imposed by the model and the subfilter stress (from MOST).
Figure 8 shows the partition of the total surface stress as a function of time
for the different cases and verifies that the wave stress and the subfilter
stress sums to the imposed non-dimensional pressure gradient for each case
considered.

Figure 7: Mean normalized a) resolved stress, b) subfilter stress, and c)
total shear stress profiles for $ak=0.06$, $c/u_{\ast}=6.57$; $ak=0.12$,
$c/u_{\ast}=3.91$; $ak=0.2$, $c/u_{\ast}=1.8$; and $ak=0.27$,
$c/u_{\ast}=1.4$. The resolved stress peaks just above the surface for each
case. The magnitude of the normalized subfilter surface stress is less than
the imposed stress, $|\mbox{ $\widetilde{\tau}$}{}_{xz}|_{wall}|<1$, with the
remaining contribution from the wave drag model. Note $H=\lambda$ for
$ak=0.06$ and $ak=0.12$, and $H=2\lambda$ for $ak=0.2$ and $ak=0.27$. Figure
8:  Normalized form stress ($\Delta_{z}F_{d,i}$) and total stress
($\langle\mbox{ $\widetilde{\tau}$}{}_{i3}^{T}\rangle$) as a function of non-
dimensional time for different wave steepnesses, and wave ages. The solid
lines represent the form stress, and the corresponding dashed lines represent
the total stress for $ak=0.2$, $c/u_{\ast}=1.8$; $ak=0.12$, $c/u_{\ast}=3.91$;
and $ak=0.06$, $c/u_{\ast}=6.57$.The contribution of the form stress increases
as a function of wave steepness.

####  4.3.2 Model Form Drag

To understand the effects of the wave drag model in the LES, the spanwise-
averaged velocity and the wave form stress calculated from the drag force
applied at the first grid point are analyzed in Figure 9 for $ak=0.12$. The
form stress for phase-resolved simulations is obtained from the resolved
pressure $\widetilde{p}\partial\widetilde{\eta}/\partial x$.

Figure 9:  The top panel depicts the scaled wave surface gradient distribution
and the spanwise averaged normalized form stress $D_{p}$ at two time
instances. Additionally shown in the top panel is the time averaged normalized
form stress. The wave drag force is correlated to the wave slope and has a
phase difference of $\pi/2$ with respect to the wave height distribution. The
bottom panel depicts the spanwise averaged normalized streamwise velocity at
the wave surface as a function of horizontal distance at two time instances.
Also shown in the bottom panel is the corresponding wave height distribution
(magenta lines). The velocity shows the characteristic acceleration on the
leeward side of the wave with respect to the wave trough followed by
deceleration. The data corresponds to the case with $ak=0.12$ and
$c/u_{*}=3.91$.

In the current LES, as the waves are unresolved in the vertical direction, the
model form stress $D_{p}$ depends on the incoming momentum flux $M_{wave}$ and
the drag coefficient for the wind wave transfer.  The stress is calculated
using the drag force per unit volume ${F}_{d}$ defined in Equation 19
multiplied by the vertical grid spacing $\Delta_{z}$:

$D_{p}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{C_{D}}{\rho}}}\frac{M_{wave}}{(\Delta_{x}\Delta_{y})}=\frac{1}{\rho}|{\bm{F}}_{d}|\Delta_{z}=\frac{P}{1+Q(ak)^{2}}(ak)\mbox{$\widetilde{u}$}{}U^{\Delta}\left(\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right){H\left\\{\widehat{n}_{u,k}\cdot\frac{\partial\widetilde{\eta}}{\partial
x_{k}}\right\\}}.$ (25)

The top panel depicts the spanwise averaged model form stress for two time
instants. The drag force is correlated to the surface gradient of the wave
height distribution and is out of phase with the wave surface (shown in the
bottom panel) by a factor of $\pi/2$. The magnitude of the peaks of the drag
force shows temporal dependence due to the variability of the turbulent
velocity at the first grid point (shown in the bottom panel). Additionally
shown in the top panel is the time and space averaged value of the normalized
form stress. There is considerable variation of the drag applied along the
wave surface with respect to the mean. This effect is also observed in
experiments and phase-resolved simulations (Hara and Sullivan, 2015; Sullivan
et al., 2018; Husain et al., 2019; Buckley et al., 2020; Yousefi et al.,
2020). The bottom panel of Figure 9 depicts the spanwise average of the
streamwise velocity at the wave surface (first grid point) for two distinct
time steps. The wave height distribution is also depicted for reference. The
velocity accelerates in the upstream side of the surface peak where
correspondingly the force applied by the model (top panel of Figure 9) is zero
and then decelerates as the flow encounters the next wave surface
corresponding to a maximum in the applied force. Phase-averaged models purely
apply the averaged drag and would not capture this variation in the velocity
and drag. This variation is important for calculating hydrodynamic forces for
offshore structures.

Figure 10a highlights the average normalized form stress defined in Equation
25 as a function of wave steepness. The LES-WDM shows good agreement with the
experimental data from Banner (1990), Banner and Peirson (1998), and Peirson
and Garcia (2008) and accurately predicts the change of form stress for
different steepness values considered here. The stresses from Buckley et al.
(2020) are higher than the other reported experimental values, but the
variation of form stress as a function of wave-steepness shows similar
qualitative trend. The experimental scatter is in part due to the
uncertainties in the measurements. Note that the form stress calculated here
is the stress due to the wave drag model acting on the first grid point. The
effect of changing wave age for a fixed wave steepness reduces the form stress
as the faster waves provide smaller resistance to the airflow above. The
normalized form stress as a function of the wave age $c/u_{\ast}$ is shown in
Figure 10b. The simulations are run at a fixed steepness, and the grid is set
such that $a=0.95\Delta z/2$. The LES-WDM accurately predicts the form stress
from the DNS datasets of Sullivan et al. (2000), Kihara et al. (2007), and
Yang and Shen (2010) as a function of wave age for $c/u_{\ast}<15$. For very
fast waves ($\mbox{\boldmath$\widetilde{u}$}{}_{r}<0$ at the first grid
point), the model predicts a negative form stress based on Equation (19),
following similar qualitative trends to the DNS, but the magnitude of the
negative form stress is too large. The transfer of momentum is from the waves
to the wind, a scenario common during conditions of swell, where fast moving
waves (with $c/u_{*}>25$) generated from distant systems transfer energy to
the slower wind. Future studies could work towards developing a model for fast
moving waves.

Figure 10:  Normalized form stress as a function of a) wave steepness and b)
wave age. In the left panel, the symbols correspond to simulation results from
the LES-WDM (current work), LES from Husain et al. (2019), and experimental
data from Banner (1990), Banner and Peirson (1998), Peirson and Garcia (2008),
and Buckley et al. (2020). In the left panel, apart from the LES-WDM
validation cases, the form stress from the LES at $ak=0.1$ and $c/u_{*}=4.2,\
7.84$ and $11.6$ are shown. In the right panel, the LES-WDM at different wave
ages $c/u_{*}=4.2,\ 7.84$ and $11.6$ and $ak=0.1$ is compared to DNS data from
Sullivan et al. (2000), Kihara et al. (2007), and Yang and Shen (2010). The
model shows good agreement with the existing experimental and numerical data.

##  5 Summary and Conclusions

Characterizing the momentum transfer between wind and waves is critical in
numerous atmospheric applications. Wave phase-resolving simulations through
immersed boundary or terrain following approaches come with a high
computational cost for simulating atmospheric flows. On the other extreme,
wave phase-averaged approaches with an effective roughness length are not
sufficiently generalizable to a broad range of conditions. In this work a
wave-based drag model applicable to wall-modeled Large Eddy Simulations is
proposed. The model quantifies the momentum transfer between the wind and the
wave-field through a canopy stress approach. The surface-gradient based
approach of Anderson and Meneveau (2010) is extended to moving boundaries and
the formulation adapted for the wind-wave problem. The drag coefficient is
found to be proportional to the wave steepness to give the correct scaling for
the form stress.

The model is applied and tested on various cases of turbulent boundary layer
flow over idealized water waves of different wave steepnesses and wave ages.
The ”phase-aware” nature of the model is well captured by the instantaneous
velocity profiles and the phase-averaged streamwise wave fluctuation velocity.
The predicted mean velocity profiles show good agreement with data from
experiments. The form stress predicted by the model is in good agreement with
existing experimental data and shows the correct trend with both wave age and
wave steepness. For faster waves, the model is capable of predicting a
negative form stress, where the momentum is transferred from the wave-field to
the wind. The model showed generalizability over the parameter space studied
when compared to phase-averaged roughness models and requires no tuning
parameters or empirical constants. The inherent spatial and temporal
variability associated with the the drag force imposed by the waves on the
airflow can also be well represented.  Vertical wave induced motions are not
accounted for explicitly in the model; however, both the wave-induced stress
and the pressure stress are captured implicitly. Incorporating the effects of
wave-induced motions to the vertical velocity is important for swell
conditions (where the waves move faster than the near-surface wind velocity)
and will be the focus of future studies.

The framework proposed here to incorporate wave effects can serve as a low-
cost accurate LES framework to simulate airflow in the marine atmospheric
boundary layer with a significant reduction in computational cost compared to
phase-resolved approaches. The coarse resolution used here allows the accurate
prediction of the mean velocity profiles and wave stress with an $O(10^{4})$
reduction in computational cost based on grid resolution. The generalizability
and ease of implementation of the wave drag model allows scalability to large
domains and could be useful in the context of offshore wind farm (fixed bottom
and floating turbines) simulations, marine atmospheric boundary layer
simulations, and ocean mixed layer simulations.

###### Acknowledgements.

The authors gratefully acknowledge financial support from the Princeton
University Andlinger Center for Energy and the Environment and High Meadows
Environmental Institute. The simulations presented in this article were
performed on computational resources supported by the Princeton Institute for
Computational Science and Engineering (PICSciE) and the Office of Information
Technology’s High Performance Computing Center and Visualization Laboratory at
Princeton University. Data can be made available upon request.

## References

  * Abkar and Moin (2017) Abkar, M., and P. Moin, 2017: Large-eddy simulation of thermally stratified atmospheric boundary-layer flow using a minimum dissipation model. Boundary-Layer Meteorology, 165 (3), 405–419, 10.1007/s10546-017-0288-4.
  * Anderson and Meneveau (2010) Anderson, W., and C. Meneveau, 2010: A Large-Eddy Simulation Model for Boundary-Layer Flow Over Surfaces with Horizontally Resolved but Vertically Unresolved Roughness Elements. Boundary-Layer Meteorology, 137 (3), 397–415, 10.1007/s10546-010-9537-5.
  * Anderson and Meneveau (2011) Anderson, W., and C. Meneveau, 2011: Dynamic roughness model for large-eddy simulation of turbulent flow over multiscale, fractal-like rough surfaces. Journal of Fluid Mechanics, 679, 288–314, 10.1017/jfm.2011.137.
  * Arthur et al. (2019) Arthur, R. S., J. D. Mirocha, K. A. Lundquist, and R. L. Street, 2019: Using a canopy model framework to improve large-eddy simulations of the neutral atmospheric boundary layer in the weather research and forecasting model. Monthly Weather Review, 147 (1), 31–52, 10.1175/MWR-D-18-0204.1.
  * Bae and Lozano-Durán (2021) Bae, H. J., and A. Lozano-Durán, 2021: Effect of wall boundary conditions on a wall-modeled large-eddy simulation in a finite-difference framework. Fluids, 6 (3), 10.3390/fluids6030112, URL https://www.mdpi.com/2311-5521/6/3/112.
  * Banner (1990) Banner, M. L., 1990: The influence of wave breaking on the surface pressure distribution in wind—wave interactions. Journal of Fluid Mechanics, 211, 463–495, 10.1017/S0022112090001653.
  * Banner and Peirson (1998) Banner, M. L., and W. L. Peirson, 1998: Tangential stress beneath wind-driven air–water interfaces. Journal of Fluid Mechanics, 364, 115–145, 10.1017/S0022112098001128.
  * Belcher (1999) Belcher, S., 1999: Wave growth by non-separated sheltering. European Journal of Mechanics - B/Fluids, 18 (3), 447–462, 10.1016/S0997-7546(99)80041-7, three-Dimensional Aspects of Air-Sea Interaction.
  * Bou-Zeid et al. (2005) Bou-Zeid, E., C. Meneveau, and M. Parlange, 2005: A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows. Physics of Fluids, 17 (2), 1–18, 10.1063/1.1839152.
  * Brown et al. (2001) Brown, A. R., J. M. Hobson, and N. Wood, 2001: Large-eddy simulation of neutral turbulent flow over rough sinusoidal ridges. Boundary-Layer Meteorology, 98 (3), 411–441, 10.1023/A:1018703209408.
  * Buckley and Veron (2016) Buckley, M. P., and F. Veron, 2016: Structure of the airflow above surface waves. Journal of Physical Oceanography, 46 (5), 1377 – 1397, 10.1175/JPO-D-15-0135.1.
  * Buckley and Veron (2017) Buckley, M. P., and F. Veron, 2017: Airflow measurements at a wavy air–water interface using piv and lif. Experiments in Fluids, 58 (11), 161.
  * Buckley and Veron (2019) Buckley, M. P., and F. Veron, 2019: The turbulent airflow over wind generated surface waves. European Journal of Mechanics - B/Fluids, 73, 132–143, https://doi.org/10.1016/j.euromechflu.2018.04.003, breaking Waves.
  * Buckley et al. (2020) Buckley, M. P., F. Veron, and K. Yousefi, 2020: Surface viscous stress over wind-driven waves with intermittent airflow separation. Journal of Fluid Mechanics, 905, 10.1017/jfm.2020.760.
  * Cao et al. (2020) Cao, T., B. Q. Deng, and L. Shen, 2020: A simulation-based mechanistic study of turbulent wind blowing over opposing water waves. Journal of Fluid Mechanics, 10.1017/jfm.2020.591.
  * Cao and Shen (2021) Cao, T., and L. Shen, 2021: A numerical and theoretical study of wind over fast-propagating water waves. Journal of Fluid Mechanics, 919, A38, 10.1017/jfm.2021.416.
  * Charnock (1955) Charnock, H., 1955: Wind stress on a water surface. Quarterly Journal of the Royal Meteorological Society, 81 (350), 639–640, https://doi.org/10.1002/qj.49708135027.
  * Desjardins et al. (2008) Desjardins, O., G. Blanquart, G. Balarac, and H. Pitsch, 2008: High order conservative finite difference scheme for variable density low mach number turbulent flows. Journal of Computational Physics, 227 (15), 7125–7159, https://doi.org/10.1016/j.jcp.2008.03.027.
  * Deskos et al. (2021) Deskos, G., J. C. Y. Lee, C. Draxl, and M. A. Sprague, 2021: Review of wind-wave coupling models for large-eddy simulation of the marine atmospheric boundary layer. Journal of the Atmospheric Sciences, 1–75, 10.1175/jas-d-21-0003.1.
  * Donelan (1990) Donelan, M. A., 1990: Air-sea interaction. Ocean Engineering Science, 9B, 239–292.
  * Donelan et al. (2006) Donelan, M. A., A. V. Babanin, I. R. Young, and M. L. Banner, 2006: Wave-follower field measurements of the wind-input spectral function. Part II: Parameterization of the wind input. Journal of Physical Oceanography, 36 (8), 1672–1689, 10.1175/JPO2933.1.
  * Drennan et al. (2005) Drennan, W. M., P. K. Taylor, and M. J. Yelland, 2005: Parameterizing the sea surface roughness. Journal of Physical Oceanography, 35 (5), 835–848, 10.1175/JPO2704.1.
  * Druzhinin et al. (2012) Druzhinin, O. A., Y. I. Troitskaya, and S. S. Zilitinkevich, 2012: Direct numerical simulation of a turbulent wind over a wavy water surface. Journal of Geophysical Research: Oceans, 117 (C11), https://doi.org/10.1029/2011JC007789.
  * Edson et al. (2013) Edson, J. B., and Coauthors, 2013: On the exchange of momentum over the open ocean. Journal of Physical Oceanography, 43 (8), 1589–1610, 10.1175/JPO-D-12-0173.1.
  * Fairall et al. (2003) Fairall, C. W., E. F. Bradley, J. E. Hare, A. A. Grachev, and J. B. Edson, 2003: Bulk Parameterization of Air Sea Fluxes: Updates and Verification for the COARE Algorithm. Journal of Climate, 16, 571–591.
  * Fairall et al. (1996) Fairall, C. W., E. F. Bradley, D. P. Rogers, J. B. Edson, and G. S. Young, 1996: Bulk parameterization of air-sea fluxes for tropical ocean-global atmosphere coupled-ocean atmosphere response experiment. Journal of Geophysical Research: Oceans, 101 (C2), 3747–3764, https://doi.org/10.1029/95JC03205.
  * Gadde et al. (2021) Gadde, S. N., A. Stieren, and R. J. A. M. Stevens, 2021: Large-eddy simulations of stratified atmospheric boundary layers: Comparison of different subgrid models. Boundary-Layer Meteorology, 178 (3), 363–382, 10.1007/s10546-020-00570-5, URL https://doi.org/10.1007/s10546-020-00570-5.
  * Garratt (1977) Garratt, J. R., 1977: Review of drag coefficients over oceans and continents. Monthly Weather Review, 105 (7), 915 – 929, 10.1175/1520-0493(1977)105¡0915:RODCOO¿2.0.CO;2.
  * Grare et al. (2013) Grare, L., W. L. Peirson, H. Branger, J. W. Walker, J. P. Giovanangeli, and V. Makin, 2013: Growth and dissipation of wind-forced, deep-water waves. Journal of Fluid Mechanics, 722, 5–50, 10.1017/jfm.2013.88.
  * Hao et al. (2018) Hao, X., T. Cao, Z. Yang, T. Li, and L. Shen, 2018: Simulation-based study of wind-wave interaction. Procedia IUTAM, 26, 162–173.
  * Hara and Sullivan (2015) Hara, T., and P. P. Sullivan, 2015: Wave boundary layer turbulence over surface waves in a strongly forced condition. Journal of Physical Oceanography, 45 (3), 868–883, 10.1175/JPO-D-14-0116.1.
  * Husain et al. (2019) Husain, N. T., T. Hara, M. P. Buckley, K. Yousefi, F. Veron, and P. P. Sullivan, 2019: Boundary layer turbulence over surface waves in a strongly forced condition: LES and observation. Journal of Physical Oceanography, 49 (8), 1997–2015, 10.1175/JPO-D-19-0070.1.
  * Janssen (1991) Janssen, P. A. E. M., 1991: Quasi-linear Theory of Wind-Wave Generation Applied to Wave Forecasting. Journal of Physical Oceanography, 21 (11), 1631–1642.
  * Jeffreys (1925) Jeffreys, H., 1925: On the formation of water waves by wind. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 107 (742), 189–206.
  * Jiménez and Dudhia (2018) Jiménez, P. A., and J. Dudhia, 2018: On the need to modify the sea surface roughness formulation over shallow waters. Journal of Applied Meteorology and Climatology, 57 (5), 1101–1110, 10.1175/JAMC-D-17-0137.1.
  * Kawai and Larsson (2012) Kawai, S., and J. Larsson, 2012: Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy. Physics of Fluids, 24 (1), 015 105, 10.1063/1.3678331, URL https://doi.org/10.1063/1.3678331, https://doi.org/10.1063/1.3678331.
  * Kihara et al. (2007) Kihara, N., H. Hanazaki, T. Mizuya, and H. Ueda, 2007: Relationship between airflow at the critical height and momentum transfer to the traveling waves. Physics of Fluids, 19 (1), 015 102, 10.1063/1.2409736, URL https://doi.org/10.1063/1.2409736, https://doi.org/10.1063/1.2409736.
  * Kudryavtsev and Makin (2004) Kudryavtsev, V. N., and V. K. Makin, 2004: Impact of swell on the marine atmospheric boundary layer. Journal of Physical Oceanography, 34 (4), 934–949, 10.1175/1520-0485(2004)034¡0934:IOSOTM¿2.0.CO;2.
  * Lee et al. (2013) Lee, J., M. Cho, and H. Choi, 2013: Large eddy simulations of turbulent channel and boundary layer flows at high reynolds number with mean wall shear stress boundary condition. Physics of Fluids, 25 (11), 110 808, 10.1063/1.4819342, URL https://doi.org/10.1063/1.4819342, https://doi.org/10.1063/1.4819342.
  * Lin et al. (2008) Lin, M. Y., C. H. Moeng, W. T. Tsai, P. P. Sullivan, and S. E. Belcher, 2008: Direct numerical simulation of wind-wave generation processes. J. of Fluid Mech., 616, 1–30, 10.1017/S0022112008004060.
  * Longuet-Higgins (1969) Longuet-Higgins, M. S., 1969: A nonlinear mechanism for the generation of sea waves. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 311 (1506), 371–389, 10.1098/rspa.1969.0123, URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1969.0123, https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1969.0123.
  * MacArt and Mueller (2016) MacArt, J. F., and M. E. Mueller, 2016: Semi-implicit iterative methods for low mach number turbulent reacting flows: Operator splitting versus approximate factorization. Journal of Computational Physics, 326, 569–595, https://doi.org/10.1016/j.jcp.2016.09.016.
  * Miles (1993) Miles, J., 1993: Surface-Wave Generation Revisited. Journal of Fluid Mechanics, 256, 427–441, 10.1017/S0022112093002836.
  * Miles (1957) Miles, J. W., 1957: On the generation of surface waves by shear flows. Journal of Fluid Mechanics, 3 (2), 185–204, 10.1017/S0022112057000567.
  * Moeng (1984) Moeng, C., 1984: A Large-Eddy-Simulation Model for the Study of Planetary Boundary-Layer Turbulence. Journal of Atmospheric Sciences, 41 (13), 2052–2062.
  * Nicoud et al. (2001) Nicoud, F., J. S. Baggett, P. Moin, and W. Cabot, 2001: Large eddy simulation wall-modeling based on suboptimal control theory and linear stochastic estimation. Physics of Fluids, 13 (10), 2968–2984, 10.1063/1.1389286, URL https://doi.org/10.1063/1.1389286, https://doi.org/10.1063/1.1389286.
  * Peirson and Garcia (2008) Peirson, W. I., and A. W. Garcia, 2008: On the wind-induced growth of slow water waves of finite steepness. Journal of Fluid Mechanics, 608, 243–274, 10.1017/S002211200800205X.
  * Piomelli and Balaras (2002) Piomelli, U., and E. Balaras, 2002: Wall-layer models for large-eddy simulations. Annual Review of Fluid Mechanics, 34 (1), 349–374, 10.1146/annurev.fluid.34.082901.144919.
  * Plant (1982) Plant, W. J., 1982: A relationship between wind stress and wave slope. Journal of Geophysical Research: Oceans, 87 (C3), 1961–1967, https://doi.org/10.1029/JC087iC03p01961, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JC087iC03p01961.
  * Rozema et al. (2015) Rozema, W., H. J. Bae, P. Moin, and R. Verstappen, 2015: Minimum-dissipation models for large-eddy simulation. Physics of Fluids, 27 (8), 085 107, 10.1063/1.4928700, https://aip.scitation.org/doi/pdf/10.1063/1.4928700.
  * Shaw and Schumann (1992) Shaw, R. H., and U. Schumann, 1992: Large-eddy simulation of turbulent flow above and within a forest. Boundary-Layer Meteorology, 61 (1-2), 47–64, 10.1007/BF02033994.
  * Shemdin and Hsu (1967) Shemdin, O. H., and E. Y. Hsu, 1967: Direct measurement of aerodynamic pressure above a simple progressive gravity wave. Journal of Fluid Mechanics, 30 (2), 403–416, 10.1017/S0022112067001508.
  * Smith et al. (1992) Smith, S. D., and Coauthors, 1992: Sea surface wind stress and drag coefficients: The hexos results. Boundary-Layer Meteorology, 60 (1-2), 109–142, 10.1007/BF00122064.
  * Snyder et al. (1981) Snyder, R. L., F. W. Dobson, J. A. Elliott, and R. B. Long, 1981: Array measurements of atmospheric pressure fluctuations above surface gravity waves. Journal of Fluid Mechanics, 102, 1–59, 10.1017/S0022112081002528.
  * Su et al. (1998) Su, H. B., R. H. Shaw, K. T. Pawu, C. H. Moeng, and P. P. Sullivan, 1998: Turbulent statistics of neutrally stratified flow within and above a sparse forest from large-eddy simulation and field observations. Boundary-Layer Meteorology, 88 (3), 363–397, 10.1023/A:1001108411184.
  * Sullivan et al. (2018) Sullivan, P. P., M. L. Banner, R. P. Morison, and W. L. Peirson, 2018: Turbulent flow over steep steady and unsteady waves under strong wind forcing. Journal of Physical Oceanography, 48 (1), 3–27, 10.1175/JPO-D-17-0118.1.
  * Sullivan et al. (2008) Sullivan, P. P., J. B. Edson, T. Hristov, and J. C. McWilliams, 2008: Large-eddy simulations and observations of atmospheric marine boundary layers above nonequilibrium surface waves. Journal of the Atmospheric Sciences, 65 (4), 1225–1245, 10.1175/2007JAS2427.1.
  * Sullivan et al. (2000) Sullivan, P. P., J. C. McWilliams, and C. H. Moeng, 2000: Simulation of turbulent flow over idealized water waves. Journal of Fluid Mechanics, 404, 47–85, 10.1017/S0022112099006965.
  * Sullivan et al. (2014) Sullivan, P. P., J. C. Mcwilliams, and E. G. Patton, 2014: Large-eddy simulation of marine atmospheric boundary layers above a spectrum of moving waves. Journal of the Atmospheric Sciences, 71 (11), 4001–4027, 10.1175/JAS-D-14-0095.1.
  * Taylor and Yelland (2001) Taylor, P. K., and M. J. Yelland, 2001: The dependence of sea surface roughness on the height and steepness of the waves. Journal of Physical Oceanography, 31 (2), 572 – 590, 10.1175/1520-0485(2001)031¡0572:TDOSSR¿2.0.CO;2.
  * Toba et al. (1990) Toba, Y., N. Iida, H. Kawamura, N. Ebuchi, and I. Jones, 1990: Wave dependence of sea-surface wind stress. Journal of Physical Oceanography, 20, 705–721.
  * Yang et al. (2013) Yang, D., C. Meneveau, and L. Shen, 2013: Dynamic modelling of sea-surface roughness for large-eddy simulation of wind over ocean wavefield. Journal of Fluid Mechanics, 726, 62–99, 10.1017/jfm.2013.215.
  * Yang et al. (2014) Yang, D., C. Meneveau, and L. Shen, 2014: Effect of downwind swells on offshore wind energy harvesting – A large-eddy simulation study. Renewable Energy, 70, 11–23, 10.1016/J.RENENE.2014.03.069.
  * Yang and Shen (2010) Yang, D., and L. Shen, 2010: Direct-simulation-based study of turbulent flow over various waving boundaries. Journal of Fluid Mechanics, 650, 131–180, 10.1017/S0022112009993557.
  * Yousefi et al. (2020) Yousefi, K., F. Veron, and M. P. Buckley, 2020: Momentum flux measurements in the airflow over wind-generated surface waves. Journal of Fluid Mechanics, 895, A15, 10.1017/jfm.2020.276.