2. Basic model and wave–current decomposition

The theory in this paper analyses wave–current–drag interactions while making the unrealistic simplifications of constant water density, a flat bottom and an absence of non-conservative forces apart from the canopy drag; such generalizations can be pursued later.

Thus, the incompressible fluid equations are

(2.1)\begin{equation} \left. \begin{aligned} \partial_t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} + f \hat{\boldsymbol{z}} \times \boldsymbol{u} &={-} \boldsymbol{\nabla} \phi + \boldsymbol{D}^h , \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} &= 0 , \end{aligned} \right\} \end{equation}

where $\boldsymbol {u} = (\boldsymbol {u}^h,w) = (u,v,w)$ is the velocity, $\phi$ is the density-normalized pressure as a deviation from the background hydrostatic profile $- g z$ that does not accelerate the flow ($g$ is the gravitational constant), $f$ is the Coriolis frequency (assumed constant) and $\boldsymbol {D}^h$ is the drag force (1.1).

For boundary conditions, all fields are periodic in the horizontal coordinates $(x,y)$, and in the vertical,

(2.2)\begin{equation} \left. \begin{aligned} w &= \partial_t \eta + \boldsymbol{u}^h \boldsymbol{\cdot} \boldsymbol{\nabla}^h \eta ,\\ \phi &= g\eta \end{aligned} \right\} \end{equation}

at the sea surface, $z = \eta (x,y,t)$, and

(2.3)\begin{equation} w = 0 \end{equation}

at the bottom, $z = - H$ (with $H$ a constant). The second condition in (2.2) is a statement that the total pressure is uniform at the surface interface. (Boundary stresses are omitted for now.)

Direct numerical solutions of (2.1)–(2.3) are difficult for general waves and currents especially because of the disparate time scales of wave propagation and current evolution and the tendency for $\eta$ to develop wave-related singularities (i.e. to break). So, to obtain a wave-averaged LES-type or circulation model for the currents, several approximations are made. These are listed in (1.4a,b) and (1.5). In particular, they imply that the wave solutions at leading order are linear and uninfluenced by currents and drag (§ 4).

The relations in this paper are all dimensional, but scaling estimates are made (§ 3) based on non-dimensionalizations using the primary wave scales: $a$, $\boldsymbol {x} \sim 1/k$ (with the current length scale $\ell$ assumed to be the same or larger), $t \sim 1/\sigma$, $\boldsymbol {u} \sim C = \sigma /k$, etc.

The velocity and other fields are decomposed into currents and waves, e.g.

(2.4)\begin{equation} \boldsymbol{u} = \bar{\boldsymbol{u}} + \boldsymbol{u}' , \end{equation}

with the distinction that $\boldsymbol {u}'$ has zero average over a wave period, $2{\rm \pi} /\sigma$. A further decomposition is made for the wave velocity,

(2.5)\begin{equation} \boldsymbol{u}' = \boldsymbol{u}'_w + \boldsymbol{u}'_c + \boldsymbol{u}'_d . \end{equation}

Here $\boldsymbol {u}'_w$ is the primary wave field with linear dynamics, its magnitude being ${O}(\epsilon C)$; $\boldsymbol {u}'_c$ is the small modification of the waves due to the current effects on the waves (CEW; § 6), its magnitude being ${O}(\epsilon \mu C)$; and $\boldsymbol {u}'_d$ is the small modification of the waves due to the transient drag force $\boldsymbol {D}^{h\prime }$ (§ 5), its magnitude being O $(\epsilon \nu C)$. The respective wave corrections due to drag and currents in $\boldsymbol {u}_d'$ and $\boldsymbol {u}_c'$ can be superimposed because their determining equations are linear (§§ 5 and 6). No decomposition is imposed on $\bar {\boldsymbol {u}}$ ($\sim \mu C$); the goal is to derive a wave-averaged dynamical balance equation for $\bar {\boldsymbol {u}}$ including all of these wave effects on currents (WEC; § 10).

3. Scaling assumptions and rationale

The primary scaling assumptions of $\epsilon,\mu,\nu \ll 1$ are stated in § 1, and the wave and current velocities have scaling estimates of

(3.1a,b)\begin{equation} \boldsymbol{u}'\sim \epsilon C , \quad \bar{\boldsymbol{u}} \sim \mu C , \end{equation}

respectively. For strict validity of a quasi-linear perturbation theory, $\epsilon$ must be small compared to $\mu$ or $\nu$, which often is not true in nature. Wind–wave equilibrium in the surface layer with Ekman currents scales as $\mu \sim \epsilon ^2$; i.e. $\mu \ll \epsilon$ with wave orbital motions ($\sim \epsilon C$) stronger than currents ($\sim \mu C$). Hanley, Belcher & Sullivan (Reference Hanley, Belcher and Sullivan2010) show the empirical distribution of the so-called turbulent Langmuir number, $La = \sqrt {u_*/|\boldsymbol {u}^{St}|(0)} \sim \mu ^{1/2}/\epsilon$, which is peaked near the equilibrium value of 0.3. Stronger currents would increase the $\mu /\epsilon$ ratio and big swells decrease it. With the Langmuir turbulence context in mind, the primary target regime is where the wave orbital motion is faster than the current speed, which itself is often comparable to the Stokes drift:

(3.2)\begin{equation} \boldsymbol{u}'_w \gg \bar{\boldsymbol{u}} \sim \boldsymbol{u}^{St} . \end{equation}

This ordering allows simpler formulae for the canopy drag on the waves and for the wave-induced canopy drag on the currents (§§ 5 and 9). Furthermore, for generality of the results, we assume that drag effects are comparable to current effects, both on the wave perturbations (§§ 5–7) and on the currents (§ 10); i.e. $\nu \sim \mu \ll 1$.

The theoretical focus in wave–current interaction is usually on the interactions associated with the spectrum-peak waves (i.e. at most a few wavenumber components that have the largest amplitudes $a$ in $\eta$). For these waves, the nonlinear evolution effects are often weak, even if the full wave spectrum has stronger nonlinear interactions. The success of quasi-linear theory and even the simplification to monochromatic primary waves as a modelling framework for Langmuir turbulence (McWilliams, Sullivan & Moeng Reference McWilliams, Sullivan and Moeng1997) provides motivation to continue this approach with the inclusion of canopy drag.

Once a commitment is made to a quasi-linear theory for wave–current interaction, in part because it is a doable theory for an otherwise very difficult problem, then the magnitude of $\epsilon$ need not be declared. However, with the nonlinear canopy drag law (1.1), the ratio $\mu /\epsilon$ is important for how the drag force is approximated and partitioned between waves and currents. In this paper the choice is made to assume that wave orbital velocities are stronger than currents – the most common natural situation – as the main line of argument, with tangential remarks about the alternatives (§ 9).

4. Primary wave

For simplicity, the primary wave field is assumed to be monochromatic, i.e. a single eigenmode for linear wave dynamics (but also see § 11 and Appendix B). We adopt a complex notation; e.g. for the surface elevation,

(4.1)\begin{equation} \eta'_w = a \frac{\mathcal{E}}{2} + {\rm c.c.} \end{equation}

Here $a$ is the amplitude, $\mathcal {E} = \exp [\textrm {i}\varTheta ]$ is the phase function, $\varTheta = \boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}^h - \sigma t$ is its horizontal-propagation argument, $\boldsymbol {k}$ is the horizontal wavenumber vector (with magnitude $k$), $\sigma$ is the frequency eigenvalue and $\textrm {c.c.}$ denotes the complex conjugate of the preceding expression. In finite-depth water ($kH < \infty$), the dispersion relation is

(4.2)\begin{equation} \sigma = \sqrt{gk\tanh[kH]} > 0 . \end{equation}

The propagation is in the direction of $\boldsymbol {k}$. The other eigenfunction fields are

(4.3)\begin{equation} \left. \begin{aligned} \phi'_w & = g a F(z) \frac{\mathcal{E}}{2} + {\rm c.c.},\\ \boldsymbol{u}^{h\prime}_w & = \frac{g a}{\sigma} F(z) \boldsymbol{k} \frac{\mathcal{E}}{2} + {\rm c.c.},\\ w'_w & ={-} \frac{{\rm i} a g k}{\sigma} G(z) \frac{\mathcal{E}}{2} + {\rm c.c.} \end{aligned} \right\} \end{equation}

The non-dimensional vertical structure functions are

(4.4)\begin{equation} \left. \begin{aligned} F(z) & = \cosh[kz] + \tanh[kH]\sinh[kz], \\ G(z) & = \sinh[kz] + \tanh[kH]\cosh[kz] . \end{aligned} \right\} \end{equation}

Functions $F$ and $G$ decay monotonically with depth. At $z = -H$, $G = 0$ while $F\ne 0$. If $a$ is real, then $\eta '_w$, $\phi '_w$ and $\boldsymbol {u}'_w$ are all $\propto \cos [\varTheta ]$, and $w'_w \propto \sin [\varTheta ]$. Furthermore, there is no vorticity in this eigenmode, $\boldsymbol {\zeta }'_w = \boldsymbol {\nabla } \times \boldsymbol {u}'_w = 0$, and this is quite helpful in solving for the drag and current wave corrections in §§ 5 and 6.

