# Generation of attached Langmuir circulations by a suspended macroalgal farm

Chao Yan  James C. McWilliams  Marcelo Chamecki<EMAIL_ADDRESS>Department
of Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA
90025, USA

###### Abstract

In this study, we focus on Langmuir turbulence in the deep ocean with the
presence of a large macroalgal farm using a Large Eddy Simulation method. The
wave-current interactions are modelled by solving the wave-averaged equations.
The hydrodynamic process over the farm is found to drive a persistent flow
pattern similar to Langmuir circulations but is locked in space across the
farm. These secondary circulations are generated because the cross-stream
shear produced by the rows of canopy elements leads to a steady vertical
vorticity field, which is then rotated to the downstream direction under the
effect of vortex force. Since the driving mechanism is similar to the Craik-
Leibovich type 2 instability theory, these secondary circulations are also
termed as attached Langmuir circulations. We then apply a triple decomposition
on the flow field to unveil the underlying kinematics and energy transfer
between the mean flow, the secondary flow resulting from the farm drag, and
the transient eddies. Flow visualizations and statistics suggest that the
attached Langmuir circulations result from the adjustment of the upper ocean
mixed layer to the macroalgal farm, and they will weaken (if not disappear)
when the flow reaches an equilibrium state within the farm. The triple-
decomposed energy budgets reveal that the energy of the secondary flow is
transferred from the mean flow under the action of canopy drag, while the
transient eddies feed on wave energy transferred by the Stokes drift and
energy conversion from the secondary flow.

##  1 Introduction

Macroalgae, also known as seaweeds, are an important component in temperate
marine ecosystems (Dayton, 1985; Schiel & Forster, 2015). Providing shelter,
food and protection for many species of marine living creatures, macroalgae
play a paramount role in preserving biodiversity and promoting sustainable
aquaculture production. Macroalgal forest harvesting also contributes
enormously to various applications, such as remediation of eutrophication
pollution, biofuel production, food and pharmaceutical processing, etc. The
desire to increase the productivity of aquaculture spurs the growing need for
aquafarm development in the ocean, where the canopy grows near the surface and
is supported by a floating structure (Troell et al., 2009; Stevens & Petersen,
2011). The macroalage canopy alters the surrounding flow conditions by
dampening the currents and wave motions (Rosman et al., 2007). These flow
modifications have profound implications for the nutrient uptake and
associated processes of sedimentation and recruitment (Duggins et al., 1990;
Plew, 2011b). Therefore, understanding and quantifying the diverse
hydrodynamic processes that occur in the presence of macroalgal farms is
essential in evaluating and designing optimal farm configurations, as well as
assessing their environmental impacts.

From a hydrodynamics perspective, aquatic vegetation can be classified as
submerged, emergent, or suspended based on its growth form. Submerged and
emergent vegetation are attached to the bottom floor, and occupy a fraction or
all of the water depth. The flow structures and mass transport over such
canopies have been well documented (Nepf, 2012a, b; Yan et al., 2017).
Particular attention has been given to the shear layer turbulence at the
canopy top (for submerged canopy), which prompts the generation of canopy-
scale coherent structures that dominate the momentum and scalar exchanges
between the canopy and the free flow above. Suspended canopies, such as the
macroalgal farm considered here, extend downward from the surface and occupy
the upper part of the water body (Plew et al., 2005, 2006; Stevens & Petersen,
2011). The flow and canopy interactions for this configuration remain less
explored as compared to the bottom-mounted counterpart (Stevens & Plew, 2019).

For suspended vegetation, the vertical discontinuity in drag beneath the
canopy also leads to a shear layer, which penetrates a finite distance into
the canopy and mediates the turbulent exchanges between the canopy and the
underlying flow (Plew, 2011a). Through laboratory experiments of suspended
canopies in shallow waters, Plew (2011a) concluded that the additional bottom
boundary layer (BBL) associated with the ocean floor affects the penetration
of the shear layer into the suspended canopy. Based on the measurements of
Plew (2011a), Huai et al. (2012) proposed a simple analytical model for the
vertical profile of streamwise velocity. While these studies focus on flow
over uniform canopies (i.e. essentially infinite size), where the flow has
been fully adjusted to the canopy, common aquaculture structures are of finite
size and the corresponding canopy flow displays distinct spatial distribution
patterns.

The finite dimensions and spatial arrangement of the suspended canopy lead to
flow patterns different from the fully developed scenario (Tseung et al.,
2016). According to Tseung et al. (2016) and Zhao et al. (2017), the flow over
a suspended canopy of finite size is similar to the terrestrial flow over
forest patches (Belcher et al., 2003), and it can be divided into four zones
of distinct mean flow behaviour in the downstream direction: (i) the upstream
adjustment zone, (ii) the transition zone, (iii) the fully developed zone, and
(iv) the wake zone. The distance over which the velocity profile reaches a
fully developed state is affected by canopy geometry (e.g. plant density and
stem diameter) (Rosman et al., 2010). Zhou & Venayagamoorthy (2019) examined
the effect of a circular patch of suspended canopy on the mean flow dynamics
in deep water, and found out that the patch geometry poses another impact on
the adjustment of flow pathways. In the light of these studies, we are
motivated by the water flow over an aquaculture farm of finite size in deep
ocean, and seek to explore how the ocean mixed layer (OML) evolves as it
approaches and flows over the farm under typical ocean conditions.

In the marine environment, ocean waves have a profound influence on the water
flow and the exchange of nutrients between kelp forests and ambient water. In
many studies, this effect is characterized in terms of the Stokes drift
(Gaylord et al., 2007; Rosman et al., 2007), which refers to the net motion of
fluid parcels in the direction of wave propagation that arises from the
unclosed orbital motions for finite amplitude waves (Monismith & Fong, 2004).
Rosman et al. (2007) explored the effects of giant kelp forests on ocean flows
through a field experiment at the coast of Santa Cruz, California. They
highlighted the importance of the Stokes drift in cross-shore transport within
the kelp canopy. Rosman et al. (2013) conducted experiments at a scaled
laboratory flume to examine the interaction of surface waves and currents with
kelp forests, and concluded that these interactions must be taken into account
when modeling flow and transport within kelp forests.

One of the distinct features widely observed in the upper ocean is the
presence of Langmuir circulations, which consists of counter-rotating vortices
near the ocean surface roughly aligned with the wind direction (Thorpe, 2004).
It is well accepted that the Langmuir circulations are generated by the
interaction between the wind-driven shear current and the Stokes drift
velocity induced by the surface gravity waves through the the Craik-Leibovich
(CL) type 2 instability (Craik, 1977; Leibovich, 1983). The associated ocean
flows are referred to as Langmuir turbulence (McWiliams et al., 1997), which
can be numerically modelled by adding a vortex force into the momentum
equation without the need to resolve the surface gravity waves (Skyllingstad &
Denbo, 1995; McWiliams et al., 1997; Yang et al., 2015; Chamecki et al.,
2019). The increased level of turbulence intensity promoted by Langmuir
circulations is expected to affect the supply and uptake of nutrients within
the marine ecosystem (Barton et al., 2014).

In this study, we use a fine-scale Large Eddy Simulation (LES) model to
explore the development of an OML in the presence of a large macroalgal farm
under typical current and wave regimes. The main goal of the present work is
to characterize the hydrodynamics around a macroalgal suspended farm and
advance our understanding of canopy flows in the ocean. We assume that the
ocean is deep enough so that the flow is free from the complexities of BBL.
Section 2 describes the numerical approach for modeling oceanic boundary layer
flow over a macroalgae canopy. A triple decomposition strategy is used to
separate the flow field into the contributions due to mean flow, secondary
flow resulting from the farm drag, and transient fluctuations. Section 3
describes the main characteristics of the flow field and the emergence of
persistent flow structures termed “attached Langmuir circulations”. Section 4
discusses the underlying mechanism of generation of attached Langmuir
circulations, and characterizes their spatial development. Section 5 describes
the energy conversion among the three components of the flow field.
Conclusions are drawn in section 6.

##  2 Methods

###  2.1 Mathematical model

For the past three decades, the LES technique has been widely adopted to study
turbulence in the OML. Detailed discussion of the LES framework and
assumptions underpinning its applicability can be found in the review paper by
Chamecki et al. (2019). In the present work, the dynamics of Langmuir
turbulence in the presence of a macroalgae canopy are captured using the LES
method by solving the wave-averaged equations described by McWiliams et al.
(1997). This mathematical model is built upon the original Craik-Leibovich
equations (Craik & Leibovich, 1976) with the inclusion of planetary rotation
and Stokes drift advection of scalar fields,

$\nabla\cdot{\widetilde{\boldsymbol{u}}}=0,$ (1)
$\frac{\partial{\widetilde{\boldsymbol{u}}}}{\partial{t}}+{\widetilde{\boldsymbol{u}}}\cdot\nabla{\widetilde{\boldsymbol{u}}}=-\nabla{\Pi}-f\boldsymbol{e_{z}}\times\left(\widetilde{\boldsymbol{u}}+\boldsymbol{u}_{s}-\boldsymbol{u}_{g}\right)+\boldsymbol{u}_{s}\times\widetilde{\boldsymbol{\zeta}}+\left(1-\frac{\widetilde{\rho}}{\rho_{0}}\right)g\boldsymbol{e_{z}}+\nabla\cdot{\mathsfbi{\boldsymbol{\tau}}^{d}}-\boldsymbol{F}_{D},$
(2)
$\frac{\partial{\widetilde{\rho}}}{\partial{t}}+\left(\widetilde{\boldsymbol{u}}+\boldsymbol{u}_{s}\right)\cdot\nabla\widetilde{\rho}=\nabla\cdot{\boldsymbol{\tau}_{\rho}},$
(3)

Here, the tilde indicates grid-filtered variables, $\widetilde{\rho}$ is the
filtered seawater density, $\rho_{0}$ is the reference density, $\Pi$ is the
generalized pressure, $f$ is the Coriolis frequency, $g=9.81\ \mathrm{m\
s^{-2}}$ is the gravitational acceleration, $\boldsymbol{e_{z}}$ is the unit
vector in the vertical direction, and
${\widetilde{\boldsymbol{u}}}=({\widetilde{u},\widetilde{v},\widetilde{w}})$
is the velocity vector represented in the Cartesian coordinate system
$\boldsymbol{x}=({x},{y},{z})$, with $x$, $y$, and $z$ being the downstream,
cross-stream, and vertical directions, respectively. The vertical coordinate
is defined positive upward with $z=0$ at the ocean surface. The geostrophic
current $\boldsymbol{u}_{g}=(u_{g},0,0)$ is driven by an external mean
pressure gradient force with magnitude $fu_{g}$ applied in the $y$-direction.
The canopy is treated as a source of flow resistance, and its effect is
accounted for by adding a drag force $\boldsymbol{F}_{D}$ to the momentum
equation.

In (2) and (3), ${\boldsymbol{\tau}^{d}}$ is the deviatoric part of the
subgrid-scale (SGS) stress tensor
$\boldsymbol{\tau}=\widetilde{\boldsymbol{u}}\widetilde{\boldsymbol{u}}-\widetilde{\boldsymbol{u}\boldsymbol{u}}$,
and
${\boldsymbol{\tau}_{\rho}}=\widetilde{\boldsymbol{u}}\widetilde{\rho}-\widetilde{\boldsymbol{u}\rho}$
is the SGS buoyancy flux. We assume that the changes in the seawater density
$\rho$ are caused by the varying potential temperature $\theta$, and these two
variables are linearly related by
$\rho=\rho_{0}[1-\alpha(\theta-\theta_{0})]$, where $\alpha=2\times 10^{-4}\
\mathrm{K}^{-1}$ is the thermal expansion coefficient, and $\theta_{0}$ is the
reference potential temperature at which $\rho_{0}$ is measured. The SGS
stress tensor is modeled using the Lagrangian scale-dependent dynamic
Smagorinsky SGS model (Bou-Zeid et al., 2005). Then, the SGS buoyancy flux is
parameterized using an eddy diffusivity closure with a prescribed value of SGS
Prandtl number $\Pran_{t}=0.4$. The viscous force is assumed to be negligible
for the high-Reynolds number flows considered in the present study.

The Stokes drift $\boldsymbol{u}_{s}$ induced by surface gravity waves is
imposed in the governing equations to reflect the time-averaged effects of the
wave field on the oceanic turbulence, since the surface wave motions are not
explicitly resolved in our simulations. The third term on the right-hand-side
(RHS) of (2) is the CL vortex force
$\boldsymbol{u}_{s}\times\widetilde{\boldsymbol{\zeta}}$ (here
$\widetilde{\boldsymbol{\zeta}}=\nabla\times\widetilde{\boldsymbol{u}}$ is the
vorticity field), which represents the interaction of wind-driven turbulence
and surface gravity waves. For simplicity, we only consider a steady
monochromatic wave. Assuming that the surface gravity wave propagates along
the mean wind direction (i.e. the $x$-direction), the Stoke drift velocity
reduces to $\boldsymbol{u}_{s}=(u_{s}(z),0,0)$, where $u_{s}$ is given by:

$u_{s}={U_{s}}e^{2kz}$ (4)

in which $k$ is the wavenumber and $U_{s}$ is the wave-induced Stokes drift at
the surface. Then, the vortex force
$\boldsymbol{u}_{s}\times\widetilde{\boldsymbol{\zeta}}$ reduces to
$(0,-u_{s}\widetilde{\zeta}_{z},u_{s}\widetilde{\zeta}_{y})$. Note that the
presence of the canopy can attenuate the waves and impact the Stokes drift
profile (Rosman et al., 2013). Based on the approach developed by Dalrymple et
al. (1984), we have estimated the effects of canopy drag on the surface waves
for the specific canopy and wave parameters used in this study and found only
a small attenuation of about $3\%$ in wave amplitude and $6\%$ in the
magnitude of the Stokes drift (see appendix B). These estimates are consistent
with those obtained in flume measurements by Rosman et al. (2013). For the
sake of simplicity, we neglect wave attenuation in this study.

Finally, there is evidence suggesting that surface waves can induce a mean
current in the direction of the wave propagation within aquatic canopies
(Luhar et al., 2010, 2013; Abdolahpour et al., 2017; Chen et al., 2019; van
Rooijen et al., 2020). This wave-induced current is caused mainly by the
reduction of the wave orbital velocity within the canopy, and inclusion in our
model would require explicitly resolving the surface waves. We used the
empirical results in the literature (see Abdolahpour et al., 2017; Chen et
al., 2019) to estimate the maximum magnitude of this wave-induced current for
the suspended farm simulated here, and found out that it is a reasonably small
fraction (about 20%) of the steady geostrophic current $u_{g}$ imposed in our
simulations. Thus, we expect that the overall effects of this wave-induced
current to be small, and we neglect them in adopting a wave-averaged approach.

