Hostname: page-component-7b9c58cd5d-9k27k Total loading time: 0 Render date: 2025-03-16T13:35:54.800Z Has data issue: false hasContentIssue false

Oceanic mixing and waves in the presence of a suspended canopy

Published online by Cambridge University Press:  14 March 2025

Tong Bo*
Affiliation:
Department of Atmospheric and Oceanic Sciences, University of California Los Angeles, Los Angeles, CA 90095-1565, USA
James C. McWilliams
Affiliation:
Department of Atmospheric and Oceanic Sciences, University of California Los Angeles, Los Angeles, CA 90095-1565, USA
Marcelo Chamecki
Affiliation:
Department of Atmospheric and Oceanic Sciences, University of California Los Angeles, Los Angeles, CA 90095-1565, USA
*
Corresponding author: Tong Bo, tbo@atmos.ucla.edu

Abstract

Large-eddy simulations are analysed to determine the influence of suspended canopies, such as those formed in macroalgal farms, on ocean mixed layer (OML) deepening and internal wave generation. In the absence of a canopy, we show that Langmuir turbulence, when compared with wind-driven shear turbulence, results in a deeper OML and more pronounced internal waves beneath the OML. Subsequently, we examine simulations with suspended canopies of varying densities located in the OML, in the presence of a background geostrophic current. Intensified turbulence occurs in the shear layer at the canopy’s bottom edge, arising from the interaction between the geostrophic current and canopy drag. Structures resembling Kelvin–Helmholtz (KH) instability emerge as the canopy shear layer interacts with the underlying stratification, radiating internal waves beneath the OML. Both intensified turbulence and lower-frequency motions associated with KH-type structures are critical factors in enhancing mixing. Consequently, the OML depth increases by up to a factor of two compared with cases without a canopy. Denser canopies and stronger geostrophic currents lead to more pronounced KH-type structures and internal waves, stronger turbulence and greater OML deepening. Additionally, vertical nutrient transport is enhanced as the OML deepens due to the presence of the canopy. Considering that the canopy density investigated in this study closely represents offshore macroalgal farms, these findings suggest a mechanism for passive nutrient entrainment conducive to sustainable farming. Overall, this study reveals the intricate interactions between the suspended canopy, turbulent mixing and stratification, underscoring their importance in reshaping OML characteristics.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2025. Published by Cambridge University Press

1. Introduction

The ocean mixed layer (OML) – also referred to as the ocean surface boundary layer – plays a pivotal role in the global climate system, as it directly influences momentum, heat and gas exchanges between the ocean and atmosphere (D’Asaro Reference D’Asaro2014; Chamecki et al. Reference Chamecki, Chor, Yang and Meneveau2019). A key factor governing OML characteristics is the intricate dynamics of turbulent vertical mixing (D’Asaro Reference D’Asaro2001; Belcher et al. Reference Belcher2012a ). Wind stress typically drives a vertically sheared flow in the upper ocean boundary layer, leading to shear turbulence and mixing (e.g. D’Asaro Reference D’Asaro2001). Moreover, surface gravity waves can interact with wind-driven currents and give rise to Langmuir turbulence in the OML (e.g. Leibovich Reference Leibovich1983; McWilliams et al. Reference McWilliams, Sullivan and Moeng1997; Li & Garrett Reference Li and Garrett1997). Additionally, convective processes also enhance turbulence and mixing when the atmosphere cools down the ocean (e.g. Shay & Gregg Reference Shay and Gregg1986).

The OML turbulence may generate and radiate internal waves when it interacts with the stratification beneath the OML (Chini & Leibovich Reference Chini and Leibovich2003; Dohan & Sutherland Reference Dohan and Sutherland2005; Polton et al. Reference Polton, Smith, MacKinnon and Tejada-Martínez2008; Munroe & Sutherland Reference Munroe and Sutherland2008). Similar internal wave generation has been studied in the bottom boundary layer, where broad-band wall-generated turbulence induces waves that propagate upward in the overlying stratified layer (Taylor & Sarkar Reference Taylor and Sarkar2007; Gibson et al. Reference Gibson, Bondur, Norris Keeler and Leung2006). These internal waves typically exhibit a narrow range of wavenumbers as a result of the frequency selection mechanism. In the OML, the frequency selection mechanism may likewise lead to downward-propagating internal waves in the underlying stratification. The prevalence of Langmuir turbulence in the OML typically plays a dominant role on internal wave generation and radiation (Polton et al. Reference Polton, Smith, MacKinnon and Tejada-Martínez2008).

Various obstacle structures situated near the ocean surface boundary, such as aquacultural farms, engineered offshore platforms and sea ice, have the potential to interact with ocean currents and waves. The existence of these structures may consequently impact turbulence generation and mixing processes, thereby influencing OML dynamics as well as its interaction with the underlying stratification. This study investigates turbulence, mixing and internal waves associated with a suspended canopy located within the OML. The investigation of the suspended canopy is motivated by offshore macroalgal farming, which emerges as a potential sustainable strategy for carbon sequestration, biofuel production, food supply and bioremediation (Frieder et al. Reference Frieder2022; Arzeno-Soltero et al. Reference Arzeno-Soltero, Saenz, Frieder, Long, DeAngelo, Davis and Davis2023).

The canopy imposes drag forces on the flow (e.g. Thom Reference Thom1971; Jackson Reference Jackson1997), resulting in attenuation of current velocity and wave motions (Rosman et al. Reference Rosman, Koseff, Monismith and Grover2007; Monismith et al. Reference Monismith, Alnajjar, Daly, Valle-Levinson, Juarez, Fagundes, Bell and Woodson2022). In the case of a submerged canopy located on the bottom boundary, shear layer turbulence typically develops at the top of the canopy due to drag discontinuities (e.g. Finnigan Reference Finnigan2000; Nepf Reference Nepf2012). Similarly, for a suspended canopy near the surface boundary, a shear layer can form at the bottom of the canopy, leading to the generation of turbulence and the exchange of momentum and scalars between the canopy and the underlying flow (Plew Reference Plew2011). In addition, the interplay between the canopy, ocean currents and surface gravity waves can induce modifications to Langmuir turbulence within the OML (Yan et al. Reference Yan, McWilliams and Chamecki2021; Bo et al. Reference Bo, McWilliams, Yan and Chamecki2024b ).

While canopy flow dynamics have been extensively studied under neutral flow conditions, relatively less attention has been paid to canopies in stably stratified flows. The presence of stratification can lead to more complex interactions between the canopy flow and vertical density gradients (Belcher et al. Reference Belcher, Harman and Finnigan2012b ). Plew et al. (Reference Plew, Spigel, Stevens, Nokes and Davidson2006) observed vertical mixing arising from the interplay between the canopy shear layer and stratification in a field study. Although the potential for Kelvin–Helmholtz (KH) instability was inferred, direct observational evidence was lacking. Previous studies have investigated the development of KH instability and the radiation of internal waves in stratified shear layers, in the absence of a canopy (Pham et al. Reference Pham, Sarkar and Brucker2009; Pham & Sarkar Reference Pham and Sarkar2010, Reference Pham and Sarkar2011; Pham et al. Reference Pham, Sarkar and Winters2012). Additionally, Sutherland & Linden (Reference Sutherland and Linden1998) analysed lee vortices and internal waves generated in the wake of a thin vertical barrier. Nevertheless, the influence of suspended canopies on stratified flow dynamics remains largely unexplored.

The suspended canopy can alter the dynamics of the OML and the underlying stratification. On the other hand, OML characteristics typically govern the vertical transport of various scalars in the upper ocean, such as nutrient, pollutant and oxygen (e.g. Polovina et al. Reference Polovina, Mitchum and Evans1995; Chamecki et al. Reference Chamecki, Chor, Yang and Meneveau2019), thereby leading to feedback on the canopy. For example, suspended macroalgal farms are often located within the OML to maximize light exposure. However, nutrient concentrations are typically lower near the sea surface and higher below the OML, and the risk of farm starvation escalates if nutrient availability within the OML remains limited (Frieder et al. Reference Frieder2022; Arzeno-Soltero et al. Reference Arzeno-Soltero, Saenz, Frieder, Long, DeAngelo, Davis and Davis2023). Therefore, investigating the interaction between the canopy, OML turbulence and the underlying stratification is also crucial for predicting nutrient supply to the OML and optimizing macroalgal farm performance (Frieder et al. Reference Frieder2022; Bo et al. Reference Bo, McWilliams, Frieder, Davis and Chamecki2024a ).

In this study, we use large-eddy simulations (LES) to explore the impact of a suspended canopy on OML dynamics, including turbulent mixing, OML depth variability and internal wave radiation beneath the OML. We first present simulations without a canopy, as background cases for comparison, to show the influence of Langmuir turbulence on OML mixing. Subsequently, we focus on simulations that incorporate a canopy, examining its importance in reshaping OML characteristics compared with Langmuir turbulence. Section 2 describes the numerical framework and simulation set-up. Section 3 compares OML characteristics with and without the canopy, and analyses turbulence generation, KH instability and internal wave propagation. Section 4 concludes the study and discusses broader implications for canopy flows and the OML dynamics.

2. Methods

2.1. Model description

The present LES framework is based on a set of wave-averaged and grid-filtered equations for mass, momentum, heat and passive tracer:

(2.1) \begin{equation} \nabla \cdot \breve{\boldsymbol{u}}=0, \end{equation}
(2.2) \begin{eqnarray} \frac {\partial \breve{\textbf{u}}}{\partial {t}} + \breve{\textbf {u}}\cdot \nabla \breve{\textbf {u}} & = & -\nabla \Pi - f\textbf {e}_z\times \left (\breve{\boldsymbol {u}} + {\textbf {u}}_s-\textbf {u}_g\right ) + \textbf {u}_s\times \breve{\boldsymbol {\zeta }} \nonumber \\ && + \left (1-\frac {\breve{\rho}}{\rho _0}\right )g\textbf {e}_z - \nabla \cdot \boldsymbol {\tau }^{sgs} - \textbf {F}_D, \end{eqnarray}
(2.3) \begin{equation} \frac {\partial {\breve{\theta }}}{\partial {t}}+\left (\breve{\textbf {u}} + \textbf { u}_s\right ) \cdot \nabla \breve{\theta } = - \nabla \cdot {\unicode{x03C0}}_\theta ^{sgs}, \end{equation}
(2.4) \begin{equation} \frac {\partial {\breve{C}_N}}{\partial {t}}+\left (\breve{\textbf {u}} + \textbf {u}_s\right ) \cdot \nabla \breve{C}_N = - \nabla \cdot {\unicode{x03C0}}_{C_N}^{sgs} + \mathcal {S}. \end{equation}

The mathematical framework presented above was first introduced in McWilliams et al. (Reference McWilliams, Sullivan and Moeng1997), which extended the original Craik–Leibovich equations in Craik & Leibovich (Reference Craik and Leibovich1976) by incorporating the effects of planetary rotation and advection of scalar by Stokes drift. The code employed in this study has been validated against Langmuir turbulence simulations by McWilliams et al. (Reference McWilliams, Sullivan and Moeng1997) and applied to prior studies on canopy flows and oceanic boundary layer flows (Yan et al. Reference Yan, McWilliams and Chamecki2021,Reference Yan, McWilliams and Chamecki2022; Bo et al. Reference Bo, McWilliams, Yan and Chamecki2024b ).

The breve notation $\breve{}$ in (2.1), (2.2), (2.3) and (2.4) indicates grid-filtered variables. In the Cartesian coordinate system $\textbf {x} = (x,y,z)$ , the velocity vector is $\breve{\textbf {u}}=(\breve{u},\breve{v},\breve{w})$ , i.e. the streamwise, cross-stream and vertical components, respectively. The filtered potential temperature is $\breve{\theta}$ in (2.3), and the filtered passive tracer concentration is $\breve{C}_N$ in (2.4) (here specified as nitrate).