In a multiscale theory, as here, $a$ can be a function of ‘slow’ time $\tau$ due to drag and current effects (§§ 5–8), in which case $a(\tau )$ can become complex due to a slow evolution in the primary wave phase. In this paper, for simplicity, $a$ is assumed to be real and positive in discussing specific canopy-related results (§§ 5 and 9–11 and Appendices A and B), while its more general complex form is retained for the current interactions (§§ 6 and 7). The latter allows for ‘wave–current resonance’ that might also incorporate canopy effects, but that is not a central focus in this paper; also see the extended discussion of wave–current resonance and several examples in WCI22, chapters 7 and 8.

An important property of surface waves is their Lagrangian mean flow, the Stokes drift velocity. For small $\epsilon$, this is defined by

(4.5)\begin{equation} \boldsymbol{u}^{St} = \left\langle \left( \int^t \boldsymbol{u}' \,{\rm d}t \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u}' \right\rangle , \end{equation}

where the angle brackets denote an average over the wave period (i.e. this has the same meaning as the overbar symbol). In this paper only the leading-order approximation to $\boldsymbol {u}^{St}$ is needed. With $\boldsymbol {u}' \approx \boldsymbol {u}'_w$,

(4.6)\begin{equation} \left. \begin{aligned} \boldsymbol{u}^{hSt}(z) & \approx \frac{|a|^2 \sigma \boldsymbol{k}}{2\sinh^2[kH]} \cosh[2k(z+H)] ,\\ w^{St}(z) & \approx 0 . \end{aligned} \right\} \end{equation}

The drift is only in the horizontal direction. Also, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}^{St} = 0$, as is usual in conservative quasi-linear wave theories (WCI22, appendix C).

For completeness, these primary wave formulae are also listed in the deep-wave limit, $kH \rightarrow \infty$ (even though that is inconsistent with bottom-rooted marine vegetation affected by surface wave motions):

(4.7a–c)\begin{equation} \sigma \rightarrow \sqrt{gk}, \quad F(z),G(z) \rightarrow {\rm e}^{kz}, \quad \frac{\cosh[2k(z+H)]}{2\sinh^2[kH]} \rightarrow {\rm e}^{2kz} . \end{equation}

Everything in this section is well known.

5. Effect of drag on the waves

With the decomposition (2.5) and the assumption that wave orbital motions are stronger than currents (3.2), the largest canopy drag term is

(5.1)\begin{equation} \boldsymbol{D}^h_w ={-} \mathcal{D} | \boldsymbol{u}^{h \prime}_w | \boldsymbol{u}^{h \prime}_w . \end{equation}

Substituting from (4.3) with real $a > 0$,

(5.2)\begin{equation} \boldsymbol{D}^{h}_w ={-} \mathcal{D} \left(\frac{(ga)^2F^2k}{\sigma^2}\right) |{\cos[\varTheta]}| \cos[\varTheta] \boldsymbol{k} . \end{equation}

The functional dependence on $\varTheta$ can be approximated as

(5.3a,b)\begin{equation} \left\langle |{\cos[\varTheta]}| \cos[\varTheta] \right\rangle = 0 , \quad |{\cos[\varTheta]}| \cos[\varTheta] \approx c_1 \cos[\varTheta] , \end{equation}

as demonstrated in figure 1. The coefficient $c_1 = 8/3{\rm \pi} \approx 0.849$ is the projection of the left-hand side of (5.2) on $\cos [\varTheta ]$, i.e. it is the value that gives the least-square error in this approximation averaged over a primary wave period; the root-mean-square (r.m.s.) relative error in the approximation is <20 %. Thus, the wave forcing by canopy drag is purely oscillatory, with $\bar {\boldsymbol {D}}^h_w = 0$, and it is approximately monochromatic. The latter allows more extensive analytic results compared to the exact form of $\boldsymbol {D}^h_w$; see Appendices A and B for generalizations.

Figure 1. Plots of $|{\cos [\varTheta ']}| \cos [\varTheta ']$ (black) and $c_1 \cos [\varTheta ']$ (red) over one wave period in $\varTheta ' = \sigma t / 2{\rm \pi}$.

The wave drag has the same phase function as the eigenmode, viz.

(5.4)\begin{equation} \boldsymbol{D}^{h\prime}_w \approx{-} \mathcal{D}c_1 \left(\frac{(ga)^2F^2k}{\sigma^2}\right) \cos[\varTheta] \boldsymbol{k} ; \end{equation}

and $\boldsymbol {u}'_d$ is needed for the WEC model (§ 10). With this approximation for $\boldsymbol {D}^{h\prime }_w$ as the forcing term for the O $(\nu )$ correction to $\boldsymbol {u}'_w$, this makes the solution for $\boldsymbol {u}'_d$ itself monochromatic. Obviously, a more complete solution with $\boldsymbol {D}^{h\prime }$ could be obtained with higher harmonics of $\varTheta$, and Appendix A shows that the errors incurred by this simplification are minor.

The solution procedure for $\boldsymbol {u}'_d$ (as for $\boldsymbol {u}'_c$ in § 6 and in WCI22, chapter 5) is decomposed into vortical and irrotational parts:

(5.5)\begin{equation} \boldsymbol{u}'_d = \boldsymbol{u}^{\zeta\prime}_d + \boldsymbol{u}^{\varPhi\prime}_d , \end{equation}

where $\boldsymbol {\nabla } \times \boldsymbol {u}^{\varPhi \prime }_d = 0$. As a consequence of $\boldsymbol {\zeta }'_w = 0$, the vortical component is relatively easy to obtain, while the irrotational component may be more difficult.

The wave vorticity equation with respect to drag at O $(\nu )$ is

(5.6)\begin{equation} \partial_t \boldsymbol{\zeta}'_d = \mathrm{curl}[\partial_t \boldsymbol{u}^{\zeta\prime}_d] = \mathrm{curl}[\boldsymbol{D}^{h\prime}_w] . \end{equation}

There is no vertical component in the drag force (1.1), and

(5.7)\begin{equation} \mathrm{curl}[\boldsymbol{D}^{h\prime}_w] = c_1 \frac{g^2a^2k}{\sigma^2} {\left(\begin{array}{c} k^y \partial_z[\mathcal{D}F^2] \\ - k^x \partial_z[\mathcal{D}F^2] \\ F^2 (\hat{\boldsymbol{z}} \boldsymbol{\cdot} \boldsymbol{k} \times \boldsymbol{\nabla}^h \mathcal{D}) \end{array}\right)} \cos[\varTheta] . \end{equation}

To invert (5.6), the time integral of the phase function is

(5.8)\begin{equation} \int^t \cos[\varTheta] \,{\rm d}t ={-}\frac{1}{\sigma}\sin[\varTheta] . \end{equation}

Therefore, using (5.4)–(5.7),

(5.9)\begin{equation} \boldsymbol{\zeta}'_d = c_1 \frac{g^2a^2k}{\sigma^3} {\left(\begin{array}{c} - k^y \partial_z[\mathcal{D}F^2] \\ + k^x \partial_z[\mathcal{D}F^2] \\ - F^2 (\hat{\boldsymbol{z}} \boldsymbol{\cdot} \boldsymbol{k} \times \boldsymbol{\nabla}^h \mathcal{D}) \end{array}\right)} \sin[\varTheta] . \end{equation}

Furthermore, the vorticity relation for $\zeta '_d = \boldsymbol {\nabla } \times \boldsymbol {u}'_d$ in (5.9) can be inverted for the vortical component of the wave velocity to obtain

(5.10)\begin{equation} \boldsymbol{u}^{h\zeta\prime}_d = \mathcal{D} c_1 \left(\frac{(ga)^2F^2k}{\sigma^3}\right) \sin[\varTheta] \boldsymbol{k} , \quad w^{\zeta\prime}_d = 0 . \end{equation}

Notice that the vertical component of $\boldsymbol {\zeta }'_d$ is zero if $\boldsymbol {\nabla }^h \mathcal {D} =0$, i.e. if the vegetation is horizontally uniform.

The irrotational component $\boldsymbol {u}_d^{\varPhi \prime }$ in (5.5) requires the solution of an elliptic boundary-value problem. It is derived by taking the divergence of the momentum equation at O $(\nu )$ and utilizing a combination of the vertical momentum equation and the conditions (2.2) and (2.3) at the boundaries to eliminate the velocity in favour of a velocity potential function $\varPhi '_d$ defined by

(5.11)\begin{equation} \boldsymbol{u}_d^{\varPhi\prime} ={-} \boldsymbol{\nabla}\!\int^t \varPhi'_d \,{\rm d}t \end{equation}

(WCI22, chapter 5). The derivation includes recognizing that $\boldsymbol {\nabla } \boldsymbol {\cdot } \partial _t \boldsymbol {u}'_d = 0$ and $\boldsymbol {\zeta }'_w = 0$. It also includes the linearized wave forms of (2.2) and (2.3) in a Taylor series expansion about $z = 0$:

(5.12)\begin{equation} \left. \begin{aligned} w' & = \partial_t \eta' + \boldsymbol{\nabla}^h \boldsymbol{\cdot} (\eta'\bar{\boldsymbol{u}}^{h}) \quad \mathrm{and} \quad \phi' = g \eta' \quad \mathrm{at} \ z = 0 ,\\ w' & = 0\quad \mathrm{at} \ z ={-} H . \end{aligned} \right\} \end{equation}

For this drag-effect derivation, the $\bar {u}^h$ boundary-condition term is not considered, but it will be included for the CEW component $\boldsymbol {u}_c^{\varPhi \prime }$ in § 6. Also, the wave vertical momentum equation has no drag component (i.e. $\hat {\boldsymbol {z}} \boldsymbol {\cdot } \boldsymbol {D}^{h\prime }_w = 0$). The combination of the relations in this paragraph is the following partial differential equation (PDE) system:

(5.13)\begin{equation} \left. \begin{aligned} \nabla^2 \varPhi'_d & = \mathrm{div}[\boldsymbol{D}^{h\prime}_w] , \\ \left[\partial_z - \frac{\sigma^2}{g} \right] \varPhi'_d & ={-} \frac{2}{g} \partial_t\partial_{\tau} \phi'_w \quad \mathrm{at} \ z = 0 ,\\ \partial_z \varPhi'_d & = 0 \quad \mathrm{at} \ z ={-} H . \end{aligned} \right\} \end{equation}