###  2.2 Numerical representation of macroalgal farm

For the cultivation of macroalgae, the aquaculture structures being deployed
in the open ocean are varied, but common practice is to suspend seeded
materials from surface buoys and mooring structures (Charrier et al., 2018).
One possible configuration for the cultivation strategy for the macroalgal of
interest (giant kelp) is shown in figure 1a. The macroalgal farm comprises
parallel lines of seeded growing ropes with a length of $W_{\mathrm{MF}}=8$ m
coiled around a backbone (or longline). Each backbone line, with a length
$L_{\mathrm{MF}}$, is anchored at each end and connected to surface buoys (not
shown). Each macroalgae consists of 8 fronds with an average length
$h_{\mathrm{MF}}=19$m, which are assumed to be in an upright posture by virtue
of the buoyancy provided by the gas-filled floats (called pneumatocysts). The
lateral spacing between two adjacent rows of canopy elements is
$S_{\mathrm{MF}}=26$m.

The frond surface area of the cultivated macroalgae species is obtained by
conversion of vertically-resolved algal biomass generated from a macroalgal
growth model (Christina Frieder, personal communication) using allometric
relationships (Fram et al., 2008). To simplify the numerical modeling, the
frond surface area is redistributed uniformly within each canopy row in the
horizontal directions, while the spatial arrangement of the row structure is
resolved in the simulation. The fraction occupied by the macroalgae elements
has a total foliage area density (FAD) profile denoted as $a(z)$, which is
shown in figure 1b. FAD is the total (one-sided) frond surface area per unit
volume of space ($\mathrm{m^{-1}}$), without explicit differentiation among
blades, fronds, and stipes, etc. Since our main focus here is to examine the
adjustment of OML as it flows over the farm, canopy parameters such as $a(z)$,
$h_{\mathrm{MF}}$, $S_{\mathrm{MF}}$ are kept constant (the only exception
being the length $L_{\mathrm{MF}}$) and a sensitivity study to farm design is
beyond the scope of this study.

Figure 1: Schematic of the spatial morphology of the suspended macroalgae
farm: (a) spatial arrangement of the macroalgal farm; (b) frond area density
profile for each macroalgae row, $a(z)$, normalized by the canopy height
$h_{\mathrm{MF}}$.

The drag per unit mass $\boldsymbol{F}_{D}$ in (2) represents the effect of
the canopy as a momentum sink for the flow field, and it is parameterized as
(Shaw & Schumann, 1992; Pan et al., 2014),

$\boldsymbol{F}_{D}=\frac{1}{2}C_{D}a(z)\mathsfbi{P}\boldsymbol{\cdot}|\widetilde{\boldsymbol{u}}|\widetilde{\boldsymbol{u}}$
(5)

in which $C_{D}$ is the drag coefficient and $|\widetilde{\boldsymbol{u}}|$ is
the magnitude of the resolved velocity vector. For the sake of simplicity, the
tilde symbols used to denote resolved variables are omitted hereafter. The
coefficient tensor
$\mathsfbi{P}=P_{x}\boldsymbol{e_{x}}\boldsymbol{e_{x}}+P_{x}\boldsymbol{e_{y}}\boldsymbol{e_{y}}+P_{z}\boldsymbol{e_{z}}\boldsymbol{e_{z}}$
is employed here to account for the projection of total foliage area onto the
orthogonal planes with normal in each one of the Cartesian directions. Note
that the expression for $\mathsfbi{P}$ involves the dyadic products of the
standard basis vectors $\boldsymbol{e_{x}}$, $\boldsymbol{e_{y}}$, and
$\boldsymbol{e_{z}}$, so that $\mathsfbi{P}$ is also a second-order tensor.
This projection operation is commonly used for terrestrial canopies (Legg &
Powell, 1979; Aylor & Flesch, 2001; Pan et al., 2014), and the coefficients
$P_{x}$, $P_{y}$, and $P_{z}$ depend on the geometry of the canopy and thus on
the specific details of each plant species (Aylor & Flesch, 2001). In the
absence of observational data to specify these coefficients, we make the
assumption of isotropic distribution of FAD (e.g. the fraction of FAD
projected towards each direction is always the same), which corresponds to
$P_{x}=P_{y}=P_{z}=1/2$.

The drag coefficient $C_{D}$ is a key input parameter in the drag model (5)
that can affect the accuracy for the prediction of turbulence statistics
(Pinard & Wilson, 2001). Generally, $C_{D}$ is estimated from the reduced
momentum balance based on experimental measurements, where large uncertainty
exists depending on the formulations of the momentum equation being used
(Cescatti & Marcolla, 2004; Pan et al., 2016) and quality of measurements
(Pinard & Wilson, 2001; Marcolla et al., 2003). Many numerical studies of
atmospheric boundary layer flows used a height-averaged $C_{D}$ of constant
value for terretrial canopies (Shaw & Schumann, 1992; Dupont & Brunet, 2008;
Finnigan et al., 2009). For flexible canopy like the macroalgae, the canopy
elements can bend back and forth with the moving water, leading to reduced
fluid drag relative to the rigid and upright vegetation (Boller & Carrington,
2006; Luhar & Nepf, 2011). Pan et al. (2014) introduced a velocity-dependent
$C_{D}$ in their LES study to account for the reconfiguration of the flexible
cornfield in response to the surrounding flow (Vogel, 1989). However, giant
kelp elements do not bend with the flowing water in the same way as many
terrestrial plants or seagrasses do, because they possess many gas-filled
floats that can keep the fronds upward to the surface via buoyancy forces
(Koehl & Wainwright, 1977; Henderson, 2019). In our LES cases, we use the
value of $C_{D}=0.0148$ reported in the experimental study of Utter & Denny
(1996), which measured the drag coefficient on Macrocystis pyrifera fronds by
towing a single plant from a boat in a field experiment. It should be noted
that Utter & Denny (1996) modeled the canopy drag by a power law of the local
velocity with an exponent of 1.6 to account for the drag reduction resulting
from plant reconfiguration, while we assume the relationship between these two
variables to be quadratic (equation (5)).

Apart from the fluid drag force, the macroalgae plants are subjected to
elastic and buoyant forces, both of which act to resist bending. The subtle
balance among these forces determines the posture of macroalgae elements
(Luhar & Nepf, 2011; Henderson, 2019). Estimates given in appendix A show that
kelp stipes remain approximately upright in the flow, except for an
oscillatory motion with amplitude comparable to the wave orbital displacement.
In fact, this assumption is implicit in the parametric model for the drag
force (5): the wave orbital velocity is not included in the drag calculation,
implying that the macroalgae plants oscillate with the wave orbital velocity
(note that this assumption is consistent with the idea that the macroalgae
canopy does not impact the waves).

###  2.3 Numerical scheme

The present LES framework employs a Cartesian grid using a vertically
staggered arrangement, with the horizontal velocity components, pressure and
potential temperature $(u,v,p,\theta)$ defined at the cell center, while the
vertical velocity component $(w)$ is stored at the cell face. Spatial
derivatives in the horizontal directions are treated with pseudo-spectral
differentiation, while the derivatives in the vertical direction are
discretized using a second-order centered-difference scheme. Aliasing errors
associated with the non-linear terms are removed via padding based on the 3/2
rule. Time advancement is performed using the fully explicit second-order
accurate Adams-Bashforth scheme. The numerical code has been validated against
the LES study of McWiliams et al. (1997) for Langmuir turbulence in deep ocean
by Yang et al. (2015).

Figure 2: Sketch of the LES computational model for Langmuir turbulence with
the presence of suspended macroalgae farm in deep ocean: (a) side view and (b)
plan view. A fringe region of length $L_{\mathrm{fr}}$ towards the end of the
domain is used to force the velocity and potential temperature back to the
inflow, so periodic conditions are satisfied in the horizontal plane.

The LES domain with dimensions of $L_{x}\times L_{y}\times L_{z}$ is shown in
figure 2. For clarity, the origin of the coordinate system is defined at the
leading edge of the farm in the central longitudinal plane, and the $z$-axis
is pointing upward. The top boundary is specified as a non-deforming surface
exposed to wind shear stress. A sponge layer is imposed within the bottom 20%
of the domain to damp out fluctuations of velocity and temperature, thus
avoiding the reflection of the internal gravity waves.

The backbone line is at a depth $h_{\mathrm{b}}=20$m below the surface while
the canopy height is $h_{\mathrm{MF}}=19$m, leaving a canopy-free layer at the
top 1m near the ocean surface to represent typical harvest practices. A domain
depth of $L_{z}=6h_{\mathrm{b}}$ is chosen to avoid the interference with the
bottom boundary condition as the flow is deflected below the canopy. The
cross-stream domain size $L_{y}=8S_{\mathrm{MF}}$ is tailored to encompass
$N=8$ parallel rows of macroalgae elements, the longitudinal axes of which are
aligned in the downstream direction. Periodic boundary conditions are imposed
in the horizontal directions, which will enable us to exclude the complexities
brought by the limited width of the farm. The inlet is positioned
$L_{u}=7.5h_{\mathrm{b}}$ upwind from the farm leading edge, and the outlet is
at a distance $L_{d}=12.5h_{\mathrm{b}}$ downstream of the farm trailing edge.
Thus, the domain size in the downstream direction is
$L_{x}=L_{\mathrm{MF}}+20h_{\mathrm{b}}$.

A fringe region of length $L_{\mathrm{fr}}=5h_{\mathrm{b}}$ is used at the end
of the domain (see figure 2) to enable simulations of spatially evolving
boundary layer flows in a periodic domain using pseudo-spectral numerics
(Stevens et al., 2014). Specifically, the inflow turbulence profile at the
inlet of the domain is provided by a precursor simulation carried out with
identical conditions in the absence of the farm. After the precursor
simulation reaches a fully developed turbulence regime, a region of length
$L_{\mathrm{fr}}$ is duplicated from the precursor simulation on the fringe
region of the actual simulation at the end of every time step. Then, any
variable $\phi$ (i.e. velocity and potential temperature) in the fringe region
is determined as a weighted average of fields in the precursor and actual
simulations (also see Stevens et al., 2014),

$\phi(x,y,z,t)=f(x)\cdot\phi_{\mathrm{pre}}(x,y,z,t)+\left[1-f(x)\right]\cdot\phi_{\mathrm{act}}(x,y,z,t),$
(6)

in which $\phi_{\mathrm{pre}}$ and $\phi_{\mathrm{act}}$ are, respectively,
the field in the precursor and actual domains, and $f(x)$ is the weighting
function expressed as,