In (2.2), $g$ is the gravitational acceleration and $\textbf {e}_z$ is the unit vector in the vertical direction. The term $f$ represents the Coriolis frequency and $\breve{\boldsymbol {\zeta }} =\nabla \times \breve{\textbf {u}}$ denotes the filtered vorticity. Here $\textbf {u}_s$ is the Stokes drift associated with surface gravity waves, and $\textbf {u}_g$ corresponds to a geostrophic current representing mesoscale ocean flows, driven by an external pressure gradient $f\textbf {e}_z\times \textbf {u}_g$ . The term $\breve{\rho}$ is the filtered density, with $\rho _0$ denoting the reference density. Density variations are assumed to be caused only by the potential temperature $\breve{\theta}$ through a linear relationship $\rho =\rho _0[1-\alpha (\theta -\theta _0)]$ , where $\alpha =2\times 10^{-4}\,\rm K^-{^1}$ is the thermal expansion coefficient and $\theta _0$ is the reference potential temperature corresponding to $\rho _0$ . In (2.2), $\Pi$ represents the modified pressure, with $\Pi =\breve{p}/\rho _0+\breve{\textbf {u}}\cdot \textbf {u}_s+|\textbf {u}_s|^2/2$ and $\breve{p}$ denoting the filtered dynamic pressure (e.g. Chamecki et al. Reference Chamecki, Chor, Yang and Meneveau2019). The term $\boldsymbol {\tau }^{sgs}$ represents the subgrid-scale (SGS) stress tensor, which is determined by using the Lagrangian scale-dependent dynamic Smagorinsky model (Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005). The term $\textbf {F}_D$ denotes the drag force imposed by the canopy onto the flow, and more detailed treatment of canopy drag will be presented in a later paragraph.

In (2.3), ${\unicode{x03C0} }_\theta ^{sgs}$ is the SGS heat flux, and in (2.4), ${\unicode{x03C0}}_{C_N}^{sgs}$ is the SGS tracer flux. The SGS heat and tracer fluxes are modelled using an eddy diffusivity closure, with diffusivity obtained from SGS viscosity and a prescribed value of turbulent Prandtl number $Pr_t = 0.4$ . Molecular viscosity and diffusivity are assumed to be negligible for the high Reynolds number flows examined in this study. It is worthwhile to note that more complex interactions between fine-scale turbulence and canopy elements may occur, where molecular viscosity and diffusivity could play a more important role. Future efforts should focus on improving the SGS model to better capture these interactions in canopy flow simulations. The term $\mathcal {S}$ in (2.4) is a sink or source term, representing nutrient uptake or exudation by macroalgae. In this study, we set $\mathcal {S}=0$ and only focus on the transport of nutrient driven by turbulence and waves.

The model employed in this study does not explicitly resolve surface gravity waves. Instead, the time-averaged effects of surface waves are incorporated by imposing the Stokes drift $\textbf {u}_s$ . We consider a scenario of monochromatic deep water waves propagating in the $x$ -direction, with amplitude $A_w$ and frequency $\Omega _w=\sqrt {gk_w}$ , where $k_w$ is the wavenumber. The Stokes drift velocity is then given by $\textbf {u}_s=(u_s,0,0)$ , where

(2.5) \begin{equation} u_s = U_s\textrm {e}^{2kz} \end{equation}

and $U_s=\Omega _w k_w A_w^2$ is the Stokes drift velocity at the surface. The impact of waves on OML turbulence is represented by the Craik–Leibovich vortex force $\textbf {u}_s\times \breve{\boldsymbol {\zeta }} = (0,-u_s \breve{\zeta }_z,u_s \breve{\zeta }_y)$ , corresponding to the third term on the right-hand side of (2.2).

The expression for the canopy drag force $F_D$ in (2.2) is (Shaw & Schumann Reference Shaw and Schumann1992; Pan et al. Reference Pan, Chamecki and Isard2014; Yan et al. Reference Yan, McWilliams and Chamecki2021)

(2.6) \begin{equation} \textbf {F}_D = \frac {1}{2}C_D a \textbf {P}\cdot \left (|\breve{\textbf {u}}|\breve{\textbf {u}}\right ), \end{equation}

with $C_D$ representing the drag coefficient. The canopy drag interacts with ocean currents based on the quadratic drag law. The canopy parameters are determined according to the specific macroalgae under investigation here (giant kelp, Macrocystis pyrifera). We use $C_D=0.0148$ based on the experimental study of Utter & Denny (Reference Utter and Denny1996), with more detailed discussion on this choice of $C_D$ provided in Yan et al. (Reference Yan, McWilliams and Chamecki2021). In (2.6), $a$ stands for the frond surface area density (m $^{-1}$ , area per volume). The coefficient tensor $\textbf {P}$ represents the projection of frond surface area into each direction, and here we use $\textbf {P}=(1/2)\textbf {I}$ , with $\textbf {I}$ denoting the identity matrix. In the subsequent sections, we will omit the breve symbols that denote grid-filtered variables for simplicity.

The present LES framework uses a pseudospectral method in the horizontal directions and a second-order central finite-difference method in the vertical direction. Zero padding is applied in each horizontal direction based on the 3/2 rule for full dealiasing. Periodic boundary conditions are employed in the horizontal directions. A constant wind shear stress $\tau _w$ is imposed at the upper boundary in the $x$ -direction. No-flux boundary conditions are used for temperature and nutrient concentration at the upper boundary. A sponge layer is applied within the lower 20 % of the domain to mitigate the reflection of internal gravity waves. The absorption rate is determined based on the buoyancy frequency in the simulation (Yan et al. Reference Yan, McWilliams and Chamecki2021), and a sensitivity test will be presented later to ensure the effectiveness of the sponge layer.

2.2. Simulation set-up

The simulation forcing conditions in this study are generally consistent with McWilliams et al. (Reference McWilliams, Sullivan and Moeng1997). A geostrophic current $\textbf {u}_g=(u_g,0,0)$ is imposed in the $x$ -direction, the same as the direction of wind and surface gravity waves. This geostrophic flow is driven by an external pressure gradient $fu_g$ in the $y$ -direction, with $u_g=0.2\,\rm m\,s^-{^1}$ . Additional cases with varying $u_g$ values are also conducted, as detailed in a following paragraph. The Coriolis frequency $f=10^{-4}$ s $^{-1}$ corresponds to a latitude of around 45 $^\circ$ N. A constant wind stress of $\tau _w=0.037\,\rm N\,m^-{^2}$ is applied at the surface boundary, corresponding to a wind speed at 10 m height above the surface of $5\,\rm m\,s^-{^1}$ and a sea surface drag coefficient of 0.001–0.0015 (e.g. Smith Reference Smith1988). The friction velocity is $u_*=\sqrt {\tau _w/\rho }=0.0061\,\rm m\,s^-{^1}$ . The amplitude of surface gravity waves is $A_w=0.8$ m. The wavelength is $\lambda _w=60$ m, with a wave period of 6.2 s. This corresponds to a surface Stokes velocity of $U_s=0.068\,\rm m\,s^-{^1}$ and a turbulent Langmuir number of $La_t=\sqrt {u_*/U_s}=0.3$ .

The simulations are conducted on a $400\times 200\times 120\,\rm m^3$ domain, with $192\times 96\times 240$ uniformly distributed grid cells. The vertical resolution is sufficient for resolving the Ozmidov length scale. A sensitivity test on grid resolution was performed by doubling the number of grid cells in the horizontal directions, which yielded consistent results. A more comprehensive analysis of grid resolution can be found in a previous study by Yan et al. (Reference Yan, McWilliams and Chamecki2021). The simulations are continued for 90 000 s, i.e. $tf=9$ , and the time step is $0.15$ s.

Figure 1. Simulation set-up. $(a)$ Initial temperature profile. The vertical coordinate $z$ is normalized by the initial OML depth $h_0$ . Dotted black lines indicate the canopy base (upper line) and the sponge layer (lower line). $(b)$ Initial nutrient profile. $(c)$ A schematic of the canopy simulation (side view). The frond surface area density $a$ is constant and uniform throughout the extent of the canopy. Values of $a$ for different simulations are detailed in table 1.

The initial mixed layer depth $h_0$ is 26 m, and a stably stratified layer is beneath it with a uniform temperature gradient $\partial \theta /\partial z|_\infty =0.01\,\rm K\,m^-{^1}$ (figure 1 a). The underlying stratification corresponds to a buoyancy frequency of $N=0.0044\,\rm m^-{^1}$ , where $N$ is defined as

(2.7) \begin{equation} N = \sqrt {\alpha g \frac {\partial {\theta }}{\partial {z}}}. \end{equation}

According to Li & Garrett (Reference Li and Garrett1997), Langmuir turbulence (in the absence of a canopy) can induce rapid deepening of the OML to a depth of approximately 10–17 $u_*/N$ . In addition, OML deepening driven by surface winds can rapidly reach a depth of around $2^{3/4}u_*/\sqrt {Nf}$ , as demonstrated in prior studies on the stratified Ekman layer (Pollard et al. Reference Pollard, Rhines and Thompson1973; Pham & Sarkar Reference Pham and Sarkar2017). We therefore set the initial OML depth $h_0$ to be greater than both $17 u_*/N$ and $2^{3/4}u_*/\sqrt {Nf}$ to prevent rapid OML deepening due to either wind-driven turbulence or Langmuir turbulence. This choice allows us to examine any further modifications to the OML due to the presence of the canopy. The initial nutrient concentration is set as a step function, with a value of 0 within the OML ( $-h_0\lt z\lt 0$ ) and 10 000 $\unicode{x03BC}$ mol m $^{-3}$ for $z\lt -h_0$ (figure 1 b). This value represents a typical upper range of nitrate concentration in the tropical, subtropical and temperate oceans (Voss et al. Reference Voss, Bange, Dippner, Middelburg, Montoya and Ward2013; Claustre et al. Reference Claustre, Johnson and Takeshita2020; Frieder et al. Reference Frieder2022).

The initial velocity profiles are uniform, with $u_0=u_g$ and $v_0=0$ . We initialize the flow with laminar conditions, introducing small perturbations to trigger turbulence, without any pre-existing OML turbulence. In this study, the canopy shear layer turbulence can be over an order of magnitude stronger than classical OML turbulence without a canopy, as demonstrated in a subsequent section. Generally, weak initial turbulence leads to only minor deviations from laminar initial conditions in the formation of KH structures (Kaminski & Smyth Reference Kaminski and Smyth2019). Therefore, pre-existing OML turbulence is expected to have a small influence on the canopy shear layer dynamics examined here, and the use of laminar initial conditions is a reasonable choice. It is also worth noting that the presence of strong turbulence in the initial conditions may inhibit the development of KH instability (Kaminski & Smyth Reference Kaminski and Smyth2019). This study focuses on one set of oceanic forcing conditions. Exploration of other conditions that generate stronger OML turbulence, as well as how such pre-existing turbulence may affect canopy flow dynamics, is beyond the scope of this research.

In the horizontal directions, the canopy extends across the entire domain with periodic boundary conditions. This set-up effectively assumes a horizontally homogeneous canopy, representing an interior region of a large farm, where no further spatial development of the flow occurs. In the vertical direction the canopy is located between the sea surface and the canopy base depth of $h_b=20\,\rm m$ (figure 1 c). The canopy base depth is shallower than the initial mixed layer depth, consistent with common design practices for suspended macroalgal farms (Frieder et al. Reference Frieder2022). The frond surface area density $a$ is constant and uniform within the extent of the canopy. The bottom edge of the canopy features a sharp interface, where the frond area density transitions from $a$ to 0 across two grid cells in the vertical direction. The sharp canopy edge employed in this study aligns with realistic farm design practices. We note that the sharpness of this transition at the canopy edge can influence the canopy shear layer dynamics (Ouwersloot et al. Reference Ouwersloot, Moene, Attema and De Arellano2017).

We conduct four simulations, namely cases C0.1 to C3, with the canopy density $a$ ranging from $0.1\,\rm m^-{^1}$ to $3\,\rm m^-{^1}$ (table 1). These values fall within the typical range of frond area density in macroalgal farms, estimated through the conversion of algal biomass (Frieder et al. Reference Frieder2022). As a comparison, we also perform a simulation without the canopy (case NC). Additionally, we investigate a scenario without surface gravity waves for two simulations (cases NC-NSW and C3-NSW) by setting the wave-induced Stokes drift velocity to zero. All the simulations mentioned above are conducted with a geostrophic current of $u_g=0.2\,\rm m\,s^-{^1}$ . We also run simulations for the highest canopy density with a weaker geostrophic current of $u_g=0.1\,\rm m\,s^-{^1}$ and no geostrophic current ( $u_g=0\,\rm m\,s^-{^1}$ ), namely cases C3-WG and C3-NG.

Table 1. Simulation parameters. Here ‘NC’ represents cases with no canopy, and ‘C’ represents cases with a canopy, with 0.1, 0.3, 1 and 3 denoting the canopy frond area density $a$ . The suffix ‘-NSW’ indicates simulations conducted without surface gravity waves, where the Stokes drift velocity is zero. The suffix ‘-WG’ represents the case with a weaker geostrophic current $u_g=0.1\,\rm m\,s^-{^1}$ , and ‘-NG’ represents the case with no geostrophic current $u_g=0\,\rm m\,s^-{^1}$ . Additionally, a simulation is conducted with a doubled domain depth, named ‘-DEEP’.