The ‘forcing’ term in (5.13) includes the possibility of a ‘slow-time’ evolution of the primary wave field (4.1)–(4.3) through its amplitude, $a(\tau )$, i.e.

(5.14)\begin{equation} - \frac{2}{g} \partial_t\partial_{\tau} \phi'_w ={-} 2 \sigma \sin[\varTheta] \partial_\tau a \end{equation}

for a slow-time coordinate, $\tau = \nu t$. This slow-time dependence is determined through a solvability condition for (5.13) that assures there is no resonant forcing of $\varPhi '_d$, i.e. that there is no projection of the forcing on the primary wave eigenmode (WCI22, appendix E).

It will be shown that canopy drag induces a slow amplitude decay for the primary wave, and this condition must be combined with an analogous constraint on $\varPhi '_c$, hence providing an amplitude dependence also on $\bar {\boldsymbol {u}}$; both of these statements will be further explained in §§ 6–8. The system (5.13) must be solved in conjunction with the evolution of the current $\bar {\boldsymbol {u}}$ also on the slow time scale $\tau$ (§ 10).

With (5.4) the drag forcing in the PDE in (5.13) is

(5.15)\begin{align} \mathrm{div}[\boldsymbol{D}^{h\prime}_w] &={-} c_1 \left(\frac{(gaF)^2k}{\sigma^2}\right) (\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{\nabla}^h)[\mathcal{D}\cos[\varTheta]] \nonumber\\ &= c_1 \left(\frac{(gaF)^2k}{\sigma^2}\right) \{k^2 \mathcal{D}\sin[\varTheta] - (\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{\nabla}^h \mathcal{D})\cos[\varTheta] \} . \end{align}

In the particular case of horizontally uniform vegetation, $\boldsymbol {\nabla }^h \mathcal {D} = 0$, and we can factor $\varPhi '_d$ for its $(\boldsymbol {x}^h,z,t)$ dependences,

(5.16)\begin{equation} \varPhi'_d = \varphi_d(z) \sin[\varTheta] , \end{equation}

where $\varphi _d$ satisfies a one-dimensional (1-D) boundary-value problem from (5.13):

(5.17)\begin{equation} \left. \begin{aligned} {}[\partial^2_z - k^2] \varphi_d & = \mathcal{D}(z) c_1 \frac{g^2a^2k^3}{\sigma^2} F^2(z), \\ \left[\partial_z - \frac{\sigma^2}{g} \right] \varphi_d & ={-} 2 \sigma \partial_\tau a \quad \mathrm{at} \ z = 0, \\ \partial_z \varphi_d & =0 \quad \mathrm{at} \ z ={-} H . \end{aligned} \right\} \end{equation}

Furthermore, this 1-D problem only needs to be solved once at each slow time $\tau$ as $a(\tau )$ evolves. Or, even more efficiently, it can be partitioned into two problems that only need to be solved once – one with only ordinary differential equation forcing by the interior factor of $a^2$, the other with only surface boundary forcing by the factor of $\partial _\tau a$ – then the two answers are superposed with the current values of $a(\tau )$ and $\partial _\tau a(\tau )$ as their factors.

6. Current effects on the waves (CEW)

In WCI22 (chapter 5) the quasi-linear theory is presented for CEW in $\boldsymbol {u}'_c$, and the relevant parts of it will be repeated briefly. The wave momentum equation is linearized with respect to wave–current advection, i.e. CEW occurs at O $(\mu )$ compared to the primary wave balance with respect to the current effects, and (2.1) is approximated as

(6.1)\begin{equation} \left. \begin{aligned} \partial_t \boldsymbol{u}'_c + f \hat{\boldsymbol{z}} \times \boldsymbol{u}'_c & ={-} \boldsymbol{\nabla} \phi'_c - \boldsymbol{\nabla} (\boldsymbol{u}'_w \boldsymbol{\cdot} \bar{\boldsymbol{u}}) + \boldsymbol{u}'_w \times \bar{\boldsymbol{\zeta}} , \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}'_c &= 0 , \end{aligned} \right\} \end{equation}

where $\bar {\boldsymbol {\zeta }}$ is the current vorticity. Besides the pressure gradient, the right-hand-side wave forces are the gradient of the wave Bernoulli head and the wave vortex force. The associated vertical boundary conditions on $\boldsymbol {u}'_c$, $\phi '_c$ and $\eta '_c$ are (5.12) now with a non-zero $\bar {\boldsymbol {u}}^h$. Because the wave components for both drag and current satisfy linear equations, they can be superposed in $\boldsymbol {u}'$.

In this and the next sections $a$ is complex (see § 4).

As before in (5.5), the wave component $\boldsymbol {u}'_c$ is decomposed into vortical and divergent components,

(6.2)\begin{equation} \boldsymbol{u}'_c = \boldsymbol{u}_c^{\zeta\prime} + \boldsymbol{u}_c^{\varPhi\prime}. \end{equation}

Taking the results from WCI22 (section 5.2), with modifications for finite $kH$, the wave vorticity associated with (6.1) is

(6.3)\begin{align} \boldsymbol{\zeta}'_c &= \boldsymbol{\nabla} \times \left(\int^t \boldsymbol{u}'_w \,{\rm d}t\right) \times (\, f\hat{\boldsymbol{z}} + \bar{\boldsymbol{\zeta}} ) \nonumber\\ &= \boldsymbol{\nabla} \times \frac{ag}{\sigma^2} \{ {\rm i} F (\boldsymbol{k} \times (\,f + \bar{\zeta}^z) \hat{\boldsymbol{z}} ) + G k (\hat{\boldsymbol{z}} \times \bar{\boldsymbol{\zeta}}^h) \} \frac{\mathcal{E}}{2} + {\rm c.c.} \end{align}

The gradient cross-product in this relation can be inverted to yield the vortical velocity component $\boldsymbol {u}_c^{\zeta \prime } = (\boldsymbol {\nabla } \times )^{-1} \boldsymbol {\zeta }'_c$, disregarding a possible contribution to the divergent component. The result is

(6.4)\begin{equation} \boldsymbol{u}_c^{\zeta\prime} = \frac{ga}{\sigma^2} {\left(\begin{array}{c} {\rm i}k^y(\,f+\bar{\zeta}^z)F(z) - k\bar{\zeta}^yG(z) \\ k\bar{\zeta}^xG(z) - {\rm i}k^x(\,f+\bar{\zeta}^z)F(z) \\ {\rm i}(k^x\bar{\zeta}^y - k^y \bar{\zeta}^x)F(z) \end{array}\right)} \frac{\mathcal{E}}{2} + {\rm c.c.} \end{equation}

As with $\boldsymbol {u}_d^{\varPhi \prime }$ in § 5, the irrotational (divergent) velocity component is the gradient of a potential function,

(6.5)\begin{equation} \boldsymbol{u}^{\varPhi\prime}_c ={-} \boldsymbol{\nabla} \!\int^t \varPhi'_c \,{\rm d}t , \end{equation}

where $\varPhi '_c = \phi '_c + \boldsymbol {u}'_w \boldsymbol {\cdot } \bar {\boldsymbol {u}}$ is a combined pressure and Bernoulli-head potential. Following the same derivation path as outlined in § 5, the associated boundary-value problem for $\varPhi '_c$ is the following:

(6.6)\begin{align} \left. \begin{aligned} \nabla^2 \varPhi'_c & ={-} \boldsymbol{\nabla} \left(\int^t \phi'_w \,{\rm d}t\right) \boldsymbol{\cdot} \nabla^2 \bar{\boldsymbol{u}} \\ & = \frac{ag}{\sigma} ( F(z)\boldsymbol{k} \boldsymbol{\cdot} \nabla^2 \bar{\boldsymbol{u}}^h - {\rm i}G(z) k \nabla^2 \bar{w} ) \frac{\mathcal{E}}{2} + {\rm c.c.},\\ \left[\partial_z - \frac{\sigma^2}{g} \right] \varPhi'_c & ={-} \frac{2}{g} \partial_t\partial_{\tau} \phi'_w - \frac{2}{g} \boldsymbol{\nabla}^h (\partial_t\phi'_w) \boldsymbol{\cdot} \bar{\boldsymbol{u}}^h - \boldsymbol{\nabla}^h \left(\int^t \phi'_w \,{\rm d}t\right) \boldsymbol{\cdot} \partial_z \bar{\boldsymbol{u}}^h\\ & \quad - \frac{1}{g} (\partial_t\phi'_w) (\boldsymbol{\nabla}^h \boldsymbol{\cdot} \bar{\boldsymbol{u}}^h ) \quad \mathrm{at} \ z = 0\\ & = \left( 2 {\rm i} \sigma \partial_\tau a + a\sigma\left({-}2 \boldsymbol{k}\boldsymbol{\cdot}\bar{\boldsymbol{u}}^h + \frac{g}{\sigma^2}\boldsymbol{k}\boldsymbol{\cdot}\partial_z\bar{\boldsymbol{u}}^h + {\rm i}\boldsymbol{\nabla}^h\boldsymbol{\cdot}\bar{\boldsymbol{u}}^h\right)\right) \frac{\mathcal{E}}{2} + {\rm c.c.}\\ & \quad \qquad \mathrm{at} \ z = 0 ,\\ \partial_z \varPhi'_c & = \hat{\boldsymbol{z}} \boldsymbol{\cdot} \boldsymbol{u}^{h\prime}_w \times \bar{\boldsymbol{\zeta}}^h = \frac{a\sigma}{2k} \frac{\cosh[2kH]}{\sinh[kH]} \boldsymbol{k} \boldsymbol{\cdot} \partial_z \bar{\boldsymbol{u}}^h \frac{\mathcal{E}}{2} + {\rm c.c.}\quad \mathrm{at} \ z ={-} H \end{aligned} \right\} \end{align}

(cf. WCI22, section 5.2). This system is functionally consistent with (5.13) for $\varPhi '_d$ except for the replacement of drag-effect forcing there with current-interaction forcing and the presumption that the slow-time coordinate is $\tau = \mu t$. As in § 5 for $\boldsymbol {u}'_d$, the elliptical system (6.6) is to be solved as $\bar {\boldsymbol {u}}$ evolves, along with the evaluation of (6.4), to obtain $\boldsymbol {u}'_c$ using (6.2) and (6.5).

If the currents themselves are three-dimensional (3-D), then (6.6) cannot be reduced to a 1-D vertical problem like (5.17), but they can be if $\bar {\boldsymbol {u}} = \bar {\boldsymbol {u}}^h(z)$ only, in which case the familiar Doppler-shifting effect on the frequency and wave vertical profile is recovered (WCI22, appendix J.1).

7. Synthesis of the wave perturbations

In §§ 5 and 6 the effects on the wave modifications due to wave drag and CEW are expressed in functionally similar ways. To make a more generally useful theory, it seems worthwhile to combine them, taking advantage of the fact that the associated $\boldsymbol {u}'$ fields are determined by linear relationships, hence are superposable. Both formulations include the possibility of a slow-time evolution of the primary wave amplitude, $a(\tau )$, and they have a solvability condition (§ 8) that brings these two effects together. Therefore, a formal synthesis is now made under the scale-ordering assumption that $\nu \sim \mu$ (§ 3), and this can be considered the distinguished limit for this multiscale asymptotic theory. If these parameters had disparate magnitudes, then either the drag or current effects would be dominant, with corresponding consequences for the wave-averaged evolution (§ 10); however, the assumption that they are comparable is the more general case that incorporates both possibilities. The combined slow-time coordinate is thus the fastest of the preceding $\tau$ coordinates.

In this section and the following one, the wave amplitude $a(\tau )$ is a complex function only of slow time (§ 4); thus, the drag-related formulae are slightly modified from § 5.

Consider a combined wave-modification velocity,

(7.1)\begin{equation} \boldsymbol{u}'_m = \boldsymbol{u}'_d + \boldsymbol{u}'_c , \end{equation}

which itself is decomposed into vortical and irrotational components, $\boldsymbol {u}_m^{\zeta \prime }$ and $\boldsymbol {u}_m^{\varPhi \prime }$. The accompanying fields $\eta '_m$, $\phi '_m$ and $\varPhi '_m$ are analogous combinations.

The goal is a wave-averaged theory, and thus the ‘fast’ dependence should be separated from the wave-averaged dependence. This is done by factoring $\mathcal {E}/2$ in the various $\boldsymbol {u}'_m$ relations and focusing on their generally complex coefficient fields; e.g. for

(7.2)\begin{equation} \varPhi'_m = Q \frac{\mathcal{E}}{2} + {\rm c.c.}, \end{equation}

the complex factor $Q(\boldsymbol {x},\tau )$ is entirely a wave-averaged quantity. For $\boldsymbol {u}'_m$, choose an analogous complex factor, $\boldsymbol {U}_m = \boldsymbol {U}^\zeta + \boldsymbol {U}^\varPhi$, so that

(7.3)\begin{equation} \boldsymbol{u}'_m = \boldsymbol{U}_m \frac{\mathcal{E}}{2} + {\rm c.c.} \end{equation}

For the vortical velocity components in (5.10) and (6.4), the combined factor is

(7.4)\begin{equation} \boldsymbol{U}^\zeta_m ={-} \frac{{\rm i} \mathcal{D}c_1(gF(z))^2 |a|a k}{\sigma^3}\boldsymbol{k} + \frac{ga}{\sigma^2} {\left(\begin{array}{c} {\rm i}k^y(\,f +\bar{\zeta}^z)F(z) - k\bar{\zeta}^yG(z) \\ k\bar{\zeta}^xG(z) - {\rm i}k^x(\,f +\bar{\zeta}^z)F(z) \\ {\rm i}\hat{\boldsymbol{z}}\boldsymbol{\cdot}\boldsymbol{k} \times \bar{\boldsymbol{\zeta}}^h F(z)\end{array}\right)} . \end{equation}

For the irrotational velocity component, the elliptical PDE system for the combined factor $Q$ from (5.13) and (6.6) is the following:

(7.5)\begin{align} \left. \begin{aligned} & [\nabla^2 + 2{\rm i} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\nabla}^h - k^2] Q ={-} c_1 \frac{(gF(z))^2 |a|a k}{\sigma^2} ( \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{\nabla}^h\mathcal{D} +{\rm i} k^2 \mathcal{D} )\\ & \qquad\qquad\qquad\qquad\quad\qquad + \frac{ag}{\sigma} \left( F(z)\boldsymbol{k} \boldsymbol{\cdot} \nabla^2 \bar{\boldsymbol{u}}^h - {\rm i}G(z) k \nabla^2 \bar{w} \right) ,\\ & \left[\partial_z - \frac{\sigma^2}{g} \right] Q = 2 {\rm i} \sigma \partial_\tau a + a\sigma\left({-}2 \boldsymbol{k}\boldsymbol{\cdot}\bar{\boldsymbol{u}}^h + \frac{g}{\sigma^2}\boldsymbol{k}\boldsymbol{\cdot}\partial_z\bar{\boldsymbol{u}}^h + {\rm i}\boldsymbol{\nabla}^h\boldsymbol{\cdot}\bar{\boldsymbol{u}}^h\right) \quad \mathrm{at} \ z = 0 ,\\ & \partial_z Q = \frac{a\sigma}{2k} \frac{\cosh[2kH]}{\sinh[kH]} \boldsymbol{k} \boldsymbol{\cdot} \partial_z \bar{\boldsymbol{u}}^h \quad \mathrm{at} \ z ={-} H . \end{aligned} \right\} \end{align}

Notice that this system has combined the forcing term $\propto \partial _\tau a$ that is present in both (5.13) and (6.6). Finally, the associated irrotational velocity factor is

(7.6)\begin{equation} \boldsymbol{U}^\varPhi_m = \frac{1}{\sigma} ( Q \boldsymbol{k} -{\rm i}\boldsymbol{\nabla} Q ). \end{equation}

In both (7.5) and (7.6) the 3-D gradient operator $\boldsymbol {\nabla }$ has been expanded to include its action both on $\mathcal {E}$ and on the wave-averaged factor $Q$.

Once $a(\tau )$ is known (§ 8), then (7.3)–(7.6) determine the $\boldsymbol {u}'_m$ field sufficiently to evaluate the WEC terms in the current evolution equations (§§ 9 and 10).

8. Solvability condition for wave amplitude $a(\tau )$

In perturbation expansions, a common technique is to enforce an orthogonality condition between a leading-order quantity – here the primary wave $\boldsymbol {u}'_w$ – and its small perturbation $\boldsymbol {u}'_m$ to avoid ‘resonant’ forcing of the perturbation that upsets the asymptotic ordering in space and/or time. Combining (5.13) and (6.6) and formally defining the inhomogeneous terms as $R_\bullet$ factors in the PDE and boundary conditions gives the following PDE system for the wave-correction velocity potential:

(8.1)\begin{equation} \left. \begin{aligned} \nabla^2 \varPhi'_m & = R_i ,\\ \left( \partial_z - \frac{\sigma^2}{g} \right) \varPhi'_m & = R_0 \quad \mathrm{at} \ z = 0 ,\\ \partial_z \varPhi'_m & = R_H \quad\mathrm{at} \ z ={-}H . \end{aligned} \right\} \end{equation}