$f(x)=\left\\{\begin{array}[]{ll}\frac{1}{2}\left[1-\cos{\left(\upi\frac{x-x_{s}}{x_{e}-x_{s}}\right)}\right],&x_{s}\leq
x\leq x_{e}\\\\[2.0pt] 1,&x>x_{e}.\end{array}\right.$ (7)

Here, $x$ represents the downstream position,
$x_{s}=L_{x}-L_{u}-L_{\mathrm{fr}}$ is the starting point of the fringe
region, $x_{e}=L_{x}-L_{u}-\frac{1}{4}L_{\mathrm{fr}}$ is the position beyond
which $\phi=\phi_{\mathrm{pre}}$. The length of the fringe region must be
large enough to enable a smooth transition of the field $\phi$ from the farm
wake flow to the inflow condition. To avoid any possible upstream influence
from the fringe region, only solutions up to $x=x_{s}-3h_{\mathrm{b}}$ are
analyzed.

###  2.4 Simulation parameters

Our major goal is to report new flow features that develop around suspended
aquafarms under realistic oceanic conditions. Therefore, instead of exploring
the vast parameter space of possible ocean states (e.g. varying degrees of
wind, waves, currents, and surface buoyancy forcing, etc.), we only focus on
one set of very typical conditions encountered in the deep ocean. The flow is
driven by two main forcings, i.e. the overlying atmospheric flow and a
geostrophic current, in a uniformly rotating environment with the Coriolis
frequency $f=1.0\times 10^{-4}\ \mathrm{s^{-1}}$ (corresponding to a latitude
of 45∘N). The simulation parameters are chosen to be the same as those used in
McWiliams et al. (1997), which serves as benchmark case in the literature on
Langmuir turbulence (Polton et al., 2008; Skitka et al., 2020). A constant
wind stress $\tau_{w}=0.37\ \mathrm{N\ m^{-2}}$ is applied at the air-sea
interface and aligned with the wave field in the downstream direction. The
corresponding wind speed at 10-m height is $U_{10}=5\ \mathrm{m\ s^{-1}}$, and
the friction velocity at the ocean surface is $u_{*}=6.1\times 10^{-3}\
\mathrm{m\ s^{-1}}$. The wave field consists of monochromatic waves with
wavelength $\lambda=60\ \mathrm{m}$ (corresponding to a wave period
$T_{w}=6.2\ \mathrm{s}$) and amplitude $a_{w}=0.8\ \mathrm{m}$, corresponding
to $U_{s}=0.068\ \mathrm{m\ s^{-1}}$. The resulting turbulent Langmuir number
$La_{t}=\sqrt{u_{*}/U_{s}}=0.3$, which is typical for wind-wave equilibrium
conditions in the open ocean (Belcher, 2012).

A geostrophic current $u_{g}=0.2\ \mathrm{m\ s^{-1}}$ in the downstream
direction is superimposed on the flow field to represent the effect of
mesoscale flow features, which are considered to behave as a constant flow on
the time and spatial scales of interest here (5 hrs and a few kilometers). The
upper mixed layer is bounded by a stably stratified layer below with a
constant temperature gradient $\mathrm{d}\theta/\mathrm{d}z=0.01\ \mathrm{K\
m^{-1}}$. Since surface heating or cooling would add another layer of
complexity associated with buoyancy effects on turbulence, we assume zero
buoyancy flux at the ocean surface for the simulations considered here.

Table 1 summarizes the simulation parameters and resolution of six different
cases considered here. In the table, $N_{x}$, $N_{y}$, and $N_{z}$ are the
number of grid points in the $x$, $y$, and $z$ directions, respectively.
Simulation cases CLT/LT and CST/ST represent the modelling of Langmuir
turbulence and pure shear-driven turbulence in the presence/absence of
macroalgal farm, respectively. These four cases are carried out to evaluate
the effects of macroalgae canopy and the role of surface gravity waves on the
flow features. The shear-driven cases CST and ST are conducted in the absence
of any surface wave forcing, i.e. the wave-induced Stokes drift velocity is
zero. For a boundary layer flow within and under a suspended canopy of finite
size, whether or not the boundary layer can reach a fully developed stage
depends on the length of the canopy (Tseung et al., 2016). Thus, Langmuir
turbulence in the presence of a longer farm ($L_{\mathrm{MF}}=800$ m),
referred to as case CLTL, is performed to explore the limit of fully developed
flow. We focus mostly on the results of the CLT simulation and use CLTL only
when investigating the downstream flow development. The mesh is uniformly
distributed, with a horizontal resolution $\Delta_{h}=2$ m and vertical
resolution of $\Delta_{z}=0.5$ m. To confirm that the resolution is
sufficient, case CLTF is performed under the same setup as CLT, but with
finer-scale resolution (twice the resolution) in all three directions.

Cases LT and ST are initialized with a converged solution based on a initial
mixed layer depth (MLD) of 20 m, from which the inertial oscillations have
been removed. The turbulence is confined to the upper mixed layer and the
water column below is stably stratified. Then, cases LT and ST serve as
precursor simulations to provide time-varying turbulent inflow conditions for
cases CLT(F/L) and CST, respectively. The simulations CLT(F/L) and CST are
first carried out for 15,000 s to allow for the adjustment of the surface
boundary layer to the macroalgae canopy. After the turbulent flow has reached
a quasi-equilibrium state, the flow field is averaged over another 9,000 s to
obtain turbulence statistics. Finally we note that even though turbulence
scales with $u_{*}/La_{T}^{2/3}$ in Langmuir turbulence (Grant & Belcher,
2009), we use the surface friction velocity $u_{*}$ as the scaling velocity
throughout the paper to facilitate the comparison between Langmuir and shear-
driven turbulence.

Case | Canopy | Wave | $La_{t}$ | $L_{\mathrm{MF}}(\mathrm{m})$ | $L_{x}(\mathrm{m})\times L_{y}(\mathrm{m})\times L_{z}(\mathrm{m})$ | $N_{x}\times N_{y}\times N_{z}$
---|---|---|---|---|---|---
CLT | Yes | Yes | 0.3 | 400 | 800 $\times$ 208 $\times$ 120 | 400 $\times$ 104 $\times$ 240
CLTF | Yes | Yes | 0.3 | 400 | 800 $\times$ 208 $\times$ 120 | 800 $\times$ 208 $\times$ 480
CLTL | Yes | Yes | 0.3 | 800 | 1200 $\times$ 208 $\times$ 120 | 600 $\times$ 104 $\times$ 240
CST | Yes | No | N/A | 400 | 800 $\times$ 208 $\times$ 120 | 400 $\times$ 104 $\times$ 240
LT | No | Yes | 0.3 | N/A | 400 $\times$ 208 $\times$ 120 | 200 $\times$ 104 $\times$ 240
ST | No | No | N/A | N/A | 400 $\times$ 208 $\times$ 120 | 200 $\times$ 104 $\times$ 240
Table 1: Parameters of the LES runs

A snapshot of the vertical velocity $w/u_{*}$ on a horizontal plane at
$z=-0.25h_{\mathrm{b}}$ for case CLTF is shown in figure 3. The elongated
streaks of downward vertical velocity readily observed upstream from the farm
leading edge are signatures of Langmuir circulations. They are oriented to the
right of the wind direction (i.e. $x$-direction), and are transient structures
that are continuously generated and dissipated. As the OML flows into the
farm, however, a persistent pattern with stronger downward and upward
velocities alternating laterally is clearly seen, roughly parallel to the
canopy rows. The magnitude of $w/u_{*}$ within the farm region can be as large
as 8.0 (the colorbar has been saturated), while the typical values for
Langmuir and shear turbulence in the absence of the farm for the same ocean
conditions are 1.6 and 0.75, respectively (e.g. see McWiliams et al., 1997).
This quasi-stationary pattern of alternating upwelling and downwelling regions
indicates the existence of counter-rotating cells, hereafter referred to as
attached Langmuir circulations (as discussed below). These secondary flow
structures extend beyond the trailing edge in the farm wake zone.

Figure 3: Snapshot of the normalized vertical velocity $w/u_{*}$ on a
horizontal plane ($z=-0.25h_{\mathrm{b}}$) for case CLTF. The black dashed
rectangles represent the region occupied by macroalgae canopy. The blue and
red color indicate downwelling and upwelling regions.

###  2.5 Flow decomposition

The statistics for cases LT and ST are obtained by averaging both temporally
and horizontally, indicated by $\left\langle\overline{\ \cdot\
}\right\rangle$. Note that the time average and spatial average are indicated
by an overbar and a pair of angled brackets, respectively. The physical
quantities for CLT and CST are first averaged in the temporal dimension.
Because of the three-dimensional spatial heterogeneity of the flow, these
time-averaged statistics are subject to larger random errors than the spatial-
temporal averaging used for cases LT and ST. Thus, either a spatial or phase
averaging operation in the cross-stream $y$ direction is also used, indicated
respectively by $\left\langle\ \right\rangle_{y}$ or $\left\langle\
\right\rangle_{p}$. Given the idealized cross-stream canopy heterogeneity, the
cross-phase average defined here, different from the wave-phase average
introduced in deriving equation (2), corresponds to averaging over equivalent
positions in cross-stream phases. For any time-averaged field
$\overline{\phi}$, the cross-phase averaging can be expressed as,

$\left\langle{\overline{\phi}}\right\rangle_{p}(x,y,z)=\frac{1}{N}\sum_{n=0}^{N-1}{\overline{\phi}}(x,y+nS_{\mathrm{MF}},z),$
(8)

where $N=8$ is the number of macroalgae rows.

Hereafter, we use the cross-stream average to define the (primary) mean field
$\langle\overline{\phi}\rangle_{y}(x,z)$, and the deviations from the mean
field are decomposed into a secondary-flow component and a transient
component. Thus, instantaneous flow quantities, such as the velocity field
$\boldsymbol{u}$, can be represented by,

$\boldsymbol{u}=\overline{\boldsymbol{u}}+\boldsymbol{u}^{\prime}=\left<\overline{\boldsymbol{u}}\right>_{y}+\overline{\boldsymbol{u}}^{c}+\boldsymbol{u}^{\prime},$
(9)

Here, $\boldsymbol{u}^{\prime}$ denotes the transient fluctuation from
$\overline{\boldsymbol{u}}$, while the secondary-flow disturbance
$\overline{\boldsymbol{u}}^{c}=\overline{\boldsymbol{u}}-\left<\overline{\boldsymbol{u}}\right>_{y}$
is stationary in time and represents the lateral structure of the time-
averaged velocity field induced by the farm geometry. As the transient
fluctuation and secondary-flow disturbance are uncorrelated, the covariance
between the velocity component $u_{i}$ and any field $\phi$ can be written as,

$\left<\overline{u_{i}\phi}\right>_{y}=\left<\overline{u_{i}}\right>_{y}\left<\overline{\phi}\right>_{y}+\left<\overline{u_{i}}^{c}\overline{\phi}^{c}\right>_{y}+\left<\overline{u_{i}^{\prime}\phi^{\prime}}\right>_{y},$
(10)

The three terms on the RHS represent the of contributions from the mean flow,
the secondary-flow part, and the transient part, respectively.

Finally, in some cases we further average results in the vertical direction
(depth-averaged), from the free surface $z=0$ to a fixed depth $z=z_{t}$ with
$z_{t}=-2h_{b}$, which are then represented by

$\langle\overline{\boldsymbol{u}}\rangle_{yz}=\frac{1}{|z_{t}|}\int_{z_{t}}^{0}\langle\overline{\boldsymbol{u}}\rangle_{y}\mathrm{d}z.$
(11)

##  3 Langmuir turbulence in the presence of canopy

###  3.1 Adjustment of the mean flow

Figure 4: (a) Hodographs of the mean velocity vector
($\langle\overline{u}\rangle_{y},\ \langle\overline{v}\rangle_{y}$) in the
vertical at four different downstream positions as noted in the legend are
also included, and downstream variation of the depth-averaged mean velocity
vector ($\langle\overline{u}\rangle_{yz},\ \langle\overline{v}\rangle_{yz}$)
(black line); (b) profiles of the resolved momentum stress
$\left<\overline{u^{\prime}w^{\prime}}\right>_{y}$ at these selected
downstream locations. Circles indicate values at the surface $z/h_{b}=0$, and
asterisks indicate the canopy bottom $z/h_{b}=-1$.

The OML undergoes significant changes as it approaches and flows over the
farm. Here, we present the mean flow for case CLTL to offer a more complete
picture of the spatial development of the upper OML. Figure 4a shows the
hodographs of the mean horizontal velocity vector
($\langle\overline{u}\rangle_{y},\ \langle\overline{v}\rangle_{y}$) at four
different downstream positions. Upstream from the canopy leading edge
($x/h_{b}=-5$, purple line), the hodograph follows a typical Stokes-Ekman
spiral in Langmuir turbulence, with the cross-stream velocity pointing to the
right of the wind stress (i.e., $\langle\overline{v}\rangle_{y}<0$) and most
of the shear located near the surface (the horizontal velocity is nearly
uniform within most of the OML depth due to strong vertical mixing). As the
flow moves into the farm ($x/h_{b}=10,\ ,20\ ,30$), the hodographs become very
distorted due to the large effect of the canopy drag. Despite the very complex
behavior of the mean flow, some features are noteworthy. At $x/h_{b}=20$ (blue
line), the cross-stream component of the flow switches direction within the
OML, and at $x/h_{b}=30$ (red line), the cross-stream flow is completely
reversed (i.e., to the left of the wind direction within the entire depth of
the OML). Also included in the figure is the downstream variation of the
depth-averaged horizontal velocity vector ($\langle\overline{u}\rangle_{yz},\
\langle\overline{v}\rangle_{yz}$) (black line). Downstream from the leading
edge, we can see that the depth-averaged mean flow direction changes sign at
$x/h_{b}\approx 18$, indicating a change in the direction of cross-stream
advection within the farm.

The overall change in the direction of the cross-stream flow can be understood
based on the differences of surface and bottom boundary layers in the presence
of the rotation. In the northern hemisphere, the horizontal transport is
oriented to the right of the wind stress in surface Ekman layers, and to the
left of the main current in bottom Ekman layers (McWilliams, 2006). In the
present case, the canopy introduces a vertically distributed drag that is more
pronounced near the bottom of the farm (where the LAD and the mean velocities
are larger). Therefore, the sign of $\langle\overline{v}\rangle_{y}$ depends
critically on the relative importance of shear stresses at the top and bottom
of the farm. Specifically, if the stress near the ocean surface dominates over
the stress around the canopy bottom, then $\langle\overline{v}\rangle_{y}$ is
aligned to the right of the wind as in the wind-stress driven mixed layer
(McWiliams et al., 1997); if the stress at the canopy bottom prevails, then
the flow behaves more like a bottom boundary layer above the canopy bottom and
$\langle\overline{v}\rangle_{y}$ is directed to the left of the wind (right of
the bottom stress) (Taylor & Sarkar, 2008). Because the former scales with
$u_{*}$ and the latter with $u_{g}$, we expect the flow behavior for a fixed
canopy configuration to depend on the ratio $u_{g}/u_{*}$. Figure 4b shows the
vertical profiles of $\left<\overline{u^{\prime}w^{\prime}}\right>_{y}$ at the
selected four downstream locations. It clearly shows that the turbulence
within the farm has not reached a fully developed state in the downstream
direction, and the complexity of hodographs from figure 4a also reflects this
fact. Along the $x-$direction, the flow transitions from a surface-stress-
dominated regime to a bottom-stress-dominated flow, which explains the switch
in mean cross-stream flow direction shown in figure 4a.

Figure 5 displays the mean vertical velocity
$\langle\overline{w}\rangle_{y}/u_{*}$ along the $x-z$ plane for case CLTL.
The region occupied by the macroalgae canopy is highlighted in a dashed
rectangle. The $\langle\overline{w}\rangle_{y}/u_{*}$ exhibits a small value
near the inlet, which implies that the macroalgae farm poses a minor impact on
the inflow. As the flow approaches the macroalgae farm, the canopy drag
obstructs the fluid. The associated pressure gradient across the leading edge
decelerates the flow within a region upstream of the canopy (termed the
“impact region” in Belcher et al., 2003) and induces a downward vertical
motion under the canopy near the leading edge by continuity. Similarly, the
pressure drop across the trailing edge causes the wake flow return to its
inflow profile, leading to an upward motion into the wake of the farm.

Figure 5: The time- and cross-stream-averaged vertical velocity
$\langle\overline{w}\rangle_{y}$, nomalized by $u_{*}$, for case CLTL along
the $x-z$ plane. The black dashed rectangle represents the location where the
macroalgae is planted. The black solid line marks the mixed layer depth, which
is defined as the position where the temperature exceeds a certain percentage
of the mixed layer value.

The solid line in figure 5 illustrates the downstream variation of the mixed
layer depth (MLD), denoted as $z_{i}$. As it develops downstream, the shear
turbulence near the bottom of the macroalgae canopy gradually erodes the
stratification by entraining denser water into the upper mixed layer. Here, we
define MLD as the location at which the potential temperature first exceeds a
certain percentage of the mixed layer temperature $\theta_{\mathrm{ML}}$. Thus

$z_{i}=\\{z:\langle\theta\rangle_{y}(x,z)-\theta_{\mathrm{ML}}=\chi\theta_{\mathrm{ML}}\\}$
(12)

where $\chi$ is a predefined constant. This definition is adapted from the
potential temperature contour method in Sullivan et al. (1998). The downstream
evolution of the MLD indicates that, for the present configuration in which
the MLD is comparable to the depth of the backbone line, the shear layer at
the bottom of the farm creates a local perturbation in the depth of the OML,
which seems to recover downstream from the farm.

###  3.2 Attached Langmuir circulations

Figure 6 shows the contours of the secondary-flow part of the vertical
velocity $\langle\overline{w}^{c}\rangle_{p}/u_{*}$ for case CLT in the cross-
sections noted in the caption. In the figure, we can observe a regular pattern
of $\langle\overline{w}^{c}\rangle_{p}$ alternating between positive and
negative values along the cross-stream direction, indicating the steady
upwelling and downwelling motions induced by the presence of the canopy. This
organized pattern is the signature of pairs of steady counter-rotating
circulatory flows with axis approximately aligned in the streamwise direction.
These upweling and downweling regions extend to the bottom of the OML. We
infer that these flows are primarily driven by the wave-current interaction
since these features are not observed in the shear-driven case CST (not
shown). We refer to these flow structures as attached Langmuir circulations
because (i) their position is determined by the spatial structure of the
canopy, and (ii) their formation depends critically on the wave-induced Stokes
drift via a mechanism that resembles the CL type 2 instability, which will be
described in section 4.

Figure 6: The normalized secondary-flow part of vertical velocity
$\langle\overline{w}^{c}\rangle_{p}/u_{*}$, averaged over time and cross-
phase, for case CLT on a $x-y$ plane at $z=0.5h_{\mathrm{b}}$ (a), and $y-z$
plane (facing upstream) at $x=2.5h_{\mathrm{b}}$ (b), $x=7.5h_{\mathrm{b}}$
(c), $x=12.5h_{\mathrm{b}}$ (d), and $x=17.5h_{\mathrm{b}}$ (e). The black
solid line marks the mixed layer depth. The black dashed rectangles represent
the location where the macroalgae is planted. The velocity has been cross-
phase-averaged and remapped to the entire plane. The extreme colors of the
colorbar are saturated to highlight the spatial variation of the strength of
the cell pattern.

While the standard Langmuir circulations appear as unsteady structures that
move around in the flow (see figure 3), the attached Langmuir cells are more
steady and regularly spaced. For the present canopy configuration, the
separation between neighboring pairs of attached Langmuir cells is determined
by the lateral spacing between consecutive rows of macroalgae elements, but
test runs suggest that this could change if the distance between canopy rows
is significantly larger (not shown). As the flow moves downstream, the
strength of the canopy-induced Langmuir circulations exhibits a non-monotonic
variation. The downwelling velocity reaches its maximum value at $x\approx
7.5h_{\mathrm{b}}$ with a magnitude of approximately $8u_{*}$ (figure 6b). The
cell pattern then gradually decays until $x\approx 12.5h_{\mathrm{b}}$ (figure
6c), and recovers at a lower level further downstream towards the trailing
edge of the farm (figure 6d). The orientation of Langmuir cells can be
identified by the elongated downwelling streaks. Owing to the non-zero
component in the mean cross-stream velocity (see figure 4a), the canopy-
attached Langmuir circulations are oblique to the downstream direction. The
upwelling and downwelling bands are mildly deflected to the right of the wind
for $x/h_{\mathrm{b}}<10.0$, and then aligned somewhat to the left of the wind
for $x/h>10.0$, in agreement with the change in cross-stream velocity
discussed in the previous section. This complex pattern is discussed further
in section 4, where results for the long farm case (case CLTL) are presented.

Figure 7: The cross-stream- and depth-averaged secondary-flow part of velocity
variances for CLT (solid lines) and CST (dashed lines), together with the
results from CLTF (dash-dotted line). The vertical dotted lines mark the
leading and trailing edge of the farm.

As clearly seen in figures 3-6, the flow field has not reached a fully
developed state at the trailing edge of the farm (true for both the short and
long farms). For canopy flows, the canopy-drag length is defined as
$L_{c}=\left(\frac{1}{2}C_{D}\overline{a}\right)^{-1}$ where
$\overline{a}=W_{\mathrm{MF}}/(S_{\mathrm{MF}}h_{\mathrm{b}})\int_{-h_{\mathrm{b}}}^{0}a(z)P_{x}\
\mathrm{d}z$ is the effective FAD. This length scale neglects the vertical and
horizontal structure of the canopy, and characterizes the distance over which
the flow adjusts to the mean drag of canopy elements (Belcher et al., 2003;
Rominger & Nepf, 2011). The values reported in section 2 yield
$\overline{a}h_{b}\approx 7.0$ and $L_{c}\approx 19.2h_{\mathrm{b}}$. Note
that the short and long farms have lengths of approximately equal to $L_{c}$
and $2L_{c}$, suggesting that the upper mixed layer flow does not fully adjust
to the canopy in these two cases.

To quantify the strength of the attached Langmuir circulations, we focus on
the three components of velocity variances due to the contribution from the
secondary flow. Figure 7 shows the downstream variation of the depth-averaged
mean velocity variances for cases CLT and CST. The results from CLTF are also
included to examine the sensitivity to grid resolution. The comparison shows
that the finer resolution simulation (CLTF) yields relatively larger variances
than CLT in all three velocity components, but the overall variations observed
in CLTF conform qualitatively to those in CLT. Thus, we consider the
simulations with moderate resolution (CLT and CST, etc.) to be a good starting
point to explore Langmuir turbulence in the presence of marine plants. It is
interesting to note that $\langle\overline{u}^{c}\overline{u}^{c}\rangle_{yz}$
shows negligible differences within the farm between CLT and CST, suggesting
that the canopy effect on the streamwise velocity component of the secondary
flow is not impacted by the surface waves. This also indicates that
$\langle\overline{u}^{c}\overline{u}^{c}\rangle_{yz}$ is dominated by the
lateral variation in mean velocity due to the spatially varying drag. For case
CST, the magnitudes of $\langle\overline{v}^{c}\overline{v}^{c}\rangle_{yz}$
and $\langle\overline{w}^{c}\overline{w}^{c}\rangle_{yz}$ within the canopy
exceed their upstream levels by roughly an order of magnitude, suggesting that
the presence of canopy rows leads to some secondary circulations driven by
adjustment to the canopy drag, which may also be impacted by spatial variation
in the turbulent stresses (i.e. Prandtl’s secondary flow of the second kind)
(Bradshaw, 1987). In the simulation with the Stokes drift (case CLT), however,
$\langle\overline{v}^{c}\overline{v}^{c}\rangle_{yz}$ and
$\langle\overline{w}^{c}\overline{w}^{c}\rangle_{yz}$ are about two orders of
magnitude greater than that in the Stokesless simulation (case CST). The
downstream enhancement and reduction of
$\langle\overline{v}^{c}\overline{v}^{c}\rangle_{yz}$ and
$\langle\overline{w}^{c}\overline{w}^{c}\rangle_{yz}$ within the canopy for
case CLT are consistent with the pattern of the vertical velocity in figure 6.
Therefore, we conclude that, for the present configuration, the presence of
Stokes drift is a key factor enabling the mean streamwise flow structure
induced by the farm drag to develop into strong secondary circulations. As
discussed above, these eddies are roughly two-dimensional with centerline
approximately aligned in the downstream direction, justifying the nomenclature
“attached Langmuir circulations”. Based on these results, hereafter we
interpret the streamwise component of the secondary flow as a product of the
spatial structure of the canopy drag, and the crosswise and vertical
components of the secondary flow in simulations with Stokes drift as attached
Langmuir circulations.

###  3.3 Langmuir turbulence intensity

Figure 8: The transient part of the vertical velocity variance
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ for case CLT
in the $x-y$ plane at $z=0.5h_{\mathrm{b}}$ (a), and $y-z$ plane at
$x=2.5h_{\mathrm{b}}$ (b), $x=7.5h_{\mathrm{b}}$ (c), $x=12.5h_{\mathrm{b}}$
(d), and $x=17.5h_{\mathrm{b}}$ (e).

Langmuir turbulence intensity is often characterized by large vertical
velocity variance. Our interest is centered on how the macroalgae farm alters
the spatial evolution of turbulence levels and associated turbulent mixing
efficiency. In the figure 8, we plot the time- and cross-phase-averaged
vertical velocity variance due to transient eddies
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ for case CLT.
Similar to that in standard Langmuir turbulence, the vertical intensity
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ peaks at a
subsurface level, even in the presence of a shear layer near the surface due
to canopy discontinuity (top 1m). In the nearfield downstream from the leading
edge ($0<x/h_{\mathrm{b}}<4$),
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ is decreased
within the canopy and increased near the canopy bottom (figure 8a and b). This
is because the canopy drag dampens the vertical kinetic energy within the
canopy, but the shear layer at the canopy bottom can inject additional energy
from the mean flow via shear production (see section 5). Further downstream,
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ first
increases, with the maximum value occurring at $9<x/h_{\mathrm{b}}<11$ (figure
8a), and then decreases towards the trailing edge. The energetics of the upper
mixed layer, which will be covered in section 5, suggest that the enhancement
and reduction of
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ are mainly
determined by two processes: (i) the energy exchanges with the attached
Langmuir circulations and (ii) the shear production associated with the
lateral/vertical shear in streamwise velocity caused by the canopy structure.
In the downstream cross-section (figure 8c-d), a clear pattern emerges with
increased $\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ at
the bottom and outer edge of the canopy rows and reduced intensity in the
lower half of the canopy row where the leaf area density is high (figure 1b).

Figure 9 shows the comparison of the RMS of the transient vertical velocity
fluctuation $\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}/u_{*}$
between Langmuir (case CLT, upper panel) and shear turbulence (case CST, lower
panel). Upstream from the leading edge,
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}$ from case CLT is
about twice as large as that from case CST. This is because Langmuir
turbulence yields significantly higher vertical velocity intensity compared to
the pure shear-driven turbulence scenario (McWiliams et al., 1997; D’Asaro,
2001; Harcourt & D’Asaro, 2008). In the absence of surface wave forcing
(figure 9b), the contour of
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}/u_{*}$ is similar to
what is expected for open-channel flow over a suspended canopy (see figure 16e
in Tseung et al., 2016). The shear layer at the canopy bottom grows
continually downstream and penetrates upward into the canopy, leading to the
augmentation of $\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}$
within the growing shear layer. Towards the end of the farm, the shear layer
penetrates over the entire canopy depth, a phenomenon that usually occurs for
sparse canopies (Nepf, 2012a). Interestingly, in the simulation that includes
the wave-induced Stokes drift (figure 9a), the shear layer turbulence seems to
merge with Langmuir turbulence at around $x/h_{\mathrm{b}}\approx 4$, and the
turbulence levels near the ocean surface are further enhanced within the
canopy (for $6<x/h_{\mathrm{b}}<12$) as compared to the Stokesless counterpart
(figure 9b). This can be attributed to the presence of attached Langmuir
circulations described above in section 3.2. This difference between the two
cases also confirms that the enhancement of
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ in figure 8 is
due to the turbulence modulation by the attached Langmuir circulations.

