1. Introduction
Mass transport by waves has been widely studied due to its crucial role in transporting material in the ocean, such as coastal sediment (Henderson & Allen Reference Henderson and Allen2004; Hsu, Elgar & Guza Reference Hsu, Elgar and Guza2006; Ruessink et al. Reference Ruessink, Michallet, Abreu, Sancho, Van der A, Van der Werf and Silva2011; Miles Reference Miles2013; Li et al. Reference Li, Zhang, Zhao, Qi, Cai and Zheng2024), spilt oil (Christensen & Terrile Reference Christensen and Terrile2009; Drivdal, Broström & Christensen Reference Drivdal, Broström and Christensen2014; Laxague et al. Reference Laxague, Özgökmen, Haus, Novelli, Shcherbina, Sutherland, Guigand, Lund, Mehta, Alday and Molemaker2018) and marine plastics (DiBenedetto et al. Reference Di benedetto, Ouellette and Koseff2018, Reference Di benedetto, Clark and Pujara2022; Calvert et al. Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019, Reference Calvert, McAllister, Whittaker, Raby, Borthwick and van den Bremer2021; van Sebille et al. Reference van Sebille, Aliani, Law, Maximenko, Alsina, Bagaev, Bergmann, Chapron, Chubarenko and Cózar2020; Núñez et al. Reference Núñez, Romano, García-Alba, Besio and Medina2023; Liao et al. Reference Liao, Zou, Vinh, Pan and Kaiser2024).
For non-breaking waves, mass transport occurs due to both the rotational wave motion near the surface and bottom boundary layers and the irrotational motion in the core region outside the boundary layers (Longuet–Higgins Reference Longuet-Higgins1953). Inside the boundary layer the horizontal and vertical velocities are no longer in quadrature, their correlation gives rise to wave radiation stress which generates a Eulerian steady streaming inside the near-bed wave boundary layer and a steady vorticity in the free-surface boundary layer (Russell & Osorio Reference Russell and Osorio1957; Longuet-Higgins Reference Longuet-Higgins1960; Zou Reference Zou2002; Zou & Hay Reference Zou and Hay2003; Zou, Bowen & Hay Reference Zou, Bowen and Hay2006). Both the steady streaming and the steady vorticity outside the boundary layer drive the mass transport in the core fluid region (Iskandarani & Liu Reference Iskandarani and Liu1991; Vittori & Blondeaux Reference Vittori and Blondeaux1996; Blondeaux, Brocchini & Vittori Reference Blondeaux, Brocchini and Vittori2002; Ng Reference Ng2004a ,Reference Vittori and Blondeaux b ).
Fluid particles in the irrotational core region experience a net forward drift known as the Stokes drift due to the non-closed trajectory over each wave cycle (Stokes Reference Stokes1847). For a monochromatic progressive wave with amplitude
$a$
, angular frequency
$\omega$
and wavenumber
$k$
at constant depth
$h$
, the first-order surface elevation and velocity field are given by
$\eta =a{\cos}(kx-\omega t)$
and
\begin{align}\left\{\begin{array}{@{}l} u = a\omega \dfrac{{\cosh \left[k\left( {z + h} \right)\right]}}{{\sinh kh}}\cos \left( {kx - \omega t} \right)\!, \\[12pt] w = a\omega \dfrac{{\sinh \left[k\left( {z + h} \right)\right]}}{{\sinh kh}}\sin \left( {kx - \omega t} \right)\!, \\ \end{array}\right. \end{align}
where z is positive upward with
$z=0$
at the still water level. The first-order trajectory of a particle can thus be obtained by integrating (1.1) with time as illustrated in figure 1. The Stokes drift velocity can be derived from the net drift after a wave cycle (Ursell Reference Ursell1953)
Figure 1. Schematic of the Eulerian and Lagrangian descriptions of mass transport induced by a progressive wave.
\begin{equation}U_{s}=\overline{\int u\textrm{d}t\frac{\partial u}{\partial x}+\int w\textrm{d}t\frac{\partial u}{\partial z}}=\frac{\left(ak\right)^{2}c}{2}\frac{\cosh \left[2k(z+h)\right]}{\sinh^{2}\!{kh}},\end{equation}
where
$c = \omega /k$
is the wave propagation speed. Integrating
$\rho U_{s}$
across the water column yields the Stokes transport
\begin{equation}\int _{-h}^{0}\rho U_{s}\textrm{d}z=\frac{\rho ga^{2}}{2c}=\frac{E}{c}\end{equation}
where
$E=\rho ga^{2}/2$
is the wave energy, ρ is the fluid density, g is the gravitational acceleration.
Besides examining the trajectory of a Lagrangian particle, Stokes transport can also be obtained in an Eulerian system. In the first-order problem, the velocity in (1.1) is always defined from the bottom
$z=-h$
up to
$z=0$
even during
$\eta \lt 0$
when physically there is no velocity between
$z=\eta$
and
$z=0$
. The velocity in the wave trough area given by (1.1) is therefore virtual and represents the analytical extension of the actual velocity from
$z=\eta$
to
$z=0$
and, by the same argument, from
$z=0$
to
$z=\eta$
in the crest area during
$\eta \gt 0$
. Effectively, the wave crest where
$\eta \gt 0$
carries ‘positive’ mass forward while the wave trough where
$\eta \lt 0$
carries ‘negative’ mass backward (Longuet-Higgins Reference Longuet-Higgins1969). Therefore, the residual net forward mass transport is approximately
$\rho \overline{u| _{z=0}\eta }$
, which leads to the same result of
$E/c$
. Physically, the horizontal velocity observed by a sensor at a fixed height
$z=z_{o}$
between wave crest and trough (
$-a\leq z_{o}\leq a$
) is intermittently zero when
$z_{o}\gt \eta$
as
\begin{align}{\left. u \right|_{z = {z_o}}} = \left\{ {\begin{array}{c@{\quad}l} {a\omega \dfrac{{\cosh \left[k({z_o} + h)\right]}}{{\sinh kh}}\cos \left( {kx - \omega t} \right)} & {{z_o} \le \eta ,} \\[6pt] 0 & {{z_o} \gt \eta ,} \\ \end{array} } \right.\end{align}
which yields a non-zero forward mean velocity
$\overline{u}(z_{o})$
between the crest and trough, as illustrated in figure 1, and hence a forward net mass transport equal to
$\int _{-a}^{a}\rho\overline{u}\textrm{d}z_{o}$
. The first-order particle orbits and thus mass transport can be significantly modified by higher-order nonlinear harmonics. Using a perturbation method, Chen, Hsu & Chen (Reference Chen, Hsu and Chen2010) derived the Lagrangian solution of particle trajectory for monochromatic progressive waves over a flat bottom up to fifth order. Chen et al. (Reference Chen, Li, Hsu and Ng2012) extended the flat bottom solution to a sloping bottom case up to third order. Recently, Li et al. (Reference Li, Chen, Zou and Hsu2023) obtained the particle trajectory solution for partial reflected waves over a sloping bottom. However, these studies adopted monochromatic waves as the leading-order wave, therefore, do not account for the effect of second-order subharmonic waves generated through difference interaction between first-order waves (primary waves) of different frequencies. As a result, subharmonic wave-induced mass transport has not been studied to authors’ knowledge.
The subharmonic wave is often regarded as an infragravity wave as its period falls into the infragravity band of 25∼250 s (Bertin et al. Reference Bertin, de Bakker, van Dongeren, Coco, André, Ardhuin, Bonneton, Bouchette, Castelle and Crawford2018). They are generated by energy transfer from short waves propagating into nearshore shallow water (Norheim, Herbers & Elgar Reference Norheim, Herbers and Elgar1998; Sheremet et al. Reference Sheremet, Guza, Elgar and Herbers2002; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015) driven by the wave group-modulated radiation stress gradient. When short waves break in the surf zone, infragravity wave height may grow up to over 1 m towards the shore from a few centimetres in the offshore. For example, the well-known surf beats observed by Munk (Reference Munk1949) and Tucker (Reference Tucker1950) are one kind of infragravity wave in the surf zone. Infragravity waves are related to various nearshore and coastal processes such as rip currents, run-up/overtopping, as well as beach and dune erosion (Russell Reference Russell1993; de Vries et al. Reference van Thiel de Vries, van Gent, Walstra and Reniers2008; Castelle et al. Reference Castelle, Scott, Brander and McCarroll2016; Bertin et al. Reference Bertin, de Bakker, van Dongeren, Coco, André, Ardhuin, Bonneton, Bouchette, Castelle and Crawford2018). When the depth is shallow relative to the wave group length but not the short wavelength, explicit theoretical solutions of the group-forced subharmonic were obtained for a flat bottom (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962; Davey & Stewartson Reference Davey and Stewartson1974; Mei et al. Reference Mei, Michael, Yue and Mei2005) and a slowly varying bathymetry (Janssen, Battjes & van Dongeren Reference Janssen, Battjes and van Dongeren2003; Zou Reference Zou2011). Calvert et al. (Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019) recently extended the solution to arbitrary depth relative to the wave group length. However, those solutions diverge when the subharmonic response to wave group forcing approaches resonance as the depth becomes shallow relative to the short wavelength. In this case, solutions in integral form were derived for a plane sloping bed (Symonds, Huntley & Bowen Reference Symonds, Huntley and Bowen1982; van Leeuwen Reference van Leeuwen1992; Schäffer Reference Schäffer1993) and for a slowly varying bathymetry (Liao et al. Reference Liao, Li, Liu and Zou2021, Reference Liao, Zou, Liu, Contardo and Li2023).
However, despite improved understanding of generation and evolution of subharmonic waves in the past years, not until recently did its effect on particle movements start to be recognised. The low-frequency infragravity wave motion causes large excursion of water particles because the particle excursion is proportional to wave period (Bondehagen, Kalisch & Roeber Reference Bondehagen, Kalisch and Roeber2024). Calvert et al. (Reference Calvert, Whittaker, Raby, Taylor, Borthwick and van den Bremer2019) showed how the set down in a transient wave group affects the particle drift throughout the water column. Bjørnestad et al.’s (Reference Bjørnestad, Buckley, Kalisch, Streßer, Horstmann, Frøysa, Ige, Cysewski and Carrasco‐Alvarez2021) Lagrangian measurements in the surf zone indicate that infragravity waves may trigger short wave breaking by enhancing the peak forward particle velocity at the surface. Through theory and field observations, Herbers & Janssen (Reference Herbers and Janssen2016) observed energetic infragravity fluctuations in the drift velocity of surface buoy in deep water that are several orders of magnitude larger and out of phase with the Eulerian infragravity motions. McAllister & van den Bremer (Reference McAllister and van den Bremer2019) found the Lagrangian measurement of surface drifter is significant in the group-bounded subharmonic band, which affects the ‘apparent’ steepness derived from time histories of drifter motion and poses a potential issue for wave buoys to measure acceleration. Flores, Williams & Horner-Devine (Reference Flores, Williams and Horner-Devine2022) show the first measurements of infragravity waves’ modulation of salinity near a river plume.
Although little is known about subharmonic wave-induced mass transport, the contribution of subharmonic waves to sediment transport are well recognised. In intermediate water, group-forced subharmonics are in anti-phase coupling with short wave groups, therefore, backward transport of high concentration sediment under large waves and subharmonic wave troughs exceeds forward transport of low concentration sediment under small waves and subharmonic wave crests, leading to net backwards sediment transport. As waves shoal and break in the nearshore, the phase relationship between subharmonic wave and wave groups evolves, which in turn alters the net sediment transport by group-forced subharmonic (Deigaard, Jakobsen & Fredsøe Reference Deigaard, Jakobsen and Fredsøe1999; Baldock et al. Reference Baldock, Manoonvoravong and Pham2010, Reference Baldock, Alsina, Caceres, Vicinanza, Contestabile, Power and Sanchez-Arcilla2011; de Bakker et al. Reference de Bakker, Brinkkemper, van der Steen, Tissier and Ruessink2016a ). On some occasions, infragravity waves may contribute to up to 60 % of onshore transport of suspended sediment in the inner surf zone over a gently sloping beach (de Bakker et al. Reference de Bakker, Brinkkemper, van der Steen, Tissier and Ruessink2016a ). Besides the sediment transport due to coupling between short waves and infragravity waves, field measurements indicated that onshore sediment transport may be induced by tidal-modulated infragravity waves in wave-dominated inlets or shallow and small estuaries (Bertin & Olabarrieta Reference Bertin and Olabarrieta2016; Williams & Stacey Reference Williams and Stacey2016; Bertin et al. Reference Bertin, Mendes, Martins, Fortunato and Lavaud2019; Mendes et al. Reference Mendes, Fortunato, Bertin, Martins, Lavaud, Nobre Silva, Pires-Silva, Coulombier and Pinto2020).
Currently, the literature on subharmonic wave-induced mass transport is scarce. In this work, theoretical solutions of the mass transport and the associated drift velocity induced by group-forced subharmonic waves in the shoaling zone are derived for the first time. The model is validated by the GLOBEX experiment (Ruessink et al. Reference Ruessink, Michallet, Bonneton, Mouaze and Lara2013). It was also validated against the mass transport velocity obtained by Lagrangian particle-tracking simulations using the flow field generated by the non-hydrostatic wave model SWASH (Zijlema, Stelling & Smit Reference Zijlema, Stelling and Smit2011). The theoretical solutions reveal two distinct kinematic mechanisms for subharmonic wave-induced mass transport: the forward transport in analogy to the Stokes transport of short waves, and the backward transport due to the anti-phase coupling between subharmonic waves and wave groups. The latter is proportional to the velocity skewness contributed by subharmonic waves. The theoretical model is combined with the SWASH model and the GLOBEX experiment (Ruessink et al. Reference Ruessink, Michallet, Bonneton, Mouaze and Lara2013) to examine the role of water depth, bed slope, wave period and bandwidth. The newly proposed subharmonic wave-induced mass transport provides the third mechanism for sediment transport in addition to the well-known boundary layer streaming (Nielsen Reference Nielsen2006) and suspension-advection mechanisms (Deigaard et al. Reference Deigaard, Jakobsen and Fredsøe1999). The subharmonic wave-induced mass transport is distributed throughout the water column hence is also expected to affect the transport of buoyant particles such as microplastics.
The paper is organised as follows. The theoretical solution of subharmonic wave surface elevation and velocity as well as the subharmonic wave-induced mass transport velocity and the bulk transport are derived and compared with the GLOBEX experiment in § 2. The set-up and postprocessing of the SWASH model and the Lagrangian particle-tracking simulation are described in § 3.1. Results of simulated particle trajectory, mass transport velocity and mass transport are demonstrated and compared with theoretical solutions in § 3.2. Discussions and conclusions are provided in § 4.
Figure 2. Definition sketch of variables and coordinate system for a short wave group propagating over an undulating bottom topography with mild slope. Here,
$\eta$
and
$\tilde{\eta }$
denote total surface elevation and the second-order subharmonic wave underneath the short wave group.
2. Theory of subharmonic wave-induced mass transport
2.1. Surface elevation and horizontal velocity of subharmonic wave
Consider shoaling unidirectional irregular waves propagating in the positive x-direction over an undulating seabed described by
$z=-h(x)$
. Let
$z=\eta (x,t)$
and
$u(x,z,t)$
be the surface elevation and horizontal velocity of the waves, respectively. The small parameter corresponding to amplitude is denoted by short wave steepness
$\varepsilon \equiv ka$
, where
$a$
and
$k$
are the characteristic wave amplitude and wavelength, respectively. The asymptotic expansion of
$\eta$
and
$u$
is given by
\begin{equation}\left.\begin{array}{c} u=u^{(1)}+u^{\left(2\right)}+\cdots \\[3pt] \eta =\eta ^{(1)}+\eta ^{\left(2\right)}+\cdots \end{array}\right\},\end{equation}
where the superscript
$(i)$
denotes the variables at the order of
$O(\varepsilon ^{i})$
.
Assume
$\eta ^{(1)}$
has frequency band
$[\omega _{p}(1-\alpha ),\omega _{p}(1+\alpha )]$
and can be written in the form of a discrete Fourier series
\begin{align}\left\{\begin{array}{@{}r@{\,}c@{\,}l}\displaystyle{\eta ^{(1)}}({x,t}) &=& \textit{Re} \bigg[ {\sum\limits_{n = {n_1}}^{{n_1} + N} {{{\hat a}_n}} ( x )\exp {\text{i}}{\theta _n}} \bigg], \\[14pt] \displaystyle{\theta _n} &=& \int_{{x_0}}^x {{k_n}} {\text{d}}x' - {\omega _n}t, \\ \end{array} \right.\end{align}
where
$\omega _{n}=n{\unicode[Arial]{x0394}} \omega$
is the radian frequency (
${\unicode[Arial]{x0394}} \omega$
= frequency resolution) for the nth component,
$k_{n}$
is the wavenumber satisfying the dispersion relationship
$\omega _{n}^{2}=gk_{n}\tanh k_{n}h$
,
$\hat{a}_{n}$
is the complex amplitude,
$n_{1}$
is the index of the lower limit frequency
$N$
is the number of frequency components and
$\textit{Re}[{\boldsymbol\cdot}]$
denotes the real part of the complex argument and will be omitted herein after for brevity.
Letting
$h_{x}\equiv \textrm{d}h/\textrm{d}x$
and
$h_{xx}\equiv \textrm{d}^{2}h/\textrm{d}x^{2}$
, we make following assumptions on the parameters:
\begin{align}\varepsilon \ll {\left( {kh} \right)^3},\,\,\,\,\,\alpha = O\left( \varepsilon \right)\!,\,\,\,\,\,\varepsilon kh \ll 1,\,\,\,\,\,\left| {{h_x}} \right| = O\big({{\varepsilon ^n}kh} \big)\left( {n \gt 1} \right)\!,\,\,\,\,\,\left| {{h_{xx}}} \right| = O\big({\varepsilon ^{2n}}{k^2}h\big),\end{align}
which is equivalent to the assumptions in Zou (Reference Zou2011).
Accordingly, the first-order horizontal velocity is given by
\begin{align} \left\{\begin{array}{c} u^{(1)}(x,z,t)=\displaystyle\sum _{n=n_{1}}^{n_{1}+N}\hat{u}_{n}(x)\dfrac{\cosh \left[k_{n}(z+h)\right]}{\cosh k_{n}h}\exp \textrm{i}\theta _{n},\\[15pt] \hat{u}_{n}(x)=\dfrac{\hat{a}_{n}(x)\omega _{n}}{\tanh k_{n}h}.\end{array}\right. \end{align}
Assuming the complex amplitude
$\hat{a}_{n}(x)$
only changes in magnitude following the conservative shoaling
\begin{equation}\frac{\hat{a}_{n}\left(x\right)}{\hat{a}_{n}\left(x_{0}\right)}=\sqrt{\frac{c_{g}\left(x_{0}\right)}{c_{g}\left(x\right)}}, c_{g}=\left.\frac{\partial \omega }{\partial k}\right| _{\omega ={\omega _{p}}}=\sqrt{\frac{g\tanh k_{p}h}{k_{p}}}\left(\frac{1}{2}+\frac{k_{p}h}{\sinh 2k_{p}h}\right)\!,\end{equation}
where
$c_{g}$
is the wave group speed corresponding to central frequency
$\omega _{p}$
, and (2.2) can be rewritten as
\begin{equation}\eta ^{(1)}(x,t)=A\left(x,t\right)\exp \left[\textrm{i}\left(\int _{x_{0}}^{x}k_{p}\textrm{d}x'\,-\omega _{p}t\right)\right]\!,\end{equation}
where the modulated complex amplitude
$A(x,t)$
is given by
\begin{align} \left\{\begin{array}{@{}c}A\left(x,t\right)=\displaystyle\sum _{n=n_{1}}^{n_{1}+N}\hat{a}_{n}\!\left(x\right)\exp \left[\textrm{i}\left(\int _{x_{0}}^{x}k_{n,p}{\textrm d}x'-\omega _{n,p}t\right)\right]\!,\\[18pt] \omega _{n,p}=\omega _{p}-\omega _{n}, \qquad k_{n,p}=k_{n}-k_{p} .\end{array} \right.\end{align}
Substituting (2.7) into the formula for wave energy
$E=\rho g| A| ^{2}/2$
yields the modulated wave energy
\begin{align}\left\{ {\begin{array}{c} {E\left( {x,t} \right) = \bar E\left( x \right) + \tilde E\left( {x,t} \right)} \\[4pt] \displaystyle{\bar E\left( x \right) = \dfrac{1}{2}\rho g\sum\limits_{n = {n_1}}^{{n_1} + N} | {{\hat a}_n}\left( x \right){|^2},} \\[10pt] \displaystyle{\tilde E\left( {x,t} \right) = \rho g\sum\limits_{m = {n_1}}^{{n_1} + N} {\sum\limits_{n = {n_1}}^{m - 1} {{{\hat a}_m}} \left( x \right)\hat a_n^*\left( x \right)\exp \left[ {{\rm i}\left( {\int _{x_{0}}^{x} {{k_{m,n}}} {\text{d}}x' - {\omega _{m,n}}t} \right)} \right]\!,} } \\[10pt] {{\omega _{m,n}} = {\omega _m} - {\omega _n},\,\,\,\,\,\,{k_{m,n}} = {k_m} - {k_n},} \\ \end{array} } \right.\end{align}
where
$\overline{E}$
and
$\tilde{E}$
denote the steady and unsteady components of the wave energy, respectively,
$n_{1}$
is the index of the lower limit frequency of
$\eta ^{(1)}$
and N is the number of frequency components of
$\eta ^{(1)}$
. The density
$\rho$
is omitted hereinafter for being constant and for brevity. Throughout this paper
$\widetilde{({\boldsymbol\cdot})}$
and
$\overline{({\boldsymbol\cdot} )}$
denote the low-frequency unsteady and steady components of a time-dependent variable, which are equivalent to averaging over the short wave and wave group phase, respectively.
Since
$\omega _{m,n}=\omega _{m}-\omega _{n}=\omega _{m-n}$
,
$\tilde{E}(x,t)$
in (2.8) can be rewritten by letting
$m-n\rightarrow l$
and
$\omega _{m,n}\rightarrow \omega _{l}$
and changing one of the double summations into the sum over new index
$l$
\begin{align} \tilde E({x,t}) & = g\sum\limits_{l = 1}^N {\sum\limits_{m = {n_1}}^{{n_1} + N} {{{\hat a}_m}}}(x)\hat a_{m - l}^*\left( x \right)\exp \left[ {{\text{i}}\left( {\int_{{x_0}}^x {{k_{m,m - l}}} {\text{d}}x' - {\omega _l}t} \right)} \right] \nonumber\\ & = g\sum\limits_{l = 1}^N {\left[ {\sum\limits_{m = {n_1}}^{{n_1} + N} {{{\hat a}_m}} \left( x \right)\hat a_{m - l}^*\left( x \right)} \right]} \exp \left[ {{\text{i}}\left( {\int_{{x_0}}^x {{k_l}} {\text{d}}x' - {\omega _l}t} \right) + O(\varepsilon)} \right]\!, \end{align}
where
\begin{align}{\omega _l} = O(\varepsilon {\omega _p}),{k_l} = \frac{{{\omega _l}}}{{{c_g}}} = O\left( {\varepsilon k} \right)\!,\end{align}
is the radian frequency and wavenumber of the wave group. The unsteady wave energy
$\tilde{E}$
can be expressed in terms of discrete Fourier series
\begin{equation}\tilde{E}(x,t)=\sum _{l=1}^{N}\hat{E}\left(x,\omega _{l}\right)\exp \left(-\textrm{i}\omega _{l}t\right)\!,\end{equation}
and comparing (2.11) with (2.9) yields
\begin{align}\left\{ {\begin{array}{c} \displaystyle{\hat E\left( {x,{\omega _l}} \right) = \hat E\left( {{x_0},{\omega _l}} \right)\frac{{{c_g}\left( {{x_0}} \right)}}{{{c_g}\left( x \right)}}\exp \left[ {{\text{i}}\int_{{x_0}}^x {{k_l}} {\text{d}}x' + O\left( \varepsilon \right)} \right]\!,} \\[14pt] \displaystyle{\hat E\left( {{x_0},{\omega _l}} \right) = g\left[ {\sum\limits_{m = {n_1}}^N {{{\hat a}_m}} \left( {{x_0}} \right)\hat a_{m - l}^*\left( {{x_0}} \right)} \right].} \\ \end{array} } \right.\end{align}
Let
$\tilde{\eta }$
,
$\tilde{u}$
and
$\tilde{Q}$
denote respectively the low-frequency surface elevation, horizontal velocity and volume flux, then, accurate to
$O(\varepsilon ^{2})$
,
$\tilde{Q}$
can be written (Longuet-Higgins Reference Longuet-Higgins1969; Monismith et al. Reference Monismith, Cowen, Nepf, Magnaudet and Thais2007)
\begin{align}\tilde Q = \tilde uh + \frac{{\tilde E}}{c} + O({\varepsilon ^4}),\end{align}
which has contributions from both the depth-uniform subharmonic orbital velocity
$\tilde{u}(x,t)$
and the modulated Stokes transport
$\tilde{E}/c$
. The linearised governing equations of
$\tilde{\eta }$
and
$\tilde{Q}$
are (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962; Zou Reference Zou2011; Liao et al. Reference Liao, Zou, Liu, Contardo and Li2023)
\begin{align}\left\{ {\begin{array}{c} \displaystyle{\frac{{\partial \tilde \eta }}{{\partial t}} + \frac{{\partial \tilde Q}}{{\partial x}} = 0,} \\[14pt] \displaystyle{\frac{{\partial \tilde Q}}{{\partial t}} + gh\frac{{\partial \tilde \eta }}{{\partial x}} = - \frac{{\partial \tilde S}}{{\partial x}},} \\ \end{array} } \right.\end{align}
where
$\tilde{S}(x,t)$
is the modulated wave radiation stress given by (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962)
\begin{equation}\tilde{S}=\tilde{E}\left(\frac{2c_{g}}{c}-\frac{1}{2}\right).\end{equation}
Similar to (2.11), we write
$\tilde{\eta }$
,
$\tilde{u}$
,
$\tilde{M}$
and
$\tilde{S}$
as discrete Fourier series
\begin{align}(\tilde \eta ,\tilde u,\tilde Q,\tilde S) = \sum\limits_{l = 1}^N {[\hat \eta (x,{\omega _l}),\hat u(x,{\omega _l}),\hat Q(x,{\omega _l}),\hat S(x,{\omega _l})]} \exp ( - {\text{i}}{\omega _l}t),\end{align}
where the argument
$\omega _{l}$
for complex amplitude will be omitted hereinafter for brevity. The complex amplitude
$\hat{Q}$
and
$\hat{S}$
can be derived by substituting (2.11)–(2.12) into (2.13) and (2.15), respectively,
\begin{align}\left\{\begin{array}{c} \displaystyle{\hat Q(x) = \hat uh + \frac{{\hat E}}{{c}} \cong {{\hat u}}h + \frac{{\hat E({x_0})}}{c}\frac{{{c_g}({x_0})}}{{{c_g}(x)}}\exp \left( {{\text{i}}\int_{{x_0}}^x {{k_l}} {\text{d}}x'} \right)\!,} \\[14pt] \displaystyle{\hat S(x) = \hat E\left( {\frac{{2{c_g}}}{c} - \frac{1}{2}} \right) \cong \left( {\frac{{2{c_g}}}{c} - \frac{1}{2}} \right)\hat E({x_0})\frac{{{c_g}({x_0})}}{{{c_g}(x)}}\exp \left( {{\text{i}}\int_{{x_0}}^x {{k_l}} {\text{d}}x'} \right).} \end{array}\right. \end{align}
Substituting (2.16) into (2.14) yields the governing equation for
$\hat{\eta }$
and
$\hat{Q}$
\begin{align} \left\{ \begin{array}{c} {\hat \eta (x) = \dfrac{1}{{{\text{i}}{\omega _l}}}\dfrac{{\partial \hat Q}}{{\partial x}},} \\[10pt] {\hat Q(x) = \dfrac{1}{{{\text{i}}{\omega _l}}}\dfrac{{\partial \hat S}}{{\partial x}} + \dfrac{{gh}}{{{\text{i}}{\omega _l}}}\dfrac{{\partial \hat \eta }}{{\partial x}}.} \end{array}\right. \end{align}
We note that the derivation from (2.1) to (2.18) essentially repeats many existing theoretical studies of the group-forced subharmonic over a mild sloping bottom (Liu Reference Liu1989; van Leeuwen Reference van Leeuwen1992; Schäffer Reference Schäffer1993; Janssen et al. Reference Janssen, Battjes and van Dongeren2003; Battjes et al. Reference Battjes, Bakkenes, Janssen and van Dongeren2004; Zou Reference Zou2011; Contardo et al. Reference Contardo, Lowe, Hansen, Rijnsdorp, Dufois and Symonds2021; Liao et al. Reference Liao, Li, Liu and Zou2021).
Over a mild bottom slope relative to the wave group length scale with the following parameter regime:
\begin{align}\beta = \frac{{{h_x}}}{{{k_l}h}} = O({\varepsilon ^{n - 1}}) \ll 1,\frac{{{h_{xx}}}}{{k_l^2h}} = O({\varepsilon ^{2n - 2}}) = O({\beta ^2}),\end{align}
the solution of
$\hat{\eta }$
accurate to
$O(\beta )$
under the off-resonance condition
\begin{equation}\frac{\beta }{\mu }=O\left(\beta \right)\!, \mu =1-\frac{c_{g}^{2}}{gh},\end{equation}
in intermediate depth is taken from the existing literature (Janssen et al. Reference Janssen, Battjes and van Dongeren2003; Zou Reference Zou2011; Liao et al. Reference Liao, Li, Liu and Zou2021)
\begin{align}\hat \eta \left( x \right) & = - \frac{{\hat S}}{{gh\mu }}\left\{1 + \frac{{{\text{i}}\beta }}{\mu }\left[ {\left( {1 - \mu } \right)\frac{{\text{d}}}{{{\text{d}}\ln h}}\ln \left[ {{{\left( {\frac{{2{c_g}}}{c} - \frac{1}{2}} \right)}^2}c_g^{ - 3}} \right] - \frac{{\text{d}}}{{{\text{d}}\ln h}}\ln ( {h{\mu ^2}} )} \right]\right.\nonumber\\[5pt]& \left.\vphantom{\left[ {{{\left( {\frac{{2{c_g}}}{c} - \frac{1}{2}} \right)}^2}c_g^{ - 3}} \right]}\quad + O( {{\beta ^2}}) \right\},\end{align}
where
$\hat{S}(x)$
is given in (2.17). The parameter μ describes the extent of deviation from resonance (Janssen et al. Reference Janssen, Battjes and van Dongeren2003) and decreases from unity in deep water to zero at shallow water as
$kh\rightarrow 0$
, and
$\beta$
measures the bed slope relative to the wave group length scale (Battjes et al. Reference Battjes, Bakkenes, Janssen and van Dongeren2004). As
$kh\rightarrow 0$
, the
$O(\beta )$
terms in (2.21) increase in magnitude as
\begin{equation}\frac{\beta }{\mu }=O\left(\frac{\varepsilon ^{n-1}}{\left(kh\right)^{2}}\right)\!,\end{equation}
hence approaching
$O(1)$
when
$kh$
approaches
$\varepsilon^{({n-1)}/{2}}$
(the so-called near-resonance condition). In this case, the off-resonant solution (2.21) starts to diverge and is no longer valid.
Without assuming the depth and resonance regime, Liao et al. (Reference Liao, Zou, Liu, Contardo and Li2023) obtained the leading-order solution of
$\hat{\eta }$
as ((B11) therein)
\begin{align} \begin{cases} \displaystyle\hat{\eta }\left(x\right)=\hat{\eta }_{f}^{+}\left(x_{0}\right)\hat{\eta }_{h}^{+}\left(x,x_{0}\right)+\hat{\eta }_{f}^{-}\left(x_{0}\right)\hat{\eta }_{h}^{-}\left(x,x_{0}\right)\\[10pt] \qquad \displaystyle+\int _{x_{0}}^{x}\,\dfrac{k_{l}^{2}}{2\textrm{i}k_{f}}\dfrac{\hat{S}}{gh}\left[\hat{\eta }_{h}^{+}(x,x')-\hat{\eta }_{h}^{-}(x,x')\right]\textrm{d}x',\\[15pt] \displaystyle \hat{\eta }_{h}^{\pm }\left(x_{2},x_{1}\right)=\left[\dfrac{h\left(x_{2}\right)}{h\left(x_{1}\right)}\right]^{-0.25}\exp \left(\textrm{i}\int _{x_{1}}^{x_{2}}\,\pm k_{f}\textrm{d}x^{\prime}\right)\!, \end{cases} \end{align}
with an error factor
$1+O(\beta )$
, where
$k_{f}=\omega _{l}/\sqrt{gh}$
is the wavenumber of shallow-water waves at wave group frequency
$\omega _{l} ; \hat{\eta }_{f}^{+}$
and
$\hat{\eta }_{f}^{-}$
are the complex amplitudes of forward and backward propagating subharmonic waves of frequency
$\omega _{l}$
, respectively;
$x_{0}$
is the boundary for integration where
$\hat{\eta }_{f}^{\pm }(x_{0})$
is known; function
$\hat{\eta }_{h}^{\pm }(x_{2},x_{1})$
describes the change in complex amplitude as a shallow-water wave propagates from
$x_{1}$
to
$x_{2}$
. Liao et al. (Reference Liao, Zou, Liu, Contardo and Li2023) derived (2.23) using a novel Green’s function. Unlike previous solutions, the Green’s function-based solution treats all subharmonic waves as free harmonics continuously radiated at shallow-water wave speed away from each point source in the forcing field of the group modulated wave radiation stress gradient, phase locked with the wave group. The solution describes the generation, propagation and interference of freely propagating subharmonic waves without distinguishing between free and bound waves. It was shown that the superposition of these free harmonics results in what are termed bound subharmonics, due to the modulation of each free harmonic’s phase by the wave group. Additionally, it was demonstrated that, during shoaling, the radiated subharmonic waves propagating with the wave groups accumulate energy along the travel distance, leading to their rapid growth and the phase lag with respect to the wave groups.
It was also shown in Liao et al. (Reference Liao, Zou, Liu, Contardo and Li2023) that (2.21) can be decomposed into
$\hat{\eta }=\hat{\eta }^{+}+\hat{\eta }^{-}$
as follows:
\begin{align}\hat{\eta }^{\pm }\left(x\right)&=\mp \frac{\hat{S}}{gh\mu }\frac{k_{l}\pm k_{f}}{2k_{f}}\left[1-\textrm{i}\beta \frac{\textrm{d}}{\textrm{dln}h}\ln \left[\left(\frac{2c_{g}}{c}-\frac{1}{2}\right)^{2}c_{g}^{-3}\right]\right .\notag\\[5pt]& \left . -\beta \frac{k_{l}\pm k_{f}}{\textrm{i}k_{l}\mu }\frac{\textrm{d}}{\textrm{dln}h}\ln \left[\left(\frac{2c_{g}}{c}-\frac{1}{2}\right)\frac{1}{c_{g}h^{0.25}\mu }\right]-\frac{\beta }{\textrm{i}k_{l}\mu }\frac{\textrm{d}}{\textrm{dln}h}\left(k_{l}\pm k_{f}\right)+O(\beta ^{2})\vphantom{\left[\left(\frac{2c_{g}}{c}-\frac{1}{2}\right)^{2}c_{g}^{-3}\right]}\right]\!,\end{align}
where the subscripts + and – denote the forward and backward propagating components, respectively. The parameter
$\hat{\eta }_{l}^{\pm }$
in (2.24) will be used as the boundary values
$\hat{\eta }_{f}^{\pm }(x_{0})$
at
$x=x_{0}$
in (2.23). It follows that the model only considers the backward propagating component generated by the forcing of radiation stress gradient and thus does not account for subharmonic waves reflected from the bathymetry or the shoreline. This can be seen from the fact that the decomposition of
$\hat{\eta }=\hat{\eta }^{+}+\hat{\eta }^{-}$
in (2.24) is always valid even for a flat bottom with
$\beta =0$
, where the sum
$\hat{\eta }^{+}+\hat{\eta }^{-}$
recovers the equilibrium solution by Longuet-Higgins & Stewart (Reference Longuet-Higgins and Stewart1962).
As shown later, the analytical consideration of subharmonic wave-induced mass transport requires its velocity
$\tilde{u}$
or
$\hat{u}$
, which was not obtained in Liao et al. (Reference Liao, Zou, Liu, Contardo and Li2023). To derive that, substituting (2.21) into the second equation in (2.18) and retaining terms to order
$O(\beta )$
yields the solution for
$\hat{Q}$
under the off-resonance condition in intermediate depth
\begin{align}\hat{Q}\left(x\right)&=\frac{-\hat{S}c_{g}}{gh\mu }\left\{1+\frac{\textrm{i}\beta }{\mu }\left(\frac{2-3\mu -\mu ^{2}}{1-\mu }\frac{\textrm{d}}{\textrm{dln}h}\ln \left(\frac{2c_{g}}{c}-\frac{1}{2}\right)-\frac{3-4\mu -\mu ^{2}}{1-\mu }\frac{\textrm{d}}{\textrm{dln}h}\ln c_{g}\right . \right .\notag\\[8pt]&\quad \left . \left . -1+\frac{\mu -2}{1-\mu }\frac{\textrm{d}}{\textrm{dln}h}\ln \mu \right)+O(\beta ^{2})\right\}.\end{align}
Similarly, substituting (2.23) into the second equation in (2.18) yields the leading-order solution in shallow water
\begin{align}\hat{Q}\left(x\right)&=\frac{\hat{S}}{c_{g}}+\sqrt{gh}\left[\hat{\eta }_{f}^{+}\left(x_{0}\right)\hat{\eta }_{h}^{+}\left(x,x_{0}\right)-\hat{\eta }_{f}^{-}\left(x_{0}\right)\hat{\eta }_{h}^{-}\left(x,x_{0}\right)+\int _{x_{0}}^{x}\,\frac{k_{l}^{2}}{2\textrm{i}k_{f}}\frac{\hat{S}}{gh}\left[\hat{\eta }_{h}^{+}\left(x,y\right)\right . \right .\notag\\&\quad \left . \left .+\;\hat{\eta }_{h}^{-}\left(x,y\right)\right]\textrm{d}y\vphantom{\frac{k_{l}^{2}}{2\textrm{i}k_{f}}\frac{\hat{S}}{gh}}\right].\end{align}
As reflected in (2.18),
$\hat{Q}$
is the combined effect of the direct response of the subharmonic volume flux to the forcing of radiation stress gradient
$\partial \hat{S}/\partial x$
and the volume flux carried by subharmonic waves (the
$\partial \hat{\eta }/\partial x$
term). At leading order, the former leads to the first term
$\hat{S}/c_{g}$
on the right-hand side of (2.26), while the latter is essentially multiplying the forward and backward propagating subharmonic waves in (2.23) with their corresponding phase speeds
$\sqrt{gh}$
and
$-\sqrt{gh}$
.
Combining (2.17) and (2.26) yields the complex amplitude of subharmonic orbital velocity
$\hat{u}$
\begin{align}\nonumber\\[-18pt]\hat{u}\left(x\right)=\frac{\hat{Q}-\hat{E}/c}{h}.\end{align}
As a special case, consider a bichromatic wave system over a flat bed at a water depth deep for the primary waves but shallow for the wave group. Substituting
$\beta =0$
,
$\hat{E}=g\hat{a}_{n}\hat{a}_{n-l}^{*}\exp [\textrm{i}(\int _{x_{0}}^{x}k_{l}\textrm{d}x'-\omega _{l}t)]$
(2.12),
$\hat{S}=\hat{E}(({2c_{g}}/{c})-({1}/{2}))$
and the deep-water relationship
$c_{g}=1/2c=g/(2\omega _{p})$
into (2.25), and then into (2.27) yields
\begin{equation}\hat{u}\left(x\right)=-\frac{\omega _{p}}{2h}\frac{\hat{a}_{n}\hat{a}_{n-l}^{*}}{1-1/\left(4k_{p}h\right)}\exp \left(\textrm{i}\int _{x_{0}}^{x}k_{l}\textrm{d}x^{\prime}\right)\!,\end{equation}
which is equivalent to (2.3) in van den Bremer et al. (Reference van den Bremer, Whittaker, Calvert, Raby and Taylor2019).
2.2. Eulerian transport
2.2.1. Generalised theory
In an Eulerian coordinate system, the wave-induced mass transport in homogeneous media with unit
$\rho$
is defined as
\begin{equation}M(x,t)=\int _{-h}^{\eta }u\textrm{d}z=\int _{-h}^{0}u\textrm{d}z+\int _{0}^{\eta }u\textrm{d}z.\end{equation}
The horizontal velocity
$u$
can be expanded vertically around
$z=0$
as
\begin{equation}u\left(x,t\right)=\left.u\right| _{z=0}+z\left.\frac{\partial u}{\partial z}\right| _{z=0}+\frac{1}{2}z^{2}\left.\frac{\partial ^{2}u}{\partial z^{2}}\right| _{z=0}+\ldots .\end{equation}
Substituting (2.30) into (2.29) yields
\begin{equation}M(x,t)=\int _{-h}^{0}u\textrm{d}z+u_{0}\eta +\frac{1}{2}\left.\frac{\partial u}{\partial z}\right| _{z=0}\eta ^{2}+\frac{1}{6}\left.\frac{\partial ^{2}u}{\partial z^{2}}\right| _{z=0}\eta ^{3}+\ldots,\end{equation}
where the subscript 0 indicates evaluation at
$z=0$
Substituting the asymptotic expansion of
$u$
and
$\eta$
in (2.1) into (2.31) yields
$M=M^{(1)}+M^{(2)}+M^{(3)}+M^{(4)}\ldots$
, where the superscript
$(i)$
denotes the variable at the order of
$O(\varepsilon ^{i})$
, i.e.
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad M^{(1)}\left(x,t\right)=\int _{-h}^{0}u^{(1)}\textrm{d}z,\\[-10pt]\nonumber \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad M^{\left(2\right)}\left(x,t\right)=\int _{-h}^{0}u^{\left(2\right)}\textrm{d}z+u_{0}^{(1)}\eta ^{(1)},\\[-10pt]\nonumber \end{align}
$$\begin{align}&\qquad M^{\left(3\right)}\left(x,t\right)=\int _{-h}^{0}u^{\left(3\right)}\textrm{d}z+u_{0}^{(1)}\eta ^{\left(2\right)}+u_{0}^{\left(2\right)}\eta ^{(1)}+\frac{1}{2}\left.\frac{\partial u^{(1)}}{\partial z}\right| _{z=0}\big[\eta ^{\left(1\right)}\big]^{2},\end {align}$$
\begin{align}& M^{\left(4\right)}\left(x,t\right)=\int _{-h}^{0}u^{\left(4\right)}\textrm{d}z+u_{0}^{(1)}\eta ^{\left(3\right)}+u_{0}^{\left(2\right)}\eta ^{\left(2\right)}+u_{0}^{\left(3\right)}\eta ^{(1)}+\frac{1}{2}\left.\frac{\partial u^{\left(2\right)}}{\partial z}\right| _{z=0}\big[\eta ^{\left(1\right)}\big]^{2}\notag\\&\qquad\qquad\quad+\left.\frac{\partial u^{(1)}}{\partial z}\right| _{z=0}\eta ^{\left(1\right)}\eta ^{\left(2\right)}+\frac{1}{6}\left.\frac{\partial ^{2}u^{(1)}}{\partial z^{2}}\right| _{z=0}\big[\eta ^{(1)}\big]^{3}.\end{align}
Our goal is to obtain the time-averaged mass transport
$\bar{M}$
at each order and seek terms involving the second-order subharmonic velocity
$\tilde{u}$
to assess its contribution. Physically, an Eulerian velocity of order
$O(\varepsilon ^{3})$
such as
$u^{(3)}$
and
$u^{(4)}$
will depend on
$\tilde{u}$
and so will the mass transport terms involving them. This dependence of mass transport on
$\tilde{u}$
, however, is dynamic rather than kinematic. The latter may be defined as the difference in the drift velocity of a Lagrangian particle after subtracting
$\tilde{u}$
from the background flow, whereby the former requires the derivation of the change in flow field at higher orders due to subtracting
$\tilde{u}$
using the momentum equations and boundary conditions. The present study focuses on the kinematic contribution of subharmonic waves to the mass transport, thereby
$u^{(3)}$
and
$u^{(4)}$
are neglected hereinafter. The subharmonic-induced mass transport solution will be similar to the kinematic Stokes drift velocity (1.2) without the accompanying second-order subharmonic return flow.
To proceed, we will first consider a more general wave velocity field without the assumptions
$\alpha =O(\varepsilon )\ll 1, k_{l}h=O(\alpha kh)=O(\varepsilon kh)\ll 1$
to derive the result of subharmonic-induced mass transport and then simplify it to the different parameter regimes. The first-order horizontal velocity
$u^{(1)}$
takes the same expression as in (2.4). The following general form is now adopted for subharmonic velocity
$\tilde{u}$
contributing to
$u^{(2)}$
as suggested by the solution of subharmonic velocity obtained in Longuet-Higgins & Stewart (Reference Longuet-Higgins and Stewart1962) ((3.10) therein):
\begin{align} \left\{\begin{array}{@{}c} u^{\left(2\right)}=\tilde{u}+\left(\text{superharmonic terms}\right)\!,\\[10pt] \displaystyle\tilde{u}(x,z,t)=\sum\limits_{n,m}\widehat{U}_{m,n}\dfrac{\cosh \left[k_{m,n}(z+h)\right]}{\cosh k_{m,n}h}\exp \textrm{i}\theta _{m,n}, \end{array}\right. \end{align}
with
\begin{align} \begin{cases} \theta _{m,n}=k_{m,n}x-\omega _{m,n}t,\\ k_{m,n}=k_{m}-k_{n},\\ \omega _{m,n}=\omega _{m}-\omega _{n} .\end{cases} \end{align}
Both indices
$n$
and
$m$
were summed over the same frequency domain of
$u^{(1)}$
.
Substituting (2.4) into (2.32
a) yields
$M^{(1)}$
\begin{equation}M^{(1)}=\sum _{m}\dfrac{\hat{u}_{m}}{k_{m}}\tanh k_{m}h\exp \textrm{i}\theta _{m}.\end{equation}
After obtaining
$M^{(i)}$
,
$\eta ^{(i)}$
is needed for deriving
$M^{(j)}$
(
$j\gt i)$
and can be obtained from the continuity equation as
\begin{equation}\eta ^{\left(i\right)}=-\int \frac{\partial M^{\left(i\right)}}{\partial x}\textrm{d}t.\end{equation}
Substituting (2.35) into (2.36) yields
$\eta ^{(1)}$
\begin{equation}\eta ^{(1)}=\sum _{m}\frac{\hat{u}_{m}}{\omega _{m}}\tanh k_{m}h\exp \textrm{i}\theta _{m},\end{equation}
which together with
$u^{(2)}$
in (2.33) and
$u^{(1)}$
in (2.4), is further substituted into (2.31b
) to yield
$M^{(2)}$
. Repeating this process, the terms in
$M^{(4)}$
involving subharmonic velocity
$\tilde{u}$
can be identified
\begin{align}&M^{\left(4\right)}\left(x,t\right)\notag\\&=\underset{u_{0}^{(1)}\eta ^{\left(3\right)}}{\underbrace{u_{0}^{(1)}\int \frac{\partial }{\partial x}\left[u_{0}^{(1)}\int \!\left(\int _{-h}^{0}\frac{\partial \tilde{u}}{\partial x}\textrm{d}z\right)\textrm{dt}-\tilde{u}_{0}\eta ^{(1)}\right]\textrm{d}t}} \underset{u_{0}^{\left(2\right)}\eta ^{\left(2\right)}}{\underbrace{-u_{0}^{\left(2\right)}\!\int \!\frac{\partial }{\partial x}\left(\int _{-h}^{0}\tilde{u}\textrm{d}z+u_{0}^{(1)}\eta ^{(1)}\!\right)\textrm{d}t}}\notag\\&\quad +\frac{1}{2}\left.\frac{\partial \tilde{u}}{\partial z}\right| _{z=0}\big[\eta ^{\left(1\right)}\big]^{2}-\left.\frac{\partial u^{(1)}}{\partial z}\right| _{z=0}\eta ^{\left(1\right)}\int \frac{\partial }{\partial x}\left(\int _{-h}^{0}\tilde{u}\textrm{d}z\right)\textrm{d}t+(\text{terms without}\ \tilde{u}).\end{align}
Notice that the first term on the right-hand side of (2.38) is part of the term
$u_{0}^{(1)}\eta ^{(3)}$
, accounting for the interaction of
$u_{0}^{(1)}$
with
$\tilde{u}$
through the component of
$\eta ^{(3)}$
. The phase average of this term turns out to be non-zero and therefore contributes to the mass transport.
Invoking (2.4) for
$u^{(1)}$
and (2.33) for
$\tilde{u}$
and (2.37) for
$\eta ^{(1)}$
, the expression of (2.38) for
$M^{(4)}$
can be worked out. Maintaining only the steady terms involving the subharmonic velocity
$\widehat{U}_{m,n}$
yields
\begin{equation}\bar{M}^{\left(4\right)}\left(x\right)=\sum _{m,n}\underset{\text{Term}\;\textrm{I}}{\underbrace{\frac{1}{2}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{\tanh k_{m,n}h}{\omega _{m,n}}}}+\underset{\text{Term II}}{\underbrace{\frac{1}{4}D_{m,n}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}h}},\end{equation}
with
\begin{align}D_{m,n}&=\left(\frac{k_{m,n}}{\omega _{m,n}}+\frac{k_{n}}{\omega _{n}}\right)\frac{k_{m}}{\omega _{m}}\frac{\tanh k_{m}h}{k_{m}h}+\left(\frac{k_{m,n}}{\omega _{m,n}}+\frac{k_{m}}{\omega _{m}}\right)\frac{k_{n}}{\omega _{n}}\frac{\tanh k_{n}h}{k_{n}h}\notag\\[5pt]&\quad +\left(\frac{k_{m}}{\omega _{m}}+\frac{k_{n}}{\omega _{n}}\right)\frac{k_{m,n}}{\omega _{m,n}}\frac{\tanh k_{m,n}h}{k_{m,n}h}\notag\\[5pt]&\quad +\left[\left(\frac{k_{n}}{\omega _{m}}+\frac{k_{m}}{\omega _{n}}\right)\frac{k_{m,n}}{\omega _{m,n}}+\frac{k_{m,n}^{2}}{{\unicode[Arial]{x03C9}} _{m}\omega _{n}}\right]\frac{\tanh k_{m}h\tanh k_{n}h\tanh k_{m,n}h}{k_{m,n}h}.\end{align}
Equation (2.39) is a generalised subharmonic-induced mass transport solution for arbitrary depth and bandwidth and can be further simplified depending on the short wave relative depth (
$k_{n}h$
,
$k_{m}h$
), subharmonic wave relative depth (
$k_{m,n}h$
) and bandwidth
$\alpha$
. Detailed derivations are provided in Appendix A.
Noteworthily, the correlation between subharmonic velocity and surface elevation
$\overline{\tilde{u}_{0}\tilde{\eta }}$
only contributes partly to
$\bar{M}^{(4)}$
through the
$u_{0}^{(2)}\eta ^{(2)}$
term in (2.31). In fact, considering
$\overline{\tilde{u}_{0}\tilde{\eta }}$
only yields
\begin{equation}\overline{\tilde{u}_{0}\tilde{\eta }}=\sum _{m,n}\frac{1}{2}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{\tanh k_{m,n}h}{\omega _{m,n}}+\frac{1}{4}\frac{k_{m,n}}{\omega _{m,n}}\left(\frac{\tanh k_{m}h}{\omega _{m}h}+\frac{\tanh k_{n}h}{\omega _{n}h}\right)\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}h,\end{equation}
where we find the entire Term I in (2.39) but only part of the Term II.
2.2.2. Narrow-banded waves with wave groups long compared with the depth
Among the simplifications of (2.39) presented in Appendix A, the scenario of narrow-banded waves with wave groups long compared with the depth, i.e.
\begin{equation}\alpha \ll 1, \left| k_{m,n}\right| h\ll 1,\end{equation}
is consistent with the theory in § 2.1 and thus of particular interest in this study. The corresponding simplified (2.39) is given by (A22). Normalising
$c$
and
$c_{g}$
with
$\sqrt{gh}$
yields an alternative form of (A22)
\begin{equation}\bar{M}^{\left(4\right)}\left(x\right)=\underset{\text{Term}\textrm{I}}{\underbrace{\frac{\overline{\tilde{u}^{2}}h}{c_{g}}}}+\underset{\text{Term II}}{\underbrace{f(k_{p}h)\frac{\overline{\left[\left.u^{(1)}\right| _{z=0}\right]^{2}\tilde{u}}}{g}}},\end{equation}
where the function
$f(k_{p}h)$
is
\begin{equation}f\left(k_{p}h\right)=\frac{3\sinh 2k_{p}h+6k_{p}h\cosh ^{2} k_{p}h+2k_{p}h\sinh ^{2} k_{p}h}{\sinh 2k_{p}h+2k_{p}h}.\end{equation}
The function
$f(k_{p}h)$
increases monotonically with
$k_{p}h$
. For short waves in deep water
$k_{p}h\gg 1$
,
$f(k_{p}h)$
approaches
$\sim 3+4\omega _p^2h/g$
asymptotically; for shallow water
$k_{p}h\ll 1, f(k_{p}h)=3+O(k_{p}^{2}h^{2})$
. Note that, in deep water
$k_{p}h\gg 1$
, the shallow-water regime for subharmonic wave (
$|k_{m,n}|h\ll 1$
) is only valid for the subharmonic component with low enough frequencies or when the bandwidth
$\alpha$
is so small that
$|k_{m,n}|h=O(\alpha k_{p}h)\ll 1$
.
2.2.3. Second definition
Alternatively, let
$z=\eta _{c}(x)$
and
$z=\eta _{t}(x)$
be the wave crest and trough, i.e. the local maximum and minimum values of η, respectively. Below
$z=\eta _{c}(x)$
the time-averaged Eulerian velocity is zero while between
$z=\eta _{c}(x)$
and
$z=\eta _{t}(x)$
the velocity observed at a fixed height is intermittently zero similar to (1.4). The subharmonic wave-induced mass transport can be defined as the contribution of
$\tilde{u}$
to the integral
\begin{equation}\bar{M}\left(x\right)=\int _{\eta _{t}}^{\eta _{c}}\overline{u}\textrm{d}z.\end{equation}
This definition, however, is more difficult to analyse than (2.29) despite the fact that both should yield the same results. As the first step, substituting the expansion of
$u$
into (2.45) yields
\begin{equation}\bar{M}=\int _{\eta _{t}}^{\eta _{c}}\overline{u^{(1)}+u^{\left(2\right)}+u^{\left(3\right)}}\textrm{d}z+O(\varepsilon ^{5}).\end{equation}
Invoking
$u^{(2)}$
in (2.33), (2.46) becomes
\begin{equation}\bar{M}=\int _{\eta _{t}}^{\eta _{c}}\overline{u^{(1)}}+\overline{\tilde{u}}\textrm{d}z+\left(\text{terms without}\; \tilde{u}\right)+O(\varepsilon ^{5}).\end{equation}
As shown later,
$\int _{\eta _{t}}^{\eta _{c}}\overline{\tilde{u}}\textrm{d}z$
can be estimated from experimental data, which implies
$\int _{\eta _{t}}^{\eta _{c}}\overline{\tilde{u}}\textrm{d}z$
is equal to
$\overline{\tilde{u}\tilde{\eta }}$
although we do not seek its proof in this work. However, it remains challenging to analyse either analytically or experimentally the contribution of
$\tilde{u}$
to
$\int _{\eta _{t}}^{\eta _{c}}\overline{u^{(1)}}\textrm{d}z$
through terms analogous to
$u_{0}^{(1)}\eta ^{(3)}$
in (2.31).
2.3. Lagrangian transport
In the Lagrangian system, the Lagrangian velocity
$\boldsymbol{u}_{{L}}(t)=[u_{{L}}(t),w_{{L}}(t)]$
describes the velocity of a moving particle passively following the background Eulerian velocity
$\boldsymbol{u}(x,z,t)=[u(x,z,t),w(x,z,t)]$
. The mass transport is obtained by integrating vertically the steady component in
$u_{{L}}$
, i.e.
$\bar{M}=\int _{-h}^{0}\overline{u_{{L}}}\textrm{d}z$
. Compared with the bulk transport
$\bar{M}$
obtained in an Eulerian system,
$\overline{u_{{L}}}$
offers more information on the vertical distribution of how fluid particles at different depths physically drift.
2.3.1. Generalised theory
The horizontal and vertical displacements of a fluid particle released from
$\boldsymbol{r}_{0}=(x_{0},z_{0})$
at
$t=0$
underneath waves
${\unicode[Arial]{x0394}} \boldsymbol{r}(t)=[{\unicode[Arial]{x0394}} x(t),{\unicode[Arial]{x0394}} z(t)]$
are given by the time integration of its Lagrangian velocity
$\boldsymbol{u}_{{L}}(t)$
\begin{equation}{\unicode[Arial]{x0394}} \boldsymbol{r}(t)=\int \boldsymbol{u}_{{L}}(t)\textrm{d}t\,,\end{equation}
where
$\boldsymbol{u}_{{L}}(t)$
is obtained by querying the Eulerian velocity field following the particle coordinate in real time, i.e.
\begin{equation}\boldsymbol{u}_{{L}}\left(t\right)=\boldsymbol{u}\left(\boldsymbol{r}_{0}+{\unicode[Arial]{x0394}} \boldsymbol{r}\left(t\right)\!,t\right).\end{equation}
The right-hand side of (2.49) can be expanded around
$\boldsymbol{r}_{0}$
using Taylor series as
\begin{equation}\left.\begin{array}{c} \boldsymbol{u}_{{L}}\left(t\right)=\left\{\left(1+{\unicode[Arial]{x0394}} \boldsymbol{r}{\boldsymbol\cdot} \nabla +\dfrac{1}{2}\left[{\unicode[Arial]{x0394}} \boldsymbol{r}{\boldsymbol\cdot} \left({\unicode[Arial]{x0394}} \boldsymbol{r}{\boldsymbol\cdot} \boldsymbol{H }\right)\right]+\ldots \right)\boldsymbol{u}\right\}_{\left(\boldsymbol{r}_{0},t\right)}\\ \nabla =\left(\dfrac{\partial }{\partial x},\dfrac{\partial }{\partial z}\right)\!,\boldsymbol{H }=\left[\begin{array}{c@{\quad}c} \dfrac{\partial ^{2}}{\partial x^{2}} & \dfrac{\partial ^{2}}{\partial x\partial z}\\[10pt] \dfrac{\partial ^{2}}{\partial x\partial z} & \dfrac{\partial ^{2}}{\partial x^{2}} \end{array}\right] \end{array}\right\},\end{equation}
where the subscript coordinates
$(\boldsymbol{r}_{0},t)$
will be omitted hereinafter.
Applying Stokes expansion in
$\varepsilon$
to
$\boldsymbol{u}$
,
$\boldsymbol{u}_{{L}}$
and
${\unicode[Arial]{x0394}} \boldsymbol{r}$
and then substituting them into (2.50) yields
\begin{align} \boldsymbol{u}_{{L}}^{(1)}(x,z,t)=\boldsymbol{u}^{(1)},\end{align}
\begin{align} \boldsymbol{u}_{{L}}^{\left(2\right)}(x,z,t)=\boldsymbol{u}^{\left(2\right)}+{\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \nabla \boldsymbol{u}^{(1)},\end{align}
\begin{align} \boldsymbol{u}_{{L}}^{\left(3\right)}(x,z,t)=\boldsymbol{u}^{\left(3\right)}+{\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \nabla \boldsymbol{u}^{\left(2\right)}+{\unicode[Arial]{x0394}} \boldsymbol{r}^{\left(2\right)}{\boldsymbol\cdot} \nabla \boldsymbol{u}^{(1)}+\frac{1}{2}\big[{\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \big({\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \boldsymbol{H }\big)\big]\boldsymbol{u}^{(1)},\end{align}
\begin{align}\boldsymbol{u}_{{L}}^{\left(4\right)}(x,z,t)&=\boldsymbol{u}^{\left(4\right)}+{\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \nabla \boldsymbol{u}^{\left(3\right)}+{\unicode[Arial]{x0394}} \boldsymbol{r}^{\left(2\right)}{\boldsymbol\cdot} \nabla \boldsymbol{u}^{\left(2\right)}+{\unicode[Arial]{x0394}} \boldsymbol{r}^{\left(3\right)}{\boldsymbol\cdot} \nabla \boldsymbol{u}^{(1)}\notag\\&\quad+\frac{1}{2}\big[{\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \big({\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \boldsymbol{H }\big)\big]\boldsymbol{u}^{\left(2\right)}+\big[{\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}{\boldsymbol\cdot} \big({\unicode[Arial]{x0394}} \boldsymbol{r}^{\left(2\right)}{\boldsymbol\cdot} \boldsymbol{H }\big)\big]\boldsymbol{u}^{(1)}\notag\\&\quad+\frac{1}{6}\left(\big[{\unicode[Arial]{x0394}} x^{(1)}\big]^{3}\frac{\partial ^{3}}{\partial x^{3}}+3\big[{\unicode[Arial]{x0394}} x^{(1)}\big]^{2}{\unicode[Arial]{x0394}} z^{(1)}\frac{\partial ^{3}}{\partial x^{2}\partial z}+3\big[{\unicode[Arial]{x0394}} z^{(1)}\big]^{2}\!{\unicode[Arial]{x0394}} x^{(1)}\frac{\partial ^{3}}{\partial x\partial z^{2}}\right .\notag\\&\quad +\left.\big[{\unicode[Arial]{x0394}} z^{(1)}\big]^{3}\frac{\partial ^{3}}{\partial z^{3}}\right)\boldsymbol{u}^{(1)}.\end{align}
The horizontal component of
$\boldsymbol{u}^{(1)}=[u^{(1)},w^{(1)}]$
can be found from (2.4), while the vertical component
$w^{(1)}$
can be obtained by integrating the continuity equation vertically
\begin{equation}w^{(1)}(x,z,t)=w^{(1)}\left(x,-h,t\right)-\int\limits_{-h}^{z}\frac{\partial }{\partial x}u^{(1)}(x,z',t)\textrm{d}z',\end{equation}
where
$w^{(1)}(x,-h,t)$
is the bottom vertical velocity, which is of
$O(\varepsilon h_{x})$
due to kinematic condition
$w=-uh_{x}$
at the bottom
$z=-h$
. Because
$| h_{x}| \ll O(\varepsilon kh)$
, for
$kh=O(1)$
we have
$w^{(1)}(x,-h,t)=O(\varepsilon h_{x})\ll O(\varepsilon ^{2})$
while for
$kh\gg 1$
the bottom effect is neglected. Therefore,
$w^{(1)}(x,-h,t)$
will be excluded in both
$w^{(1)}$
and
$w^{(2)}$
. Substituting (2.4) into (2.52) and taking
$w^{(1)}(x,-h,t)=0$
yields
\begin{equation}w^{(1)}(x,z,t)=\sum _{n=n_{1}}^{n_{1}+N}-\textrm{i}\hat{u}_{n}(x)\frac{\sinh \left[k_{n}(z+h)\right]}{\cosh k_{n}h}\exp \textrm{i}\theta _{n}.\end{equation}
Substituting (2.4) and (2.53) into
$\boldsymbol{u}^{(1)}=[u^{(1)},w^{(1)}]$
and then into (2.51a
) yields
$\boldsymbol{u}_{{L}}^{(1)}$
. Then the
${\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}$
can be obtained following
${\unicode[Arial]{x0394}} \boldsymbol{r}^{(i)}=\int \boldsymbol{u}_{{L}}^{(i)}(t)\textrm{d}t$
. After
${\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}$
is obtained,
$\boldsymbol{u}^{(2)}$
is still needed to evaluate
$\boldsymbol{u}_{{L}}^{(2)}$
. The horizontal component of
$\boldsymbol{u}^{(2)}=[u^{(2)},w^{(2)}]$
can be found from (2.33), while the vertical component
$w^{(2)}$
is obtained similar to
$w^{(1)}$
as
\begin{align} \begin{cases} w^{\left(2\right)}=\tilde{w}+\left(\text{superharmonic terms}\right)\\[5pt] \tilde{w}(x,z,t)=\sum\limits_{n,m}-\textrm{i}\widehat{U}_{m,n}\dfrac{\sinh \left[k_{m,n}(z+h)\right]}{\cosh k_{m,n}h}\exp \textrm{i}\theta _{m,n} \end{cases} \end{align}
Substituting
$\boldsymbol{u}^{(2)}=[u^{(2)},w^{(2)}]$
,
${\unicode[Arial]{x0394}} \boldsymbol{r}^{(1)}$
and
$\boldsymbol{u}^{(1)}$
into (2.50) yields
$\boldsymbol{u}_{{L}}^{(2)}$
and, by
${\unicode[Arial]{x0394}} \boldsymbol{r}^{(i)}=\int \boldsymbol{u}_{{L}}^{(i)}(t)\textrm{d}t$
,
${\unicode[Arial]{x0394}} \boldsymbol{r}^{(2)}$
. Repeating this process, the expression of
$\boldsymbol{u}_{{L}}^{(4)}$
can be worked out after lengthy algebra. Since we are interested in the resultant
$\overline{u_{{L}}^{(4)}}$
involving
$\widehat{U}_{m,n}$
, it will be useful to write out the expression of
$u_{{L}}^{(4)}$
\begin{align}u_{{L}}^{\left(4\right)}(x,z,t)&={\unicode[Arial]{x0394}} x^{\left(2\right)}\frac{\partial \tilde{u}}{\partial x}+{\unicode[Arial]{x0394}} z^{\left(2\right)}\frac{\partial \tilde{u}}{\partial z}+{\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} x^{\left(2\right)}\frac{\partial ^{2}u^{(1)}}{\partial x^{2}}\notag\\[5pt]&\quad +\big[{\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} z^{\left(2\right)}+{\unicode[Arial]{x0394}} z^{(1)}{\unicode[Arial]{x0394}} x^{\left(2\right)}\big]\frac{\partial ^{2}u^{(1)}}{\partial x\partial z}+{\unicode[Arial]{x0394}} z^{(1)}{\unicode[Arial]{x0394}} z^{\left(2\right)}\frac{\partial ^{2}u^{(1)}}{\partial z^{2}}+{\unicode[Arial]{x0394}} x^{\left(3\right)}\frac{\partial u^{(1)}}{\partial x}\notag\\[5pt]&\quad +{\unicode[Arial]{x0394}} z^{\left(3\right)}\frac{\partial u^{(1)}}{\partial z}+\frac{1}{2}\big[{\unicode[Arial]{x0394}} x^{(1)}\big]^{2}\frac{\partial ^{2}\tilde{u}}{\partial x^{2}}+{\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} z^{(1)}\frac{\partial ^{2}\tilde{u}}{\partial x\partial z}+\frac{1}{2}\big[{\unicode[Arial]{x0394}} z^{(1)}\big]^{2}\frac{\partial ^{2}\tilde{u}}{\partial z^{2}}\notag\\[5pt]&\quad +(\text{terms without }\tilde{u}\text{ or }\tilde{w}),\end{align}
and evaluate each part separately. Note that in (2.55) the exclusion of terms without
$\tilde{u}\text{ or }\tilde{w}$
were not exhausted (e.g. expanding
${\unicode[Arial]{x0394}} x^{(2)}$
in
${\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} x^{(2)}({\partial ^{2}u^{(1)}})/({\partial x^{2}})$
also generates those terms). The steady terms in each part are
\begin{align}\overline{{\unicode[Arial]{x0394}} x^{\left(2\right)}\frac{\partial \tilde{u}}{\partial x}+{\unicode[Arial]{x0394}} z^{\left(2\right)}\frac{\partial \tilde{u}}{\partial z}}&=\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{k_{m,n}}{\omega _{m,n}}\frac{\cosh \left[2k_{m,n}(z+h)\right]}{2\cosh ^{2} k_{m,n}h}\notag\\&\quad +\frac{B_{m,n}^{\textrm{I}}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{4\cosh k_{m}h\cosh k_{n}h\cosh k_{m,n}h}+(\text{terms without }\tilde{u}\text{ or }\tilde{w}),\end{align}
\begin{align}\overline{{\unicode[Arial]{x0394}} x^{\left(3\right)}\frac{\partial u^{(1)}}{\partial x}+{\unicode[Arial]{x0394}} z^{\left(3\right)}\frac{\partial u^{(1)}}{\partial z}}&=\sum _{m,n}\frac{B_{m,n}^{\textrm{II}}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{4\cosh k_{m}h\cosh k_{n}h\cosh k_{m,n}h}\notag\\&\quad +(\text{terms without }\tilde{u}\text{ or }\tilde{w}),\end{align}
\begin{align}&\overline{{\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} x^{\left(2\right)}\frac{\partial ^{2}u^{(1)}}{\partial x^{2}}+\left[{\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} z^{\left(2\right)}+{\unicode[Arial]{x0394}} z^{(1)}{\unicode[Arial]{x0394}} x^{\left(2\right)}\right]\frac{\partial ^{2}u^{(1)}}{\partial x\partial z}+{\unicode[Arial]{x0394}} z^{(1)}{\unicode[Arial]{x0394}} z^{\left(2\right)}\frac{\partial ^{2}u^{(1)}}{\partial z^{2}}}\notag\\&\quad=\sum _{m,n}\frac{B_{m,n}^{\textrm{III}}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{4\cosh k_{m}h\cosh k_{n}h\cosh k_{m,n}h}+(\text{terms without }\tilde{u}\text{ or }\tilde{w}),\end{align}
\begin{align}&\overline{\frac{1}{2}\left[{\unicode[Arial]{x0394}} x^{(1)}\right]^{2}\frac{\partial ^{2}\tilde{u}}{\partial x^{2}}+{\unicode[Arial]{x0394}} x^{(1)}{\unicode[Arial]{x0394}} z^{(1)}\frac{\partial ^{2}\tilde{u}}{\partial x\partial z}+\frac{1}{2}\left[{\unicode[Arial]{x0394}} z^{(1)}\right]^{2}\frac{\partial ^{2}\tilde{u}}{\partial z^{2}}}\notag\\&\quad =\sum _{m,n}\frac{B_{m,n}^{\text{IIII}}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{4\cosh k_{m}h\cosh k_{n}h\cosh k_{m,n}h}+(\text{terms without }\tilde{u}\text{ or }\tilde{w}),\end{align}
with the coefficients
$B_{m,n}^{\textrm{I}}$
to
$B_{m,n}^{\text{IIII}}$
\begin{align}\left[\begin{array}{c} B_{m,n}^{\textrm{I}}\\[8pt]B_{m,n}^{\textrm{II}}\\[8pt] B_{m,n}^{\textrm{III}}\\[8pt] B_{m,n}^{\text{IIII}} \end{array}\right]=\left[\begin{array}{c} \dfrac{k_{m,n}}{\omega _{m,n}}\left\{\dfrac{k_{n}}{\omega _{m}}\cosh \left[2k_{n}(z+h)\right]+\dfrac{k_{m}}{\omega _{n}}\cosh \left[2k_{m}(z+h)\right]\right\}\\[16pt] \begin{array}{c} \dfrac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\dfrac{\left\{\omega \cosh \left[2k(z+h)\right]\right\}_{m,n}}{\omega _{m,n}}-\dfrac{k_{m,n}^{2}}{\omega _{m}\omega _{n}}\cosh \left[2k_{m,n}(z+h)\right]\\[18pt] -\dfrac{k_{m,n}^{2}}{\omega _{m}\omega _{n}}\cosh \left[2k_{m,n}(z+h)\right]\\[18pt] \dfrac{\left\{\omega k^{2}\cosh \left[2k(z+h)\right]\right\}_{m,n}}{\omega _{m}\omega _{n}\omega _{m,n}} \end{array} \end{array}\right]\!,\end{align}
where all the subscripts
$m,n$
on the right-hand side denote the difference between the mth and nth frequency components. Collecting (2.56)–(2.59) yields
$\overline{u_{{L}}^{(4)}}$
\begin{align}\overline{u_{{L}}^{\left(4\right)}}(x,z)&=\sum _{m,n}\underset{\text{Term}\;\textrm{I}}{\underbrace{\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{k_{m,n}}{\omega _{m,n}}\frac{\cosh \left[2k_{m,n}(z+h)\right]}{2\cosh ^{2} k_{m,n}h}}}\notag\\[8pt]&\quad +\underset{\text{Term II}}{\underbrace{\frac{\left(B_{m,n}^{\textrm{I}}+B_{m,n}^{\textrm{II}}+B_{m,n}^{\textrm{III}}+B_{m,n}^{\text{IIII}}\right)\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{4\cosh k_{m}h\cosh k_{n}h\cosh k_{m,n}h}}},\end{align}
and integrating
$\overline{u_{{L}}^{(4)}}$
vertically yields the mass transport
$\bar{M}^{(4)}(x)=\int _{-h}^{0}\overline{u_{{L}}^{(4)}}(x,z)\textrm{d}z$
\begin{equation}\bar{M}^{\left(4\right)}\left(x\right)=\sum _{m,n}\frac{1}{2}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{\tanh k_{m,n}h}{\omega _{m,n}}+\frac{1}{4}D_{m,n}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}h,\end{equation}
with
\begin{align}D_{m,n}&=\left(\frac{1}{\omega _{m}}\frac{\tanh k_{m}h-\tanh k_{m,n}h}{h}+\frac{1}{\omega _{n}}\frac{\tanh k_{n}h+\tanh k_{m,n}h}{h}\right)\frac{k_{m,n}}{\omega _{m,n}}\notag\\[6pt]&\quad +\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{n}h+\tanh k_{m,n}h}{k_{m}h}-\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{m}h-\tanh k_{m,n}h}{k_{n}h}\right)\notag\\[6pt]&\quad -\frac{2k_{m,n}^{2}}{\omega _{m}\omega _{n}}\frac{\tanh k_{m}h-\tanh k_{n}h}{k_{m,n}h}\notag\\[6pt]&\quad +\frac{\omega _{m}k_{m}\left(\tanh k_{n}h+\tanh k_{m,n}h\right)-\omega _{n}k_{n}\left(\tanh k_{m}h-\tanh k_{m,n}h\right)}{\omega _{m}\omega _{n}\omega _{m,n}h}.\end{align}
Term I on the right-hand side of (2.62) is the same as in (2.39) obtained in Eulerian system. Also, both equations are equivalent for the coefficient
$D_{m,n}$
(see Appendix B for a proof). Therefore, the result of
$\bar{M}^{(4)}$
derived in the Lagrangian system (2.62) is consistent with that in the Eulerian system (2.39).
At the leading order, the mass transport is the sum of Term I and II, as labelled in (2.39) and (2.61). Term I originates from the term
$\overline{\int \boldsymbol{u}^{(2)}\textrm{d}t{\boldsymbol\cdot} \nabla \boldsymbol{u}^{(2)}}$
and accounts for the net displacement solely due to the subharmonic motion (self-interaction), in analogy to the Stokes drift velocity due to the primary wave orbital motion. Term II takes the form of the weighted sum of all bispectral components of horizontal velocity skewness
$\textrm{Re}[\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}]$
arising from the triad-difference interaction among two primary and one subharmonic components. Therefore, Term II may be regarded as the effect of subharmonic wave-induced skewness on the mass transport velocity. Note that these two terms do not necessarily balance with each other, and the residual of their sum has to be compensated at least partly by a steady component of
$u^{(4)}$
depending on the boundary conditions.
2.3.2. Narrow-banded waves with wave groups long compared with the depth
As shown in Appendix C, (2.61) for
$\overline{u_{L}^{(4)}}$
can be further simplified depending on the short wave relative depth (
$k_{n}h$
,
$k_{m}h$
), subharmonic wave relative depth (
$k_{m,n}h$
) and bandwidth
$\alpha$
similar to the simplifications to (2.39) in Appendix A. In the case of narrow-banded waves with wave groups that are long compared with the depth, equation (2.61) reduces to (C16), i.e.
\begin{equation}\overline{u_{L}^{\left(4\right)}}(x,z)=\frac{\overline{\tilde{u}^{2}}}{c_{g}}+\zeta \left(k_{p}h,\frac{z}{h}\right)\frac{\overline{\left[\left.u^{(1)}\right| _{z=0}\right]^{2}\tilde{u}}}{cc_{g}},\end{equation}
where
\begin{equation}\zeta \left(k_{p}h,\frac{z}{h}\right)=\left[\frac{5}{2}+\frac{k_{p}h}{\sinh 2k_{p}h}+2k_{p}(z+h)\tanh \left[2k_{p}(z+h)\right]\right]\frac{\cosh \left[2k_{p}(z+h)\right]}{\cosh ^{2} k_{p}h}.\end{equation}
The second term of (2.64) is the product of the characteristic transport velocity
$\overline{[u^{(1)}| _{z=0}]^{2}\tilde{u}}/(cc_{g})$
and the depth-dependent coefficient
$\zeta$
which is a function of
$k_{p}h$
and the
$z/h$
, as shown in figure 3(a). For short waves in deep water
$k_{p}h\gg 1$
,
$\zeta$
approaches
$({4\omega _{p}^{2}h})/({g})\exp({2\omega _{p}^{2}z}/{g})$
asymptotically; for shallow water
$k_{p}h\ll 1, \zeta$
becomes
$3+O(k_{p}^{2}h^{2})$
.
Figure 3. (a) Contour line of
$\zeta (k_{p}h,z/h)$
in (2.64) which describes the vertical distribution of the skewness-induced transport velocity due to subharmonics proportional to
$\overline{[u^{(1)}| _{z=0}]^{2}\tilde{u}}/(cc_{g})$
at different short wave relative depths
$k_{p}h$
. (b) Theoretical prediction of (2.66), i.e. the ratio of Term I (‘Stokes drift’ type contribution) over Term II (velocity skewness contribution) for each pair of bichromatic waves in the subharmonic wave-induced mass transport given by (2.43). Dashed lines: flat bottom solution of Longuet-Higgins & Stewart (Reference Longuet-Higgins and Stewart1962); solid lines: Green’s function-based solution in § 2.1 based on Liao et al. (Reference Liao, Zou, Liu, Contardo and Li2023).
$\beta =h_{x}/(k_{l}h)$
is the relative bottom slope at wave group length scale,
$\mu =1-c_{g}^{2}/(gh)$
is the parameter of departure from resonance.
In (2.64) the weighting coefficient for all bispectral components of velocity skewness
$Re[\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}]$
in (2.61) is the same and therefore Term II is simply proportional to
$\sum _{m,n}Re[\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}]\propto \overline{[u^{(1)}| _{z=0}]^{2}\tilde{u}}$
, which is non-zero due to phase coupling between the subharmonic wave and wave groups. As shown by Fiedler et al. (Reference Fiedler, Smit, Brodie, McNinch and Guza2019), the correlation
$\overline{\tilde{u}[u^{(1)}]^{2}}$
can be found through expanding the velocity skewness
$\overline{(u^{(1)}+\tilde{u})^{3}}$
as the contribution of subharmonic waves to the total velocity skewness.
For each pair of m and n, the ratio |Term I/Term II| in
$\bar{M}^{(4)}$
by (2.43) is given by
\begin{equation}\frac{2gh}{f(k_{p}h)c_{g}}\frac{\big| \widehat{U}_{m,n}\big| ^{2}}{\big| Re\big[\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}\big]\big| },\end{equation}
which relies on a function of
$k_{p}h$
and the dependence of
$\widehat{U}_{m,n}$
on
$\hat{u}_{m}^{*}\hat{u}_{n}$
. This dependence can be evaluated based on the equilibrium subharmonic wave solution on a flat bottom by Longuet-Higgins & Stewart (Reference Longuet-Higgins and Stewart1962) and also the Green’s function-based solution in § 2.1 and shown in figure 3(b). Both models indicate that Term I exceeds Term II for
$\mu \lt 0.36$
(
$k_{p}h\lt$
0.72), while the Green’s function-based solution indicates decreasing critical
$\mu$
or
$k_{p}h$
as |β| increases and exceeds ∼0.01. The relative importance of Term I increases with decreasing μ due to increasing subharmonic wave energy during shoaling and with decreasing β. As β decreases, the predictions of the Green’s function-based solution asymptoticly approaches that of the equilibrium solution for a flat bottom. Notice that at given depth for equilibrium solution or small enough relative slope
$| \beta |$
this ratio becomes independent of
$\beta$
and hence subharmonic wavenumber
$k_{l}$
, suggesting that it also measures the ratio between the entire Term I over Term II.
2.4. Comparisons with GLOBEX experiment
2.4.1. Experimental set-up and data processing
The test series A of GLOBEX experiment (Ruessink et al. Reference Ruessink, Michallet, Bonneton, Mouaze and Lara2013) of unidirectional irregular waves normally incident onto a beach of 1/80 bed slope is employed to validate the model in § 2.1 for subharmonic surface elevation and velocity. Then, the theoretical model in §§ 2.2 and 2.3 is evaluated both analytically and experimentally to analyse the subharmonic wave-induced mass transport. Note that the experiment did not provide direct measurement of the kinematic effect of subharmonic wave-induced mass transport and therefore does not allow validation of the mass transport model proposed. This will be addressed in § 3 through comparative Lagrangian particle-tracking simulations.
Because the model in § 2.1 ignores the reflected subharmonic waves, the incident subharmonic surface elevation
$\tilde{\eta }$
and orbital velocity
$\tilde{u}$
were obtained using the weakly nonlinear method detailed in Appendix D. The significant wave height of the incoming subharmonic wave was then obtained as
$H_{lf}=4\sqrt{\overline{(\tilde{\eta })^{2}}}$
. The phase coupling between subharmonic waves and short wave groups (biphase) was calculated following Fiedler et al., (Reference Fiedler, Smit, Brodie, McNinch and Guza2019) as
\begin{equation} \text{Biphase}=\tan ^{-1} \left(\frac{\overline{{\mathcal{H}}\left\{\tilde{\eta }\right\}{\mathcal{H}}\left\{\eta ^{(1)}\right\}^{2}}}{\overline{\tilde{\eta }\left[\eta ^{(1)}\right]^{2}}}\right)\!, \end{equation}
where
$\mathcal{H}\{{\boldsymbol\cdot} \}$
denotes the Hilbert transform and
$\eta ^{(1)}$
the short wave surface elevation with frequency >
$0.5/T_{p}$
.
To calculate the subharmonic wave-induced time-averaged velocity
$\overline{\tilde{u}}$
between wave crest
$z=\eta _{c}$
and trough
$z=\eta _{t}$
described in § 2.2.3, the incident subharmonic orbital velocity
$\tilde{u}$
measured at a fixed height above the bed was assumed to be uniform throughout the water column. At each cross-shore location with velocity measurements, 21 equidistant points from
$z=\eta _{t}$
to
$z=\eta _{c}$
were defined and at each point
$\tilde{u}$
is assigned with the simultaneous value of the
$\tilde{u}$
measured for the timepoints when it is below the free surface and otherwise zero. Time averaging the resultant
$\tilde{u}$
at the defined points yields the
$\overline{\tilde{u}}$
between the wave crest and trough.
To experimentally evaluate (2.43),
$\overline{\tilde{u}^{2}}$
is calculated from the incident subharmonic
$\tilde{u}$
while
$\overline{\tilde{u}[u^{(1)}| _{z=0}]^{2}}$
is evaluated using the measured
$\overline{\tilde{u}[u^{(1)}]^{2}}$
at
$z_{{p}}+h$
above the bed according to
\begin{equation} \overline{\tilde{u}\left[\left.u^{(1)}\right| _{z=0}\right]^{2}}=\left.\left\{\overline{\tilde{u}\left[u^{(1)}\right]^{2}}\right\}\right| _{z={z_{p}}}\frac{\cosh ^{2} k_{p}h}{\cosh ^{2} k_{p}\left(z_{{p}}+h\right)}.\end{equation}
2.4.2. Validation of theoretical model for subharmonic wave
Figure 4(a–c) shows the flume set-up, significant wave height of the primary wave and the incoming subharmonic wave of A3 case of GLOBEX experiment. In figure 4(c), also shown is the theoretical wave height of subharmonic waves forced by wave groups,
$H_{lf}=4\sqrt{\sum (\hat{\eta }_{l}\hat{\eta }_{l}^{*}/2)}$
, where the complex amplitude of subharmonic wave
$\hat{\eta }_{l}$
is given by (2.23). The theory predicted the observed growth of group-forced subharmonic waves in the shoaling zone and the outer surf zone very well, but overpredicted the observation in the inner surf zone. The overprediction in the surf zone is likely due to drastic modification of the phase of group forcing by a breaking process not included in the model (figure 4
d), which is crucial for capturing the wave height decay in the inner surf zone (Liao et al. Reference Liao, Zou, Liu, Contardo and Li2023).
Figure 4. (a) Numerical flume set-up and wave conditions of Case A3 in the GLOBEX experiment. Significant wave height of (b) short wave (c) incident subharmonic waves (frequency <
$0.5/T_{p}$
). (d) Biphase for the triad-interaction between incident subharmonic waves and short wave groups. (e) Measured time-averaged Eulerian velocity induced by subharmonic waves
$\overline{\tilde{u}}$
between the maximum (crest,
$z=\eta _{c}$
) and minimum (trough,
$z=\eta _{t}$
) surface elevation level.
The subharmonic wave-induced time-averaged velocity
$\overline{\tilde{u}}$
between wave crest
$z=\eta _{c}$
and trough
$z=\eta _{t}$
is illustrated in figure 4(e). Below the wave trough
$z=\eta _{t}$
, the oscillatory subharmonic wave leads to a zero time-averaged velocity
$\overline{\tilde{u}}$
. In the shoaling zone and outer surf zone,
$\overline{\tilde{u}}$
is negative near the wave crest
$z=\eta _{c}$
but positive near the wave trough
$z=\eta _{t}$
. This is because only the largest waves in the wave train can achieve the crest, when the subharmonic wave troughs and backward orbital velocity occurs since the group-forced subharmonic wave is in anti-phase with the wave group. The dominant forward velocity near the subharmonic wave trough leads to a net positive contribution of the term
$\int _{\eta _{t}}^{\eta _{c}}\overline{\tilde{u}}\textrm{d}z$
to the mass transport (2.47), which increases with decreasing depth towards the shore. In the inner surf zone, the negative velocity near the surface becomes positive, due to the transition from anti-phase to in-phase coupling between the subharmonic wave and wave groups (figure 4
d). Although the model is no longer accurate for x> 70 m in this region, the
$\overline{\tilde{u}}$
value evaluated from experimental measurements still offers insights into the mass transport, especially the effect of growing subharmonic waves and decaying short waves.
Figure 5 compares the measured and theoretically predicted mass transport component
$\rho \overline{\tilde{u}\tilde{\eta }}$
in the shoaling zone. In all three cases, the theoretical predictions agree with the observations very well, validating the model in § 2.1 for subharmonic surface elevation and velocity. Noteworthily, the experimental results of the term
$\int _{\eta _{t}}^{\eta _{c}}\overline{\tilde{u}}\textrm{d}z$
are nearly identical to
$\overline{\tilde{u}\tilde{\eta }}$
, suggesting that
$\int _{\eta _{t}}^{\eta _{c}}\overline{\tilde{u}}\textrm{d}z=\overline{\tilde{u}\tilde{\eta }}$
although a theoretical demonstration is desired.
Figure 5. Subharmonic wave-induced mass transport components
$\rho \overline{\tilde{u}\tilde{\eta }}$
and
$\rho \int _{\eta _{t}}^{\eta _{c}}\overline{\tilde{u}}\textrm{d}z$
in the shoaling zone for (a) A1, (b) A2 and (c) A3 irregular wave cases of GLOBEX experiment. Dots: experiment; lines: theoretical predictions. Here,
$\rho =1000 (\textrm{kg}\,\textrm{m}^{-3})$
was adopted.
2.4.3. Subharmonic wave-induced mass transport
Figure 6 shows the mass transport
$\bar{M}^{(4)}$
(2.43) evaluated in the shoaling zone using measured and theoretical velocity statistics (
$\overline{\tilde{u}^{2}}$
and
$\overline{\tilde{u}[u^{(1)}| _{z=0}]^{2}}$
). As mentioned above, the experiment did not actually measure
$\bar{M}^{(4)}$
and thus does not allow comparison with (2.43). The comparison in figure 6 is essentially between the experimental and theoretical subharmonic velocity instead of
$\bar{M}^{(4)}$
.
Figure 6. Subharmonic wave-induced mass transport
$\bar{M}^{(4)}$
(2.43) in the shoaling zone for (a) A1, (b) A2 and (c) A3 irregular wave cases of GLOBEX experiment. Dots: experiment; lines: theoretical predictions. Red and blue markers and lines denote the Term I (‘Stokes drift’ effect) and Term II (skewness effect) components in (2.43), respectively. Green crosses denote the same result as the dots in figure 5, i.e. the component
$\overline{\tilde{u}\tilde{\eta }}$
. Here,
$\rho =1000\;(\textrm{kg}\,\textrm{m}^{-3})$
was adopted.
The theory agrees with the experiment very well for case A1. For cases A2 and A3 the theory slightly underestimated the magnitudes of both term I (‘Stokes drift’ effect) and term II (skewness effect), indicating an underestimation of subharmonic velocity magnitude. Nevertheless, the theory aligns with experimental results regarding total transport, primarily due to the opposite signs of the two terms that cancel out the underestimations.
Term I contributes to a forward drift as it is proportional to the energy of the subharmonic velocity, whereas Term II accounts for a backward drift because the near anti-phase coupling between subharmonic wave and wave group results in a negative
$\overline{\tilde{u}[u^{(1)}]^{2}}$
. In all three cases as depth decreases towards the shore the total transport is initially negative (seaward) and turns positive (shoreward) near the end of the shoaling zone due to the positive Term I growing faster than the negative Term II and becoming predominant prior to wave breaking. The results also suggest that the transition of total transport occurs at shallower relative depth
$k_p{h}$
for larger wave period and smaller bandwidth.
The measured Eulerian mass transports (crosses) in all three cases were lower than term I (red dots) alone but greater than the total transport (black dots), especially for a smaller period (figure 6
a) and greater depth. This aligns with (2.41) that
$\overline{\tilde{u}\tilde{\eta }}$
(the subscript of
$\tilde{u}_{0}$
for
$z=0$
omitted since subharmonic waves are in shallow water in this case) accounts for the entire Term I contribution but only a fraction of the Term II contribution that is generally negative. It also stresses the necessity to include those negative velocity skewness contributions not accounted for by
$\overline{\tilde{u}\tilde{\eta }}$
to correct its overprediction.
3. Numerical model of subharmonic wave-induced mass transport
3.1. Non-hydrostatic wave model SWASH
To validate the theoretical model for subharmonic wave-induced mass transport, a Lagrangian particle-tracking simulation is conducted, which requires a high-resolution flow field in space and time not available from the experiment. Thus, the non-hydrostatic wave model SWASH (Zijlema et al. Reference Zijlema, Stelling and Smit2011) was used to investigate the subharmonic-induced mass transport for bichromatic wave groups propagating over a plane sloping ramp. The multilayer non-hydrostatic model SWASH is essentially a RANS (Reynolds-Averaged Navier–Stokes) solver without VOF (Volume Of Fluid) and level set type of free-surface capturing scheme, capable of resolving strong nonlinear wave motion. It solves the momentum and mass conservation equations
$$\begin{align}& \frac{\partial u}{\partial t}+\frac{\partial uu}{\partial x}+\frac{\partial wu}{\partial z}=-\frac{1}{\rho }\frac{\partial (p_{h}+p_{nh})}{\partial x}+\frac{\partial \tau _{xz}}{\partial z}+\frac{\partial \tau _{xx}}{\partial x},\end{align}$$
$$\begin{align}&\quad\,\, \frac{\partial w}{\partial t}+\frac{\partial uw}{\partial x}+\frac{\partial ww}{\partial z}=-\frac{1}{\rho }\frac{\partial p_{nh}}{\partial z}+\frac{\partial \tau _{zz}}{\partial z}+\frac{\partial \tau _{zx}}{\partial x},\end{align}$$
$$\begin{align}&\qquad\qquad\qquad\qquad\,\, \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0,\end{align}$$
where p
h
is the hydrostatic pressure, p
nh
is the non-hydrostatic pressure and τ is the turbulent stress described by the standard K-ε model, where
$K$
is the turbulent kinetic energy and
$\epsilon$
is its dissipation rate in the present work. The model has been extensively validated against both laboratory (de Bakker, Tissier & Ruessink Reference de Bakker, Tissier and Ruessink2016b
) and field data (Rijnsdorp, Ruessink & Zijlema Reference Rijnsdorp, Ruessink and Zijlema2015) of the nearshore nonlinear evolution of subharmonic waves. Detailed descriptions of the model can be found in Rijnsdorp, Smit & Zijlema (Reference Rijnsdorp, Smit and Zijlema2014) and Smit et al. (Reference Smit, Janssen, Holthuijsen and Smith2014).
3.1.1. Model set-up and data processing
The numerical flume has a bed profile of a plane sloping ramp with bottom slope ranging between 1/30 and 1/120 connected with two flat bottom plateaux (figure 7). The toe of the ramp is 300 m away from the incident boundary and the plateau at the shallow end is 480 m long equipped with a 140 m sponge layer near the outgoing boundary to minimise reflection of short waves. The depth at the shallow end was carefully chosen to prevent wave breaking, thereby no surf zone occurs in the model domain. Eight sigma layers with thicknesses 6.06, 7.27, 8.73, 10.47, 12.57, 15.08, 18.10, 21.71 per cent of full depth were adopted. The bottom boundary layer and the associated mass transport caused by boundary layer streaming was not considered in the numerical model by adopting a zero bottom friction factor. The bathymetry and wave condition parameters are summarised in table 1.
Figure 7. Schematic configuration of numerical flume set-up used for the SWASH model;
$h_{x}=\textrm{d}h/\textrm{d}x$
is the bottom slope.
Table 1. Bathymetry parameters and bichromatic wave conditions for SWASH model.
Preliminary analysis of incoming and outgoing subharmonic waves using the same method as in § 2.4 indicated that the reflection is less than 5 % for the subharmonic wave amplitude. The following analysis is based on the result before decomposition into incoming/outgoing components as there is no fundamental difference.
The SWASH model outputs the horizontal velocities at the centre of each vertical layer, where the absolute z-coordinate varies in time with the undulating surface elevation. To derive the Eulerian velocity at fixed heights, the velocity output at a time-dependent height was first vertically mapped onto fixed heights, including ten equidistant layers between the wave crest
$z=\eta _{t}$
and wave trough
$z=\eta _{c}$
and eight sigma layers between the bottom
$z=-h$
and wave trough
$z=\eta _{t}$
, using spline interpolation.
The horizontal Eulerian velocity
$u$
and surface elevation
$\eta$
were first demeaned and then low-pass filtered with cutoff frequency of
$0.5\omega$
, where
$\omega =(\omega _{1}+\omega _{2})/2$
is the central frequency of bichromatic waves, to obtain the surface elevation
$\tilde{\eta }$
and Eulerian velocity
$\tilde{u}$
of the subharmonic wave. The first-order velocity
$u^{(1)}$
was obtained by filtering through the band
$[0.5\omega ,1.5\omega ]$
.
To make best use of the vertical distribution of velocity resolved by SWASH, the following relationship:
\begin{equation}\frac{\cosh \left[2k_{p}(z+h)\right]}{\cosh ^{2} k_{p}h}\left[\left.u^{(1)}\right| _{z=0}\right]^{2}=\big[1+\tanh ^{2} k_{p}(z+h)\big]\big[u^{(1)}\big]^{2},\end{equation}
is used to change (2.64) into an alternative form
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\frac{\overline{\tilde{u}^{2}}}{c_{g}}+\zeta \left(k_{p}h,\frac{z}{h}\right)\big[1+\tanh ^{2} k_{p}(z+h)\big]\frac{\overline{\left[u^{(1)}\right]^{2}\tilde{u}}}{cc_{g}}.\end{equation}
Equation (3.5) is evaluated both theoretically and numerically for comparison with the result from the Lagrangian particle-tracking simulation. The numerical evaluation calculates
$\overline{\tilde{u}^{2}}$
and
$\overline{[u^{(1)}]^{2}\tilde{u}}$
first and then substitutes them into (3.5). The theoretical evaluation calculates
$\overline{\tilde{u}^{2}}$
and
$\overline{[u^{(1)}]^{2}\tilde{u}}$
as
\begin{equation}\overline{\tilde{u}^{2}}=\frac{1}{2}\left| \hat{u}_{l}^{2}\right| , \overline{\left[u^{(1)}\right]^{2}\!\tilde{u}}=\frac{1}{\sqrt{2}}\textrm{std}\big\{\widetilde{\big[u^{(1)}\big]^{2}}\big\}\left| \hat{u}_{l}\right| \cos \phi,\end{equation}
where
$| \hat{u}_{l}| \cos \phi$
(
$\phi$
= biphase) is theoretically predicted from (2.27) while
$\textrm{std}\{\widetilde{[u^{(1)}]^{2}}\}$
(i.e. the standard deviation of low-pass filtered
$[u^{(1)}]^{2}$
in time dimension) is calculated from simulated
$u^{(1)}$
.
3.1.2. Lagrangian particle tracking and subharmonic wave-induced mass transport velocity
The Lagrangian particle tracking was conducted using the velocity field by the SWASH model. Let
$x_{{L}}(t)$
and
$z_{{L}}(t)$
be the horizontal and vertical coordinates of a tagged particle,
$u_{{L}}(t)$
and
$w_{{L}}(t)$
be the horizontal and vertical velocities, the particle displacement in one time step from
$t$
to
$t+{\unicode[Arial]{x0394}} t$
was given by
\begin{align} \begin{cases} x_{{L}}\left(t+{\unicode[Arial]{x0394}} t\right)-x_{{L}}\left(t\right)=u_{{L}}\left(t\right){\unicode[Arial]{x0394}} t,\\[4pt] z_{{L}}\left(t+{\unicode[Arial]{x0394}} t\right)-z_{{L}}\left(t\right)=w_{{L}}\left(t\right){\unicode[Arial]{x0394}} t, \end{cases} \end{align}
where the velocities
$u_{{L}}(t)$
and
$w_{{L}}(t)$
were interpolated from the simultaneous SWASH flow field to the particle location
$x_{{L}}(t)$
and
$z_{{L}}(t)$
using the square inverse distance weighted interpolation (Shepard Reference Shepard1968).
To derive the mass transport velocity at a given location (x, z), a particle is released from (x, z) and tracked for one wave group cycle of 35.72 s. This procedure is repeated 40 times with the release timepoint being increased by 1/20 wave group period successively to establish an ensemble average. A linear regression with respect to the tracking time
$t$
is applied to the resultant 40 records of
$x_{{L}}(t)$
to obtain the mass transport velocity at (x, z).
3.2. Model results
3.2.1. Trajectories of Lagrangian particles
Figure 8(a) shows the numerical results of time-averaged total horizontal Eulerian velocity profiles over a 1/120 bed slope. Similar to the illustration in figure 1, the velocity is positive between the wave trough and crest, leading to a forward mass transport. Below the wave trough, a return current is formed and offsets the forward mass transport, resulting in a zero net mass transport in the numerical flume. The return current is skewed towards the surface and becomes more uniform vertically as the water depth decreases.
Figure 8. (a) Time-averaged total Eulerian horizontal velocity for the
$h_{x}=-1/120$
case simulated by the SWASH model. (b)–(d) Trajectories of fluid particles initially released at the timepoint of wave crest along the vertical profile of local depths (b)
$k_{0}h=0.33$
, (c)
$k_{0}h=0.25$
and (b)
$k_{0}h=0.17$
and tracked for one wave group cycle. The magenta dashed lines in (a) denote the wave crest and trough heights. Here,
$k_{0}=\omega ^{2}/g$
is the deep-water wavenumber. Green and red circles denote the starting and ending positions of the tagged particles, respectively.
Figure 8(b–d) shows the trajectories of particles released at various heights for three water depths shown in figure 8(a). Each particle was tracked over one wave group cycle, after which the particles near the surface moved backward while those near the bottom moved forward. This phenomenon results from the stronger return Eulerian current near the surface, as shown in figure 8(a), overpowering the forward Stokes drift near the surface, while the opposite is true near the bottom. At greater depths (figure 8 b), a small proportion of particles near the surface moved backward while, as depth decreases, more particles moved backwards as the return current becomes more uniformly distributed in shallow water (figure 8 d). Figure 8 also shows that particles move downward after a wave group cycle, more so at a shallow depth (figure 8 d) than at a deeper depth (figure 8 b). Near the bed, however, particles no longer move downward and can only move forward along the bed.
3.2.2. Validation of subharmonic wave-induced mass transport solution
Figure 9 compares the profile of the subharmonic wave-induced mass transport velocity obtained by using Lagrangian particle tracking (dots), the value of
$\overline{u_{{L}}^{(4)}}$
by (3.5) evaluated numerically (SWASH, solid lines) and theoretically (theory, dashed lines). Both predictions using numerical and theoretical velocities agree well with the particle-tracking results, validating the theoretical model of subharmonic wave-induced mass transport.
Figure 9. Profiles of mass transport velocity induced by subharmonic wave forced by bichromatic wave groups with amplitudes
$[a_{1},a_{2}]=[0.27,0.09]$
m and radian frequencies
$[\omega _{1},\omega _{2}]=[0.9676,0.7917]$
rad s−1 over bottom slopes of 1/30 (a, d), 1/60 (b, e) and 1/120 (c, f). Upper:
$kh=0.54$
; lower:
$kh=0.44$
. Black dots: difference between the mass transport velocity evaluated using Lagrangian particle-tracking model driven by SWASH flow field (cf. § 3.1.2) with and without the subharmonic velocity. Black, red and blue lines denote the subharmonic wave-induced mass transport velocity
$\overline{u_{{L}}^{(4)}}$
, the Term I (‘Stokes drift’ effect), and Term II (skewness effect) components predicted by (3.5). Equation (3.5) was evaluated using both simulated velocity (SWASH, solid lines) and theoretical subharmonic velocity (theory, dashed lines).
The total mass transport velocity is positive at the bottom and decreases upwards, approaching zero or even a negative value towards the surface. This is well captured by the theoretical profile of
$\overline{u_{{L}}^{(4)}}$
, especially when evaluated using the simulated profile of
$\overline{[u^{(1)}]^{2}\tilde{u}}$
. The theoretical profiles indicate that this is due to the greater negative Term II near the surface whereas profiles of Term I were almost vertically uniform, aligning with the theoretical assumption of shallow water for a subharmonic wave. Term I increases significantly in magnitude with decreasing bed slope and depth while Term II does not, resulting in the enhanced dominance of Term I over Term II over a mild slope and in shallow water, consistent with figure 3(b). Over the 1/30 bed slope at kh = 0.44 (figure 9
d), Term I marginally surpasses Term II while the advantage becomes apparent for 1/60 bed slope and (figure 9
e). As the bed slope decreases to 1/120, Term I increases significantly while Term II does not, making Term I predominant and hence close to the total velocity (figure 9
f).
3.2.3. Influences of bed slope on subharmonic wave-induced mass transport
Figure 10(a) presents the numerical results of
$\bar{M}^{(4)}$
(2.43) for bed slopes 1/30–1/120, showing that the mass transport increases with decreasing depth at a greater rate for gentler bed slopes. However, as the bed slope decreases from 1/90 to 1/120, the increase in mass transport becomes marginal. For bed slopes 1/30 and 1/60, the transport is negative (seaward) at deeper depths while it is positive near the end of the shoaling zone, consistent with the results derived from experimental data (figure 6). The results also indicate that, when the bed slope is not too mild (seemingly greater than 1/90 in this case), the transition from negative to positive transport occurs at deeper depths for milder bed slopes, which also aligns with figure 2(b).
Figure 10. (a) Subharmonic wave induced mass transport
$\bar{M}^{(4)}$
(40) evaluated with SWASH-simulated velocity. (b) Percentage ratio between
$\bar{M}^{(4)}$
and Stokes transport
$\overline{E}/c$
. Here,
$\rho =1000\,(\textrm{kg}\,\textrm{m}^{-3})$
was adopted.
Figure 11. Fluid particle trajectories over ten wave group cycles (fifty wave cycles). Upper: particles driven by flow field simulated by SWASH model initially released at 5 different heights at water depth (a)
$k_{0}h=0.27$
and (b)
$k_{0}h=0.17$
; lower: particles driven by flow field after removing the subharmonic orbital velocity. Releasing profiles in (c) and (d) are same as in (a) and (b), respectively. The red circles in (c)–(d) are duplicates of those in (a)–(b) for comparison.
$k_{0}=\omega ^{2}/g$
is the deep-water wavenumber. Here,
${\unicode[Arial]{x0394}} x$
denotes the horizontal displacement of particles. Green and red/yellow circles denote the starting and ending position of the tagged particles with/without subharmonic orbital velocity, respectively.
The percentage ratio of
$\bar{M}^{(4)}$
to the Stokes transport of short waves
$\overline{E}/c$
(figure 10
b) shows a similar trend to
$\bar{M}^{(4)}$
in figure 10(a), indicating little effect of the bed slope on the Stokes transport. Specifically, as the bed slope decreases from 1/30 to 1/90, the maximum subharmonic wave-induced mass transport at the shallow-water end increased from approximately 0.5 % to 2 % of the Stokes transport but remained almost unchanged as the bed slope decreases further from 1/90 to 1/120. It is also noticed that for bed slope 1/30 the negative
$\bar{M}^{(4)}$
around kh = 0.57 is nearly 0.5 % of the Stokes transport, which is comparable to the ratio of the maximum
$\bar{M}^{(4)}$
to the Stokes transport at the end of shoaling zone. Although the ratio between
$\bar{M}^{(4)}$
and the Stokes transport is generally low, below 2 %, according to (2.13) and (2.26),
$\tilde{u}$
generally scales with
$\tilde{E}/(hc)$
, thus this ratio relates to the short wave amplitude approximately as
\begin{equation}\frac{\bar{M}^{\left(4\right)}}{\overline{E}/c}\propto \frac{\overline{\left(\frac{\tilde{E}}{hc}\right)^{2}\frac{h}{c_{g}}}}{\overline{E}/c}\sim \frac{\left(a_{1}a_{2}/h\right)^{2}}{a_{1}^{2}+a_{2}^{2}}=\left(\frac{h^{2}}{a_{1}^{2}}+\frac{h^{2}}{a_{2}^{2}}\right)^{-1},\end{equation}
therefore, higher relative importance of
$\bar{M}^{(4)}$
is expected for larger waves. This is natural since
$\bar{M}^{(4)}$
scales with
$a^{4}$
while
$\overline{E}$
scales with
$a^{2}$
. More results are needed to exam this relationship.
3.2.4. Long-term effect of subharmonic wave-induced mass transport
Despite being secondary to Stokes transport of short waves, the effect of subharmonic wave-induced mass transport on particle movements accumulates with time and travel distance and becomes important in the long term. Figure 11(a–b) illustrates the trajectories of particles across 10 wave group cycles released at 5 heights at two water depths over a 1/120 bed slope. At the greater depth in figure 11(a), the particle movement maintained a similar trend as that over one wave group cycle (figure 8 c). At a shallow depth (figure 11 b), the particle movement near the bed maintained a similar trend along the bed as depicted in figure 8(d), while those released near the surface initially moved backward after one wave group cycle (figure 8 d) but eventually moved forward after ten wave group cycles. This is due to the marked downward movement of particles initially near the surface (figure 8 d), ultimately aligning them with the forward-moving trajectories of particles in the middle and lower water column.
Figure 11(c,d) demonstrates the same results as (figure 11 a,b) but without the subharmonic velocity. A comparison between figure 11(a) and c and the difference between the red and yellow circles in figure 11(c,d) reveal that the subharmonic orbital motion induces a weak net forward drift. This forward drift becomes more significant at shallow depths (figures 11 b and 11 d) especially for those released near the surface that drifted slowly at the leading order, which by contrast highlights the importance of the contribution of subharmonic wave-induced drift.
4. Conclusion and discussion
This study develops and validates a theoretical solution for the first time for mass transport as well as the associated drift velocity induced by group-forced subharmonic (infragravity) waves over an undulating bottom topography. The fourth-order solution for Lagrangian mass transport velocity incorporates the kinematic effects driven by both the subharmonic wave-induced velocity variance and skewness. The theory was validated by GLOBEX experiment and Lagrangian particle-tracking simulations driven by the SWASH model. Additionally, the analytical solution was applied to the experimental data from the GLOBEX campaign to demonstrate its applicability under field conditions.
The theory identifies two key kinematic mechanisms contributing to subharmonic wave-induced mass transport: (i) a forward drift analogous to Stokes drift, arising from the self-interaction of subharmonic orbital motion that is proportional to the variance (energy) of the subharmonic velocity, and (ii) a depth-decaying backward drift arising from the negative velocity skewness induced by the near anti-phase coupling between subharmonic waves and short wave groups. For relative bed slope |β| less than ∼0.01, these two mechanisms are comparable at around
$kh\approx 0.72$
, and this critical value decreases with increasing bed slope, wave period and decreasing bandwidth. In deeper water or over steeper bed slopes, the backward drift can exceed the forward drift near the surface, producing a net seaward transport in the upper water column. As water depth decreases, forward drift increases and dominates net positive mass transport, especially over mild slopes.
Although subharmonic-induced transport is typically smaller than the classical Stokes transport (e.g. reaching up to ∼2 % of Stokes transport at
$kh\approx 0.42$
over a 1/120 slope), its influence on material transport accumulates over time and distance and is vertically distributed across the entire water column. These characteristics, therefore, make subharmonic-induced transport particularly relevant in shallow coastal regions, especially where Stokes drift is nearly cancelled by the Eulerian return flow, leaving subharmonic contributions the key driver of particle motion. Based on the expression of mass transport obtained, it is evident that the relative importance of subharmonic-induced transport to Stokes transport increases with short wave amplitude.
The subharmonic wave-induced mass transport provides the new third mechanism for sediment transport. Despite the skewness-induced sediment transport (Deigaard et al. Reference Deigaard, Jakobsen and Fredsøe1999; Yu, Hsu & Hanes Reference Yu, Hsu and Hanes2010; Baldock et al. Reference Baldock, Alsina, Caceres, Vicinanza, Contestabile, Power and Sanchez-Arcilla2011) arising from group-modulated sediment concentration near the bottom, however, the present skewness-induced mass transport is a purely kinematic effect of subharmonic wave motion so that it affects the whole water column from surface to bottom, therefore, influences the movement of non-inertial particles at the wave surface (e.g. buoyant microplastics). The latter was previously thought to be affected solely by Stokes drift from short waves.
Acknowledgements
We thank Q. Liu from Tianjin University for checking the calculation of fourth-order Lagrangian velocity in the shallow-water limit. We are grateful to the anonymous reviewers whose input greatly improved the manuscript.
Funding
Z.L. and Q.Z. have been supported by the Natural Environment Research Council of the UK (Grant No. NE/V006088/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Subharmonic wave-induced mass transport in Eulerian system
1. Broad-banded waves with
$\alpha =O(1)$
then
\begin{equation} k_{m,n}h=O\left(\alpha k_{m}h\right)=O\left(k_{m}h\right) .\end{equation}
-
1.1. Deep water
$| k_{m,n}| h=O(k_{m}h)\gg 1$
In this case
$\tanh k_{m,n}h=\text{sign}(k_{m,n})[1+O(e^{-{k_{m}}h})]=\text{sign}(\omega _{m,n})[1+O(e^{-{k_{m}}h})]$
, equation (2.39) becomes
\begin{equation}\bar{M}^{\left(4\right)}=\sum _{m,n}\frac{1}{2}\frac{\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}}{\left| \omega _{m,n}\right| }+\frac{\max(k_{m}\omega _{m},k_{n}\omega _{n})}{\omega _{m}\omega _{n}\left| \omega _{m,n}\right| }\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n},\end{equation}
with an error factor
$1+O(e^{-{k_{m}}h})$
.
Invoking the deep-water dispersion relationship
\begin{equation}k_{m}=\frac{\omega _{m}^{2}}{g}\big[1+O\big(e^{-2{k_{m}}h}\big)\big],\end{equation}
equation (A2) reduces to
\begin{equation}\bar{M}^{\left(4\right)}=\sum _{m,n}\frac{1}{2}\frac{\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}}{\left| \omega _{m,n}\right| }+\frac{\max \left(\omega _{m}^{3},\omega _{n}^{3}\right)}{\omega _{m}\omega _{n}\left| \omega _{m,n}\right| }\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{g},\end{equation}
with an error factor
$[1+O(e^{-2{k_{m}}h})]$
.
-
1.2. Intermediate depth
$| k_{m,n}| h=O(k_{m}h)=O(1)$
In this case the generalised solution (2.39) must be invoked.
-
1.3. Shallow water
$| k_{m,n}| h=O(k_{m}h)\ll 1$
In this case, the dispersion relationship becomes
\begin{equation}\frac{\omega _{m}}{k_{m}}=\frac{\omega _{n}}{k_{n}}=\frac{\omega _{m,n}}{k_{m,n}}=\sqrt{gh}\big[1+O\big(k_{m}^{2}h^{2}\big)\big],\end{equation}
and (2.39) reduces to
\begin{equation}\bar{M}^{\left(4\right)}=\sqrt{\frac{h}{g}}\sum _{m,n}\frac{1}{2}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}+\frac{3}{g}\sum _{m,n}\frac{1}{2}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n},\end{equation}
with an error factor
$1+O(k_{m}^{2}h^{2})$
.
Furthermore, recall
$\tilde{u}$
and
$u^{(1)}$
in (2.33) and (2.4), respectively, we have the relationships
\begin{equation}\left.\overline{\tilde{u}^{2}}\right| _{z=0}=\overline{Re\left[\sum _{m,n}\widehat{U}_{m,n}e^{\textrm{i}{\theta _{m,n}}}\right]Re\left[\sum _{m,n}\widehat{U}_{m,n}e^{\textrm{i}{\theta _{m,n}}}\right]}=\frac{1}{2}\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*},\end{equation}
\begin{equation}\left.\overline{\tilde{u}\left[u^{(1)}\right]^{2}}\right| _{z=0}=\overline{Re\left[\sum _{m}\hat{u}_{m}e^{\textrm{i}{\theta _{\textrm{m}}}}\!\right]Re\left[\sum _{n}\hat{u}_{n}e^{\textrm{i}{\theta _{n}}}\!\right]Re\left[\sum _{m,n}\widehat{U}_{m,n}e^{\textrm{i}{\theta _{m,n}}}\!\right]}=\frac{1}{2}\sum _{m,n}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}.\end{equation}
And therefore (A6) can be recovered in the time domain as
\begin{equation}\bar{M}^{\left(4\right)}=\sqrt{\frac{h}{g}}\overline{\tilde{u}^{2}}+\frac{3}{g}\overline{\left[u^{(1)}\right]^{2}\tilde{u}},\end{equation}
where the velocities at arbitrary depth instead of at
$z=0$
are used with the same error factor
$1+O(k_{m}^{2}h^{2})$
.
2. Narrow-banded waves with
$\alpha \ll 1$
with central frequency/wavenumber
$\omega _{p}$
/
$k_{p}$
then
\begin{equation} \left| k_{m,n}\right| h=O\left(\alpha k_{m}h\right)=O\left(\alpha k_{p}h\right)\ll O\left(k_{p}h\right) .\end{equation}
-
2.1. Deep water for both short and subharmonic waves
$1\ll | k_{m,n}| h\ll O(k_{p}h)$
In this case, equation (A4) applies but with the coefficient of the second term approximated as
\begin{align}\frac{{\max (\omega _m^3,\omega _n^3)}}{{{\omega _m}{\omega _n}|{\omega _{m,n}}|}} = \frac{{{\omega _p}}}{{|{\omega _{m,n}}|}}\left[ {1 + O\left( \alpha \right)} \right]\!,\end{align}
hence it becomes
\begin{equation}\bar{M}^{\left(4\right)}=\sum _{m,n}\frac{1}{2}\frac{\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}}{\left| \omega _{m,n}\right| }+\frac{\omega _{p}}{\left| \omega _{m,n}\right| }\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{g}\end{equation}
with an error factor
$1+O(\alpha ,e^{-2| k_{m,n}| h})$
.
-
2.2. Deep water for short wave, intermediate depth for subharmonic waves
$| k_{m,n}| h=O(1)\ll O(k_{p}h)$
Let
$(k_{m}h,k_{n}h)=O(k_{p}h)\gg 1$
but keep
$k_{m,n}h=O(1)$
and invoke the deep-water dispersion relationship (A3), equation (2.39) becomes approximately
\begin{equation}\bar{M}^{\left(4\right)}=\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{\tanh k_{m,n}h}{2\omega _{m,n}}+\frac{\omega _{p}\tanh k_{m,n}h}{\omega _{m,n}}\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{g},\end{equation}
with an error factor
$1+O(\alpha ,k_{p}^{-1}h^{-1})$
for the second term.
-
2.3. Deep water for short wave, shallow water for subharmonic waves
$| k_{m,n}| h\ll 1\ll O(k_{p}h)$
Let
$|k_{m,n}|h\ll 1$
in (2.39), it becomes
\begin{equation}\bar{M}^{\left(4\right)}=\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{k_{m,n}h}{2\omega _{m,n}}+\frac{\omega _{p}k_{m,n}h}{\omega _{m,n}}\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{g},\end{equation}
with an error factor
$1+O(k_{m,n}^{3}h^{3})$
. For narrow-banded waves introduce the approximation
\begin{align} \left\{ \begin{array}{c} {k_n} = \dfrac{{{\omega _p}}}{c} + \dfrac{{{\omega _n} - {\omega _p}}}{{{c_g}}} + O\big[ {{{\big( {{\omega _n} - {\omega _p}} \big)}^2}} \big] \\[9pt] {k_m} = \dfrac{{{\omega _p}}}{c} + \dfrac{{{\omega _m} - {\omega _p}}}{{{c_g}}} + O\big[ {{{\big( {{\omega _m} - {\omega _p}} \big)}^2}} \big] \\ \end{array}\right. \end{align}
where
\begin{equation} c=\frac{\omega _{p}}{k_{p}}=\frac{g}{\omega _{p}}\big[1+O\big(e^{-{k_{p}}h}\big)\big], c_{g}=\left.\frac{\partial \omega }{\partial k}\right| _{\omega ={\omega _{p}}}=\frac{g}{2\omega _{p}}\big[1+O\big(e^{-{k_{p}}h}\big)\big] ,\end{equation}
are the deep-water wave phase and group speed corresponding to central frequency
$\omega _{p}$
, respectively. The coefficient of Term I in (A14) is thus approximately
\begin{equation} \frac{k_{m,n}}{\omega _{m,n}}=\frac{1}{c_{g}}\left[1+O\left(\alpha \right)\right]=\frac{2\omega _{p}}{g}\big[1+O\big(\alpha ,e^{-{k_{p}}h}\big)\big] .\end{equation}
Therefore, (A14) reduces to
\begin{equation}\bar{M}^{\left(4\right)}=\frac{\omega _{p}h}{g}\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}+\frac{2\omega _{p}^{2}h}{g^{2}}\sum _{m,n}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n},\end{equation}
with an error factor
$1+O(k_{m,n}^{3}h^{3},\alpha ,e^{-2{k_{p}}h})$
.
Like (A6), (A18) can be recovered in time domain as
\begin{equation}\bar{M}^{\left(4\right)}=\frac{2\omega _{p}}{g}\overline{\tilde{u}^{2}}h+\frac{4\omega _{p}^{2}}{g^{2}}\overline{\left[\left.u^{(1)}\right| _{z=0}\right]^{2}\tilde{u}}h.\end{equation}
-
2.4. Intermediate depth for short waves, shallow water for subharmonic waves
$| k_{m,n}| h\ll O(k_{p}h)=O(1)$
Let
$(k_{m}h,k_{n}h)=k_{p}h[1+O(\alpha )]=O(1)$
and
$k_{m,n}h\ll O(1)$
, equation (2.39) becomes
\begin{equation}\bar{M}^{\left(4\right)}=\frac{h}{2c_{g}}\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}+\left[\!\left(\frac{1}{cc_{g}}+\frac{1}{c^{2}}\!\right)\!\frac{\tanh k_{p}h}{k_{p}h}+\frac{1+\tanh ^{2} k_{p}h}{cc_{g}}\!\right]\!\sum _{m,n}\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}h}{2},\end{equation}
with an error factor
$[1+O(k_{m,n}^{2}h^{2},\alpha )]$
and
$c$
and
$c_{g}$
take the general expressions
\begin{align}c = \frac{{{\omega _p}}}{{{k_p}}} = \sqrt {\frac{{g\tanh {k_p}h}}{{{k_p}}}} ,\quad {c_g} = {\left. {\frac{{\partial \omega }}{{\partial k}}} \right|_{\omega = {\omega _p}}} = c\left( {\frac{1}{2} + \frac{{{k_p}h}}{{\sinh 2{k_p}h}}} \right).\end{align}
Equation (A20) can be recovered in time domain as
\begin{align}{{\bar M}^{\left( 4 \right)}} = \frac{{\overline {{{\tilde u}^2}} h}}{{{c_g}}} + \left( {\frac{3}{2}\frac{{\tanh {k_p}h}}{{{k_p}h}} + \frac{{2\cosh 2{k_p}h + 1}}{{2{{\cosh }^2}{k_p}h}}} \right)\frac{{\overline {{{[{u^{\left( 1 \right)}}{|_{z = 0}}]}^2}\tilde u} h}}{{c{c_g}}}.\end{align}
And for
$k_{p}h\gg 1$
(A22) reduces to (A19).
-
2.5. Shallow water for both short and subharmonic waves
$| k_{m,n}| h\ll O(k_{p}h)\ll 1$
In this case the result is given by equations (A6) and (A9), which can also be obtained for
$k_{p}h\ll 1$
in (A22).
Appendix B. Proof of equivalence between
$\boldsymbol{D}_{\boldsymbol{m},\boldsymbol{n}}$
in Eulerian and Lagrangian systems
The coefficient
$D_{m,n}$
obtained in Lagrangian system (2.63) writes as
\begin{align}{D_{m,n}} & = \left( {\frac{1}{{{\omega _m}}}\frac{{\tanh {k_m}h - \tanh {k_{m,n}}h}}{h} + \frac{1}{{{\omega _n}}}\frac{{\tanh {k_n}h + \tanh {k_{m,n}}h}}{h}} \right)\frac{{{k_{m,n}}}}{{{\omega _{m,n}}}} \nonumber\\[3pt] &\quad + \frac{{{k_m}{k_n}}}{{{\omega _m}{\omega _n}}}\left( {\frac{{{\omega _m}}}{{{\omega _{m,n}}}}\frac{{\tanh {k_n}h + \tanh {k_{m,n}}h}}{{{k_m}h}}} \right. \nonumber\\[3pt] &\quad \left. { - \frac{{{\omega _n}}}{{{\omega _{m,n}}}}\frac{{\tanh {k_m}h - \tanh {k_{m,n}}h}}{{{k_n}h}}} \right) \nonumber\\[3pt] &\quad + \frac{{{\omega _m}{k_m}\left( {\tanh {k_n}h + \tanh {k_{m,n}}h} \right) - {\omega _n}{k_n}\left( {\tanh {k_m}h - \tanh {k_{m,n}}h} \right)}}{{{\omega _m}{\omega _n}{\omega _{m,n}}h}} \nonumber\\[3pt] &\quad - \frac{{2k_{m,n}^2}}{{{\omega _m}{\omega _n}}}\frac{{\tanh {k_m}h - \tanh {k_n}h}}{{{k_{m,n}}h}}. \end{align}
Using the following relationship of the tanh function:
\begin{equation} \tanh k_{m,n}h-\tanh k_{m}h+\tanh k_{n}h=\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h ,\end{equation}
the first term in (B1) may be rewritten as
\begin{align}&\left(\frac{1}{\omega _{m}}\frac{\tanh k_{m}h-\tanh k_{m,n}h}{h}+\frac{1}{\omega _{n}}\frac{\tanh k_{n}h+\tanh k_{m,n}h}{h}\right)\frac{k_{m,n}}{\omega _{m,n}}\notag\\[3pt]&\quad =\frac{k_{m,n}}{\omega _{m,n}}\left(\frac{k_{m}}{\omega _{m}}\frac{\tanh k_{m}h}{k_{m}h}+\frac{k_{n}}{\omega _{n}}\frac{\tanh k_{n}h}{k_{n}h}\right)+\frac{k_{m,n}^{2}}{\omega _{m}\omega _{n}}\frac{\tanh k_{m,n}h}{k_{m,n}h}.\end{align}
The second and third terms can be combined as
\begin{align}&\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{n}h+\tanh k_{m,n}h}{k_{m}h}-\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{m}h-\tanh k_{m,n}h}{k_{n}h}\right)\notag\\[3pt]&\qquad +\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{n}h+\tanh k_{m,n}h}{k_{n}h}-\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{m}h-\tanh k_{m,n}h}{k_{m}h}\right)\notag\\[3pt]&\quad =\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{m}h+\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h}{k_{m}h}\right .\notag\\&\qquad\left. -\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{n}h-\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h}{k_{n}h}\right)\notag\\[3pt]&\qquad+\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{n}h+\tanh k_{m,n}h}{k_{n}h}-\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{m}h-\tanh k_{m,n}h}{k_{m}h}\right)\notag\\[3pt]&\quad=\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\tanh k_{m}h}{k_{m}h}+\frac{\tanh k_{n}h}{k_{n}h}+\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h}{k_{m}h}\right .\notag\\[3pt]&\qquad\left. +\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h}{k_{n}h}+\frac{\omega _{m}}{\omega _{m,n}}\frac{\tanh k_{n}h}{k_{n}h}+\frac{\omega _{n}}{\omega _{m,n}}\frac{\tanh k_{m}h}{k_{m}h}\right) \notag\\[3pt]&\quad =\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\tanh k_{m}h}{k_{m}h}+\frac{\tanh k_{n}h}{k_{n}h}\right .\notag\\[3pt]&\qquad\left. +\frac{1}{\omega _{m,n}}\frac{\left(\omega _{n}+\omega _{m,n}\right)\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h+\left(\omega _{m}-\omega _{m,n}\right)\tanh k_{m,n}h}{k_{m}h}\right .\notag\\[3pt]&\qquad\left. +\frac{1}{\omega _{m,n}}\frac{\left(\omega _{m}-\omega _{m,n}\right)\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h+\left(\omega _{n}+\omega _{m,n}\right)\tanh k_{m,n}h}{k_{n}h}\right) \notag\\[3pt]&\quad =\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\tanh k_{m}h}{k_{m}h}+\frac{\tanh k_{n}h}{k_{n}h}\right)+\left(\frac{\omega _{m}}{k_{m}}+\frac{\omega _{n}}{k_{n}}\right)\frac{k_{m,n}}{\omega _{m,n}}\frac{\tanh k_{m,n}h}{k_{m,n}h}\notag\\[3pt]&\qquad +\left(\frac{\omega _{n}}{k_{m}}+\frac{\omega _{m}}{k_{n}}\right)\frac{k_{m,n}}{\omega _{m,n}}\frac{\tanh k_{m,n}h\tanh k_{m}h\tanh k_{n}h}{k_{m,n}h}+\frac{k_{m,n}^{2}}{k_{m}k_{n}}\frac{\tanh k_{m}h-\tanh k_{n}h}{k_{m,n}h}.\end{align}
Substituting (B3) back into (B1) and using the relationship
\begin{equation} \frac{k_{m,n}^{2}}{\omega _{m}\omega _{n}}\frac{\tanh k_{m,n}h}{k_{m,n}h}-\frac{k_{m,n}^{2}}{\omega _{m}\omega _{n}}\frac{\tanh k_{m}h-\tanh k_{n}h}{k_{m,n}h}=\frac{k_{m,n}^{2}}{\omega _{m}\omega _{n}}\frac{\tanh k_{m}h\tanh k_{n}h\tanh k_{m,n}h}{k_{m,n}h} ,\end{equation}
equation (B1) becomes
\begin{align}D_{m,n}&=\frac{k_{m,n}}{\omega _{m,n}}\left(\frac{\tanh k_{m}h}{\omega _{m}h}+\frac{\tanh k_{n}h}{\omega _{n}h}\right)+\frac{k_{m}k_{n}}{\omega _{m}\omega _{n}}\left(\frac{\tanh k_{m}h}{k_{m}h}+\frac{\tanh k_{n}h}{k_{n}h}\right)\notag\\&\quad+\left(\frac{\omega _{m}}{k_{m}}+\frac{\omega _{n}}{k_{n}}\right)\frac{k_{m,n}}{\omega _{m,n}}\frac{\tanh k_{m,n}h}{k_{m,n}h}\notag\\&\quad+\left[\left(\frac{k_{n}}{\omega _{m}}+\frac{k_{m}}{\omega _{n}}\right)\frac{k_{m,n}}{\omega _{m,n}}+\frac{k_{m,n}^{2}}{{\unicode[Arial]{x03C9}} _{m}\omega _{n}}\right]\frac{\tanh k_{m}h\tanh k_{n}h\tanh k_{m,n}h}{k_{m,n}h},\end{align}
which is same as that obtained in Eulerian system (2.40).
Appendix C. Subharmonic wave-induced Lagrangian mass transport velocity
-
1. Broad-banded waves with
$\alpha =O(1)$
then(C1)
\begin{equation} k_{m,n}h=O\left(\alpha k_{m}h\right)=O\left(k_{m}h\right) .\end{equation}
-
1.1. Deep water
$| k_{m,n}| h=O(k_{m}h)\gg 1$
In this case we have the following approximations:
\begin{equation} \cosh \left[2k_{m}(z+h)\right]=\frac{1}{2}e^{2{k_{m}}(z+h)}\big[1+O\big(e^{-4{k_{m}}h}\big)\big] ,\end{equation}
\begin{equation} \cosh \left[2k_{m,n}(z+h)\right]=\frac{1}{2}e^{2\left| k_{m,n}\right| (z+h)}\big[1+O\big(e^{-4{k_{m}}h}\big)\big] ,\end{equation}
\begin{equation} \cosh ^{2} k_{m,n}h=\frac{1}{4}e^{2\left| k_{m,n}\right| h}\big[1+O\big(e^{-2{k_{m}}h}\big)\big] ,\end{equation}
\begin{equation} \cosh k_{m}h\cosh k_{n}h\cosh k_{m,n}h=\frac{1}{8}e^{2\ max(k_{m},k_{n})h}\big[1+O\big(e^{-2{k_{m}}h}\big)\big] ,\end{equation}
and equation (2.61) becomes approximately
\begin{align}\overline{u_{L}^{\left(4\right)}}&=\sum _{m,n}\frac{k_{m,n}}{\omega _{m,n}}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}e^{2\left| k_{m,n}\right| z}\notag\\&\quad +\frac{k_{m,n}\max \left(\omega _{m}k_{m},\omega _{n}k_{n}\right)+k_{m}k_{n}\max \left(\omega _{m},\omega _{n}\right)+\max \left(\omega _{m}k_{m}^{2},\omega _{n}k_{n}^{2}\right)}{\omega _{m}\omega _{n}\left| \omega _{m,n}\right| }\notag\\&\qquad\times \hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}e^{2\max \left(k_{m},k_{n}\right)z}\notag\\& =\sum _{m,n}\frac{k_{m,n}}{\omega _{m,n}}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}e^{2\left| k_{m,n}\right| z}+\frac{2\max \left(k_{m}^{2}\omega _{m},k_{n}^{2}\omega _{n}\right)}{\omega _{m}\omega _{n}\left| \omega _{m,n}\right| }\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}e^{2\max \left(k_{m},k_{n}\right)z},\end{align}
with an error factor
$1+O(e^{-2{k_{m}}h})$
. Invoking the deep-water dispersion relationship (A3), equation (C6) becomes
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\sum _{m,n}\frac{\omega _{m}+\omega _{n}}{g}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}e^{\frac{2\left| \omega _{m}^{2}-\omega _{n}^{2}\right| z}{g}}+\frac{2\max \left(\omega _{m}^{5},\omega _{n}^{5}\right)}{\omega _{m}\omega _{n}\left| \omega _{m,n}\right| }\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{g^{2}}e^{2\max \left(\omega _{m}^{2},\omega _{n}^{2}\right)\frac{z}{g}},\end{equation}
with an error factor
$[1+O(e^{-2{k_{m}}h})]$
. Integrating (C7) from
$z=-\infty$
to
$z=0$
recovers (A4).
-
1.2. Intermediate depth
$| k_{m,n}| h=O(k_{m}h)=O(1)$
In this case the generalised solution (2.61) must be invoked.
-
1.3. Shallow water
$| k_{m,n}| h=O(k_{m}h)\ll 1$
In this case we substitute the shallow-water dispersion relationship (A5)
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\frac{1}{\sqrt{gh}}\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}+\frac{3}{gh}\sum _{m,n}{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}},\end{equation}
which can be recovered in time domain as
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\frac{\overline{\tilde{u}^{2}}}{\sqrt{gh}}+\frac{3}{gh}\overline{\left[u^{(1)}\right]^{2}\tilde{u}}.\end{equation}
Multiplying (C9) with h yields (A9).
-
2. Narrow-banded waves with
$\alpha \ll 1$
then(C10)
\begin{equation} \big|k_{m,n}\big| h=O\big(\alpha k_{m}h\big)=O\big(\alpha k_{p}h\big)\ll O\big(k_{p}h\big) .\end{equation}
-
2.1. Deep water for both short and subharmonic waves
$1\ll | k_{m,n}| h\ll O(k_{p}h)$
In this case equation (C7) applies but with
$\omega _{m}$
and
$\omega _{n}$
replaced by
$\omega _{p}[1+O(\alpha )]$
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\sum _{m,n}\frac{2\omega _{p}}{g}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}e^{\frac{4\omega _{p}\left| \omega _{m,n}\right| }{g}z}+\frac{2\omega _{p}^{3}}{\left| \omega _{m,n}\right| }\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{g^{2}}e^{\frac{2\omega _{p}^{2}}{g}z},\end{equation}
with an error factor
$[1+O(\alpha ,e^{-2| k_{m,n}| h})]$
. Integrating (C11) vertically from
$z=-\infty$
to
$z=0$
recovers (A12).
-
2.2. Deep water for short wave, intermediate depth for subharmonic waves
$| k_{m,n}| h=O(1)\ll O(k_{p}h)$
Let
$(k_{m}h,k_{n}h)=O(k_{p}h)\gg 1$
but keep
$k_{m,n}h=O(1)$
and invoke the deep-water dispersion relationship (A3), then equation (2.61) becomes approximately
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\sum _{m,n}\frac{\omega _{p}}{g}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}\frac{\cosh \left[2k_{m,n}(z+h)\right]}{\cosh ^{2} k_{m,n}h}+\frac{2\omega _{p}^{3}}{g^{2}}\frac{\tanh k_{m,n}h}{\omega _{m,n}}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}e^{2{k_{p}}z},\end{equation}
with an error factor
$[1+O(\alpha ,k_{p}^{-1}h^{-1})]$
. For
$| k_{m,n}| h\gg 1$
, (C12) reduces to (C11).
-
2.3. Deep water for short wave, shallow water for subharmonic waves
$| k_{m,n}| h\ll 1\ll O(k_{p}h)$
In this case the following result can be obtained either from (2.61) similar to the derivation of (A14) or considering
$| k_{m,n}| h\ll 1$
in (C12)
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\frac{\omega _{p}}{g}\sum _{m,n}\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}+\frac{4\omega _{p}^{4}h}{g^{3}}\sum _{m,n}\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}e^{\frac{2\omega _{p}^{2}z}{g}}.\end{equation}
which becomes in time domain
\begin{equation}\overline{u_{L}^{\left(4\right)}}=\frac{2\omega _{p}}{g}\overline{\tilde{u}^{2}}+\frac{8\omega _{p}^{4}h}{g^{3}}\overline{\left[u^{(1)}\right]^{2}\tilde{u}},\end{equation}
with an error factor
$[1+O(\alpha ,k_{p}^{-1}h^{-1},k_{m,n}^{2}h^{2})]$
.
-
2.4. Intermediate depth for short waves, shallow water for subharmonic waves
$| k_{m,n}| h\ll O(k_{p}h)=O(1)$
Let
$(k_{m}h,k_{n}h)=k_{p}h[1+O(\alpha )]=O(1)$
and
$k_{m,n}h\ll O(1)$
, equation (2.61) becomes approximately
\begin{align}\overline{u_{L}^{\left(4\right)}}&=\sum _{m,n}\frac{\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}}{2c_{g}}+\left(\frac{2\cosh 2k_{p}(z+h)}{cc_{g}}+\frac{1}{c^{2}}\frac{\partial }{\partial \omega }\left\{\omega \cosh \left[2k(z+h)\right]\right\}\right .\notag\\&\quad\left . +\frac{1}{\omega _{p}^{2}}\frac{\partial }{\partial \omega }\left\{\omega k^{2}\cosh \left[2k(z+h)\right]\right\}\right)\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{4{\cosh ^{2} k_{p}h}}h\notag\\&=\sum _{m,n}\frac{\widehat{U}_{m,n}\widehat{U}_{m,n}^{*}}{2c_{g}}+\left(\frac{1}{cc_{g}}\left[2+2k_{p}(z+h)\tanh \left[2k_{p}(z+h)\right]\right]+\frac{1}{c^{2}}\right)\notag\\&\quad\times\frac{\hat{u}_{m}^{*}\hat{u}_{n}\widehat{U}_{m,n}}{2}\frac{\cosh 2k_{p}(z+h)}{\cosh ^{2} k_{p}h},\end{align}
with an error factor
$[1+O(k_{m,n}^{2}h^{2},\alpha )]$
and
$c$
and
$c_{g}$
take the general expressions as in (A16). Equivalently in time domain (C15) becomes
\begin{align}\overline{u_{L}^{\left(4\right)}}&=\frac{\overline{\tilde{u}^{2}}}{c_{g}}+\left[\frac{5}{2}+\frac{k_{p}h}{\sinh 2k_{p}h}+2k_{p}(z+h)\tanh \left[2k_{p}(z+h)\right]\right]\notag\\&\quad\times\frac{\cosh 2k_{p}(z+h)}{\cosh ^{2} k_{p}h}\frac{\overline{\left[\left.u^{(1)}\right| _{z=0}\right]^{2}\tilde{u}}}{cc_{g}}.\end{align}
For
$k_{p}h\gg 1$
(C16) reduces to (C14). Integrating (C16) vertically from
$z=-h$
to
$z=0$
recovers (A22).
-
2.5. Shallow water for both short and subharmonic waves
$| k_{m,n}| h\ll O(k_{p}h)\ll 1$
In this case the result is given by equations (C8)–(C9), which can also be obtained for
$k_{p}h\ll 1$
in (C16).
Appendix D. Separation of incident and reflected subharmonic waves
The subharmonic surface elevation was decomposed into the incident (
$\tilde{\eta }^{+}$
) and reflected (
$\tilde{\eta }^{-}$
) components in the time domain using the weakly nonlinear method proposed by van Dongeren (Reference van Dongeren1997) as per Ruju, Lara & Losada (Reference Ruju, Lara and Losada2012) and Liu et al. (Reference Liu, Yao, Liao, Li, Zhang and Zou2023)
\begin{equation}\tilde{\eta }^{+}=\frac{\sqrt{gh}\tilde{\eta }+\tilde{Q}}{c_{g}+\sqrt{gh}},\tilde{\eta }^{-}=\frac{c_{g}\tilde{\eta }-\tilde{Q}}{c_{g}+\sqrt{gh}},\end{equation}
where
$c_{g}$
is the group velocity corresponding to the peak frequency;
$\tilde{Q}$
is the low-frequency volume flux per unit obtained by filtering the volume flux Q through in the infragravity frequency band
$[0,0.5/T_{p}]$
. The volume flux Q was estimated using the surface elevation h(x, t) and the horizontal velocity u(x, z, t) measured at a fixed height above the bed as follows.
The volume flux per unit width Q is defined as
\begin{equation}Q\left(t\right)=\int _{-h}^{\eta }u\left(z,t\right){\textrm d}z,\end{equation}
where we assume the velocity was measured at a fixed height of the probe
$z=z_{{p}}$
, which can be expressed as the discrete Fourier series
\begin{equation}u(z_{{p}},t)=\sum _{f=0}^{1.5f_{p}}\hat{u}\left(z_{{p}},f\right)\exp \left(-2\pi ift\right)\!,\end{equation}
where
$\hat{u}$
denotes the complex amplitude. The upper limit of frequency
$1.5f_{p}$
is used here to remove the contributions from the flume resonance and superharmonics.
The profile of
$\hat{u}(z,f)$
was then calculated from the single point value of
$\hat{u}(z_{{p}},f)$
using the linear wave theory as
$\hat{u}(z,f)=\hat{u}(z_{{p}},f)\cosh [k(h+z)]/\cosh [k(h+z_{{p}})]$
, with k being the wavenumber at frequency f. Therefore, the horizontal velocity in the time domain is given by
\begin{equation}u(z,t)=\sum _{f=0}^{1.5f_{p}}\hat{u}\left(z_{{p}},f\right)\frac{\cosh \left[k\left(h+z\right)\right]}{\cosh \left[k\left(h+z_{{p}}\right)\right]}\exp \left(-2\pi ift\right).\end{equation}
Substituting (D4) into (D2) yields
\begin{equation}Q(t)=\sum _{f=0}^{1.5f_{p}}\frac{\hat{u}\left(z_{{p}},f\right)}{\cosh \left[k\left(h+z_{{p}}\right)\right]}\frac{\sinh \left[k\left(h+\eta \right)\right]}{k}\exp \left(-2\pi ift\right).\end{equation}
Accurate to the second order in wave steepness, the Taylor expansion
$\sinh [k(h+\eta )]=\sinh kh+k\eta \cosh kh+O(\eta ^{2})$
is adopted to further reduce (D5) to
\begin{equation} Q(t)=\!\sum _{f=0}^{1.5f_{p}}\frac{\hat{u}\left(z_{{p}},f\right)\sinh kh}{k\cosh \left[k\left(h+z_{{p}}\right)\right]}\exp \left(-2{\unicode[Arial]{x03C0}} ift\right)+\eta\! \sum _{f=0}^{1.5f_{p}}\frac{\hat{u}\left(z_{{p}},f\right)\cosh kh}{\cosh \left[k\left(h+z_{{p}}\right)\right]}\exp \left(-2{\unicode[Arial]{x03C0}} ift\right)\! .\end{equation}
The orbital velocity of incident subharmonic wave is obtained by subtracting the reflected component from the total velocity as
\begin{equation} \tilde{u}^{+}=\tilde{u}+\tilde{\eta }^{-}\sqrt{\frac{g}{h}} .\end{equation}
Open access