The homogeneous solution to this system is the $\phi '$ eigenmode, $\mathcal {G} = F(z)\mathcal {E}$, with the eigenvalue $\sigma$ and wavenumber vector $\boldsymbol {k}$, i.e. the primary wave (§ 4). There is an analogous $\mathcal {G}^{\dagger} = F(z)\mathcal {E}^{\dagger}$ for any other $\boldsymbol {k}^{\dagger}$ with the same wavenumber magnitude $k$ (i.e. for a circle in $\boldsymbol {k}$-space of radius $k$).

The solvability condition that suppresses excitation of a homogeneous solution in $\varPhi '$ is obtained as follows (WCI22, appendix E). Multiply the PDE in (8.1) by the complex conjugate of $\mathcal {G}^{\dagger}$ (denoted by a $*$ superscript), integrate over the domain volume and then integrate by parts to insert the boundary conditions from (8.1). The result is the solvability condition:

(8.2)\begin{equation} \iiint {\rm d}\boldsymbol{x} \,\mathcal{G}^{{{\dagger}} *} R_i(\boldsymbol{x}) = \iint {\rm d}\boldsymbol{x}^h \,(\mathcal{G}^{{{\dagger}} *}(0)R_0(\boldsymbol{x}^h) - \mathcal{G}^{{{\dagger}} *}({-}H)R_H(\boldsymbol{x}^h) ) . \end{equation}

Define a horizontal wavenumber projection operator,

(8.3)\begin{equation} \mathcal{P}^{\dagger}[ {\,{\cdot}\,} ] \equiv \frac{1}{\mathcal{A}} \iint{\rm d}\boldsymbol{x}^h\,\mathcal{E}^{{{\dagger}}*}[\boldsymbol{\,{\cdot}\,}], \end{equation}

with the property that

(8.4)\begin{equation} \mathcal{P}^{\dagger}[ \mathcal{E}' ] = \delta_{\boldsymbol{k}',\boldsymbol{k}^{\dagger}} . \end{equation}

In the above, $\mathcal {A}$ is the horizontal area of the domain.

The solvability condition (8.2) can be more compactly expressed as an integral only over $z$:

(8.5)\begin{equation} \int^0_{{-}H} F(z) \mathcal{P}^{\dagger}[ R_i] \,{\rm d}z = F(0) \mathcal{P}^{\dagger}[ R_0]- F({-}H)\mathcal{P}^{\dagger}[ R_H] . \end{equation}

This condition must be satisfied for all $\boldsymbol {k}^{\dagger}$. Because $R_0$ contains a term with $(\partial _\tau a) \mathcal {E}$ and all the other $R_\bullet$ terms depend only on $a \mathcal {E}$, this condition gives a first-order differential equation for $a(\tau )$. This assures that $Q(\boldsymbol {x})$ will not have a component proportional to $F(z)$ without any associated $\boldsymbol {x}^h$ structure. Once the $a(\tau )$ are determined from (8.5), then (7.5) is a well-behaved PDE system for $\varPhi '_m$.

Notice that the solvability condition (8.5) for $a(\tau )$ depends only on $\mathcal {D}$, the primary wave properties, and the current $\bar {\boldsymbol {u}}$, i.e. it does not require knowing the solution for $\boldsymbol {u}'_m$. So $a(\tau )$ can be determined independently.

In advance of more complete solutions with both drag and current effects on the waves, the simple situation of drag without any current or Coriolis effect (i.e. $\bar {\boldsymbol {u}} = f = 0$) and with a horizontally uniform canopy (i.e. $\boldsymbol {\nabla }^h \mathcal {D} = 0$) can be evaluated using (8.5). In this case, the projection operations are trivial except for $\boldsymbol {k}^{\dagger} = \boldsymbol {k}$, and the solvability condition becomes

(8.6)\begin{equation} \int_{{-}H}^0 F(z) \left(-{\rm i} \frac{(g^2 |a|a k^3)}{\sigma^2} F^2 \mathcal{D} c_1 \right) {\rm d}z = F(0) 2{\rm i}\sigma\partial_\tau a . \end{equation}

After rearrangement, one obtains

(8.7a,b)\begin{equation} \partial_\tau a ={-} J |a|a , \quad J = c_1 \frac{g^2k^3}{2\sigma^3} \int_{{-}H}^0 F^3(z)\mathcal{D}(z)\,{\rm d}z . \end{equation}

With a slow-time initial condition of $a(\tau _0) = a_0 > 0$, the solution is a real, positive amplitude function:

(8.8a,b)\begin{equation} a(\tau) = \frac{a_0}{1 + r(\tau-\tau_0)} , \quad r = a_0 J = c_1 \frac{g^2k^3a_0}{2\sigma^3} \int_{{-}H}^0 F^3(z)\mathcal{D}(z) \,{\rm d}z . \end{equation}

This shows a slow-time, inverse-linear decay in the primary wave amplitude at a rate $r$ that depends on the primary wave properties, including its initial amplitude, the canopy density and the drag coefficient. (A similar linear-decay dependence for wave amplitude is found for canopy-drag effects in Dalrymple et al. (Reference Dalrymple, Kirby and Hwang1984), but with somewhat different functional dependences than in $r$ associated with how $\mathcal {D}$ is specified.)

With non-zero currents ($\bar {\boldsymbol {u}} \ne 0$), the implications of (7.4)–(7.6) and (8.5) must mostly wait for combined drag and current solutions in a computational code built on top of the equations derived here. However, an important issue is the possibility of wave–current resonance, which is discussed extensively in WCI22 (chapters 7 and 8 and appendices J and K). If the current field has a wavenumber component $\boldsymbol {p}$, and if $|\boldsymbol {k}^{\dagger} | = | \boldsymbol {k} \pm \boldsymbol {p}| = k$, then there can be a slow-time energy exchange between the amplitudes $a(\tau ; \boldsymbol {k})$ and $a(\tau ; \boldsymbol {k}^{\dagger} )$. This occurs through the quadratic product of primary wave and current fields in the quasi-linear advection terms in (6.1). That is, new primary waves can arise through the wave–current interaction, and wave energy can be exchanged among these resonant components. In principle, their amplitudes can become as large as that of the original primary wave. Several idealized examples of wave–current resonance for particular $\bar {\boldsymbol {u}}(\boldsymbol {x})$ fields are analysed in WCI22, and in all such cases the slow-time amplitude behaviour is oscillatory between primary wave modes under the condition of steady currents. For most simple currents $\bar {\boldsymbol {u}}(\boldsymbol {x})$ and primary wave $\boldsymbol {u}'_w$ fields, however, there is no resonance and no current-induced change in $a$. If there are multiple primary waves present (Appendix B), then the solvability conditions are still valid, although this would compound the resonant possibilities.

As the currents evolve with WEC and other influences (§ 10), it is presently unknown how often wave–current resonance might be an important phenomenon. (All of the papers cited in § 1 ignore this possibility except for WCI22.) A computational simulation strategy for dealing with it is also discussed in WCI22 (chapter 9), but it has not yet been implemented. For now, it remains an open issue for solving a quasi-linear wave model like the one here.

9. Wave-averaged drag force

The wave-averaged drag force is

(9.1)\begin{equation} \bar{\boldsymbol{D}}^h ={-} \mathcal{D} \langle |\boldsymbol{u}^h| \boldsymbol{u}^h \rangle , \end{equation}

and the total velocity is defined in (2.4) and (2.5). In support of the WEC model in § 10, the goal is to obtain a wave-averaged formula for $\bar {D}$ without having to do wave-resolved evaluations, nor even having to solve the elliptic problems for $\boldsymbol {u}^{\varPhi \prime }_d$ or $\boldsymbol {u}^{\varPhi \prime }_c$.

For the results in this and the following section, $a$ is real and positive (§ 4), and the approximation (5.4) is made for $\boldsymbol {D}_w^{h\prime }$.

As a first step, consider an approximate form for the quadratic term in the drag law for a sum of velocities, $\boldsymbol {u} + \boldsymbol {v}$, when $|\boldsymbol {v}| \ll |\boldsymbol {u}|$:

(9.2)\begin{align} | \boldsymbol{u} + \boldsymbol{v} | (\boldsymbol{u} + \boldsymbol{v}) &\approx |\boldsymbol{u} + \boldsymbol{v}| \boldsymbol{u} + |\boldsymbol{u}|\boldsymbol{v} \nonumber\\ &\approx |\boldsymbol{u}|\boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{v})\frac{\boldsymbol{u}}{|\boldsymbol{u}|} + |\boldsymbol{u}|\boldsymbol{v} , \end{align}

neglecting second-order terms in $\boldsymbol {v}$.

Assume, as in (3.2) and §§ 3–5, that the largest velocity is $\boldsymbol {u}^{h\prime }_w$. Inserting (2.4) and (2.5) into (9.1) – with the preceding approximation and $\nu \sim \mu$ ordering for $\boldsymbol {u}^{h\prime }_m$ and $\bar {\boldsymbol {u}}^h$ – gives

(9.3)\begin{equation} \bar{\boldsymbol{D}}^h ={-}\mathcal{D} \left\langle \left( |\boldsymbol{u}^{h\prime}_w|\boldsymbol{u}^{h\prime}_w + \boldsymbol{u}^{h\prime}_w\boldsymbol{\cdot}(\bar{\boldsymbol{u}}^h+\boldsymbol{u}^{h\prime}_m) \frac{\boldsymbol{u}^{h\prime}_w}{|\boldsymbol{u}^{h\prime}_w|} + |\boldsymbol{u}^{h\prime}_w| (\bar{\boldsymbol{u}}^h+\boldsymbol{u}^{h\prime}_m) \right)\right\rangle . \end{equation}