Figure 9: The transient part of the vertical velocity standard deviation
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}/u_{*}$ for CLT (a)
and CST (b) in the $x-z$ plane. The black solid line marks the mixed layer
depth.

Since transient eddies and attached Langmuir circulations coexist as the fluid
impinges upon and flows over the farm (figure 3), it is desirable to compare
the energy associated with transient eddies from that of attached Langmuir
circulations. In figure 10, we plot the vertical velocity variances due to the
contribution from the transient eddies and attached Langmuir circulations as
noted in the caption. Again, only some minor differences exist between CLT and
CLTF within the farm region, building confidence on the use of the coarser
simulations to analyze the flow. To evaluate if the flow has fully adjusted to
the canopy towards the end of the farm in cases CLT, the results from CLTL are
also shown. The discrepancies between cases CLT and CLTL (black and red lines)
are mainly located near the end of the farm in CLT ($x/h_{\mathrm{b}}=20$) due
to the trailing edge effect. As the farm extends further downstream (case
CLTL, $L_{\mathrm{MF}}=40h_{\mathrm{b}}$),
$\langle\overline{w}^{c}\overline{w}^{c}\rangle^{\mathrm{CLT}}_{yz}$ does not
become uniform but still evolves in the streamwise direction within the farm
(black solid line). It is observed that the attached Langmuir circulations
gradually attenuate in strength from $x/h_{\mathrm{b}}\approx 20$ and
eventually fade away at $x/h_{\mathrm{b}}\approx 32$ (black solid line). This
suggests that their existence is a result of flow adjustment to the suspended
farm of finite size rather than a fully developed state. While the attached
Langmuir circulations disappear, the vertical velocity variance of transient
eddies for case CLTL
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{\mathrm{CLTL}}_{yz}$ is
increasing from $x/h_{\mathrm{b}}\approx 30$ towards the end of the farm
(black dashed line). The enhanced
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{\mathrm{CLTL}}_{yz}$ of
transient eddies is mainly attributed to the canopy shear in the horizontal
direction, which no longer assists the generation of attached Langmuir
circulations as the flow has reached an equilibrium state. Except in the
nearfield downstream of the leading edge,
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{\mathrm{CLT}}_{yz}$ is much
larger than
$\langle\overline{w}^{c}\overline{w}^{c}\rangle^{\mathrm{CLT}}_{yz}$
throughout the remaining part of the farm.

Figure 10: Downstream variations of
$\langle\overline{w^{\prime}w^{\prime}}\rangle_{yz}/u_{*}^{2}$ (dashed line)
and $\langle{\overline{w}^{c}\overline{w}^{c}}\rangle_{yz}/u_{*}^{2}$ (solid
line) for CLT (red), CST (blue), CLTF (green), and CLTL (black). The end of
the farm is located at $x/h_{\mathrm{b}}=20$ for CLT/CST/CLTF, and
$x/h_{\mathrm{b}}=40$ for CLTL. The grey dotted lines mark the values of the
depth-averaged vertical velocity variance for normal Langmuir turbulence (case
LT, upper line) and pure shear-driven turbulence (case ST, lower line).

###  3.4 Comparison against standard Langmuir circulations

To compare the attached Langmuir cells with the traditional Langmuir cells
that appear in the absence of the farm, we employ a conditional sampling
approach for the LES solutions to educe the coherent structure of both fields
(also see McWiliams et al., 1997; Kukulka et al., 2010; Van Roekel et al.,
2012; Shrestha & Anderson, 2019). Based on the preconception of the form of
cell structure, we identify the Langmuir cells by searching for the strong
downwelling motion. The conditioning event $\mathscr{E}$ is defined as all
$(x_{s},y_{s},t)$ instances that satisfy
$w(x_{s},y_{s},z_{*},t)\leq-{\sigma_{w}}\big{|}_{max}$, where $\sigma_{w}$ is
the RMS of transient vertical velocity and $z_{*}$ is the depth at which
$\sigma_{w}$ attains its maximum value, denoted as
${\sigma_{w}}\big{|}_{max}$. The ordered pair $(x_{s},y_{s})$ represents a set
of grid points in the horizontal space. For case LT,
$\sigma_{w}=\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}$ and
$(x_{s},y_{s})$ enumerates the entire horizontal domain; while for case CLT,
$\sigma_{w}=\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}$ is a
function of $x_{s}$, and $(x_{s},y_{s})$ only contains grid points at the
center of the canopy spacing along the $x$-direction. Thus, the conditional
average for any quantity, denoted as $\widehat{\phi}$, is obtained with,

$\widehat{\phi}(x_{s},y_{s},x^{\prime},y^{\prime},z,t)=\left\langle\phi(x_{s}+x^{\prime},y_{s}+y^{\prime},z,t)\big{|}\mathscr{E}\right\rangle,\
$ (13)

It should be noted that $(x,y)$ is the absolute coordinate in the horizontal
plane based on the Cartesian system defined in figure 2, while $(x_{s},y_{s})$
denotes the reference point with $(x^{\prime},y^{\prime})$ being the distance
from $(x_{s},y_{s})$ in the horizontal direction. Only when the flow is
horizontally inhomogeneous should $(x_{s},y_{s})$ be equal to $(x,y)$. To
reduce the sampling error, the sampled flow field for case CLT is then further
smoothed by moving average with window size in the streamwise direction given
by $x_{s}-h_{\mathrm{b}}/2<x<x_{s}+h_{\mathrm{b}}/2$.