Another simulation (case C3-DEEP, $u_g=0.2\,\rm m\,s^-{^1}$ ) is performed with a doubled domain depth of 240 m and a vertical grid cell number of 480. This allows us to explore potential internal wave propagation in deeper water and confirm that internal wave reflection from the bottom boundary does not pose an issue in the other simulations with a shallower domain. Additionally, a simulation with a finite canopy length of 400 m is conducted to examine horizontal edge effects and the potential streamwise development of canopy flow (simulation parameters provided in appendix B).

3. Results and analyses

3.1. An overview of the results

In this section, we first provide an overview of the OML dynamics and the distinct flow features associated with the canopy. Subsequently, we present more detailed analyses of turbulence, internal waves and vertical mixing. The velocity field is decomposed as

(3.1) \begin{equation} \textbf {u} = \overline {\textbf {u}} + \textbf {u}', \end{equation}

where $\textbf {u}'$ represents high-frequency fluctuations (the turbulence component), while $\overline {\textbf {u}}$ corresponds to other lower frequency processes such as internal waves. The separation between the two components is based on a time average over every 100 time steps, i.e. 15 s, and this averaging period is generally shorter than the internal wave period and the OML evolution time scale. In the following results, the time-averaged horizontal velocity $\overline {u}$ and $\overline {v}$ will be normalized by the geostrophic velocity $u_g$ , while the time-averaged vertical velocity $\overline {w}$ and all turbulence components $\textbf {u}'$ will be normalized by the friction velocity $u_*$ . Similarly, the covariance between velocity and any variable $\phi$ can be decomposed as

(3.2) \begin{equation} \overline {\textbf {u}\phi } = \overline {\textbf {u}}\overline {\phi } + \overline {\textbf {u}'\phi '}. \end{equation}

The first term on the right-hand side represents the contribution from the low-frequency flow, and the second term denotes the turbulent flux.

Figure 2. The mixed layer depth $h_{OML}$ , normalized by its initial value $h_0$ , as a function of time. $(a)$ Simulations with varying canopy densities (details in the legend). Generally, higher canopy density (cases C0.1–C3) leads to more pronounced OML deepening compared with cases without a canopy (cases NC and NC-NSW). $(b)$ Simulations with different geostrophic currents $u_g$ . Weaker geostrophic currents lead to less OML deepening.

The evolution of the mixed layer depth in different simulations is compared in figure 2. The OML depth $h_{OML}$ is defined as the minimum depth where the vertical temperature gradient equals to the gradient of the underlying stratified flow (Taylor & Sarkar Reference Taylor and Sarkar2007), i.e.

(3.3) \begin{equation} h_{OML} = min (|z|) \Big |_ {\partial \langle \overline \theta \rangle _{xy}/\partial z \ge \partial \theta /\partial z|_\infty }, \end{equation}

where $\left \langle \ \right \rangle _{xy}$ represents averaging over the horizontal plane. The OML depth remains relatively unchanged in the two cases without the canopy (NC and NC-NSW, with and without the Stokes drift, respectively), with a slight increase of around 10 % compared with the initial depth $h_0$ in case NC and less than 5 % increase in case NC-NSW. By contrast, the OML is significantly deepened due to the presence of the canopy (figure 2 a). This increase in OML depth exhibits a positive correlation with the canopy density $a$ , with an increase to more than twice the original value in the case with the highest canopy density. In addition, for a constant canopy density, weaker geostrophic currents $u_g$ lead to less OML deepening (figure 2 b). Notably, in the absence of a geostrophic current (case C3-NG), the OML depth increases by only around 10%, comparable to the case without a canopy (case NC). These results suggest that OML deepening arises from the interaction between the geostrophic current and canopy, which involves the generation of shear layer turbulence as analysed below.

Figure 3. Side views of the cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$ . Snapshots at around $tf=1$ are shown for four representative cases: $(a)$ NC, no canopy, with the Stokes drift; $(b)$ NC-NSW, no canopy, no Stokes drift; $(c)$ C0.3, canopy density of 0.3 m $^{-1}$ ; $(d)$ C3, canopy density of 3 m $^{-1}$ ; $(e)$ C3-WG, canopy density of 3 m $^{-1}$ and a weaker geostrophic current; and $(f)$ C3-NG, canopy density of 3 m $^{-1}$ and no geostrophic current. Black dotted lines represent the canopy base depth (lines also included in the no-canopy cases for comparison purposes). Note that the sponge layer is not shown in the plots.

Along with variations in the OML depth, distinct flow patterns are observed across different simulations. In the side-view snapshots of the two cases without the canopy (figures 3 a and 3 b), more pronounced internal wave patterns are observed in the case with Stokes drift (NC) compared with the case without Stokes drift (NC-NSW). This comparison indicates the importance of Langmuir turbulence, arising from the interaction between surface waves and wind-driven currents, for internal wave generation (Polton et al. Reference Polton, Smith, MacKinnon and Tejada-Martínez2008) (details examined later).

In the presence of a canopy, patterns resembling KH instability (e.g. Kundu et al. Reference Kundu, Cohen and Dowling2015) emerge near the canopy bottom edge. The dominant wavelength of these patterns positively depends on the canopy density (comparing figures 3 c and 3 d). Additionally, the KH-type patterns also lead to the radiation of internal waves in the deeper water. In the case with a weaker $u_g$ , the dominant wavelength of the KH structures decreases (figure 3 e). In the absence of a geostrophic current ( $u_g = 0$ ), the KH structures do not occur (figure 3 f), and the internal wave patterns resemble those observed without a canopy.

The impact of Stokes drift on the canopy shear layer is minimal. The two canopy cases, with and without Stokes drift (C3 and C3-NSW, figure 2 a), exhibit similar changes in OML depth, and the KH structures in these two cases are also consistent (case C3-NSW not shown). This result aligns with expectations, given that the KH instability primarily relates to the vertical shear generated by the canopy, and Stokes drift is weak at that depth. In this study, we focus on the canopy simulation with Stokes drift (case C3) to facilitate comparison with the Langmuir turbulence simulation (case NC).

Overall, the presence of a canopy significantly impacts OML dynamics, including the emergence of KH-type patterns, the deepening of the OML, and the radiation of internal waves in the underlying stratification. Further analysis of these effects will be conducted in the following sections. It is also worth mentioning that velocity and temperature profiles, as well as the KH-type patterns, evolve over time as the OML depth gradually increases. The temporal evolution of velocity and temperature dynamics will be discussed in § 3.3.

Figure 4. Velocity and temperature profiles in four representative simulations at the final simulation time $tf=9$ . $(a,b)$ Horizontally averaged streamwise velocity $\langle \overline {u}\rangle _{xy}$ and cross-stream velocity $\langle \overline {v}\rangle _{xy}$ , both normalized by the geostrophic velocity $u_g=0.2\,\rm m\,s^-{^1}$ . $(c,d)$ Horizontally averaged temperature $\langle \overline {\theta }\rangle _{xy}$ and the vertical gradient of temperature. Horizontal dotted lines represent the canopy base depth.

3.2. Final velocity and temperature profiles

Figure 4 compares the velocity and temperature profiles for four representative cases (the same as the ones depicted in figure 3), namely NC, NC-NSW, C0.3 and C3, at the final simulation time $tf=9$ , when the OML depth reaches its maximum. In the two cases without a canopy, the horizontal velocity components $u$ and $v$ exhibit a more uniform distribution with depth in the presence of Stokes drift and Langmuir turbulence compared with their absence (cases NC and NC-NSW, represented by the two black lines in figures 4 a and 4 b). This aligns with the findings by McWilliams et al. (Reference McWilliams, Sullivan and Moeng1997), indicating that Langmuir turbulence leads to more efficient vertical momentum transport than wind-driven shear turbulence.

A peak in the vertical temperature gradient $\partial \langle \theta \rangle _{xy}/\partial z$ forms as the OML deepens (figure 4 $d$ ), consistent with results from previous studies (e.g. Li & Garrett Reference Li and Garrett1997). This local maximum of the vertical temperature gradient is commonly referred to as a pycnocline (here a thermocline), located between the OML and the underlying uniform stratification (e.g. Taylor & Sarkar Reference Taylor and Sarkar2007). The case with the Stokes drift demonstrates a deeper thermocline and a more notable peak in $\partial \langle \theta \rangle _{xy}/\partial z$ , because Langmuir turbulence can enhance vertical mixing within the OML (Li & Garrett Reference Li and Garrett1997; Kukulka et al. Reference Kukulka, Plueddemann, Trowbridge and Sullivan2010). Nevertheless, the initial OML depth $h_0$ is greater than the threshold value of 10–17 $u_*/N$ , where Langmuir turbulence is expected to cause rapid OML deepening (Li & Garrett Reference Li and Garrett1997). As a result, the observed increase in the OML depth is relatively small, e.g. around 10%, in the absence of a canopy.

When a canopy is present in the OML, the streamwise velocity is significantly decelerated near the sea surface due to the canopy drag force, leading to the formation of a shear layer below the canopy (figure 4 $a$ ). The extent of deceleration and the thickness of the shear layer exhibit a positive dependence on the canopy density (comparing cases C1 and C3). Additionally, a positive cross-stream velocity arises (to the left of the streamwise current) due to the combined effects of the Ekman spiral and canopy drag (Yan et al. Reference Yan, McWilliams and Chamecki2021) (figure 4 $b$ ).

The shear layer attached to the canopy results in the emergence of KH-type patterns and significant deepening of the OML. By the end of the simulation, the local maximum of the vertical temperature gradient $\partial \langle \theta \rangle _{xy}/\partial z$ reaches close to three times its initial value in the case with the highest canopy density (figure 4 $d$ ). The increased OML depth generally corresponds with the thickness of the canopy shear layer (comparing figures 4 $a$ and 4 $d$ ), suggesting that canopy-generated shear is the predominant mechanism driving mixing and OML deepening. In addition, a secondary peak of $\partial \langle \theta \rangle _{xy}/\partial z$ emerges within the canopy in case C3 (e.g. at around $z/h_b=-0.3$ in figure 4 $d$ ; also see figure 3 $d$ ). The formation of this secondary peak is attributed to the upward advection of gradients from the thermocline by the KH-type structures, as well as the localized inhibition of turbulence and mixing inside the canopy.

3.3. Kelvin–Helmholtz instability

Patterns resembling KH instability emerge as the canopy shear layer interacts with the underlying stratification (figure 3). In this section, we further investigate characteristics of the observed KH-type patterns. The KH instability occurs when the destabilizing effect of velocity shear overcomes the stabilizing effect of stratification (e.g. Kundu et al. Reference Kundu, Cohen and Dowling2015). The gradient Richardson number $Ri_g$ compares the two competing effects and is defined as

(3.4) \begin{equation} Ri_g = \frac {N^2}{S^2}, \end{equation}

where $N$ is the buoyancy frequency as described by (2.7) and $S$ is the vertical shear given by

(3.5) \begin{equation} S = \sqrt {\left (\frac {\partial {u}}{\partial {z}}\right )^2+\left (\frac {\partial {v}}{\partial {z}}\right )^2}. \end{equation}

The $Ri_g$ is calculated from the horizontally averaged velocity and temperature profiles, and $Ri_g\lt 0.25$ is found in the canopy shear layer when the KH-type structures occur (please refer to the Supplementary material for the movie of $Ri_g$ ). This is consistent with the condition favourable for the growth of KH instability and mixing obtained from linear stability analysis (Miles Reference Miles1961; Howard Reference Howard1961).

Figure 5. $(a)$ The dominant $x$ -direction wavenumber $k_{x,peak}$ versus canopy shear layer thickness $\delta _s$ . Different colours represent simulations with different canopy density and $u_g$ (C0.1–C3 and C3-WG, see the legend), and each point corresponds to a time interval of every $tf=0.0015$ . $(b)$ The dominant phase speed $c_{x,peak}$ of the KH-type structures versus the current speed of the inflection point in the shear layer. Note that the results are normalized by $u_g=0.1\,\rm m\,s^-{^1}$ in case C3-WG and by $u_g=0.2\,\rm m\,s^-{^1}$ in all the other cases.

Moreover, the wavelength of the most unstable mode of KH instability is typically 7–8 times the shear layer thickness (Sutherland Reference Sutherland and Rahman2005). The canopy shear layer thickness $\delta _s$ is calculated as (e.g. Finnigan Reference Finnigan2000)