The first right-hand-side term has already been evaluated to have a zero average in (5.3a,b). Because both $\boldsymbol {u}'_w$ and $\boldsymbol {u}'_m$ are $\propto \mathcal {E}$ in their fast-time dependence, here two other wave quantities that are either linear or cubic in their signed dependence on $\mathcal {E}$ also have zero average because of the odd sign symmetry in different phases of $\varTheta$. Thus,

(9.4)\begin{equation} \left. \begin{aligned} \langle |\boldsymbol{u}^{h\prime}_w| \boldsymbol{u}^{h\prime}_w \rangle & = 0, \\ \left\langle (\boldsymbol{u}^{h\prime}_w\boldsymbol{\cdot}\boldsymbol{u}^{h\prime}_m) \frac{\boldsymbol{u}^{h\prime}_w}{|\boldsymbol{u}^{h\prime}_w|} \right\rangle & = 0, \\ \langle |\boldsymbol{u}^{h\prime}_w| \boldsymbol{u}^{h\prime}_m \rangle & = 0 . \end{aligned} \right\} \end{equation}

This is because they are proportional to the trigonometric functions of $\varTheta$ in a way that has equal magnitude and opposite sign over the interval of $[0,2{\rm \pi} ]$.

In contrast, $(\boldsymbol {u}^{h\prime }_w \boldsymbol {\cdot } \bar {\boldsymbol {u}}^h )\boldsymbol {u}^{h\prime }_w / |\boldsymbol {u}^{h\prime }_w|$ and $|\boldsymbol {u}^{h\prime }_w|$ are single-signed and have non-zero averages. (These cancellations also extend to the relatively weak higher harmonics in $\boldsymbol {u}'_w$ (figure 1); see Appendix A.) Furthermore, for terms linear in the magnitude of $|\cos [\varTheta ]|$, the averaging operation gives

(9.5)\begin{equation} \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} |{\cos[\varTheta]}| \,{\rm d}\varTheta = \frac{2}{\rm \pi} . \end{equation}

Thus, with (4.3) the wave-averaged drag force (9.1) is

(9.6)\begin{align} \bar{\boldsymbol{D}}^h &={-}\mathcal{D} \left\langle (\boldsymbol{u}^{h\prime}_w \boldsymbol{\cdot} \bar{\boldsymbol{u}}^h)\frac{\boldsymbol{u}^{h\prime}_w}{|\boldsymbol{u}^{h\prime}_w|} + (|\boldsymbol{u}^{h\prime}_w|+\bar{\boldsymbol{u}}^h)\bar{\boldsymbol{u}}^h \right\rangle \nonumber\\ &={-} \mathcal{D} \left( \frac{2ag}{{\rm \pi} \sigma} F(z) \left(k\bar{\boldsymbol{u}}^h + \frac{\boldsymbol{k}\boldsymbol{\cdot}\bar{\boldsymbol{u}}^h}{k} \boldsymbol{k} \right) +|\bar{\boldsymbol{u}}^h|\bar{\boldsymbol{u}}^h \right) . \end{align}

In an ad hoc way (i.e. beyond leading order in $\mu$), the last term is added to have a well-behaved limit in the absence of a primary wave when the drag is entirely due to the current.

Thus, the canopy exerts a horizontal drag force on the currents that is enhanced by the primary wave. The spatial structure of this drag is the product of the canopy distribution $\mathcal {D}(\boldsymbol {x})$, the wave eigenmode profile $F(z)$ and the horizontal current itself $\boldsymbol {u}^h(\boldsymbol {x})$. It is retarding in the direction of the current in all three of the terms in (9.6); i.e. the drag work, $\bar {\boldsymbol {D}}^h \boldsymbol {\cdot } \bar {\boldsymbol {u}}^h$, is negatively proportional to each of $(\bar {\boldsymbol {u}}^h)^2$, $(\boldsymbol {k}\boldsymbol {\cdot }\bar {\boldsymbol {u}}^h)^2$ and $|\bar {\boldsymbol {u}}^h|(\bar {\boldsymbol {u}}^h)^2$. The first two terms in (9.6) have a scaling magnitude of $(\nu \mu ) C^2/\ell$, while the last term is $(\nu \mu ^2/\epsilon ) C^2/\ell$, where $\ell$ is a length scale for both the waves and the currents. As desired, (9.6) is expressed entirely in terms of wave-averaged quantities. Notice that there is no dependence on $\boldsymbol {u}'_m$, so the wave perturbations need not be solved for, apart from whatever $a(\tau )$ dependence is caused by the solvability condition (§ 8).

It would be straightforward to extend the representation of $\bar {\boldsymbol {D}}^h$ to include a bottom drag stress and again have a wave enhancement over the dependence on only $\bar {\boldsymbol {u}}^h$. The form is likely to be similar to (9.6), which is functionally different from presently used parametrization formulae for wave-enhanced bottom stress in a turbulent bottom boundary layer (e.g. Feddersen et al. Reference Feddersen, Guza, Elgar and Herbers2000). This issue is worth further consideration.

10. Wave effects on the currents (WEC)

The approximate wave-averaged current evolution equations for the general model (2.1) are

(10.1)\begin{equation} \left. \begin{aligned} {}[\partial_\tau + \bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla}] \bar{\boldsymbol{u}} + f\hat{\boldsymbol{z}} \times \boldsymbol{u}^h + \boldsymbol{\nabla} \bar{\phi} &={-} \boldsymbol{\nabla} \bar{B} + \langle \boldsymbol{u}' \times \boldsymbol{\zeta}' \rangle + \bar{\boldsymbol{D}}^h,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \bar{\boldsymbol{u}} &= 0 , \end{aligned} \right\} \end{equation}

with vertical boundary conditions of $\bar {w} = 0$ at $z = 0$ and $-H$ (i.e. a rigid flat top and bottom; see WCI22, chapter 6, but with $\bar {\boldsymbol {D}}^h$ added). The time coordinate is again $\tau = \mu t$, based on the advective time scale of current evolution, although one can also view this as valid for $\tau = \nu t$ when $\nu \sim \mu$ (§ 7).

The added forces due to the presence of waves are on the right-hand side of (10.1). The wave-averaged drag $\bar {D}^h$ is in (9.6). The wave-averaged Bernoulli head is

(10.2)\begin{equation} \bar{B} = \tfrac{1}{2} \langle \boldsymbol{u}^{\prime 2} \rangle / 2 \end{equation}

(WCI22, chapter 6), and it does not have to be evaluated explicitly because the ‘pressure Poisson’ equation from the incompressibility condition in (10.1) gives a combined expression for the potential force, $\bar {\phi } + \bar {B}$, as part of the time integration of (10.1) for $\bar {\boldsymbol {u}}$. This assumes that the lateral boundary conditions are consistent with their combination, as certainly is true for horizontal periodicity. The approximate rigid-lid condition on $\bar {w}$ is associated with a higher-order kinematic surface boundary condition that is a diagnostic relation for wave-averaged sea level, $\bar {\eta }$, which is much smaller than $\eta _w'$ (WCI22, appendix B).

The total wave-averaged vortical force is a combination of drag and current effects on the waves:

(10.3)\begin{equation} \boldsymbol{VF} = \langle \boldsymbol{u}' \times \boldsymbol{\zeta}' \rangle \approx \langle\boldsymbol{u}'_w \times (\boldsymbol{\zeta}'_d + \boldsymbol{\zeta}'_c )\rangle , \end{equation}

where the wave vorticities are in (5.9) and (6.3). In a way that is familiar in the papers referenced in the introduction, the wave-averaged wave vortex force associated with the currents is