Figure 11: Contour plots of the conditional-averaged transient vertical
velocity $\widehat{w^{\prime}}/u_{*}$ in $y-z$ planes for case LT (a), and
case CLT at different downstream locations (b) $x_{s}=2.5h_{\mathrm{b}}$, (c)
$x_{s}=7.5h_{\mathrm{b}}$, (d) $x_{s}=12.5h_{\mathrm{b}}$, and (e)
$x_{s}=17.5h_{\mathrm{b}}$. The black solid line marks the mixed layer depth.

Figure 11 shows the contour plots of $\widehat{w}/u_{*}$ in $y^{\prime}-z$
planes for cases LT and CLT as noted in the caption. Note that the mean
vertical velocity $\langle\overline{w}\rangle_{y}$ has been removed for the
case CLT before conditional averaging operations to better compare the
distinct attached Langmuir circulations against standard Langmuir circulations
(e.g. $\langle\overline{w}\rangle_{y}$ is identically zero for LT but not for
CLT). In both cases (LT and CLT), the Langmuir cells extend down to the bottom
of the OML. The Langmuir cell pattern for LT (figure 11a) appears asymmetrical
about the longitudinal plane because of the Ekman shear. The row spacing
happens to be very close in width to the natural lateral size of the
downwelling region in standard Langmuir circulations, and this may be related
to the geometric characteristics of the attached Langmuir cells presented
here. This canopy row spacing also plays a role in determining the separation
between neighboring attached circulations as described above, and the effects
of varying row spacing should be explored in the future. The downwelling
velocity is greater than the upwelling velocity for both cases, but the
upwelling motions increase by an order of magnitude in the presence of the
canopy. This is partly caused by the fact that the obstruction of farm rows
constrains the lateral extension of upwelling regions compared to standard
Langmuir turbulence regime, producing stronger upwelling to conserve mass.

##  4 Mechanism for attached Langmuir circulations

The standard Langmuir cells in a horizontally uniform OML (e.g. case LT) are
generated through the CL2 instability, which is triggered by the wave-induced
Stokes drift velocity acting upon a cross-stream perturbation in an otherwise
horizontally uniform current (Craik, 1977; Leibovich, 1983; Suzuki & Fox-
Kemper, 2016). The instability arises from the torques produced by the
variations of vortex force
$\boldsymbol{u}_{s}\times\widetilde{\boldsymbol{\zeta}}$ that appears in (2),
which leads to overturning cellular motions with downstream vorticity
(Leibovich, 1977, 1983). This flow pattern drives the well-known Langmuir
circulations that are transient in nature in the sense that they can survive
for long periods of time but they also occasionally merge and disappear
(McWiliams et al., 1997).

Figure 12: Sketch illustrating the mechanism for attached Langmuir
circulations generated due to the presence of a farm in the upper ocean. The
cross-varying current excited by the farm is rotated by the Stoke drift,
producing the attached Langmuir circulations (black solid curves) that persist
across the farm.

In the presence of a suspended farm with row structure (cases CLT and CLTL),
the canopy drag acts within the fraction of volume occupied by the canopy
elements, thus decelerating the fluid within the farm rows and accelerating
the fluid in the spacing between rows due to continuity (figure 12). The
cross-stream variation of the current produced by the farm generates a
persistent vertical vorticity field $\zeta_{z}$ that interacts with the wave-
induced Stokes drift in a way similar to the CL2 instability. Specifically,
the vertical component of vorticity $\zeta_{z}$ associated with the cross-
stream anomaly introduces a cross-stream vortex force $-u_{s}\zeta_{z}$ that
carries fluid parcels towards the planes of local maximum $u$ where fluid
sinks due to continuity (Leibovich, 1983; Thorpe, 2004). Because the
horizontal shear is persistent within the farm and the Stoke drift associated
with the waves is horizontally uniform, such interaction gives rise to the
formation of attached Langmuir circulations that are stationary and stable
within the farm. This leads to downwelling regions in the high velocity
regions between the canopy rows and upwelling regions within the rows of
canopy elements. A schematic diagram illustrating the generation of such
circulations is shown in figure 12. The black closed curves provide an
illustration of the swirling streamlines in the plane perpendicular to the
canopy axis.

Figures 13a and b show the cross-stream and vertical components of vortex
force, i.e. $-u_{s}\langle\overline{\zeta_{z}}\rangle_{p}$ and
$u_{s}\langle\overline{\zeta_{y}}\rangle_{p}$ respectively, in the $x-y$ plane
at $z=-0.5h_{\mathrm{b}}$ for the CLTL case. In terms of magnitude, the cross-
stream component $-u_{s}\langle\overline{\zeta_{z}}\rangle_{p}$ dominates over
the vertical component $u_{s}\langle\overline{\zeta_{y}}\rangle_{p}$ down to
about $x/h_{\mathrm{b}}\leq 30$, while they are both negligibly small towards
the end of the longer farm. Consistent with the pattern of the coherent part
of vertical velocity $\langle\overline{w}^{c}\rangle_{p}/u_{*}$ (figure 13c),
the vortex force alternates in sign periodically across the farm, forming
pairs of equal magnitude, oppositely directed forces in the cross-stream
direction. Very close to the leading edge ($0<x/h_{\mathrm{b}}<2$), as the
flow just enters the farm, $-u_{s}\langle\overline{\zeta_{z}}\rangle_{p}$ is
positive (pointing in the positive $y-$direction) and negative (pointing in
the negative $y-$direction) near the left and right edges of the canopy rows,
respectively. In consequence, the action of
$-u_{s}\langle\overline{\zeta_{z}}\rangle_{p}$ drives upwelling motions within
the farm rows and downwelling motions in the spacing (see figure 13c), as
illustrated in figure 12. This pattern is clearly disrupted downstream from
the leading edge, as discussed below.

Figure 13: The time- and cross-phase-averaged vortex force: (a) cross-stream
component $-u_{s}\langle\overline{\zeta_{z}}\rangle_{p}\cdot
h_{\mathrm{b}}/u_{*}^{2}$ and (b) vertical component
$u_{s}\langle\overline{\zeta_{y}}\rangle_{p}\cdot h_{\mathrm{b}}/u_{*}^{2}$;
and (c) secondary-flow part of vertical velocity
$\langle\overline{w}^{c}\rangle_{p}/u_{*}$ for case CLTL at
$z=-0.5h_{\mathrm{b}}$

.

To further characterize the flow structure associated with the attached
Langmuir circulation, we look at the streamwise vorticity $\zeta_{x}$. Figure
14 plots the contours of $\langle\overline{\zeta_{x}}\rangle_{p}$ for case
CLTL at several cross-sections as noted in the caption. As described above, it
is the Stokes drift rotation of vertical vorticity $\zeta_{z}$ that produces
downstream vorticity $\zeta_{x}$ of alternating signs in the cross-stream
direction. Although the heterogeneous canopy in the absence of the Stokes
drift (case CST) also generates turbulence-driven secondary flows (because of
spatial variability of the turbulent stresses), it fails to yield any regular
patterns in the streamwise vorticity as those shown in Figure 14 (not shown).

To better visualize the overturning circulations, we plot streamlines on $y-z$
cross-sections in figure 14b-d (Akselsen & Ellingsen, 2019, 2020). We
determine the streamlines as isolines of the non-divergent two-dimensional
streamfunction $\psi$ computed from

$\frac{\partial^{2}\psi}{\partial y^{2}}+\frac{\partial^{2}\psi}{\partial
z^{2}}=-\zeta_{x},$ (14)

The streamlines in figure 14 portray pairs of counter-rotating vortices, with
the axes aligned to the right of the wind for $0<x/h_{\mathrm{b}}<10$ and
tilted to the left of the wind after $x/h_{\mathrm{b}}\approx 10$. Since the
attached Langmuir cells are not strictly aligned with the $x-$direction, the
use of $\langle\zeta_{x}\rangle_{p}$ only captures the largest downstream
component of the three-dimensional vortices, and thus documents weaker
overturning motions relative to the full form of coherent circulations. The
variations of $\langle\overline{\zeta_{x}}\rangle_{p}$ resemble that of the
secondary-flow part of vertical velocity in figure 6a, with the maximum
magnitude appearing in the nearfield downstream from the leading edge
($x=2.5h_{\mathrm{b}}\sim 7.5h_{\mathrm{b}}$). Towards the end of farm, the
negative downstream vortices vanish and only the weak positive vortices are
left. This is mainly because the cross-stream vortex force therewith is not
strong enough (figure 13a) to sustain a downstream counter-rotating vortex
pair.

In an idealized configuration in which the incoming mean flow is perfectly
parallel to the farm rows, we would expect an organized flow structure similar
to that shown in figure 12. However, as it is clearly seen in figures 6, 13,
and 14, the patterns that emerge from the simulation are far more complex. The
attached Langmuir cells meander in the cross-stream direction and their
amplitude changes in a non-monotonic way as a function of distance from the
leading edge of the farm. These departures from the idealized scenario are
mostly caused by the cross-stream advection, as seen by the superposition of
horizontal velocity vectors onto the streamwise vorticity in figure 14a. In
particular, the shift in cross-stream velocity from negative to positive
around $x/h_{\mathrm{b}}\approx 15$ discussed in section 3.2 produces a
similar change in the effect of advection, causing the upwelling motions to be
displaced to the right of the farm row in the region near the leading edge
(i.e., up to $x/h_{\mathrm{b}}\approx 10$) and to the left of the row for
$x/h_{\mathrm{b}}>18$.

The same pattern observed in the upwelling/downwelling regions is clearly seen
in the streamwise vorticity, as the two quantities are tied together by the
overturning structure of the flow. However, the advection of the vertical and
cross-stream components of vorticity is less effective, as clearly seen in the
patterns of the vortex force (which reflect the patterns of
$\langle\overline{\zeta_{z}}\rangle$ and
$\langle\overline{\zeta_{y}}\rangle$). This is mostly because the canopy drag
continues to generate lateral shear at the canopy edges, strongly influencing
the position of $\langle\overline{\zeta_{z}}\rangle$ and
$\langle\overline{\zeta_{y}}\rangle$. As a consequence, in the region between
$10<x/h_{\mathrm{b}}<15$, the upwelling/downwelling branches of the attached
Langmuir cells no longer coincide with the divergence/convergence of the
cross-stream vortex force (compare figures 13a and 13c), leading to the
weakening of the attached Langmuir cells around $x/h_{\mathrm{b}}=12$ followed
by a restrengthening at the more favorable position with the upwelling within
the canopy row. This process appears mostly as an abrupt left shift of the
flow structure at $x/h_{\mathrm{b}}\approx 12$. Towards the end of the farm,
$-u_{s}\langle\overline{\zeta_{z}}\rangle_{p}$ is significantly reduced, and
is no longer capable of driving clear attached Langmuir circulations (see
figures 13c and 14), which is also consistent with the decay of the vertical
variance for the secondary-flow component of the flow seen in figure 10.

Figure 14: The time- and cross-phase-averaged downstream vorticity
$\langle\overline{\zeta_{x}}\rangle_{p}h_{\mathrm{b}}/u_{*}$ with overlaid
horizontal velocity vector (a scale factor of 1/5 is applied to
$\langle\overline{u}\rangle_{p}$ for better visualization) for case CLTL at
$z=-0.5h_{\mathrm{b}}$ (a), and in the $y-z$ plane at four different
downstream locations (b) $x=2.5h_{\mathrm{b}}$, (c) $x=7.5h_{\mathrm{b}}$, (d)
$x=12.5h_{\mathrm{b}}$, and (e) $x=17.5h_{\mathrm{b}}$, overlying the two-
dimensional streamfunction $\langle\overline{\psi}\rangle_{p}$ (grey lines)
computed from $\zeta_{x}$.

##  5 Mixed layer energetics

In this section, we examine the budget of the kinetic energy in the mixed
layer, which will reveal the energy source for the secondary flow in our LES
solutions. Following the decomposition strategy described in section 2.5, the
total kinetic energy ($K=\left<\overline{u_{i}u_{i}}\right>_{y}/2$) is
composed of contributions due to the mean flow, secondary flow, and transient
eddies as,

$K=\frac{1}{2}\left<\overline{u_{i}u_{i}}\right>_{y}=\underbrace{\vphantom{\left<\overline{u_{i}^{\prime}u_{i}^{\prime}}\right>_{y}}\frac{1}{2}\left<\overline{u_{i}}\right>_{y}\left<\overline{u_{i}}\right>_{y}}_{K_{\mathrm{M}}}+\underbrace{\vphantom{\left<\overline{u_{i}^{\prime}u_{i}^{\prime}}\right>_{y}}\frac{1}{2}\left<\overline{u_{i}}^{c}\overline{u_{i}}^{c}\right>_{y}}_{K_{\mathrm{SE}}}+\underbrace{\frac{1}{2}\left<\overline{u_{i}^{\prime}u_{i}^{\prime}}\right>_{y}}_{K_{\mathrm{TE}}}$
(15)

Here, $K_{\mathrm{M}}$ represents the mean kinetic energy, $K_{\mathrm{SE}}$
is the kinetic energy of the secondary mean flow (which includes lateral
variations in the flow produced by the spatial structure of the farm and the
attached Langmuir circulations), and $K_{\mathrm{TE}}$ is the turbulent
kinetic energy. By manipulating the governing equations (1) and (2), the
transport equations for $K_{\mathrm{M}}$, $K_{\mathrm{SE}}$ and
$K_{\mathrm{TE}}$ can be obtained as follows,

$\displaystyle\frac{\mathrm{D}{K_{\mathrm{M}}}}{\mathrm{D}{t}}$
$\displaystyle=-C_{\mathrm{M-SE}}-C_{\mathrm{M-TE}}+S_{\mathrm{M}}+B_{\mathrm{M}}+\varepsilon_{\mathrm{M}}+D_{\mathrm{M}}+T_{\mathrm{M}}+R_{\mathrm{M}},$
(16a) $\displaystyle\frac{\mathrm{D}{K_{\mathrm{SE}}}}{\mathrm{D}{t}}$
$\displaystyle=C_{\mathrm{M-SE}}-C_{\mathrm{SE-
TE}}+S_{\mathrm{SE}}+B_{\mathrm{SE}}+\varepsilon_{\mathrm{SE}}+D_{\mathrm{SE}}+T_{\mathrm{SE}},$
(16b) $\displaystyle\frac{\mathrm{D}{K_{\mathrm{TE}}}}{\mathrm{D}{t}}$
$\displaystyle=C_{\mathrm{M-TE}}+C_{\mathrm{SE-
TE}}+S_{\mathrm{TE}}+B_{\mathrm{TE}}+\varepsilon_{\mathrm{TE}}+D_{\mathrm{TE}}+T_{\mathrm{SE}}.$
(16c)