(3.6) \begin{equation} \delta _s = \frac {max (\langle u\rangle _{xy})-min (\langle u\rangle _{xy})}{max (\left |\frac {\partial {\langle u\rangle _{xy}}}{\partial {z}}\right |)}, \end{equation}

i.e. the maximum velocity difference in the horizontally averaged profile divided by the maximum shear. The dominant wavenumber $k_{x,peak}$ is estimated from the peak in the $x$ -direction wavenumber spectrum. Here we investigate the KH wavelength scaling (Sutherland Reference Sutherland and Rahman2005) by examining the following relationship:

(3.7) \begin{equation} k_{x,peak}\delta _s \approx 1. \end{equation}

The scaling in (3.7) generally applies across all simulations spanning a range of canopy densities (figure 5 $a$ ). Increased canopy density typically results in increased shear layer thickness and consequently longer wavelengths in the induced KH-type patterns (cases C0.1 to C3). This dependence on canopy density is also apparent when comparing the side views in figures 3 $(c)$ and 3 $(d)$ . Additionally, a weaker geostrophic current results in a thinner shear layer, leading to shorter KH wavelengths (cases C3 and C3-WG, figure 5 $a$ ; also see figures 3 $d$ and 3 $e$ ). Moreover, each simulation with a specific canopy density exhibits variations in KH wavelength (figure 5 $a$ ), which relates to the temporal evolution of the canopy shear layer and KH-type patterns, as further discussed in the following paragraph.

While a representative snapshot has been shown in a previous section in figure 3, the KH-type patterns also evolve over time. Here we take case C3 as an example to illustrate the evolving features in the canopy shear layer (figure 6). The simulations are initiated with a uniform velocity profile of $u=u_g$ . As the canopy drag decelerates the ocean current, a shear layer rapidly develops beneath the canopy (e.g. for $tf\lt 0.5$ ) and the shear layer thickness remains relatively stable after $tf=0.5$ . Note that the growth of the shear layer thickness occurs within a relatively short time (figure 6 $a$ ) compared with the deepening of the OML (figure 6 $b$ ). At an earlier stage when the shear layer thickness is smaller (e.g. $tf=0.25$ in figure 6 $a$ ), the dominant KH wavelength is correspondingly shorter (comparing figure 6 $c$ with figure 3 $d$ ), as predicted by the scaling relationship in (3.7). As time progresses, the KH wavelength increases with the thickening of the shear layer (see the movie in the Supplementary material). Eventually, the shear layer becomes mismatched with the stratification due to the continuous deepening of the thermocline, resulting in a regime that deviates from classical KH instability.

Figure 6. Time evolution of the velocity and temperature field in case C3. $(a)$ Vertical profiles of horizontally averaged streamwise velocity $\langle \overline {u}\rangle _{xy}$ , at $tf=0.3$ , 0.5 and 9. $(b)$ Vertical profiles of horizontally averaged temperature $\langle \overline {\theta }\rangle _{xy}$ , at every $tf=1$ . $(c,d)$ Side views of the cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$ at $tf=0.3$ and $tf=3$ , representing stages earlier and later than figure 3 $(d)$ , respectively.

It is worth noting that the canopy-generated structures may comprise of a combination of KH waves and Holmboe waves (Carpenter et al. Reference Carpenter, Balmforth and Lawrence2010; Zagvozkin et al. Reference Zagvozkin, Vorobev and Lyubimova2019), since their distinction is often less clear in an asymmetric stratified profile (such as the one used in this study). The KH instability usually occurs in the presence of relatively strong velocity shear that satisfies $Ri_g\lt 0.25$ , whereas Holmboe instability can arise from the interaction between shear and stratification even when $Ri_g\gt 0.25$ . The KH waves typically travel at a phase speed close to the current speed of the inflection point ( $u_{infl}$ , here defined as the point where $|\partial u/\partial z|$ reaches its maximum). However, the propagating phase speed of Holmboe waves can be significantly different from $u_{infl}$ (e.g. Alexakis Reference Alexakis2007; Carpenter et al. Reference Carpenter, Balmforth and Lawrence2010).

The relationship between $u_{\textit{infl}}$ and the dominant phase speed $c_{p,peak}=\omega _{peak}/k_{x,peak}$ is examined in figure 5 $(b)$ . Although $c_{p,peak}$ is in general close to $u_{infl}$ across all simulations, deviations are evident, particularly for the highest canopy density (case C3), implying the presence of Holmboe waves. As time progresses, while the shear layer remains attached to the canopy bottom, the thermocline continues to deepen as the OML depth increases. Consequently, at later stages, the thermocline no longer coincides with the shear layer, which may lead to patterns resembling Holmboe waves (figure 6 $d$ ). The magnitude of these canopy-generated structures eventually diminishes as the mixed layer continues deepening and the thermocline becomes displaced from the shear layer.

Distinguishing between KH modes and Holmboe modes is challenging in our non-idealized scenario, and the cases presented here likely exhibit characteristics of both. We therefore refer to all these canopy-induced patterns as KH-type. In the following analyses, we primarily focus on a representative time around $tf=1$ when the KH-type structures and internal waves are most prominent (the same time as shown in figure 3). Cases without a canopy will also be analysed at the same representative time to serve as comparisons.

Figure 7. $(a)$ The TKE profiles in four representative cases, normalized by the square of the friction velocity $u_*$ . $(b)$ The vertical component of TKE. $(c)$ Ratio of the vertical component to the total TKE. $(d)$ Skewness of the vertical turbulent velocity. The results correspond to around $tf=1$ , same as the time when the side views in figure 3 are shown. Note that the ratio and skewness are not shown for regions where TKE $/u_*^2\lt 1$ . Horizontal dotted lines represent the canopy base depth.

3.4. Turbulence

We first calculate the turbulence kinetic energy (TKE) for four representative cases (NC, NC-NSW, C0.3 and C3), where