(10.4)\begin{equation} \left. \begin{aligned} \boldsymbol{SVF} & = \langle \boldsymbol{u}'_w \times \boldsymbol{\zeta}'_c \rangle = \boldsymbol{u}^{St} \times (\,f\hat{\boldsymbol{z}} + \bar{\boldsymbol{\zeta}}) - \boldsymbol{\nabla} \bar{C} ,\\ \bar{C} & = \hat{\boldsymbol{z}} \boldsymbol{\cdot} \int_{{-}H}^0 \boldsymbol{u}^{St}(z')\,{\rm d}z \times (\,f\hat{\boldsymbol{z}} + \bar{\boldsymbol{\zeta}}) \end{aligned} \right\} \end{equation}

(WCI22, section 6.1 and appendix G). The right-hand-side terms are evaluated using (4.3) and (6.4). The first term is the familiar Stokes–Coriolis and Stokes–vorticity forces. The second term $\boldsymbol {\nabla } \bar {C}$ does not have to be evaluated explicitly because it, too, combines into the total potential force for the $\bar {\boldsymbol {u}}$ evolution.

The WEC vortical force associated with the wave-drag vortex force is evaluated using (4.3) and (5.9):

(10.5)\begin{align} \boldsymbol{CVF} &= \langle \boldsymbol{u}'_w \times \boldsymbol{\zeta}'_d \rangle ={-}c_1\frac{g^3a^3k^2}{2\sigma^4} G(z)\partial_z[\mathcal{D}F^2] \boldsymbol{k}. \end{align}

Because $\mathcal {D}$, $F$, $G$ and $\partial _zF^2$ are all positive from (1.2) and (4.4), this force can be either retarding or accelerating in the direction of primary wave propagation, depending on the depth profile of the canopy structure, $\partial _z\mathcal {D} \propto \partial _z \alpha (\boldsymbol {x})$ in (1.2). A scaling estimate is $\boldsymbol {CVF} \sim \nu \epsilon ^2 C^2/\ell$, compared to $\boldsymbol {SVF} \sim \mu \epsilon ^2 C^2/\ell$; i.e. they are of similar magnitude if $\nu \sim \mu$. As remarked following (9.6), the wave-enhanced drag terms are $\bar {\boldsymbol {D}}^h \sim \nu \mu C^2/\ell$; therefore, all three wave forces are comparable if $\epsilon ^2 \sim \mu \sim \nu$.

Thus, the final form of the current momentum balance is

(10.6)\begin{align} & [\partial_\tau + \bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla}] \bar{\boldsymbol{u}} + f\hat{\boldsymbol{z}} \times \boldsymbol{u}^h + \boldsymbol{\nabla}\bar{\varPi} \nonumber\\ & \quad = \boldsymbol{VF} + \bar{\boldsymbol{D}}^h \nonumber\\ & \quad ={-}c_1\frac{g^3a^3k^2}{2\sigma^4}G(z)\partial_z[\mathcal{D}F^2] \boldsymbol{k} + \boldsymbol{u}^{St} \times (\,f\hat{\boldsymbol{z}} + \bar{\boldsymbol{\zeta}}) + \bar{\boldsymbol{D}}^h , \end{align}

where $\bar {\varPi } = \bar {\phi } + \bar {B} + \bar {C}$ is the combined force potential function, $\boldsymbol {VF} = \boldsymbol {CVF}+\boldsymbol {SVF}$ is the combined vortex force and $\bar {D}^h$ is the combined drag force on the currents (9.6). The right-hand-side terms are, respectively, the vortex force due to canopy drag on the waves $\boldsymbol {CVF}$, the Stokes vortex forces $\boldsymbol {SVF}$ and the wave-enhanced current drag. This equation is the principal result in this paper.

Notice that (10.6) is entirely specified in terms of $\mathcal {D}$, $\bar {\boldsymbol {u}}$ and the primary wave properties. Among the latter is the slow-time evolution of the primary wave amplitude $a(\tau )$ as determined from the solvability condition (8.5) that includes canopy wave damping and possible wave–current resonances. This system is readily implemented in an LES or circulation model, with the most challenging new element being the solvability condition (§ 8) that has heretofore been neglected in wave–current interaction studies.

The vortex force is a product of the quasi-linear CEW theory (§ 6) at leading order in $\mu$ (WCI22, section 5.1); the higher-order CEW and WEC terms presented in WCI22 (sections 5.2 and 6.6) are not included here, although they could be, as their forms are known. But this would require explicit solutions of the elliptic problems in §§ 5–7 that are not required for (10.6). Recall that another approximation is indicated in (5.3a,b), where the weak higher harmonics in $\boldsymbol {D}^{h\prime }_w$ are neglected; this simplifies the theory, and Appendix A shows that this approximation is fairly accurate for representing the further consequences of $\boldsymbol {D}^{h\prime }_w$ in § 5, hence in $\bar {\boldsymbol {D}}^h$ as well.

11. Weaker waves

For the preceding $\boldsymbol {u}_c'$, $\boldsymbol {u}^{St}$ and $\boldsymbol {SVF}$ results, it does not matter how small $\epsilon \ll 1$ is. However, this is not true for the canopy effects. They were derived for the scaling relation (3.2), i.e. $\mu \ll \epsilon$, which yields the analytic relations in § 5 that carry through to the current effects in § 10. If the wave orbital velocities are weaker, with $\mu \sim \epsilon$, then the total drag force (1.1) at leading order, viz.

(11.1)\begin{equation} \boldsymbol{D}^h ={-} \mathcal{D} | \boldsymbol{u}^{h\prime}_w + \bar{\boldsymbol{u}}^h | (\boldsymbol{u}^{h\prime}_w + \bar{\boldsymbol{u}}^h) \equiv \boldsymbol{D}^{h\prime}_w + \bar{\boldsymbol{D}}^h , \end{equation}

does not have explicit analytic expressions for its fluctuation and wave-averaged drag-force components for the primary wave eigenmode (4.3) and a general current field $\bar {\boldsymbol {u}}(\boldsymbol {x})$. Thus, at every location $\boldsymbol {x}$ and slow time $\tau$, these forces must be computed by taking a numerical average for $\bar {\boldsymbol {D}}^h$ over the primary wave period as an auxiliary calculation for the current evolution equation.

The resulting $\boldsymbol {D}^{h\prime } = \boldsymbol {D}^h - \bar {\boldsymbol {D}}^h$ can then be used to evaluate (5.6) as

(11.2)\begin{align} \partial_t\boldsymbol{\zeta}_d' &=\mathrm{curl}[\boldsymbol{D}^{h\prime}]\nonumber\\ &= {\left(\begin{array}{@{}c@{}} -\partial_z D^{y \prime} \\ +\partial_z D^{x\prime} \\ - \hat{\boldsymbol{z}}\boldsymbol{\cdot} \boldsymbol{\nabla}^h \times \boldsymbol{D}^{h\prime}\end{array}\right)} , \end{align}

with discretized spatial derivatives in the right-hand-side quantities. Similarly, a spatially discretized $\boldsymbol {\nabla }^h \boldsymbol {\cdot } \boldsymbol {D}^{h\prime }$ can be used to evaluate the solvability condition (8.5), also by numerical averaging with the projection operator $\mathcal {P}$. There is no need to solve for $\boldsymbol {u}_m'$ with weaker waves.

For the current evolution equation analogous to (10.6), the current drag $\bar {\boldsymbol {D}}^h$ can be inserted directly, and the canopy vortex force $\boldsymbol {CVF}$ is calculated as follows:

(11.3) \begin{align} \boldsymbol{CVF} &= \langle \boldsymbol{u}'_w \times \boldsymbol{\zeta}'_d \rangle = \left\langle\boldsymbol{u}'_w \times \int^t\mathrm{curl}[\boldsymbol{D}^{h\prime}]\,{\rm d}t\right\rangle ={-}\left\langle \left(\int^t \boldsymbol{u}'_w \,{\rm d}t\right) \times \mathrm{curl}[\boldsymbol{D}^{h\prime}] \right\rangle \nonumber\\ &= \left\langle\left(\begin{array}{@{}c@{}} \left(\int^t v'_w \,{\rm d}t\right) (\hat{\boldsymbol{z}} \boldsymbol{\cdot} \boldsymbol{\nabla}^h \times \boldsymbol{D}^{h\prime}) + \left(\int^t w'_w \,{\rm d}t\right)\partial_z\boldsymbol{D}^{x\prime} \\ \left(\int^t w'_w \,{\rm d}t\right)\partial_z\boldsymbol{D}^{y\prime} - \left(\int^t u'_w \,{\rm d}t\right)(\hat{\boldsymbol{z}} \boldsymbol{\cdot} \boldsymbol{\nabla}^h \times\boldsymbol{D}^{h\prime})\\ - \left(\int^t \boldsymbol{u}^{h\prime}_w \,{\rm d}t\right) \boldsymbol{\cdot} \partial_z \boldsymbol{D}^{h\prime} \end{array}\right)\right\rangle . \end{align}

With the primary wave (4.3) and real $a$, this becomes

(11.4)\begin{align} \boldsymbol{CVF} &= \frac{g a}{\sigma^2} \langle F \sin[\varTheta] (\hat{\boldsymbol{z}}\boldsymbol{\cdot} \boldsymbol{\nabla}^h \times \boldsymbol{D}^{h\prime}) (\hat{\boldsymbol{z}} \times \boldsymbol{k}) \nonumber\\ & \quad + k G \cos[\varTheta] \partial_z\boldsymbol{D}^{h\prime} \nonumber\\ & \quad + F \sin[\varTheta] (\boldsymbol{k} \boldsymbol{\cdot} \partial_z\boldsymbol{D}^{h\prime}) \hat{\boldsymbol{z}} \rangle . \end{align}

Again a fast-time numerical average over a primary wave period is needed, as indicated by the angle brackets, where $\varTheta$ and $\boldsymbol {D}^{h\prime }$ are the fast-time variables. In (11.4) the first two lines are horizontal components of $\boldsymbol {CVF}$, and the final line is a vertical component. Thus, in this weaker wave case, the wave-averaged effects on the currents (WEC) are still well determined, but at the cost of added numerical averages of the wave-induced drag and canopy vortex force.

If the primary wave is even weaker, with $\boldsymbol {u}'_w \ll \bar {\boldsymbol {u}}$ and $\epsilon \ll \mu$, then $\bar {\boldsymbol {D}}^h \approx - \mathcal {D} | \bar {\boldsymbol {u}}^h | \bar {\boldsymbol {u}}^h$, and $\boldsymbol {CVF}$ is much smaller than, for example, the advective tendencies in the current evolution equation. That is, wave effects on currents (WEC) are negligible, and the most important interaction is current effects on waves (CEW) (e.g. a Doppler shift in $\sigma$) and a fluctuating canopy drag force on the waves.