in which the material derivative
$\mathrm{D}/\mathrm{D}{t}=\partial/\partial{t}+\left<\overline{u_{j}}\right>_{y}\partial/\partial{x_{j}}+u_{s}\partial/\partial{x}$.
Note that the prescribed wave and current conditions, namely
$\boldsymbol{u}_{g}=(u_{g},0,0)$ and $\boldsymbol{u}_{s}=(u_{s}(z),0,0)$, have
been invoked in deriving these equations. The first two terms on the RHS of
(16) represent the magnitude of energy conversion between $K_{\mathrm{M}}$,
$K_{\mathrm{SE}}$, and $K_{\mathrm{TE}}$ as implied in the subscripts, and are
given by,

$\left.\begin{array}[]{l}\displaystyle
C_{\mathrm{M-SE}}=-\left<\overline{u_{i}}^{c}\overline{u_{j}}^{c}\right>_{y}\frac{\partial{\left<\overline{u_{i}}\right>_{y}}}{\partial{x_{j}}},\\\\[16.0pt]
\displaystyle
C_{\mathrm{M-TE}}=-\left<\overline{u_{i}^{\prime}u_{j}^{\prime}}\right>_{y}\frac{\partial{\left<\overline{u_{i}}\right>_{y}}}{\partial{x_{j}}},\\\\[16.0pt]
\displaystyle C_{\mathrm{SE-
TE}}=-\left<\overline{u_{i}^{\prime}u_{j}^{\prime}}\frac{\partial{\overline{u_{i}}^{c}}}{\partial{x_{j}}}\right>_{y}.\end{array}\right\\}$
(17)

Note that the Einstein summations convention is used. As an example,
$C_{\mathrm{M-SE}}>0$ represents the rate of production of $K_{\mathrm{SE}}$
at the expense of $K_{\mathrm{M}}$, as this term appears as a source in the
equation for $K_{\mathrm{SE}}$ (16b) and a sink in the equation for
$K_{\mathrm{M}}$ (16a). Thus, it represents the energy transfer rate from the
mean flow to the secondary flow.

The third terms on the RHS of equations (16) are the Stokes production terms
that reflect the energy conversion between the waves and the decomposed field,

$S_{\mathrm{M}}=-\left<\overline{u}\right>_{y}\left<\overline{w}\right>_{y}\partial{u_{s}}/\partial{z},\quad
S_{\mathrm{SE}}=-\left<\overline{u}^{c}\overline{w}^{c}\right>_{y}\partial{u_{s}}/\partial{z},\quad
S_{\mathrm{TE}}=-\left<\overline{u^{\prime}w^{\prime}}\right>_{y}\partial{u_{s}}/\partial{z}$
(18)

Interestingly, the Stokes production, which only makes contribution to the
turbulent kinetic energy in a horizontally homogeneous OML, now also appears
in the budget equation of mean kinetic energy in our LES experiments because
of a non-zero and spatially evolving mean vertical velocity
$\left<\overline{w}\right>_{y}$ field. The fourth term is the buoyancy
production term,

$B_{\mathrm{M}}=\alpha
g\left<\overline{w}\right>_{y}\left(\left<\overline{\theta}\right>_{y}-\theta_{0}\right),\quad
B_{\mathrm{SE}}=\alpha
g\left<\overline{w}^{c}\overline{\theta}^{c}\right>_{y},\quad
B_{\mathrm{TE}}=\alpha g\left<\overline{w^{\prime}\theta^{\prime}}\right>_{y}$
(19)

Here, $B_{\mathrm{M}}$ represents an exchange of mean kinetic energy
$K_{\mathrm{M}}$ with the potential energy. The fifth term in (16) is the SGS
dissipation term,

$\varepsilon_{\mathrm{M}}=-\left<\overline{\tau_{ij}}\right>_{y}\partial{\left<\overline{u_{i}}\right>_{y}}/\partial{x_{j}},\quad\varepsilon_{\mathrm{SE}}=-\left<\overline{\tau_{ij}}^{c}\partial{\overline{u_{i}}^{c}}/\partial{x_{j}}\right>_{y},\quad\varepsilon_{\mathrm{TE}}=-\left<\overline{\tau_{ij}^{\prime}\partial{u_{i}^{\prime}}/\partial{x_{j}}}\right>_{y}$
(20)

In light of the energy cascade phenomenology (Pope, 2000), we expect most of
the energy dissipation occurs at the small-scale transient eddies, while the
energy loss of the large-scale mean flow and secondary flow to direct SGS
dissipation effects is negligible, i.e. $\varepsilon_{\mathrm{M}},\
\varepsilon_{\mathrm{SE}}\ll\varepsilon_{\mathrm{TE}}$. Thus, we will assume
$\varepsilon\approx\varepsilon_{\mathrm{TE}}$ in interpreting the LES
solutions, and do not partition the total dissipation $\varepsilon$ into 3
components as in (20). The sixth term in (16) is the canopy destruction term,

$D_{\mathrm{M}}=-\left<\overline{u_{i}}\right>_{y}\left<\overline{F_{D,i}}\right>_{y},\quad
D_{\mathrm{SE}}=-\left<\overline{u_{i}}^{c}\overline{F_{D,i}}^{c}\right>_{y},\quad
D_{\mathrm{TE}}=-\left<\overline{u_{i}^{\prime}F_{D,i}^{\prime}}\right>_{y}$
(21)

which represents the energy gain/loss of each component of the flow field
(i.e. mean flow, secondary flow, and transient eddies) due to the action of
canopy drag. The terms in flux form are grouped together as a transport term
in (16),

$\displaystyle\begin{split}T_{\mathrm{M}}&=\frac{\partial{}}{\partial{x_{j}}}\left[\left<\overline{u_{i}}\right>_{y}\left<\overline{\tau_{ij}}\right>_{y}+\left<\overline{u_{j}}\right>_{y}\left<\overline{u}\right>_{y}u_{s}-\left<\overline{u_{j}}\right>_{y}\left<\overline{\Pi}\right>_{y}\right.\\\
&\quad\quad\quad\quad{}\left.-\left<\overline{u_{i}}\right>_{y}\left<\overline{u_{i}}^{c}\overline{u_{j}}^{c}\right>_{y}-\left<\overline{u_{i}}\right>_{y}\left<\overline{u_{i}^{\prime}u_{j}^{\prime}}\right>_{y}\right],\\\
\end{split}$ (22a) $\displaystyle T_{\mathrm{SE}}$
$\displaystyle=\frac{\partial{}}{\partial{x_{j}}}\left<\overline{u_{i}}^{c}\overline{\tau_{ij}}^{c}-\frac{1}{2}\overline{u_{i}}^{c}\overline{u_{i}}^{c}\overline{u_{j}}^{c}-\overline{u_{j}}^{c}\overline{\Pi}^{c}-\overline{u_{i}^{\prime}u_{j}^{\prime}}\overline{u_{i}}^{c}+\overline{u_{j}}^{c}\overline{u}^{c}u_{s}\right>_{y},$
(22b) $\displaystyle T_{\mathrm{TE}}$
$\displaystyle=\frac{\partial{}}{\partial{x_{j}}}\left<\overline{u_{i}^{\prime}\tau_{ij}^{\prime}}-\overline{u_{i}^{\prime}u_{i}^{\prime}u_{j}^{\prime}}/2-\overline{u_{j}^{\prime}\Pi^{\prime}}+\overline{u_{j}^{\prime}u^{\prime}}u_{s}-\frac{1}{2}\overline{u_{j}}^{c}\overline{u_{i}^{\prime}u_{i}^{\prime}}\right>_{y}.$
(22c)

which represents the transport of kinetic energy ($K_{\mathrm{M}}$,
$K_{\mathrm{SE}}$, or $K_{\mathrm{TE}}$) through resolved momentum stresses,
SGS stresses, and pressure. The last term on the RHS of (16a) represents the
effect of Coriolis force associated with the Stokes drift and geostrophic
current,

$R_{\mathrm{M}}=f\left<\overline{v}\right>_{y}(u_{g}-u_{s})$ (23)

which transfers energy from surface waves and external larger-scale field to
the mean flow (Suzuki & Fox-Kemper, 2016).

Figure 15: Budget terms of the depth-averaged kinetic energy in the upper
surface layer for case CLTL: (a) downstream variation of the triply-decomposed
kinetic energy; and partition of conversion, Stokes production, buoyancy
production and canopy destruction for (b) $K_{\mathrm{m}}$, (c)
$K_{\mathrm{SE}}$, and (d) $K_{\mathrm{TE}}$. The terms are normalized by
$h_{\mathrm{b}}/u_{*}^{3}$.

Figure 15a shows the downstream variation of the depth-averaged kinetic energy
for the triply-decomposed field. Within the canopy region
($0<x/h_{\mathrm{b}}<40$), $K_{\mathrm{M}}$ decreases because the farm drains
the mean kinetic energy by decelerating the time-mean flow. As the OML flow
impinges upon the farm, both $K_{\mathrm{SE}}$ and $K_{\mathrm{TE}}$ increase
in the near-field downstream of the leading edge. While $K_{\mathrm{TE}}$
maintains at a high level after that, $K_{\mathrm{SE}}$ gradually decreases
towards the end of the farm. This suggests that, in the presence of a
suspended farm, the flow within the canopy region is in a highly turbulent
state but the organized secondary circulations become less intense as the
fluid moves further downstream. The downstream variations of the various
production and destruction terms in the kinetic energy budget equation (16)
for the mean flow, secondary flow, and transient eddies are depicted in figure
15b-d, respectively, using $u_{*}$ and $h_{\mathrm{b}}$ as the scaling
parameters (transport terms are not shown). To facilitate interpretation, the
curves are color coded according to the diagram depicting energy exchanges
shown in figure 16, which provides a summary of the energy budget for the
three components of the flow integrated over the entire farm.

The Stokes production $S_{\mathrm{M}}$ is the main source for $K_{\mathrm{M}}$
(figure 15b), except it is negative after $x/h_{\mathrm{b}}\approx 32$, mainly
because of the upward deflection near the trailing edge, i.e.
$\left<\overline{w}\right>_{y}>0$ that makes
$S_{\mathrm{M}}=-\left<\overline{u}\right>_{y}\left<\overline{w}\right>_{y}\partial{u_{s}}/\partial{z}<0$.
Contrary to expectations, the energy conversion term $-C_{\mathrm{M-TE}}$ is
mostly positive along the farm (green dashed line in figure 15b and d),
indicating that the transient eddies lose kinetic energy to the mean flow. The
canopy destruction term $D_{\mathrm{M}}$ is the primary sink term for
$K_{\mathrm{M}}$ as the hydrodynamic drag imparted by the farm consistently
removes the momentum from the flow (e.g.
${\partial{\left<\overline{u}\right>_{y}}}/{\partial{x}}<0$). The energy
conversion term $C_{\mathrm{M-SE}}$ (red solid line in figure 15b and c)
constitutes the secondary energy sink for $K_{\mathrm{M}}$, i.e. energy is
transferred from the mean flow to the secondary flow. This is mainly because
the leading order term of $C_{\mathrm{M-SE}}$ in (17) is
$-\left<\overline{u}^{c}\overline{u}^{c}\right>_{y}{\partial{\left<\overline{u}\right>_{y}}}/{\partial{x}}>0$.
Since the geostrophic current and Stokes drift velocity are prescribed, the
sign of Coriolis-related term $R_{\mathrm{M}}$ in (23) is directly determined
by the cross-stream velocity $\left<\overline{v}\right>_{y}$, which goes to
the right of the wind (i.e. $\left<\overline{v}\right>_{y}<0$) as in standard
Langmuir turbulence before $x/h_{\mathrm{b}}\approx 18$ and then turn to the
left of the wind (i.e. $\left<\overline{v}\right>_{y}>0$) after that (not
shown). The flow veering is largely caused by the modification of the
suspended farm on the vertical momentum transfer, given that
$f\left<\overline{v}\right>_{y}\sim\partial{\left<\overline{u^{\prime}w^{\prime}}\right>_{y}}/\partial{z}$
as yielded from a reduced form of (2).

In terms of the secondary flow, the canopy-related term $D_{\mathrm{SE}}$ is a
major source term for $K_{\mathrm{SE}}$ (black dotted line in figure 15c),
mainly because it is the spatial arrangement of the farm that leads to
persistent variations in the streamwise flow across the farm. Apart from the
energy conversion from the mean flow $C_{\mathrm{M-SE}}$, another important
source term for the secondary mean flow is the Stokes production
$S_{\mathrm{SE}}$, which is the main source of energy to the attached Langmuir
circulations. This is true everywhere except for the region
$9<x/h_{\mathrm{b}}<12$ where $S_{\mathrm{SE}}$ is negative. In this local
range, $S_{\mathrm{SE}}$ serves as an sink of $K_{\mathrm{SE}}$ and the energy
transferred from the mean flow $C_{\mathrm{M-SE}}$ is also decreasing (red
solid line), which to some extent explains the local attenuation of attached
Langmuir circulations at $x/h_{\mathrm{b}}=12.5$ (figure 6d). For
$x/h_{\mathrm{b}}>32$, $S_{\mathrm{SE}}$ is approximately zero because the
coherent vertical velocity $\overline{w}^{c}$ almost vanishes (figure 10) and
hence the momentum stress due to the secondary flow
$\left<\overline{u}^{c}\overline{w}^{c}\right>_{y}$ in (18) is negligibly
small. These three source terms ($D_{\mathrm{SE}}$, $S_{\mathrm{SE}}$, and
$C_{\mathrm{M-SE}}$) are responsible for the maintenance of secondary flow
(including the attached Langmuir circulations) in the adjustment region
downstream of the leading edge, whereas the exchange with the transient eddies
$C_{\mathrm{SE-TE}}$ constantly extracts energy from the secondary flow to
support the turbulence level (purple dashed line in figure 15c).

As shown in figure 15d, the transient eddies feed on wave energy transferred
by the Stokes drift shear (blue dash-dotted line) and energy conversion from
the secondary flow. The transient eddies lose energy mostly via three
processes: (i) energy transfer to the mean flow; (ii) energy removal due to
the canopy drag; and (iii) energy dissipation at the small scales (represented
by the SGS dissipation). As Langmuir turbulence in the presence of canopy
features strong shear layers and wave forcing, and we assume no incoming or
outgoing buoyancy flux at the surface, the buoyant production terms for the
secondary flow and transient eddies ($B_{\mathrm{SE}}$ and $B_{\mathrm{TE}}$)
are negligibly small in comparison.

Figure 16: Schematic diagram of the depth-averaged energy budget for the mean
flow, secondary flow, and transient eddies. The arrow-lines represent the
transfer of energy integrated over the entire farm length, with the direction
of net energy flow indicated by heavy arrowheads. The number alongside each
arrow-line is the farm-averaged value of the corresponding term, normalized by
$h_{\mathrm{b}}/u_{*}^{3}$. Note that the transport terms are not included
here, thus the energy budget for each component is not closed.

##  6 Conclusions