(3.8) \begin{equation} \textrm {TKE} = \frac {1}{2}\left ( \overline {u'^2}+\overline {v'^2}+\overline {w'^2} \right ). \end{equation}

The cases with a canopy display a noticeable increase in TKE at the bottom edge of the canopy (figure 7 $a$ ), corresponding to the canopy shear layer turbulence. Note that while the calculated TKE accounts for high-frequency turbulent motions induced by shear and KH instability, it does not include the energy directly associated with the KH-type structures, as these are classified as lower-frequency motions according to the definition in (3.1). The lower-frequency motions will be examined in the following section.

Langmuir turbulence (case NC, no canopy, with the Stokes drift) is characterized by a relatively higher energy content in the vertical component (figure 7 $c$ ), and similar turbulence anisotropy has been reported by McWilliams et al. (Reference McWilliams, Sullivan and Moeng1997). In comparison, canopy shear layer turbulence typically exhibits less turbulence anisotropy (cases C0.3 and C3, e.g. at around $z/h_0=-0.8$ to $-1.5$ ).

The skewness of the vertical turbulence component, defined as $\overline {{w'}^3}/\overline {{w'}^2}^{3/2}$ , is examined to illustrate the various turbulence characteristics (figure 7 $d$ ). Langmuir turbulence (case NC) is featured by its distinctively negative skewness of vertical velocity (McWilliams et al. Reference McWilliams, Sullivan and Moeng1997). This negative skewness is attributed to stronger downward motions confined within narrower regions of Langmuir circulations, compared with broader and weaker upward motions. For shear layer turbulence generated at the canopy bottom (cases C0.3 and C3), the vertical velocity skewness is positive, e.g. in the lower-half of the canopy. This positive skewness indicates the dominance of sweep events that brings high-momentum fluid into the canopy, consistent with those found in classical canopy flows (Raupach & Thom Reference Raupach and Thom1981; Katul et al. Reference Katul, Kuhn, Schieldge and Hsieh1997; Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004). It is also worth noting that the skewness sign in the suspended canopy is opposite to that of submerged benthic canopies; here, the sweep events are characterized by intensified positive vertical velocities into the canopy. Additionally, in the upper-half of the canopy where the influence of Stokes drift is greater, negative skewness is observed, which represents a characteristic of canopy-modulated Langmuir turbulence (Bo et al. Reference Bo, McWilliams, Yan and Chamecki2024b ).

Figure 8. Terms in the TKE budget (3.9) for case NC (no canopy, with the Stokes drift) and case NC-NSW (no canopy, no Stokes drift) at around $tf=1$ : $(a)$ shear production associated with the streamwise velocity; $(b)$ Stokes production; $(c)$ dissipation; $(d)$ buoyancy production; $(e)$ pressure work (vertical divergence of pressure flux); $(e)$ pressure flux. Only the dominant terms are shown. Note the differences in the horizontal axis range.

Further analysis of the TKE budget is conducted to understand the mechanisms of turbulence generation,

(3.9) \begin{eqnarray} \frac {\text {D}}{\text {D}t} \textrm {TKE} = - \overline {u_i' u_j'}\frac {\partial {\overline {u_i}}}{\partial {x_j}} - \overline {u'w'}\frac {\partial {u_s}}{\partial {z}} - \epsilon - \overline {u'_{i} F'_{D,i}} + \alpha g\overline {w'\theta '} - \frac {1}{\rho }\frac {\partial {}}{\partial {x_j}}\overline {u'_j p'} \nonumber \\ - \frac {1}{2}\frac {\partial {}}{\partial {x_j}}\overline {u'_iu'_iu'_j} + \frac {\partial {}}{\partial {x_j}}\overline {u'_i\tau '_{ij}}, \\\nonumber\end{eqnarray}

where the material derivative is defined as $\text {D}/\text {D}t = \partial /\partial t + \overline {u_j}\partial /\partial x_j + u_s\partial /\partial x$ . Note that we have used Einstein notation. The first and second terms on the right-hand side represent shear production and Stokes production, respectively. The third term $\epsilon$ denotes SGS dissipation, and the fourth term stands for canopy drag dissipation. The fifth term represents buoyancy production, followed by the sixth term that denotes pressure work. The last two terms represent TKE transport associated with the resolved stress and SGS stress.

In the two cases without a canopy (figures 8 $a$ and 8 $b$ ), shear production is the primary source of wind-driven shear turbulence (case NC-NSW), while Stokes production predominates for Langmuir turbulence (case NC). The shear production and Stokes production terms are generally balanced by the dissipation term (figure 8 $c$ ), with minimal influences from other source or sink terms or TKE transport.

Figure 9. Terms in the TKE budget (3.9) for case C3 (with the highest canopy density): red lines, at $tf=1$ ; blue lines, at $tf=0.3$ ; black lines, at $tf=9$ ; $(a,b)$ shear production associated with the streamwise velocity and cross-stream velocity, respectively; $(c)$ Stokes production; $(d)$ dissipation; $(e)$ canopy drag dissipation; $(f)$ TKE transport; $(g)$ buoyancy production; $(h)$ pressure work (vertical divergence of pressure flux); $(i)$ pressure flux. Note the differences in the horizontal axis range. Horizontal dotted lines represent the canopy base depth.

Negative buoyancy production is evident near the thermocline (figure 8 $d$ ), which indicates the energy transfer from TKE to potential energy, contributing to mixing and OML deepening. Langmuir turbulence induces stronger buoyancy fluxes than wind-driven shear turbulence, consistent with the greater increase in OML depth in case NC compared with case NC-NSW (figure 2 $a$ ). Moreover, negative pressure fluxes occur below the thermocline in the simulation with Langmuir turbulence (figure 8 $f$ ), corresponding to the observed downward-propagating internal waves (figure 3 $a$ ). The vertical divergence of pressure fluxes (the pressure work term) (figure 8 $e$ ) is in general much smaller than the Stokes production and dissipation terms, suggesting that internal wave radiation has a small influence on TKE. Nevertheless, pressure work has a comparable magnitude to buoyancy production, and this underscores the potential importance of internal waves for mixed layer deepening (Taylor & Sarkar Reference Taylor and Sarkar2007).

In the presence of a canopy (case C3, $tf=1$ , red lines in figure 9), a notable peak in shear production occurs near the canopy bottom edge where KH-type patterns emerge. Stokes production has a small magnitude compared with shear production, and this explains why removing the Stokes drift leads to minimal changes to the OML depth when a canopy is present (comparing cases C3 and C3-NSW). The TKE transport term indicates a loss of TKE in the canopy shear layer where production peaks and a gain within the canopy (figure 9 $e$ ), consistent with findings from classical canopy flow studies (e.g. Finnigan Reference Finnigan2000). In addition, a negative peak in pressure work is observed at the canopy base, corresponding to the canopy-induced KH-type patterns and radiated internal waves.

Negative buoyancy production is evident in the canopy shear layer, contributing to mixing at the thermocline (figure 9 $g$ ). The buoyancy production term is small compared with the shear production term. This buoyancy mixing associated with the canopy and KH-type patterns is consistent with the principles of the classical Osborn–Cox model (Osborn & Cox Reference Osborn and Cox1972), wherein sharp gradients of oceanic microstructures usually play a pivotal role in vertical mixing. It is also worth noting that buoyancy production here can be around an order of magnitude larger than that without a canopy (case NC), thus leading to more pronounced OML deepening in the presence of a canopy (figure 2 $a$ ).

While the above analysis focuses on $tf=1$ , the TKE budget at earlier ( $tf=0.3$ ) and later ( $tf=9$ ) stages is also presented in figure 9 (blue and black lines). In general, the vertical distribution of the shear production, dissipation and TKE transport terms is similar across different times. However, the magnitudes of shear production and dissipation are greater at later times, as turbulence develops with time. Additionally, at the later stage, the peaks of buoyancy production and pressure flux occur at greater depths (black lines in figures 9 $g$ and 9 $i$ ), because the thermocline deepens. Moreover, negative buoyancy production is not observed at the earlier stage (red line in figure 9 $g$ ). Instead, positive values of buoyancy production occur, likely indicating convectively unstable regimes, as dense filaments are advected upward in the initial stages of KH instability.

3.5. Internal waves beneath the OML

Alongside the deepening of the OML, internal waves can be generated in the underlying stratification as flow processes in the OML interact with the thermocline. As outlined in a previous section, Langmuir turbulence (case NC, no canopy) and the canopy-induced KH-type structures (cases C0.1–C3) can both lead to internal wave radiation beneath the OML (figure 3). In this section, we further analyse the properties of internal waves. Here internal wave computations are based on vertical velocity gradients $\partial w/\partial z$ following Taylor & Sarkar (Reference Taylor and Sarkar2007), as $\partial w/\partial z$ reveals internal wave patterns more clearly. The spectra of the vertical velocity $w$ generally provide consistent information with those of $\partial w/\partial z$ .

Figure 10. Side-view and cross-sectional plots of vertical velocity gradients at around $tf=1$ in case NC. $(a)$ The cross-stream averaged vertical velocity gradient $\partial \langle w\rangle _{y}/\partial z$ in the $x-z$ plane. $(b)$ The vertical velocity gradient $\partial w/\partial z$ at $y=0$ in the $x-z$ plane. $(c)$ The $x$ -direction averaged vertical velocity gradient $\partial \langle w\rangle _{x}/\partial z$ in the $y-z$ plane. $(d)$ The vertical velocity gradient $\partial w/\partial z$ at $x=0$ in the $y-z$ plane.

Figure 11. Spectra of vertical velocity gradients $\partial w/\partial z$ in case NC, versus $x$ - or $y$ -wavenumber or frequency. $(a{-}c)$ Spectra of the $x$ - or $y$ -direction averaged results. $(e,f)$ The $x$ - or $y$ -direction average of the spectra at every $x$ - or $y$ -location. Dashed lines represent the average between $z/h_0=-1$ to $0$ (within the OML), and solid lines represent the average between $z/h_0=-3$ and $-2$ (below the thermocline). Vertical red lines indicate the buoyancy frequency $N$ . Note that the frequency in the spectral plot is generally much higher than $N$ due to the Doppler shift caused by the geostrophic current. In the text, intrinsic internal wave frequencies are analysed after removing the Doppler shift effect.

In case NC, the side-view plots (figures 10 $a$ and 10 $b$ ) indicate propagating internal waves in the $x-z$ plane. Internal wave patterns are more discernible when averaging in the $y$ -direction, whereas examining a slice at any specific $y$ -location tends to obscure these patterns. It is worth noting that internal waves also propagate in the $y$ -direction, as illustrated in the cross-sectional plots (figures 10 $c$ and 10 $d$ ). Plotting $\partial u/\partial z$ instead of $\partial w/\partial z$ will enhance the clarity of internal wave patterns in the $x-z$ plane, with less discrepancy between the $y$ -direction average and a specific $y$ -location. However, $\partial u/\partial z$ does not reveal any wave patterns in the $y$ -direction, and we therefore use $\partial w/\partial z$ in our analysis.

The dispersion relation for non-rotating linear internal waves is given by

(3.10) \begin{equation} \omega _i = N\cos \Theta , \end{equation}

where $\omega _{i}$ is the intrinsic frequency, $N$ is the buoyancy frequency as defined in (2.7) and $\Theta$ is the angle between the wavenumber vector and the horizontal axis. While near-inertial oscillations can also propagate downwards below the OML (Bell Reference Bell1978), the observed wave frequency is close to $N$ and significantly higher than the Coriolis frequency $f$ , and we have thus neglected the influence of planetary rotation.

In the $x$ -direction wavenumber spectrum of case NC (figure 11 $a$ ), a peak at $k_xh_0=3.6$ is evident below the thermocline, representing the dominant wavenumber of the downward propagating internal waves, and this contrasts with the broadband turbulence observed in the surface boundary layer. Correspondingly, a peak frequency of around $\omega _{peak}/f=300$ is observed in the time spectrum (figure 11 $c$ ). The intrinsic internal wave frequency $\omega _{i}$ is estimated by accounting for the Doppler shift:

(3.11) \begin{equation} \omega _i = \omega _{peak} - k_{x,peak} u_g. \end{equation}

Here the intrinsic frequency is $\omega _{i}/f\approx 20$ , which is smaller than the buoyancy frequency $N/f=44$ , suggesting a wave propagation angle $\Theta$ of around $60^\circ$ . This in general agrees with the angle estimated from

(3.12) \begin{equation} \Theta = \cos ^{-1}\left (\sqrt {\frac {k_{x,peak}^2+k_{y,peak}^2}{k_{x,peak}^2+k_{y,peak}^2+k_{z,peak}^2}}\right ), \end{equation}

so the dispersion relation of freely propagating internal waves is satisfied. Nevertheless, the non-monotonic nature of the radiated internal waves may lead to uncertainty when applying the dispersion relation in (3.10). Additionally, we note that while the wavenumber and frequency peaks of internal waves are evident in the $x$ - or $y$ -direction averaged field, they are less clear when we calculate spectra based on slices at any specific $x$ - or $y$ -location (figure 11 $d{-}f$ ).

The presence of a narrow range of dominant wavenumbers and frequencies in internal waves is consistent with the selection mechanism proposed by Gibson et al. (Reference Gibson, Bondur, Norris Keeler and Leung2006) and Taylor & Sarkar (Reference Taylor and Sarkar2007). While energy is distributed across a wide range of frequencies in boundary layer turbulence, internal waves with either low or high frequencies exhibit slow vertical propagation speeds or rapid viscous decay rates. Consequently, an intermediate dominant frequency emerges as waves propagate beneath the OML, although they are excited by broadband turbulence.

The analysis above has focused on internal waves associated with Langmuir turbulence (case NC, no canopy, with the Stokes drift). Wind-driven shear turbulence (case NC-NSW, no canopy, no Stokes drift) is expected to generate internal waves through a similar mechanism. However, internal wave energy is much weaker in case NC-NSW compared with case NC, with less discernible wave patterns (comparing figures 3 $a$ and 3 $b$ ) and a less distinctive peak in the wavenumber spectrum (not shown). We attribute this difference to the fact that Langmuir turbulence, compared with wind-driven shear turbulence, results in stronger mixing in the OML (McWilliams et al. Reference McWilliams, Sullivan and Moeng1997; Li & Garrett Reference Li and Garrett1997) and more deeply penetrating vertical velocities that interact with the stratification (Polton & Belcher Reference Polton and Belcher2007), which thus leads to more pronounced internal waves (Chini & Leibovich Reference Chini and Leibovich2003; Polton et al. Reference Polton, Smith, MacKinnon and Tejada-Martínez2008).

Figure 12. Side views of vertical velocity gradients at around $tf=1$ in case C3. $(a)$ The cross-stream average vertical velocity gradient $\partial \langle w\rangle _{y}/\partial z$ . $(b)$ The vertical velocity gradient $\partial w/\partial z$ at $y=0$ .

In the presence of a canopy, internal wave patterns below the thermocline are primarily associated with the KH-type structures in the canopy shear layer (figure 12). A distinct peak at $k_x h_0=1.5$ is noticeable in the $x$ -direction wavenumber spectrum, both within the surface boundary layer and below the thermocline (figures 13 $a$ and 13 $d$ ). This corresponds to an observed peak frequency of around $\omega _{peak}/f=50$ (figures 13 $c$ and 13 $f$ ). In contrast to the Langmuir turbulence simulation where internal waves propagate in both the $x$ - and $y$ -directions, the simulation with a canopy exhibits minimal variations in the $y$ -direction, with no apparent discrepancy between the cross-stream averaged results and any specific $y$ -location (e.g. comparing figures 12 $a$ and 12 $b$ ). Multiple peaks in the $y$ -direction wavenumber spectrum are observed (figure 13 $e$ ), which likely represent cross-stream structures associated with secondary instability in KH and Holmboe waves (e.g. Smyth & Peltier Reference Smyth and Peltier1991; Mashayek & Peltier Reference Mashayek and Peltier2013; Mashayek et al. Reference Mashayek, Caulfield and Peltier2013; VanDine et al. Reference VanDine, Pham and Sarkar2021). Our model is unable to fully resolve these secondary structures. However, it is important to note that the energy associated with cross-stream structures is generally much lower than that of the streamwise structures. We therefore primarily focus on the streamwise structures that are more dominant in OML mixing. In the simulation with doubled horizontal resolution, consistent peaks are observed in the $y$ -direction wavenumber spectrum.

Figure 13. Spectra of vertical velocity gradients $\partial w/\partial z$ in case C3, versus $x$ - or $y$ -wavenumber or frequency (a similar set of panels to figure 11). $(a{-}c)$ Spectra of the $x$ - or $y$ -direction averaged results. $(e,f)$ The $x$ - or $y$ -direction average of the spectra at every $x$ - or $y$ -location. Dashed lines represent the average between $z/h_0=-1$ to $0$ (within the OML), and solid lines represent the average between $z/h_0=-3$ and $-2$ (below the thermocline).

The dominant intrinsic frequency of internal waves is estimated to be $\omega _{i}/f\approx -50$ by considering the Doppler shift in (3.11). The negative value indicates that internal waves would propagate in the $-x$ -direction when observed in a coordinate frame moving with the geostrophic current $u_g$ . This is consistent with analysis in the previous section (figure 5 $b$ ), where the phase speed of KH-type structures is close to the current speed at the inflection point ( $u_{infl}\lt u_g$ ). The magnitude of the intrinsic frequency $\omega _i$ is close to the buoyancy frequency $N$ , suggesting horizontal wave propagation according to the dispersion relation (3.10). This aligns with observations from the side-view plots (figure 12), where the wave phase lines are nearly vertical. In addition to the dominant internal wave patterns that primarily propagate horizontally (with a peak wavenumber and frequency associated with the KH-type structures), internal waves with other frequencies are also excited in the canopy simulation (figure 12). These waves propagate in the vertical direction, as indicated by their inclined phase lines.

Note that a sponge layer is applied at the bottom boundary to avoid internal wave reflection. We also conduct a simulation with doubled domain depth to examine any potential reflection (case C3-DEEP, details provided in appendix A). The case with the highest canopy density is chosen because it results in the deepest OML and KH-type structures closest to the bottom boundary, making it most susceptible to the influence of internal wave reflection. Overall, the wave spectra show negligible differences between cases C3 and C3-DEEP, indicating that the sponge layer effectively prevents wave reflection (appendix A).

3.6. The OML deepening and nutrient entertainment

The presence of a canopy can significantly deepen the OML, and both canopy-induced turbulence and lower frequency motions associated with KH-type structures have the potential to enhance mixing. Moreover, since nutrient transport in the ocean is typically influenced by OML dynamics, the presence of a canopy may thus substantially impact the vertical distribution of nutrients. In this section, we quantify the vertical transport processes associated with the canopy. The velocity has been decomposed into turbulence ( $\textbf {u}'$ ) and low-frequency motions ( $\overline {\textbf {u}}$ , averaged every 75 s) based on (3.1). We further decompose the low-frequency velocity as

(3.13) \begin{equation} \overline {\textbf {u}} = \widetilde {\overline {\textbf {u}}} + \overline {\textbf {u}}^l, \end{equation}

where the tilde denotes the overall time average, and the superscript $l$ denotes low-frequency fluctuations around the time average. The complete flow decomposition is thus given by

(3.14) \begin{equation} \textbf {u} = \widetilde {\overline {\textbf {u}}} + \overline {\textbf {u}}^l + \textbf {u}'. \end{equation}

The term $\textbf {u}'$ denotes fluctuations associated with turbulence, while other lower-frequency fluctuations, such as the KH-type structures and internal waves, are represented by $\overline {\textbf {u}}^l$ . Likewise, any scalar $\phi$ , e.g. temperature and nutrient concentration, can be decomposed into

(3.15) \begin{equation} \phi = \widetilde {\overline {\phi }} + \overline {\phi }^l + \phi '. \end{equation}

The total time-averaged vertical flux is thus written as

(3.16) \begin{equation} \widetilde {\overline {w\phi }} = \widetilde {\overline {w}^l\overline {\phi }^l} + \widetilde {\overline {w'\phi '}}. \end{equation}

Note that there is no vertical flux driven by the mean flow because $\widetilde {\overline {w}}=0$ .

Figure 14. Vertical fluxes of momentum and temperature for case NC (black) and case C3 (red), time-averaged from $tf=0$ to $tf=9$ and horizontally averaged. $(a,b)$ Turbulent fluxes of momentum and temperature. $(c,d)$ Fluxes driven by lower-frequency motions. Fluxes are normalized by $u_*$ and $\Delta \theta =1\,\rm K$ .

The momentum and temperature fluxes of cases NC and C3 are compared in figure 14. The intensified turbulence within the canopy shear layer considerably enhances vertical temperature mixing compared with the case without a canopy. Moreover, lower-frequency motions also substantially contribute to temperature mixing, with a comparable magnitude to the contribution from turbulence. This underscores the dual role of the canopy, which generates turbulence and enhances turbulent mixing, while also directly leading to mixing through the lower-frequency KH-type structures. The momentum flux is primarily dominated by turbulence in the canopy shear layer, with lower-frequency velocities exhibiting weaker correlation than turbulence. As a contrast, in the absence of a canopy (case NC), although internal waves are generated, these lower-frequency motions have a negligible influence on temperature and momentum mixing. This indicates that the deepening of the OML is primarily driven by Langmuir turbulence rather than the radiated internal waves.

Figure 15. $(a)$ Nutrient profiles for case NC (black) and case C3 (red), at $tf=9$ . The thin grey line shows the initial nutrient profile at $tf=0$ . $(b,c)$ Vertical nutrient fluxes driven by turbulence and lower-frequency motions, respectively. These fluxes are time-averaged from $tf=0$ to $tf=9$ and horizontally averaged. Fluxes are normalized by $u_*$ and $\Delta C_N=10000$ $\unicode{x03BC}\rm mol\,m^-{^3}$ .

Furthermore, the presence of a canopy also significantly enhances vertical nutrient transport (figure 15). Similar to temperature fluxes, both turbulence and lower-frequency KH structures are important factors in enhancing nutrient entrainment into the OML. The nutrient concentration near the sea surface is increased by more than twofold compared with the case without a canopy. Considering that the present study on the canopy flow is motivated by macroalgal farming, the increased nutrient concentration fosters a favourable condition for macroalgal growth. This provides a self-sustaining solution for nutrient supply and farm maintenance through continuous passive entrainment. Note that a secondary peak in nutrient flux occurs near the sea surface in the case with the canopy (at $z/h_0=-0.3$ in figure 15 $b$ ). A similar secondary peak is found in the temperature flux (figure 14 $b$ ). This is because the KH structures advect strong gradients from the thermocline to near the sea surface (see figure 3 $d$ ), resulting in locally increased fluxes as these gradients are mixed by turbulence within the canopy.

4. Discussion and conclusion

In this study, we investigate the impact of a suspended canopy on OML dynamics by using LES. A series of simulations are compared, with and without a suspended canopy, as well as with and without surface gravity waves (Stokes drift). Overall, the presence of a canopy generates intensified turbulence and KH-type structures in the canopy shear layer, leading to significant OML deepening and the radiation of internal waves beneath the OML.

A first comparison is made between two cases without a canopy, with and without surface waves (Stokes drift), respectively. In the absence of the Stokes drift, wind-driven shear turbulence is generated by the shear production term in the TKE budget. With the inclusion of the Stokes drift, Langmuir turbulence is produced due to the interaction between surface waves and wind-driven currents, with the Stokes production term dominating over shear production. Langmuir turbulence in general leads to stronger mixing and a greater OML depth compared with wind-driven shear turbulence.

In the presence of a suspended canopy, the depth of the OML is substantially increased compared with cases without a canopy. Higher canopy density typically results in more pronounced deepening of the OML. Turbulence is intensified in the shear layer attached to the bottom edge of the canopy. Moreover, KH-type patterns emerge in the canopy shear layer, with their dominant wavelength scaled with the thickness of the canopy shear layer. The travelling phase speed of the KH structures is generally close to the current speed of the inflection point, although deviations may occur, implying the concurrent presence of Holmboe waves.

Additionally, internal waves are generated as turbulent motions in the OML interact with the underlying stratification. Langmuir turbulence (no canopy, with the Stokes drift) induces internal waves that propagate downward with an inclined angle. With a canopy located in the OML, internal waves are radiated beneath the thermocline, irrespective of the presence of surface waves (Stokes drift). The dominant wavenumber and frequency align with those of the KH-type structures, while other wavenumbers and frequencies are also discernible in the internal wave patterns.

Along with the increased OML depth, vertical nutrient transport is substantially enhanced due to the presence of a canopy. Both high-frequency turbulence and low-frequency motions associated with KH-type patterns are crucial for driving nutrient fluxes. The enhanced vertical mixing due to the canopy suggests a mechanism for passive nutrient entrainment into the OML, which provides implications for sustainable offshore macroalgal farm development. Future field observations are needed to validate the present simulation results. However, it is important to recognize the challenges in obtaining field data from macroalgal farms, e.g. the limited availability of high-resolution measurements at relevant temporal and spatial scales. Our study reveals the key mechanisms governing canopy flow dynamics and highlights their significance for macroalgal growth, which could serve as a reference framework to guide future fieldwork.

This study focuses on simulations with a horizontally homogeneous suspended canopy, which effectively represents the interior region of a large farm. These findings can also be extended to a finite-length canopy, where horizontal edge effects are present (details in appendix B). Similar KH-type structures and internal waves are generated by the finite-length canopy. It is worth noting that in the finite-canopy simulation, fully developed OML turbulence is present in the inflow condition. This indicates that pre-existing OML turbulence does not inhibit the formation of KH structures in the canopy shear layer. In the finite-canopy simulation, the dominant wavenumber of KH structures changes downstream, because the shear layer thickness evolves as the current adjusts to the presence of the canopy. Subsequent studies could further investigate the streamwise development of the canopy shear layer and its impacts on KH structures in the finite-length canopy.

We have employed a constant and uniform surface wave field, while it is worth noting that canopies can cause wave attenuation and reduce Stokes drift. Wave amplitude decay is estimated to be less than 10 % over a distance of 2 km, based on the oceanic and canopy parameters relevant to this study (Yan et al. Reference Yan, McWilliams and Chamecki2021). Our findings apply to macroalgal farms that are sufficiently long to allow for fully developed canopy flow (with no significant spatial evolution, as seen in the finite-length canopy simulation above) but not longer than the distance over which significant wave attenuation occurs. Simulations of larger farms should include more complex canopy-wave interactions (e.g. McWilliams Reference McWilliams2023) to account for the downstream change in Stokes drift.

In this study, both the wave and current directions are aligned with the $x$ -direction. The results are expected to be similar for waves and currents approaching from different directions, as long as they remain aligned, due to the horizontally uniform configuration of the canopy (Bo et al. Reference Bo, McWilliams, Yan and Chamecki2024b ). While cases involving misaligned wave and current directions have been examined in the absence of a canopy, such conditions still warrant further investigation when a suspended canopy is present. In addition, this study focuses on monochromatic swell waves, but surface waves in reality are often broadband. Broadband waves exhibit a different vertical distribution of Stokes drift, typically being more surface-intensified, as shorter waves decay more quickly with depth compared with swell waves. Moreover, shorter waves may experience greater attenuation when they propagate through the canopy, requiring a more rigorous treatment of wave attenuation as discussed above. Nevertheless, the influences of surface waves on the canopy shear layer are in general small, because here the shear layer is deeper than the regions affected by surface waves. The inclusion of Stokes drift results in only minor variations in KH structures and OML deepening (comparing case C3 and C3-NSW in § 3.1). Therefore, incorporating more complex wave conditions in the model should not alter the key conclusions of this study.

Supplementary movie.

Supplementary movie is available at https://doi.org/10.1017/jfm.2025.68.

Acknowledgements.

The authors thank K. Davis and G. Pawlak for helpful discussions, and thank three anonymous reviewers for their insightful suggestions. The authors acknowledge high-performance computing support from Derecho (10.5065/qx9a-pg09), provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation.

Funding.

This project is funded by the U.S. Department of Energy ARPA-E MARINER (Macroalgae Research Inspiring Novel Energy Resources) program, grant number DE-AR0000920.

Declaration of interests.

The authors report no conflict of interest.

Appendix A. Dependence on domain depth

In this appendix, we examine a simulation (case C3-DEEP) with a doubled domain depth of 240 m to explore internal wave propagation in deeper water. The KH-type patterns and the OML depth are consistent with those in the shallower domain of case C3 (figure 16). Spectral analysis also reveals minimal discrepancies between the shallower and deeper simulations for $z/h_0\gt -3.5$ (figure 17). This further confirms that the results in the shallower domain are not affected by any internal wave reflection from the bottom boundary.

Figure 16. Side views of case C3-DEEP at around $tf=1$ . $(a)$ The cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$ . $(b)$ The cross-stream averaged vertical velocity gradients $\partial \langle w\rangle _{y}/\partial z$ . Black dotted lines represent the canopy base depth. Note that the sponge layer is not shown.

In the deeper water, e.g. $z/h_0\lt -3.5$ , the dominant wavenumber of internal waves decreases, and the phase lines exhibit an inclined angle (figures 16 and 17). This wavenumber difference between the shallower and deeper regions is likely associated with a transition in wave characteristics. In the shallower region, internal waves are primarily forced by turbulence in the canopy shear layer, whereas in deeper water, they evolve into freely radiating internal waves that are unaffected by turbulence. The altered dominant wavenumber, frequency and angle in the deeper region comply with the dispersion relation in (3.10). In addition, we note that the domain length must be an integer multiple of the wavelength due to periodic boundary conditions of the LES model, which may constrain the transition of internal wave characteristics in deeper water. Nevertheless, internal waves in deeper water are expected to have negligible influence on OML deepening, because there is no reflection from the bottom boundary.

Figure 17. Spectra of the $x$ - or $y$ -direction averaged vertical velocity gradients in case C3-DEEP, versus $(a)$ $x$ -wavenumber, $(b)$ $y$ -wavenumber or $(c)$ frequency. Solid lines represent the average between $z/h_0=-3$ to $-2$ , and dash–dotted lines represent the average between $z/h_0=-7$ to $-6$ .

Appendix B. A finite-length canopy

A simulation incorporating a finite-length canopy is conducted on a $800\times 200\times 120\,\rm m^3$ domain, which allows for the downstream development of canopy flow. A suspended canopy with a length of $400\,\rm m$ is positioned in the centre of the domain, with a base depth of $h_b=20\,\rm m$ . In the cross-stream direction, the canopy extends across the domain with the periodic boundary condition, effectively resulting in an infinite canopy width. The forcing conditions are in general identical to those in the simulations with a horizontally homogeneous interior canopy, with $u_g=0.2\,\rm m\,s^-{^1}$ .

We employ a precursor inflow method to simulate the spatially evolving flow associated with the finite-length canopy (Yan et al. Reference Yan, McWilliams and Chamecki2021; Bo et al. Reference Bo, McWilliams, Yan and Chamecki2024b ). Flow velocity and temperature fields at the upstream boundary are provided by a separate precursor simulation, which is conducted under identical conditions but without the canopy. A fringe region with a length of 100 m is applied at the downstream end of the domain of the canopy simulation. Within the fringe region, the flow field smoothly transitions towards the inflow conditions obtained from the precursor simulation at the end of each time step. This precursor inflow method prevents canopy wakes from cycling back through the periodic boundary conditions.

Figure 18. A side-view plot of the cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$ in the finite-length canopy simulation. The dotted rectangle shows the extent of the canopy.

A side view of the finite-length canopy is depicted in figure 18. On the upstream boundary, the OML depth of the inflow is $h_i = 16\,\rm m$ , with an underlying stratification of uniform $\partial \theta /\partial z|_\infty =0.01\,\rm K\,m^-{^1}$ . Note that in the simulation presented here, the inflow OML depth is shallower than the canopy base depth $h_b=20\,\rm m$ . The streamwise flow is decelerated within the canopy due to the drag force, resulting in downwelling near the canopy leading edge as a result of mass conservation. The thermocline tilts downward due to the canopy-driven downwelling. The KH-type structures emerge in the shear layer below the finite-length canopy.

A key difference between the finite-length canopy and the horizontally homogeneous canopy lies in the downstream evolution of KH structures. Below the finite-length canopy, the characteristic wavelength of KH structures is shorter in the upstream region and elongates downstream. This occurs because the shear layer thickness increases gradually from the leading edge of the canopy, and the KH wavelength typically scales with the canopy shear layer thickness, as discussed in the main text. The OML deepens to around 40 m as a result of the combined effects of downwelling and shear mixing induced by the canopy. The OML depth then recovers to around 30 m due to upwelling downstream from the canopy trailing edge. Additionally, internal wave patterns are observed in the underlying stratification, e.g. at $z/h_0=-3$ , $x/h_0=15$ to 25.

Another distinction between the finite-length canopy and the horizontally homogeneous canopy is the presence of leading-edge downwelling. This downwelling, resulting from the flow adapting to the finite-length canopy, can vertically advect the underlying stratification and lead to mismatches with the canopy shear layer. In other simulations where the inflow OML depth is $h_i = 26\,\rm m$ , no apparent KH structures are observed, because the downwelling displaces the stratification to a depth well below the canopy shear layer. The KH structures only emerge when the OML depth of the inflow is less than the canopy base depth (e.g. in figure 18 the inflow OML depth is $h_i = 16\,\rm m$ ).

References

Alexakis, A. 2007 Marginally unstable Holmboe modes. Phys. Fluids 19 (5), 054105.CrossRefGoogle Scholar
Arzeno-Soltero, I.B., Saenz, B.T., Frieder, C.A., Long, M.C., DeAngelo, J., Davis, S.J. & Davis, K.A. 2023 Large global variations in the carbon dioxide removal potential of seaweed farming due to biophysical constraints. Commun. Earth Environ. 4 (1), 185.CrossRefGoogle Scholar
Belcher, S.E., et al. 2012 a A global perspective on Langmuir turbulence in the ocean surface boundary layer. Geophys. Res. Lett. 39 (18), L18605.CrossRefGoogle Scholar
Belcher, S.E., Harman, I.N. & Finnigan, J.J. 2012 b The wind in the willows: flows in forest canopies in complex terrain. Annu. Rev. Fluid Mech. 44 (1), 479504.CrossRefGoogle Scholar
Bell, T.H. 1978 Radiation damping of inertial oscillations in the upper ocean. J. Fluid Mech. 88 (2), 289308.CrossRefGoogle Scholar
Bo, T., McWilliams, J.C., Frieder, C.A., Davis, K.A. & Chamecki, M. 2024 a Nutrient replenishment by turbulent mixing in suspended macroalgal farms. Geophys. Res. Lett. 51 (14), e2024GL109128.CrossRefGoogle Scholar
Bo, T., McWilliams, J.C., Yan, C. & Chamecki, M. 2024 b Langmuir turbulence in suspended kelp farms. J. Fluid Mech. 985, A11.CrossRefGoogle Scholar
Bou-Zeid, E., Meneveau, C. & Parlange, M. 2005 A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows. Phys. Fluids 17 (2), 025105.CrossRefGoogle Scholar
Carpenter, J., Balmforth, N. & Lawrence, G. 2010 Identifying unstable modes in stratified shear layers. Phys. Fluids 22 (5), 054104.CrossRefGoogle Scholar
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), 13381371.CrossRefGoogle Scholar
Chini, G. & Leibovich, S. 2003 Resonant Langmuir-circulation–internal-wave interaction. part 1. internal wave reflection. J. Fluid Mech. 495, 3555.CrossRefGoogle Scholar
Claustre, H., Johnson, K.S. & Takeshita, Y. 2020 Observing the global ocean with biogeochemical-argo. Annu. Rev. Mar. Sci. 12 (1), 2348.CrossRefGoogle ScholarPubMed
Craik, A.D. & Leibovich, S. 1976 A rational model for Langmuir circulations. J. Fluid Mech. 73 (3), 401426.CrossRefGoogle Scholar
D’Asaro, E.A. 2001 Turbulent vertical kinetic energy in the ocean mixed layer. J. Phys. Oceanogr. 31 (12), 35303537.2.0.CO;2>CrossRefGoogle Scholar
D’Asaro, E.A. 2014 Turbulence in the upper-ocean mixed layer. Annu. Rev. Mar. Sci. 6 (1), 101115.CrossRefGoogle ScholarPubMed
Dohan, K. & Sutherland, B. 2005 Numerical and laboratory generation of internal waves from turbulence. Dyn. Atmos. Oceans 40 (1-2), 4356.CrossRefGoogle Scholar
Finnigan, J. 2000 Turbulence in plant canopies. Annu. Rev. Fluid Mech. 32 (1), 519571.CrossRefGoogle Scholar
Frieder, C.A. et al. 2022 A macroalgal cultivation modeling system (macmods): evaluating the role of physical-biological coupling on nutrients and farm yield. Front. Mar. Sci. 9, 752951.CrossRefGoogle Scholar
Gibson, C.H., Bondur, V.G., Norris Keeler, R. & Leung, P.T. 2006 Energetics of the beamed zombie turbulence maser action mechanism for remote detection of submerged oceanic turbulence. J. Appl. Fluid Mech. 1 (1), 1142.Google Scholar
Howard, L.N. 1961 Note on a paper of John W. Miles. J. Fluid Mech. 10 (4), 509512.CrossRefGoogle Scholar
Jackson, G.A. 1997 Currents in the high drag environment of a coastal kelp stand off California. Cont. Shelf Res. 17 (15), 19131928.CrossRefGoogle Scholar
Kaminski, A. & Smyth, W. 2019 Stratified shear instability in a field of pre-existing turbulence. J. Fluid Mech. 862, 639658.CrossRefGoogle Scholar
Katul, G., Kuhn, G., Schieldge, J. & Hsieh, C.-I. 1997 The ejection-sweep character of scalar fluxes in the unstable surface layer. Boundary-Layer Meteorol. 83 (1), 126.CrossRefGoogle Scholar
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), 23812400.CrossRefGoogle Scholar
Kundu, P.K., Cohen, I.M. & Dowling, D.R. 2015 Fluid Mechanics. Academic press.Google Scholar
Leibovich, S. 1983 The form and dynamics of Langmuir circulations. Annu. Rev. Fluid Mech. 15 (1), 391427.CrossRefGoogle Scholar
Li, M. & Garrett, C. 1997 Mixed layer deepening due to Langmuir circulation. J. Phys. Oceanogr. 27 (1), 121132.2.0.CO;2>CrossRefGoogle Scholar
Mashayek, A., Caulfield, C. & Peltier, W. 2013 Time-dependent, non-monotonic mixing in stratified turbulent shear flows: implications for oceanographic estimates of buoyancy flux. J. Fluid Mech. 736, 570593.CrossRefGoogle Scholar
Mashayek, A. & Peltier, W. 2013 Shear-induced mixing in geophysical flows: does the route to turbulence matter to its efficiency? J. Fluid Mech. 725, 216261.CrossRefGoogle Scholar
McWilliams, J.C. 2023 Surface waves and currents in aquatic vegetation. J. Fluid Mech. 958, A14.CrossRefGoogle Scholar
McWilliams, J.C., Sullivan, P.P. & Moeng, C.-H. 1997 Langmuir turbulence in the ocean. J. Fluid Mech. 334, 130.CrossRefGoogle Scholar
Miles, J.W. 1961 On the stability of heterogeneous shear flows. J. Fluid Mech. 10 (4), 496508.CrossRefGoogle Scholar
Monismith, S., Alnajjar, M., Daly, M., Valle-Levinson, A., Juarez, B., Fagundes, M., Bell, T. & Woodson, C.B. 2022 Kelp forest drag coefficients derived from tidal flow data. Estuar. Coast. 45 (8), 24922503.CrossRefGoogle Scholar
Munroe, J.R. & Sutherland, B.R. 2008 Generation of internal waves by sheared turbulence: experiments. Environ. Fluid Mech. 8 (5-6), 527534.CrossRefGoogle Scholar
Nepf, H.M. 2012 Flow and transport in regions with aquatic vegetation. Annu. Rev. Fluid Mech. 44 (1), 123142.CrossRefGoogle Scholar
Osborn, T.R. & Cox, C.S. 1972 Oceanic fine structure. Geophys. Fluid Dyn. 3 (4), 321345.CrossRefGoogle Scholar
Ouwersloot, H., Moene, A., Attema, J. & De Arellano, J.V.-G. 2017 Large-eddy simulation comparison of neutral flow over a canopy: sensitivities to physical and numerical conditions, and similarity to other representations. Boundary-Layer Meteorol. 162 (1), 7189.CrossRefGoogle Scholar
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, 499534.CrossRefGoogle Scholar
Pham, H. & Sarkar, S. 2017 Turbulent entrainment in a strongly stratified barrier layer. J. Geophys. Res. Oceans 122 (6), 50755087.CrossRefGoogle Scholar
Pham, H.T. & Sarkar, S. 2010 Internal waves and turbulence in a stable stratified jet. J. Fluid Mech. 648, 297324.CrossRefGoogle Scholar
Pham, H.T. & Sarkar, S. 2011 Mixing events in a stratified jet subject to surface wind and buoyancy forcing. J. Fluid Mech. 685, 5482.CrossRefGoogle Scholar
Pham, H.T., Sarkar, S. & Brucker, K.A. 2009 Dynamics of a stratified shear layer above a region of uniform stratification. J. Fluid Mech. 630, 191223.CrossRefGoogle Scholar
Pham, H.T., Sarkar, S. & Winters, K.B. 2012 Near-n oscillations and deep-cycle turbulence in an upper-equatorial undercurrent model. J. Phys. Oceanogr. 42 (12), 21692184.CrossRefGoogle Scholar
Plew, D.R. 2011 Depth-averaged drag coefficient for modeling flow through suspended canopies. J. Hydraul. Engng 137 (2), 234247.CrossRefGoogle Scholar
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 (6), 519539.CrossRefGoogle Scholar
Poggi, D., Porporato, A., Ridolfi, L., Albertson, J. & Katul, G. 2004 The effect of vegetation density on canopy sub-layer turbulence. Boundary-Layer Meteorol. 111 (3), 565587.CrossRefGoogle Scholar
Pollard, R.T., Rhines, P.B. & Thompson, R.O. 1973 The deepening of the wind-mixed layer. Geophys. Fluid Dyn. 4 (4), 381404.CrossRefGoogle Scholar
Polovina, J.J., Mitchum, G.T. & Evans, G.T. 1995 Decadal and basin-scale variation in mixed layer depth and the impact on biological production in the central and north pacific, 1960-88. Deep-Sea Res. I: Oceanogr. Res. Pap. 42 (10), 17011716.CrossRefGoogle Scholar
Polton, J.A. & Belcher, S.E. 2007 Langmuir turbulence and deeply penetrating jets in an unstratified mixed layer. J. Geophys. Res. Oceans 112 (C9), C09020.CrossRefGoogle Scholar
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.CrossRefGoogle Scholar
Raupach, M. & Thom, A.S. 1981 Turbulence in and above plant canopies. Annu. Rev. Fluid Mech. 13 (1), 97129.CrossRefGoogle Scholar
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. Oceans 112 (C2), C02016.CrossRefGoogle Scholar
Shaw, R.H. & Schumann, U. 1992 Large-eddy simulation of turbulent flow above and within a forest. Boundary-Layer Meteorol. 61 (1-2), 4764.CrossRefGoogle Scholar
Shay, T.J. & Gregg, M. 1986 Convectively driven turbulent mixing in the upper ocean. J. Phys. Oceanogr. 16 (11), 17771798.2.0.CO;2>CrossRefGoogle Scholar
Smith, S.D. 1988 Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature. J. Geophys. Res. Oceans 93 (C12), 1546715472.CrossRefGoogle Scholar
Smyth, W. & Peltier, W. 1991 Instability and transition in finite-amplitude Kelvin–Helmholtz and Holmboe waves. J. Fluid Mech. 228, 387415.Google Scholar
Sutherland, B.R. 2005 Stratified shear flow: instability and wave radiation. In Instability of Flows, Adv. in Fluid Mech. Ser. (ed. Rahman, M.), vol. 41, pp. 79104. WIT Press.Google Scholar
Sutherland, B.R. & Linden, P.F. 1998 Internal wave excitation from stratified flow over a thin barrier. J. Fluid Mech. 377, 223252.CrossRefGoogle Scholar
Taylor, J.R. & Sarkar, S. 2007 Internal gravity waves generated by a turbulent bottom ekman layer. J. Fluid Mech. 590, 331354.CrossRefGoogle Scholar
Thom, A. 1971 Momentum absorption by vegetation. Q. J. R. Meteorol. Soc. 97 (414), 414428.CrossRefGoogle Scholar
Utter, B.D. & Denny, M.W. 1996 Wave-induced forces on the giant kelp macrocystis pyrifera (agardh): field test of a computational model. J. Expl Biol. 199 (12), 26452654.CrossRefGoogle ScholarPubMed
VanDine, A., Pham, H.T. & Sarkar, S. 2021 Turbulent shear layers in a uniformly stratified background: dns at high reynolds number. J. Fluid Mech. 916, A42.CrossRefGoogle Scholar
Voss, M., Bange, H.W., Dippner, J.W., Middelburg, J.J., Montoya, J.P. & Ward, B. 2013 The marine nitrogen cycle: recent discoveries, uncertainties and the potential relevance of climate change. Philos. Trans. R. Soc. B: Biol. Sci. 368(1621), 20130121.CrossRefGoogle ScholarPubMed
Yan, C., McWilliams, J.C. & Chamecki, M. 2021 Generation of attached Langmuir circulations by a suspended macroalgal farm. J. Fluid Mech. 915, A76.CrossRefGoogle Scholar
Yan, C., McWilliams, J.C. & Chamecki, M. 2022 Overlapping boundary layers in coastal oceans. J. Phys. Oceanogr. 52 (4), 627646.CrossRefGoogle Scholar
Zagvozkin, T., Vorobev, A. & Lyubimova, T. 2019 Kelvin–Helmholtz and Holmboe instabilities of a diffusive interface between miscible phases. Phys. Rev. E 100 (2), 023103.CrossRefGoogle ScholarPubMed
Figure 0

Figure 1. Simulation set-up. $(a)$ Initial temperature profile. The vertical coordinate $z$ is normalized by the initial OML depth $h_0$. Dotted black lines indicate the canopy base (upper line) and the sponge layer (lower line). $(b)$ Initial nutrient profile. $(c)$ A schematic of the canopy simulation (side view). The frond surface area density $a$ is constant and uniform throughout the extent of the canopy. Values of $a$ for different simulations are detailed in table 1.

Figure 1

Table 1. Simulation parameters. Here ‘NC’ represents cases with no canopy, and ‘C’ represents cases with a canopy, with 0.1, 0.3, 1 and 3 denoting the canopy frond area density $a$. The suffix ‘-NSW’ indicates simulations conducted without surface gravity waves, where the Stokes drift velocity is zero. The suffix ‘-WG’ represents the case with a weaker geostrophic current $u_g=0.1\,\rm m\,s^-{^1}$, and ‘-NG’ represents the case with no geostrophic current $u_g=0\,\rm m\,s^-{^1}$. Additionally, a simulation is conducted with a doubled domain depth, named ‘-DEEP’.

Figure 2

Figure 2. The mixed layer depth $h_{OML}$, normalized by its initial value $h_0$, as a function of time. $(a)$ Simulations with varying canopy densities (details in the legend). Generally, higher canopy density (cases C0.1–C3) leads to more pronounced OML deepening compared with cases without a canopy (cases NC and NC-NSW). $(b)$ Simulations with different geostrophic currents $u_g$. Weaker geostrophic currents lead to less OML deepening.

Figure 3

Figure 3. Side views of the cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$. Snapshots at around $tf=1$ are shown for four representative cases: $(a)$ NC, no canopy, with the Stokes drift; $(b)$ NC-NSW, no canopy, no Stokes drift; $(c)$ C0.3, canopy density of 0.3 m$^{-1}$; $(d)$ C3, canopy density of 3 m$^{-1}$; $(e)$ C3-WG, canopy density of 3 m$^{-1}$ and a weaker geostrophic current; and $(f)$ C3-NG, canopy density of 3 m$^{-1}$ and no geostrophic current. Black dotted lines represent the canopy base depth (lines also included in the no-canopy cases for comparison purposes). Note that the sponge layer is not shown in the plots.

Figure 4

Figure 4. Velocity and temperature profiles in four representative simulations at the final simulation time $tf=9$. $(a,b)$ Horizontally averaged streamwise velocity $\langle \overline {u}\rangle _{xy}$ and cross-stream velocity $\langle \overline {v}\rangle _{xy}$, both normalized by the geostrophic velocity $u_g=0.2\,\rm m\,s^-{^1}$. $(c,d)$ Horizontally averaged temperature $\langle \overline {\theta }\rangle _{xy}$ and the vertical gradient of temperature. Horizontal dotted lines represent the canopy base depth.

Figure 5

Figure 5. $(a)$ The dominant $x$-direction wavenumber $k_{x,peak}$ versus canopy shear layer thickness $\delta _s$. Different colours represent simulations with different canopy density and $u_g$ (C0.1–C3 and C3-WG, see the legend), and each point corresponds to a time interval of every $tf=0.0015$. $(b)$ The dominant phase speed $c_{x,peak}$ of the KH-type structures versus the current speed of the inflection point in the shear layer. Note that the results are normalized by $u_g=0.1\,\rm m\,s^-{^1}$ in case C3-WG and by $u_g=0.2\,\rm m\,s^-{^1}$ in all the other cases.

Figure 6

Figure 6. Time evolution of the velocity and temperature field in case C3. $(a)$ Vertical profiles of horizontally averaged streamwise velocity $\langle \overline {u}\rangle _{xy}$, at $tf=0.3$, 0.5 and 9. $(b)$ Vertical profiles of horizontally averaged temperature $\langle \overline {\theta }\rangle _{xy}$, at every $tf=1$. $(c,d)$ Side views of the cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$ at $tf=0.3$ and $tf=3$, representing stages earlier and later than figure 3$(d)$, respectively.

Figure 7

Figure 7. $(a)$ The TKE profiles in four representative cases, normalized by the square of the friction velocity $u_*$. $(b)$ The vertical component of TKE. $(c)$ Ratio of the vertical component to the total TKE. $(d)$ Skewness of the vertical turbulent velocity. The results correspond to around $tf=1$, same as the time when the side views in figure 3 are shown. Note that the ratio and skewness are not shown for regions where TKE$/u_*^2\lt 1$. Horizontal dotted lines represent the canopy base depth.

Figure 8

Figure 8. Terms in the TKE budget (3.9) for case NC (no canopy, with the Stokes drift) and case NC-NSW (no canopy, no Stokes drift) at around $tf=1$: $(a)$ shear production associated with the streamwise velocity; $(b)$ Stokes production; $(c)$ dissipation; $(d)$ buoyancy production; $(e)$ pressure work (vertical divergence of pressure flux); $(e)$ pressure flux. Only the dominant terms are shown. Note the differences in the horizontal axis range.

Figure 9

Figure 9. Terms in the TKE budget (3.9) for case C3 (with the highest canopy density): red lines, at $tf=1$; blue lines, at $tf=0.3$; black lines, at $tf=9$; $(a,b)$ shear production associated with the streamwise velocity and cross-stream velocity, respectively; $(c)$ Stokes production; $(d)$ dissipation; $(e)$ canopy drag dissipation; $(f)$ TKE transport; $(g)$ buoyancy production; $(h)$ pressure work (vertical divergence of pressure flux); $(i)$ pressure flux. Note the differences in the horizontal axis range. Horizontal dotted lines represent the canopy base depth.

Figure 10

Figure 10. Side-view and cross-sectional plots of vertical velocity gradients at around $tf=1$ in case NC. $(a)$ The cross-stream averaged vertical velocity gradient $\partial \langle w\rangle _{y}/\partial z$ in the $x-z$ plane. $(b)$ The vertical velocity gradient $\partial w/\partial z$ at $y=0$ in the $x-z$ plane. $(c)$ The $x$-direction averaged vertical velocity gradient $\partial \langle w\rangle _{x}/\partial z$ in the $y-z$ plane. $(d)$ The vertical velocity gradient $\partial w/\partial z$ at $x=0$ in the $y-z$ plane.

Figure 11

Figure 11. Spectra of vertical velocity gradients $\partial w/\partial z$ in case NC, versus $x$- or $y$-wavenumber or frequency. $(a{-}c)$ Spectra of the $x$- or $y$-direction averaged results. $(e,f)$ The $x$- or $y$-direction average of the spectra at every $x$- or $y$-location. Dashed lines represent the average between $z/h_0=-1$ to $0$ (within the OML), and solid lines represent the average between $z/h_0=-3$ and $-2$ (below the thermocline). Vertical red lines indicate the buoyancy frequency $N$. Note that the frequency in the spectral plot is generally much higher than $N$ due to the Doppler shift caused by the geostrophic current. In the text, intrinsic internal wave frequencies are analysed after removing the Doppler shift effect.

Figure 12

Figure 12. Side views of vertical velocity gradients at around $tf=1$ in case C3. $(a)$ The cross-stream average vertical velocity gradient $\partial \langle w\rangle _{y}/\partial z$. $(b)$ The vertical velocity gradient $\partial w/\partial z$ at $y=0$.

Figure 13

Figure 13. Spectra of vertical velocity gradients $\partial w/\partial z$ in case C3, versus $x$- or $y$-wavenumber or frequency (a similar set of panels to figure 11). $(a{-}c)$ Spectra of the $x$- or $y$-direction averaged results. $(e,f)$ The $x$- or $y$-direction average of the spectra at every $x$- or $y$-location. Dashed lines represent the average between $z/h_0=-1$ to $0$ (within the OML), and solid lines represent the average between $z/h_0=-3$ and $-2$ (below the thermocline).

Figure 14

Figure 14. Vertical fluxes of momentum and temperature for case NC (black) and case C3 (red), time-averaged from $tf=0$ to $tf=9$ and horizontally averaged. $(a,b)$ Turbulent fluxes of momentum and temperature. $(c,d)$ Fluxes driven by lower-frequency motions. Fluxes are normalized by $u_*$ and $\Delta \theta =1\,\rm K$.

Figure 15

Figure 15. $(a)$ Nutrient profiles for case NC (black) and case C3 (red), at $tf=9$. The thin grey line shows the initial nutrient profile at $tf=0$. $(b,c)$ Vertical nutrient fluxes driven by turbulence and lower-frequency motions, respectively. These fluxes are time-averaged from $tf=0$ to $tf=9$ and horizontally averaged. Fluxes are normalized by $u_*$ and $\Delta C_N=10000$$\unicode{x03BC}\rm mol\,m^-{^3}$.

Figure 16

Figure 16. Side views of case C3-DEEP at around $tf=1$. $(a)$ The cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$. $(b)$ The cross-stream averaged vertical velocity gradients $\partial \langle w\rangle _{y}/\partial z$. Black dotted lines represent the canopy base depth. Note that the sponge layer is not shown.

Figure 17

Figure 17. Spectra of the $x$- or $y$-direction averaged vertical velocity gradients in case C3-DEEP, versus $(a)$$x$-wavenumber, $(b)$$y$-wavenumber or $(c)$ frequency. Solid lines represent the average between $z/h_0=-3$ to $-2$, and dash–dotted lines represent the average between $z/h_0=-7$ to $-6$.

Figure 18

Figure 18. A side-view plot of the cross-stream averaged temperature gradients $\partial \langle \theta \rangle _{y}/\partial z$ in the finite-length canopy simulation. The dotted rectangle shows the extent of the canopy.

Supplementary material: File

Bo et al. supplementary material movie

Time evolution of canopy shear layer dynamics in case C3. Upper panel: A side view of the cross-stream averaged temperature gradients. Lower panels: Vertical profiles of the horizontally averaged squared vertical shear, squared buoyancy frequency, gradient Richardson number, and canopy drag force.
Download Bo et al. supplementary material movie(File)
File 16.4 MB