12. Multiple primary waves

Now reverting to the principal scaling assumption (3.2), a wave-averaged theory is still achievable for multiple components in the primary wave field, i.e. for $N$ components,

(12.1)\begin{equation} \eta_w' = \sum_{n=1}^{N} a_n \frac{\mathcal{E}_n}{2} + {\rm c.c.} , \end{equation}

with amplitudes $a_n$ and exponential phase factors $\mathcal {E}_n = \exp [\textrm {i} \varTheta _n]$ with $\varTheta _n = \boldsymbol {k}_n\boldsymbol {\cdot }\boldsymbol {x} - \sigma _n t$. Again, these waves are eigenmodes with $\sigma _n$ as in (4.2) for $\boldsymbol {k}_n$, accompanying fields as in (4.3), vertical structure functions $F_n$ and $G_n$ as in (4.4), and a Stokes drift $\boldsymbol {u}^{St}(z)$ as in (4.6) that is a sum over the $N$ separate contributions.

However, such a theory has added complexity in dealing with the drag effects on both the waves and currents. In particular, for use in the drag force on the waves (§ 5), the formula (5.1) is evaluated with the primary wave horizontal velocity,

(12.2)\begin{equation} \boldsymbol{u}^{h\prime}_w = \sum_n \frac{ga_n}{\sigma_n} F_n(z)\boldsymbol{k}_n \frac{\mathcal{E}_n}{2} + {\rm c.c.}, \end{equation}

which then propagates through the equivalents of §§ 5–10. The source of the complexity is that analytic expressions are no longer available for $D_w^{h\prime }$, $\bar {D}^h$ and $\boldsymbol {CVF}$. Besides necessary numerical averages over the primary wave periodic interval for these quantities and in the solvability condition (8.5) (as in § 11), the simplifications in (9.4) are no longer strictly valid, so that $\boldsymbol {u}'_m =\boldsymbol {u}'_d + \boldsymbol {u}'_c$ contributions to $\bar {D}^h$ would also need to be solved for computationally with the elliptic PDE system (7.5). This adds up to a considerable computational burden beyond what is involved with the analytic expressions in §§ 9 and 10.

As discussed in Appendix B, some quantitative exploration of $\boldsymbol {D}^h$ with multiple waves gives support for thinking that there may be useful approximate shortcuts to the full programme of numerical evaluations in the preceding paragraph. In particular, solving for $\boldsymbol {u}'_m$ to evaluate $\bar {\boldsymbol {D}}^h$ in (9.3) might not be necessary if their associated quantities in (9.4) are not too large compared to the $\boldsymbol {u}'_w$ contributions. Of course, this would need to be tested in particular problems of interest.

On the other hand, the CEW corrections in §§ 6–8 are easily generalized with multiple primary waves by simple superposition (as discussed in WCI22), as are their associated WEC terms in § 10. This superposition can even be extended to broadband wave spectra, e.g. for $\boldsymbol {u}^{St}$ in $\boldsymbol {SVF}$.

For weaker waves (§ 11), multiple primary waves add only a modest burden to the indicated numerical evaluations, with $\boldsymbol {u}_w'$ in the formulae (11.1) and (11.3) understood to be the summation over the multiple primary wave components. The wave perturbations $\boldsymbol {u}_d'$ and $\boldsymbol {u}_c'$ continue to be negligible in this case.

13. Summary and prospects

This paper presents a theory for the wave-averaged effects of surface gravity waves and currents in the presence of a vegetative canopy, suitable for a modified LES or fine-scale circulation model. The principal new effects – compared to previous wave-averaged models with only a mean-current drag force and the Stokes vortex forces associated with the currents (e.g. Yan et al. Reference Yan, McWilliams and Chameki2021) – are a wave-amplitude damping (8.5), a wave-enhanced current drag (9.6) and a wave-drag vortex force (10.5). The model includes WEC for current evolution (§ 10) plus wave damping and possible resonances for the primary wave amplitude(s), $a(\tau ; \boldsymbol {k})$. For a single primary wave field, the set of resonance amplitudes is limited to a circle in $\boldsymbol {k}$ space, which should not be burdensome to calculate (versus a broadband wave field). Interestingly, apart from its role in setting up the wave amplitude equations, the irrotational part of the wave perturbations, $\boldsymbol {u}^{\varPhi \prime }_d + \boldsymbol {u}^{\varPhi \prime }_c$, need not be solved for this model. This feature is in common with the original Stokes vortex-force models (e.g. Craik & Leibovich Reference Craik and Leibovich1976) at leading order in $\mu$, where knowledge of $\boldsymbol {u}^{\zeta \prime }_c$ suffices for determining the WEC expression in $\boldsymbol {SVF}$.

For weaker waves (§ 11) and a fortiori for multiple primary waves (§ 12 and Appendix B), a wave-averaged theory for the current evolution is feasible in principle but probably laborious in practice because of the extra computational evaluations to obtain the final WEC equations. From this perspective, the monochromatic primary wave theory presented here, with its analytic wave-averaged expressions, is useful for clarifying wave-enhanced drag and vortex force effects and for supporting a rather straightforward computational implementation for current evolution. However, the generalization to multiple primary waves or a realistic broadband wave spectrum is much more daunting and might eventually require some form of wave-resolved modelling. If this were done, some ‘long-wave’ approximation to the free-surface dynamics could avoid the complexities of resolving spectrum broadening and wave breaking (versus its parametrization). Also, it would likely be advantageous to do a split-explicit time stepping (as commonly done for circulation models by separating the barotropic and baroclinic vertical modes (Shchepetkin & McWilliams Reference Shchepetkin and McWilliams2005)) with fast time stepping of the surface gravity wave component and subsequent numerical averaging on the slow time step of the currents. This type of hybrid wave–LES model would be a pioneering task. It could bypass the need for a formal wave-averaged theory, although it is likely that such theories as here would be useful for interpreting the computational results. Nevertheless, the community experience with Langmuir turbulence and depth-induced nearshore wave breaking indicates that important lessons can be learned even with only a monochromatic wave field representing the spectrum-peak component.

It seems worthwhile to proceed to a computation implementation of (10.6), e.g. for a kelp canopy, while also including the necessary non-conservative current effects (including some form of wave-enhanced bottom drag for the currents; see § 9), and, for now, without a formal extension of the theory to include buoyancy effects. The latter has been done previously for wave-averaged CEW and WEC theories (McWilliams et al. Reference McWilliams, Restrepo and Lane2004; WCI22, chapter 10), and it probably would contain no surprises with a canopy. The waves are not influenced by buoyancy because the surface-interface gravitational force is so much stronger than the force associated with interior stratification, and the main additional WEC influence is Stokes-drift buoyancy advection. An important implementation ingredient is sufficient experimental information about $\mathcal {D}(\boldsymbol {x})$, which often is not easily obtained. Nevertheless, many nearshore and estuarine sites have marine vegetation that is likely to be influential on the local currents and material transports.

Another necessary extension for nearshore realism in the wave-averaged theory and computational model built upon it is inclusion of both bottom drag and variable bottom depth. The former is briefly alluded to at the end of § 9 and is not inherently difficult. The latter is naturally a part of the resolved current dynamics, and small-scale roughness elements are a part of a Monin–Obukhov specification of the bottom drag coefficient. However, further thought is needed for how to efficiently include finite-height, finite-width (versus horizontally slowly varying, where a ray theory would be appropriate) depth variations in the wave component of the flow. Perhaps this, too, would be an argument for developing a long-wave free-surface LES model as mentioned in the second paragraph of this section. At the present time, oceanic LES examples for flow over such topography are rare, even in the absence of surface wave effects – but see Calhoun & Street (Reference Calhoun and Street2001), Calhoun, Street & Kosoff (Reference Calhoun, Street and Kosoff2001) and Skyllingstad & Wijesekera (Reference Skyllingstad and Wijesekera2004).

Beyond the derivation of wave-averaged equations, a measured programme of idealized computational experiments would be the right next step. Among the references in § 1 are various experimental results about wave damping in the presence of currents, and a computational implementation of this theory could take advantage of their various determinations of $\mathcal {D} = \alpha (\boldsymbol {x}) C_D$. In an oceanically targeted LES model with finite depth and bottom-rooted vegetation such as kelp, simple geometries of different plant profiles and different horizontal patches could be explored with various surface stress and buoyancy flux and primary waves to assess how the combined wave–current–canopy interactions manifest. Our recent papers partly set this framework (Yan et al. Reference Yan, McWilliams and Chameki2021; Yan, McWilliams & Chamecki Reference Yan, McWilliams and Chamecki2022), but as yet without the wave-enhanced drag effects. Beyond these idealized studies, some field conditions could be targeted for experimental testing. For various motivations, it would be timely to focus on kelp examples, both the heterogeneous natural patches found in shallow water outside the surf zone and the present and future suspended-platform aqua-farm configurations.

It is a delicate balancing act between do-able theories that enable informative complex-flow simulations and their underlying approximations that are only marginally defensible compared to nature's true complexities. Formulating, solving and experimentally testing feasible, relevant problems involving waves, currents, topography and canopies still have a long way to go.