In this study, a fine-scale LES model is used to explore how Langmuir
turbulence in deep ocean evolves as it flows over and through a row-structured
macroalgae farm. The ocean flow is driven by a constant wind stress and a
geostrophic current, under the influences of surface gravity waves, planetary
rotation, and stable interior stratification. The effects of Langmuir
turbulence are accounted for by adding the CL vortex force to the momentum
equation without explicitly resolving the surface waves. For the case studied
here, the drag force at the bottom of the farm becomes dominant over the wind
forcing at the surface with increasing distance downstream from the leading
edge. As a result, the mean horizontal flow switches from the canonical
surface forced Ekman layer to a regime that resembles a bottom Ekman layer.
This transition is evident in the change of direction of the mean current
perpendicular to the wind.

Following a triple-decomposition technique, the turbulent transport is divided
into contributions from the mean flow, secondary flow, and transient
fluctuations. We find out that the row structure of the farm causes the cross-
stream variation of the current that ultimately leads to the formation of
coherent circulations via a mechanism similar to CL2 instability theory.
Specifically, the vertical vortex lines associated with this cross-varying
current are tilted by the Stokes drift, driving the formation of downstream
vortices that are stationary in time, phase locked in space, and periodically
alternating in sign across the lateral direction. Thus, we also refer to these
coherent structures as attached Langmuir circulations.

The attached Langmuir circulations are unique to the upper OML in the presence
of aquacultural farms (or other distributed roughness elements) since the
cross-stream variation of the current is excited by the canopy. They are
roughly oriented along the rows of canopy elements, which are aligned with the
wind direction within the present numerical framework. The vertical extent of
attached Langmuir circulations can occupy the entire OML, with the lateral
scale of the associated downwelling regions comparable to the row spacing in
the farm. The potential impact of varying farm geometry and ocean conditions
on these circulations is out of the scope here, but should be explored in the
future.

Because the associated upwelling motions are concentrated in regions occupied
by macroalgae elements, these attached Langmuir circulations are conducive to
vertical mixing and could increase nutrient availability within macroalgae
farm environments. The strength of transient eddies, characterized by cross-
stream and vertical turbulence intensity (i.e.
$\langle\overline{v^{\prime}v^{\prime}}\rangle^{1/2}_{y}$ and
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{y}$), is much larger
under the effect of Stokes drift associated with the surface waves (case CLT)
compared to the pure shear-driven scenario (case CST), which is also
consistent with previous studies in the absence of the canopy (McWiliams et
al., 1997; D’Asaro, 2001). For both cases, the suspended farm prompts a shear
layer development near the canopy bottom and deepens the mixed layer as the
flow moves downstream. For the simulation with Stokes drift (case CLT), in the
nearfield downstream of the leading edge, the canopy drag dampens the
turbulence, leading to the reduction of
$\langle\overline{w^{\prime}w^{\prime}}\rangle^{1/2}_{p}/u_{*}$ for
$0<x/h_{\mathrm{b}}<4$. Further downstream, the attached Langmuir circulations
promote strong enhancement of turbulence. This enhancement slowly fades as the
flow adjusts to the canopy and the strength in the secondary flow decays
(figure 10). The presence of the canopy leads to the formation of the attached
Langmuir circulations and to local enhancement of the turbulence. Both flow
modifications are expected to enhance vertical mixing within the OML and
possibly help the entrainment of nutrients from the pycnocline.

Analysis of kinetic energy budget shows that, as the flow moves downstream of
the canopy leading edge, the canopy drag acts as an energy sink for the mean
flow and transient fluctuations, while serving as a major source for the
kinetic energy of the secondary mean flow. If the canopy is long enough, the
secondary flow pattern vanishes when the oceanic turbulence is fully adjusted
to the macroalgal farm. Therefore, this flow feature arises from the
adjustment of the upper mixed layer to the aquafarm.

The conclusions drawn here are valid for conditions in which the effect of
Stokes drift dominates over that of wind stress and external pressure gradient
forcing (i.e. the solutions are posed in the Langmuir turbulence regime).
Despite the simplification made here (e.g. plant reconfiguration,
monochromatic waves, etc.), we are optimistic that the findings presented
above are relevant to realistic practice, and could serve as guidance for the
design of large scale macroalgae systems. Still, the attached Langmuir
circulations from our LES solutions and their potential implication on
nutrient uptake by aquaculture farms await field observations to confirm their
veracity.

From a fluid dynamics perspective, the physical flow presented here
encompasses a variety of processes (stratification, Coriolis acceleration,
wave-driven transport, and a canopy, etc). One of our main goals is to make it
clear that these flow features are important in practice, in conditions under
which macroalgal farms are deployed. As it turns out, most of the complexity
involved in our setup is essential for the attached Langmuir eddies to develop
(waves, mean current, non-uniform canopy, and downstream flow development).
There are some possible simplifications that would allow one to reduce the
parameter space and simplify the problem, bringing it to a more manageable
fundamental configuration (e.g., removing the effects of planetary rotation
and stratification). The results in this paper warrant further investigation
of a more fundamental nature in simplified conditions, which could help
reconcile a bit the complexity of the flow features we discovered with a more
traditional fluid dynamical investigation of the parameter space.

## Acknowledgements

This work is supported by the ARPA-E MARINER Program (DE-AR0000920). We thank
the three anonymous reviewers for their constructive comments which led to
improvements of the manuscript.

## Declaration of Interests

The authors report no conflict of interest.

##  Appendix A Motion of buoyant, flexible macroalgae in upper OML

The stipe reconfiguration in response to the flowing water depends on the
ocean parameters (wave amplitude, wave period, and current) and the mechanical
properties of macroalgae (stipe length, Young’s modulus, density, and
buoyancy). We decompose the upper OML flow into two parts, i.e. the steady
flow (geostrophic current) and oscillatory flow (wave orbital velocity), and
analyze the motion of buoyant, flexible macroalgae with respect to flow
components separately. For each plant, the stipe bundles are simplified to
have a circular cross-section, with length $l_{s}=20$ m, radius $r_{s}=0.1$ m
(corresponding second moment of area $I=\upi r_{s}^{4}/4$), Young’s modulus
$E=1\times 10^{7}$ Pa, and density $\rho_{\mathrm{s}}=595\ \mathrm{kg\
m^{-3}}$ (properties taken from Utter & Denny, 1996; Henderson, 2019).

In a unidirectional steady current (e.g. $u_{g}=0.2\ \mathrm{m\ s^{-1}}$), the
key parameters determining the form of macroalgae elements in sustained flow
conditions are the dimensionless Cauchy number Ca (fluid drag/elastic force)
and buoyancy number B (buoyancy force/elastic force) defined as (Luhar & Nepf,
2011),

$Ca=\frac{1}{2}\frac{\rho C_{D}r_{s}l_{s}u_{g}^{2}}{EI/l_{s}^{2}},$ (24)
$B=\frac{\left(\rho-\rho_{s}\right)g\upi r_{s}^{2}l_{s}}{EI/l_{s}^{2}},$ (25)

in which $\rho=1010\ \mathrm{kg\ m^{-3}}$ is the density of water. The ratio
$Ca/B$ measures the relative importance of fluid drag and buoyancy force. Note
that the flexibility of blades is neglected, because the blades can fold and
rotate in the water, while the stipe bundles constitute the essential part
governing the bending of macroalage elements. As such, only the fluid drag on
the stipe bundles (denominator in (24)) is considered instead of (5). From the
values of parameters given above, the resulting $Ca/B=2.3\times 10^{-3}\ll 1\
(Ca=304.5,\ B=1.3\times 10^{5})$, and the bending angle of the stipe bundles
$\xi=0.13\degree$ (estimated by equation (12) in Luhar & Nepf, 2011),
suggesting the buoyancy force dominates over the fluid drag and the stipe
bundles deform very little relative to its vertical position.

For wave-induced oscillatory flows, such as a sinusoidal wave with surface
elevation $\eta=a_{w}\cos{\left(kx-\sigma t_{w}\right)}$, Henderson (2019)
introduced a new dimensionless buoyancy number $\beta$ and stiffness number
$S$,

$\beta=\frac{\left(\rho-\rho_{s}\right)gr_{s}t_{w}}{\rho C_{D}l_{s}u_{w}},$
(26)

$S=\frac{EIt_{w}}{\rho C_{D}r_{s}l_{s}^{4}u_{w}},$ (27)

Here, $\sigma=2\upi/a_{w}$ is the angular frequency, $t_{w}=2\upi/\sqrt{gk}$
and $u_{w}=\sigma a_{w}$ are the wave period and orbital velocity scale,
respectively. Based on the monochromatic wave parameters reported in section
2.4 ($t_{w}=6.2$ s, $u_{w}=0.81\ \mathrm{m\ s^{-1}}$), the resulting
$\beta=1.06$ and $S=3.2\times 10^{-5}$. The relative magnitude of buoyancy and
elasticity scales with $\gamma=\beta/S^{1/2}=184\gg 1$, which suggests that
the elasticity plays a negligible role here. As $\beta$ is of order unity, the
stipe displacement and the wave-induced water motion are comparable, i.e. the
stipe bend with the waves. Note that our estimates of $S$ and $\beta$ are
different from those in Henderson (2019) because different values of wave
(e.g. period and amplitude) and canopy parameters (e.g. length and drag
coefficient) are used here.

##  Appendix B Deep-water wave attenuation by suspended canopies

Surface waves propagating through marine plants lose energy due to the drag
exerted by the canopy, leading to attenuation in wave heights (Dalrymple et
al., 1984). Canopy configuration (subemerged, emergent, suspended) and the
associated spatial distribution patterns exert a major impact on wave
attenuation (Chen et al., 2019). The following mathematical derivation is
based on the work of Dalrymple et al. (1984) for damping by rigid cylinders in
coastal regions, and considers suspended macroalgal farms in deep water
(described previously in the main text).

Assuming that energy dissipation is dominated by the canopy drag force, the
conservation of wave energy equation is,

$\frac{\partial\left(E_{w}c_{g}\right)}{\partial
x}=-\alpha_{D}\varepsilon_{D},$ (28)

in which $E_{w}=\frac{1}{2}\rho ga_{w}^{2}$ is the energy density per unit
area of sea surface waves, $a_{w}$ is the wave amplitude, and
$c_{g}=\frac{1}{2}\sqrt{g/k}$ is the wave group velocity. The prefactor
$\alpha_{D}$ accounts for the reduction in dissipation arising from the motion
of buoyant, flexible macroalgae, and it is a function of $\beta$ and $S$
defined in appendix A expressed as (Henderson, 2019),

$\alpha_{D}=\left[\frac{C_{S}S+C_{\beta}\beta^{2}}{1+C_{S}S+C_{\beta}\beta^{2}}\right]^{1/4},$
(29)

in which $C_{S}=1/4$ and $C_{\beta}=1/16$. For the highly flexible macroalgae
($\beta=1.06$ and $S=3.2\times 10^{-5}$ in appendix A), the value of
$\alpha_{D}$ is 0.51. $\varepsilon_{D}$ is the mean depth-integrated wave
dissipation due to canopy drag force,

$\varepsilon_{D}=\overline{\int_{-h_{b}}^{0}D_{x}u_{x}\ \mathrm{d}z},$ (30)

in which the overline denotes averaging over a complete wave period,
$u_{x}=\sigma ae^{kz}\cos{\left(kx-\sigma t\right)}$ is the horizontal
velocity due to wave orbital motions, $\sigma=\sqrt{gk}$ is the angular
frequency, and $D_{x}=\frac{1}{2}\rho C_{D}\langle
a\rangle_{y}P_{x}|u_{x}|u_{x}$ is the wave drag force on the canopy with
$\langle a\rangle_{y}$ being the lateral-averaged FAD. Substituting equation
(30) into equation (28) yields,

$\frac{1}{2}gc_{g}\frac{\partial a_{w}^{2}}{\partial x}=-Ga_{w}^{3},$ (31)

in which

$G=\frac{4}{3}\alpha_{D}C_{D}P_{x}\sigma^{3}\int_{-h_{b}}^{0}ae^{3kz}\
\mathrm{d}z,$ (32)

The solution of (31) is

$\frac{a_{w}}{a_{w}^{0}}=\left({1+\frac{Ga_{w}^{0}}{gc_{g}}x}\right)^{-1},$
(33)

in which $a_{w}^{0}$ (=0.8 m here) is the incident wave amplitude before
entering the macroalage canopy. From the values of parameters reported above,
the wave height decay over a 800-m (400-m) long farm is about 2.9% (1.4%), and
the corresponding decay in Stokes drift velocity is 5.7% (2.8%).

## References

  * Abdolahpour et al. (2017) Abdolahpour, M., Hambleton, M. & Ghisalberti, M. 2017 The wave‐driven current in coastal canopies. J. Geophys. Res. Oceans 122, 3660–3674.
  * Akselsen & Ellingsen (2019) Akselsen, A. H. & Ellingsen, S. Å. 2019 Weakly nonlinear transient waves on a shear current: ring waves and skewed langmuir rolls. Journal of Fluid Mechanics 863, 114–149.
  * Akselsen & Ellingsen (2020) Akselsen, A. H. & Ellingsen, S. Å. 2020 Langmuir-type vortices in wall-bounded flows driven by a criss-cross wavy wall topography. J. Fluid Mech. 900, A19.
  * Aylor & Flesch (2001) Aylor, D. E. & Flesch, T. K. 2001 Estimating spore release rates using a lagrangian stochastic simulation model. J. Appl. Meteorol. 40, 1196–1208.
  * Barton et al. (2014) Barton, A. D., Ward, B. A., Williams, R. G. & Follows, M. J. 2014 The impact of fine-scale turbulence on phytoplankton community structure. Limnol. Oceanogr. Fluids Environ. 4 (1), 34–49.
  * Belcher et al. (2003) Belcher, S. E., Jerram, N. & Hunt, J. C. R. 2003 Adjustment of a turbulent boundary layer to a canopy of roughness elements. J. Fluid Mech. 488, 369–398.
  * Belcher (2012) Belcher, S. E. et al. 2012 A global perspective on langmuir turbulence in the ocean surface boundary layer. Geophys. Res. Lett. 39 (18), L18605.
  * Boller & Carrington (2006) Boller, M. L. & Carrington, E. 2006 The hydrodynamic effects of shape and size change during reconfiguration of a flexible macroalga. J. Exp. Biol. 209 (10), 1894–1903.
  * Bou-Zeid et al. (2005) Bou-Zeid, E., Meneveau, C. & Parlange, M. B. 2005 A scale-dependent lagrangian dynamic model for large eddy simulation of complex turbulent flows. Phys. Fluids 17, 025105.
  * Bradshaw (1987) Bradshaw, P 1987 Turbulent secondary flows. Annu. Rev. Fluid Mech. 19 (1), 53–74.
  * Cescatti & Marcolla (2004) Cescatti, A. & Marcolla, B. 2004 Drag coefficient and turbulence intensity in conifer canopies. Agric. For. Meteorol. 121 (3), 197–206.
  * Chamecki et al. (2019) Chamecki, M., Chor, T., Yang, D. & Meneveau, C. 2019 Material transport in the ocean mixed layer: Recent developments enabled by large eddy simulations. Rev. Geophys. 57 (4), 1338–1371.
  * Charrier et al. (2018) Charrier, B., Wichard, T. & Reddy, C. R. K. 2018 Protocols for Macroalgae Research. CRC Press.
  * Chen et al. (2019) Chen, H., Liu, X. & Zou, Q. 2019 Wave-driven flow induced by suspended and submerged canopies. Adv. Water Resour. 123, 160–172.
  * Craik (1977) Craik, A. D. D. 1977 The generation of langmuir circulations by an instability mechanism. J. Fluid Mech. 81 (2), 209–223.
  * Craik & Leibovich (1976) Craik, A. D. D. & Leibovich, S. 1976 A rational model for langmuir circulations. J. Fluid Mech. 73 (3), 401–426.
  * Dalrymple et al. (1984) Dalrymple, R. A., Kirby, J. T. & Hwang, P. A. 1984 Wave diffraction due to areas of energy dissipation. J. Water. Port, Coastal, Ocean Eng. 110 (1), 67–79.
  * D’Asaro (2001) D’Asaro, E. A. 2001 Turbulent vertical kinetic energy in the ocean mixed layer. J. Phys. Oceanogr. 31 (12), 3530–3537.
  * Dayton (1985) Dayton, P. K. 1985 Ecology of kelp communities. Annu. Rev. Ecol. Syst. 16, 215–245.
  * Duggins et al. (1990) Duggins, D. O., Eckman, J. E. & Sewell, A. T. 1990 Ecology of understory kelp environments. II. Effects of kelps on recruitment of benthic invertebrates. J. Exp. Mar. Biol. Ecol. 143, 27–45.
  * Dupont & Brunet (2008) Dupont, S. & Brunet, Y. 2008 Influence of foliar density profile on canopy flow: A large-eddy simulation study. Agric. For. Meteorol. 148, 976–990.
  * Finnigan et al. (2009) Finnigan, J. J., Shaw, R. H. & Patton, E 2009 Turbulence structure above a vegetation canopy. J. Fluid Mech. 637, 387–424.
  * Fram et al. (2008) Fram, J. P., Stewart, H. L., Brzezinski, M. A., Gaylord, B., Reed, D. C., Williams, S. L. & MacIntyre, S. 2008 Physical pathways and utilization of nitrate supply to the giant kelp, macrocystis pyrifera. Limnol. Oceanogr. 53 (4), 1589–1603.
  * Gaylord et al. (2007) Gaylord, B., Rosman, J. H., Reed, D. C., Koseff, J. R.and Fram, J., MacIntyre, S., Arkema, K., McDonald, C., Brzezinski, M. A., Largier, J. L., Monismith, S. G., Raimondi, P. T. & Mardian, B. 2007 Spatial patterns of flow and their modification within and around a giant kelp forest. Limnol. Oceanogr. 52, 1838–1852.
  * Grant & Belcher (2009) Grant, A. L. M. & Belcher, S. E. 2009 Characteristics of langmuir turbulence in the ocean mixed layer. J. Phys. Oceanogr. 39 (8), 1871–1887.
  * Harcourt & D’Asaro (2008) Harcourt, R. R. & D’Asaro, E. A. 2008 Large-eddy simulation of langmuir turbulence in pure wind seas. J. Phys. Oceanogr. 38 (7), 1542–1562.
  * Henderson (2019) Henderson, S. M. 2019 Motion of buoyant, flexible aquatic vegetation under waves: Simple theoretical models and parameterization of wave dissipation. Coast. Eng. 152, 103497.
  * Huai et al. (2012) Huai, W. X., Hu, Y., Zeng, Y. H. & Han, J. 2012 Velocity distribution for open channel flows with suspended vegetation. Adv. Water Resour. 49, 56–61.
  * Koehl & Wainwright (1977) Koehl, M. A. R. & Wainwright, S. A. 1977 Mechanical adaptations of a giant kelp. Limnol. Oceanogr. 22, 1067–1071.
  * Kukulka et al. (2010) Kukulka, T., Plueddemann, A. J., Trowbridge, J. H. & Sullivan, P. P. 2010 Rapid mixed layer deepening by the combination of langmuir and shear instabilities: A case study. J. Phys. Oceanogr. 40 (11), 2381–2400.
  * Legg & Powell (1979) Legg, B. J. & Powell, F. A. 1979 Spore dispersal in a barley crop: A mathematical model. Agric. Meteorol. 20 (1), 47–67.
  * Leibovich (1977) Leibovich, S. 1977 Convective instability of stably stratified water in the ocean. J. Fluid Mech. 82 (3), 561–581.
  * Leibovich (1983) Leibovich, S. 1983 The form and dynamics of langmuir circulations. Annu. Rev. Fluid Mech. 15 (1), 391–427.
  * Luhar et al. (2010) Luhar, M., Coutu, S., Infantes, E., Fox, S. & Nepf, H. 2010 Wave-induced velocities inside a model seagrass bed. J. Geophys. Res. Oceans 115 (C12), C12005.
  * Luhar et al. (2013) Luhar, M., Infantes, E., Orfila, A., Terrados, J. & Nepf, H. 2013 Field observations of wave‐induced streaming through a submerged seagrass (posidonia oceanica) meadow. J. Geophys. Res. Oceans 118, 1955–1968.
  * Luhar & Nepf (2011) Luhar, M. & Nepf, H. M. 2011 Flow‐induced reconfiguration of buoyant and flexible aquatic vegetation. Limnol. Oceanogr. 56 (6), 2003–2017.
  * Marcolla et al. (2003) Marcolla, B., Pitacco, A. & Cescatti, A. 2003 Canopy architecture and turbulence structure in a coniferous forest. Boundary-Layer Meteorol. 108, 39–59.
  * McWiliams et al. (1997) McWiliams, J., Sullivan, P. & Moeng, C. 1997 Langmuir turbulence in the ocean. J. Fluid Mech. 334, 1–30.
  * McWilliams (2006) McWilliams, J. 2006 Fundamentals of Geophysical Fluid Dynamics. Cambridge University Press.
  * Monismith & Fong (2004) Monismith, S. G. & Fong, D. A. 2004 A note on the potential transport of scalars and organisms by surface waves. Limnol. Oceanogr. 49, 1214–1217.
  * Nepf (2012a) Nepf, H. M. 2012a Flow and transport in regions with aquatic vegetation. Annu. Rev. Fluid Mech. 44, 123–142.
  * Nepf (2012b) Nepf, H. M. 2012b Hydrodynamics of vegetated channels. J. Hydraul. Res. 50, 262–279.
  * Pan et al. (2014) Pan, Y., Chamecki, M. & Isard, S. A. 2014 Large-eddy simulation of turbulence and particle dispersion inside the canopy roughness sublayer. J. Fluid Mech. 753, 499–534.
  * Pan et al. (2016) Pan, Y., Chamecki, M. & Nepf, H. M. 2016 Estimating the instantaneous drag–wind relationship for a horizontally homogeneous canopy. Boundary-Layer Meteorol. 160, 63–82.
  * Pinard & Wilson (2001) Pinard, J. D. & Wilson, J. D. 2001 First- and second-order closure models for wind in a plant canopy. J. Appl. Meteorol. 40 (10), 1762–1768.
  * Plew (2011a) Plew, D. R. 2011a Depth-averaged drag coefficient for modeling flow through suspended canopies. J. Hydraul. Eng. 137, 234–247.
  * Plew (2011b) Plew, D. R. 2011b Shellfish farm-induced changes to tidal circulation in an embayment, and implications for seston depletion. Aquacult. Env. Interac. 1, 201–214.
  * Plew et al. (2006) Plew, D. R., Spigel, R. H., Stevens, C. L., Nokes, R. I. & Davidson, M. J. 2006 Stratified flow interactions with a suspended canopy. Environ. Fluid Mech. 6, 519–539.
  * Plew et al. (2005) Plew, D. R., Stevens, C. L., Spigel, R. H. & D., Hartstein N. 2005 Hydrodynamic implications of large offshore mussel farms. IEEE J. Ocean Eng. 30, 95–108.
  * Polton et al. (2008) Polton, J. A., Smith, J. A., MacKinnon, J. A. & Tejada-Martínez, A. E. 2008 Rapid generation of high-frequency internal waves beneath a wind and wave forced oceanic surface mixed layer. Geophys. Res. Lett. 35 (13), L13602.
  * Pope (2000) Pope, Stephen B. 2000 Turbulent Flows. Cambridge University Press.
  * Rominger & Nepf (2011) Rominger, J. T. & Nepf, H. M. 2011 Flow adjustment and interior flow associated with a rectangular porous obstruction. J. Fluid Mech. 680, 636–659.
  * van Rooijen et al. (2020) van Rooijen, A., Lowe, R., Rijnsdorp, D., Ghisalberti, M., Jacobsen, N. G. & McCall, R. 2020 Wave‐driven mean flow dynamics in submerged canopies. J. Geophys. Res. Oceans 125, e2019JC015935.
  * Rosman et al. (2013) Rosman, J. H., Denny, M. W., Zeller, R. B., Monismith, S. G. & Koseff, J. R. 2013 Interaction of waves and currents with kelp forests (macrocystis pyrifera): Insights from a dynamically scaled laboratory model. Limnol. Oceanogr. 58, 790–802.
  * Rosman et al. (2007) Rosman, J. H., Koseff, J. R., Monismith, S. G. & Grover, J. 2007 A field investigation into the effects of a kelp forest (macrocystis pyrifera) on coastal hydrodynamics and transport. J. Geophys. Res. 112, C02016.
  * Rosman et al. (2010) Rosman, Johanna H., Monismith, Stephen G., Denny, Mark W. & Koseff, Jeffrey R. 2010 Currents and turbulence within a kelp forest (macrocystis pyrifera): Insights from a dynamically scaled laboratory model. Limnol. Oceanogr. 55 (3), 1145–1158.
  * Schiel & Forster (2015) Schiel, D. R. & Forster, M. S. 2015 The biology and ecology of giant kelp forests. University of California Press.
  * Shaw & Schumann (1992) Shaw, R. H. & Schumann, U. 1992 Large-eddy simulation of turbulent flow above and within a forest. Boundary-Layer Meteorol. 61, 47–64.
  * Shrestha & Anderson (2019) Shrestha, K. & Anderson, W. 2019 Coastal langmuir circulations induce phase-locked modulation of bathymetric stress. Environ. Fluid Mech. 20, 873–884.
  * Skitka et al. (2020) Skitka, J., Marston, J. B. & Fox-Kemper, B. 2020 Reduced-Order Quasilinear Model of Ocean Boundary-Layer Turbulence. J. Phys. Oceanogr. 50 (3), 537–558.
  * Skyllingstad & Denbo (1995) Skyllingstad, E. D. & Denbo, D. W. 1995 An ocean large‐eddy simulation of langmuir circulations and convection in the surface mixed layer. J. Geophys. Res. Oceans 100, 8501–8522.
  * Stevens & Plew (2019) Stevens, C. & Plew, D. 2019 Bridging the separation between studies of the biophysics of natural and built marine canopies. Front. Mar. Sci. 6, 217.
  * Stevens & Petersen (2011) Stevens, C. L. & Petersen, J. K. 2011 Turbulent, stratified flow through a suspended shellfish canopy: implications for mussel farm design. Aquacult. Env. Interac. 36, 87–104.
  * Stevens et al. (2014) Stevens, R. J. A. M., Graham, J. & Meneveau, C. 2014 A concurrent precursor inflow method for large eddy simulations and applications to finite length wind farms. Renew. Energy 68, 46 – 50.
  * Sullivan et al. (1998) Sullivan, P. P., Moeng, C. H., Stevens, B., Lenschow, D. H. & Mayor, S. D. 1998 Structure of the entrainment zone capping the convective atmospheric boundary layer. J. Atmos. Sci. 55, 3042–3064.
  * Suzuki & Fox-Kemper (2016) Suzuki, N. & Fox-Kemper, B. 2016 Understanding stokes forces in the wave-averaged equations. J. Geophys. Res. Oceans 121 (5), 3579–3596.
  * Taylor & Sarkar (2008) Taylor, J. R. & Sarkar, S. 2008 Stratification effects in a bottom ekman layer. J. Phys. Oceanogr. 38 (11), 2535–2555.
  * Thorpe (2004) Thorpe, S. A. 2004 Langmuir circulation. Annu. Rev. Fluid Mech. 36 (1), 55–79.
  * Troell et al. (2009) Troell, M., Joyce, A., Chopin, T., Neori, A., Buschmann, A. H. & Fang, J. G. 2009 Ecological engineering in aquaculture - potential for integrated multi-trophic aquaculture (IMTA) in marine offshore systems. Aquaculture 297, 1–9.
  * Tseung et al. (2016) Tseung, H. L., Kikkert, G. A. & Plew, D. 2016 Hydrodynamics of suspended canopies with limited length and width. Environ. Fluid Mech. 16, 145–166.
  * Utter & Denny (1996) Utter, B. & Denny, M. 1996 Wave-induced forces on the giant kelp Macrocystis pyrifera (Agardh): field test of a computational model. J. Exp. Biol. 199, 2645–2654.
  * Van Roekel et al. (2012) Van Roekel, L. P., Fox-Kemper, B., Sullivan, P. P., Hamlington, P. E. & Haney, S. R. 2012 The form and orientation of langmuir cells for misaligned winds and waves. J. Geophys. Res. Oceans 117 (C5).
  * Vogel (1989) Vogel, S. 1989 Drag and reconfiguration of broad leaves in high winds. J. Exp. Bot. 40, 941–948.
  * Yan et al. (2017) Yan, C., Nepf, H. M., Huang, W. & Cui, G. 2017 Large eddy simulation of flow and scalar transport in a vegetated channel. Environ. Fluid Mech. 17, 497–519.
  * Yang et al. (2015) Yang, D., Chen, B., Chamecki, M. & Meneveau, C. 2015 Oil plumes and dispersion in langmuir, upper‐ocean turbulence: Large‐eddy simulations and k‐profile parameterization. J. Geophys. Res. Oceans 120, 4729–4759.
  * Zhao et al. (2017) Zhao, F., X., Huai W. & Li, D. 2017 Numerical modeling of open channel flow with suspended canopy. Adv. Water Resour. 105, 132–143.
  * Zhou & Venayagamoorthy (2019) Zhou, J. & Venayagamoorthy, S. 2019 Near-field mean flow dynamics of a cylindrical canopy patch suspended in deep water. J. Fluid Mech. 858, 634–655.