A nonlinear, non-dispersive energy balance for surfzone waves: infragravity wave dynamics on a sloping beach

Abstract

A fully nonlinear non-dispersive energy balance for surfzone waves is derived based on the nonlinear shallow water equations to study the nearshore dynamics of infragravity (IG) waves. Based on simulations of waves on a relatively moderate and mild beach slope with a non-hydrostatic wave-flow model (SWASH), the new theory shows that spatial gradients in IG energy flux are nearly completely balanced by the combined effect of bottom stresses and predominantly nonlinear triad interactions. The new balance confirms many features of existing weakly nonlinear theories, and yields an improved description in the inner surfzone where waves become highly nonlinear. A gain of IG energy flux throughout the shoaling and outer surfzones is driven by triad interactions between IG waves and pairs of sea-swell (SS) waves. The IG energy flux decreased in the inner surfzone, primarily through an energy cascade to the swell-band and superharmonic frequencies where wave energy is ultimately dissipated. Dissipation by bottom friction was weak on both slopes. The IG wave breaking, characterized by triads between three IG or two IG waves and one SS wave, was significant only deep inside the surfzone of the mild slope. Even though IG waves broke on the mild slope, nonlinear interactions between IG waves and pairs of SS waves were responsible for at least half of the net IG flux loss.

1. Introduction

Linear ocean surface waves shoaling from deep to shallow water refract towards shore and steepen (e.g. wavelength decreases and amplitude increases (Lamb 1932)). Nonlinear processes become increasingly important in shallow water, where waves pitch forward and break (Stoker 1957). During shoaling, nonlinear interactions cause a significant transfer of energy from the spectral peak frequency towards superharmonic and subharmonic frequencies (Phillips 1960; Hasselmann 1962; Longuet-Higgins & Stewart 1962; Hasselmann et al. 1976). These triad interactions are off-resonant in deep and intermediate depths, $kd\ge 1$ (where $k$ is the wave number and $d$ the local water depth), resulting in relatively small energy transfers. In shallow water, nonlinear triads approach resonance and energy transfer towards superharmonics causes sea-swell (SS) waves to become skewed and asymmetric (e.g. Elgar & Guza 1985; Elgar et al. 1997), until they eventually break in the surfzone. Energy transfers to subharmonics excite so-called infragravity (IG) waves. The IG wave height is generally $O$(cm) in deep water (e.g. Webb, Zhang & Crawford 1991; Aucan & Ardhuin 2013) but can reach $O$(m) in shallow water during storms (e.g. Sheremet et al. 2014; Matsuba, Shimozono & Sato 2020).

The nearshore dynamics of IG waves have been widely studied during the past decades through theoretical, laboratory, field and numerical efforts (see Bertin et al. (2018) for a recent review). Various theories have been developed to explain the substantial growth of IG waves in the nearshore (e.g. Symonds, Huntley & Bowen 1982; Schäffer 1993; Janssen, Battjes & Van Dongeren 2003; Nielsen & Baldock 2010; Contardo et al. 2021; Liao et al. 2021, and many others). Such theoretical background combined with numerical modelling (e.g. Reniers et al. 2002; Van Dongeren et al. 2003; Lara, Ruju & Losada 2011), laboratory experiments (e.g. Boers 1996; Baldock et al. 2000; Baldock & Huntley 2002) and field campaigns (e.g. Herbers, Elgar & Guza 1994; Herbers et al. 1995b; Herbers, Elgar & Guza 1995a; Okihiro, Guza & Seymour 1992), has significantly advanced our understanding of the IG wave patterns, growth rates, dissipation and phase relationship with SS wave groups (e.g. Battjes et al. 2004; Van Dongeren et al. 2007; Baldock 2012; de Bakker et al. 2015). Several studies used an IG wave energy balance to understand and quantify the growth and decay of IG waves in the nearshore (e.g. Henderson & Bowen 2002; Henderson et al. 2006; Ruju, Lara & Losada 2012; de Bakker, Tissier & Ruessink 2016). Such theories have generally been able to explain the IG dynamics in the shoaling region and outer surfzone, but failed to provide a complete description of the dynamics in the highly nonlinear inner surfzone due to inherent assumptions of weak wave nonlinearity. In this paper, we derive a new fully nonlinear IG energy balance that aims to describe the IG wave dynamics throughout the surfzone up to the mean waterline.

Nonlinear nearshore wave dynamics have often been studied using a wave energy balance that assumes cross-shore propagation of normally incident waves (e.g. Phillips 1970),

(1.1)\begin{equation} \frac{\partial E_f}{\partial t} + \frac{\partial F_f}{\partial x} = S_f^{NL} + D_f, \end{equation}

where $t$ is time, $x$ the cross-shore coordinate, $f$ frequency, $E_f=E(x,f,t)$ the energy density spectrum, $F_f=F(x,f,t)$ the energy flux, $S_f^{NL}=S(x,f,t)$ represents the nonlinear interactions that conservatively distribute energy over frequencies, and $D_f=D(x,f,t)$ accounts for the dissipation of wave energy (e.g. due to wave breaking and bottom friction). Terms in the wave energy balance have often been evaluated using a perturbation approach where wave nonlinearity ($\delta ={a}/{d}$, where $a$ is the wave amplitude) and dispersive effects ($\mu ={d}/{L}$, where $L$ is the wavelength) are accounted for up to some order (e.g. Freilich & Guza 1984; Agnon & Sheremet 1997; Herbers & Burton 1997). For example, Herbers, Russnogle & Elgar (2000) studied the spectral energy balance of breaking SS waves by evaluating the flux from linear theory $F_f=\rho g c_g E_f$ (with $c_g$ the wave group velocity, $g$ the gravitational acceleration and $\rho$ the water density) and with a nonlinear interaction term $S^{NL}$ from classical Boussinesq scaling, $\delta =O(\mu ^2)$. A perturbation approach necessarily truncates solutions at finite order and is formally invalid (no series convergence) when $\delta =O(1)$. Furthermore, energy balances (or stochastic wave models) derived using the perturbation approach depend on an infinite hierarchy of higher-order moments (bispectrum, trispectrum, etc.) without natural closure (e.g. Smit & Janssen 2016).

Nearshore IG wave dynamics have been examined with various approximate wave energy balances. Many studies assumed small nonlinear contributions, and estimated the IG flux from linear theory (e.g. Thomson et al. 2006; Van Dongeren et al. 2007; Torres-Freyermuth, Lara & Losada 2010; de Bakker et al. 2015; Liao et al. 2021). A subset of nonlinear contributions to the energy flux can be included by using a fully nonlinear low-frequency energy balance (e.g. Phillips 1970; Schäffer 1993), which assumes that IG frequencies are much lower than SS frequencies. Nonlinear contributions from correlations between three IG wave components are included, but (typically stronger) correlations between combinations of IG and SS components are not accounted for. Henderson et al. (2006) derived a weakly nonlinear energy balance that includes correlations between one IG and two SS components, the interaction that drives the off-resonant ‘bound’ IG wave observed well seaward of the surfzone (Longuet-Higgins & Stewart 1962). These triads become near-resonant during SS shoaling and breaking, and have been shown to explain much of the observed cross-shore variation of IG energy flux (e.g. Henderson et al. 2006; Ruju et al. 2012; Guedes, Bryan & Coco 2013; Rijnsdorp, Ruessink & Zijlema 2015; Mendes et al. 2018). de Bakker et al. (2015, 2016) used $S^{NL}$ from Boussinesq scaling to show that nonlinear interactions between two IG and one SS component can become significant and must be included to explain the loss of IG flux near the shoreline where the IG and SS wave heights are similar. Interactions between three IG components have been detected on mild sloping beaches, and have been associated with IG waves pitching forward and breaking (Van Dongeren et al. 2007; de Bakker et al. 2015, 2016).

Weakly nonlinear energy balances have significantly improved our understanding of nearshore IG wave dynamics but typically break down (i.e. the balance does not close) in the strongly nonlinear inner surfzone (e.g. Ruju et al. 2012; Guedes et al. 2013; Rijnsdorp et al. 2015; Mendes et al. 2018). Neglected or inaccurately estimated nonlinear terms distort our understanding of the dissipation of wave energy $D_f$ when it is estimated as the residual of the energy balance (1.1). Here, we consider a fully nonlinear wave energy balance for shallow water, where $a/d= \delta =O(1)$ and dispersive effects are cooperatively small $O(\mu ) \ll O(\delta )$ (i.e. Ursell numbers $Ur= \delta /\mu ^2 \gg 1$). We derive a frequency-resolved energy balance following the methodology of Henderson & Bowen (2002) in § 2 based on the nonlinear shallow water equations (NLSWE) that are expected to well describe the macroproperties of nonlinear and non-dispersive wave motions. Using simulations with a fully nonlinear and dispersive wave model (SWASH, § 3) of irregular waves propagating over a moderate (1/30) and mild (1/100) planar slope (§ 4), we use the energy balance to study nearshore IG wave dynamics (§ 5). This includes a description of the spatial variations of the nonlinear interactions contributing to the IG energy balance, and an analysis of the terms responsible for nearshore IG flux losses. In § 6, the new NLSWE energy balance is compared with existing theories, the role of IG wave breaking is described, and limitations of the non-dispersive energy balance for SS waves are discussed. Results are summarized in § 7.

2. Theory

2.1. Bulk energy balance

We consider shallow water waves in a horizontal ($x,y$) plane with a single-valued free-surface $z=\eta (x,y,t)$, bottom $z=-d(x,y)$, with components of the depth-averaged flow $u_m$, and continuous space and time derivatives. Assuming a very small characteristic depth over wavelength $\mu$ (or equivalently, assuming negligible vertical accelerations) yields the two-dimensional NLSWE,

(2.1)\begin{gather} \frac{\partial \eta}{\partial t} + \frac{\partial D u_n}{\partial x_n} = 0, \end{gather}

(2.2)\begin{gather}\frac{\partial u_m}{\partial t} + u_n\frac{\partial u_m}{\partial x_n} + g\frac{\partial \eta}{\partial x_m} ={-}\frac{\tau_m}{D}. \end{gather}

The Einstein summation convention is used, with the instantaneous water depth $D=d+\eta$, gravity acceleration $g$ and seabed stress $\tau _m$. A momentum balance in conservative form is retrieved by adding (2.2) $\times D$ and (2.1) $\times u_m$,

(2.3)\begin{equation} \frac{\partial Q_m}{\partial t} + \frac{\partial }{\partial x_n} \left( Du_n u_m + \frac{1}{2}gD^2 \right) = gD\frac{\partial d}{\partial x_m} - \tau_m, \end{equation}

with $Q_m=Du_m$.

An average bulk energy balance follows from multiplying (2.3) with $u_m/2$; (2.2) with $Du_m/2$; (2.1) with $g\eta$, summing, and taking the expected value,

(2.4)\begin{equation} \frac{\partial E}{\partial t} + \frac{\partial F_n}{\partial x_n} = S^\tau. \end{equation}

The energy density $E$, energy flux $F$ and a frictional dissipation source $S^\tau$ are defined as

(2.5)\begin{gather} E = \tfrac{1}{2}g \langle\eta^2\rangle + \tfrac{1}{2}\left\langle Du_m u_m \right\rangle, \end{gather}

(2.6)\begin{gather}F_n = \left\langle Q_n \left[ \frac{1}{2}u_m u_m + g\eta \right]\right\rangle, \end{gather}

(2.7)\begin{gather}S^\tau ={-}\left\langle\tau_m u_m\right\rangle, \end{gather}

where $\langle \ldots \rangle$ denotes a moving average over the fast time scale of the longest considered wave motions. With $S^\tau$ neglected and a smooth solution, $E$ is conserved to leading order. For stationary one-dimensional conditions, ${\partial E}/{\partial t}= {\partial F}/{\partial x}={\partial Q}/{\partial x}=0$, the conservation of hydraulic head is recovered,

(2.8)\begin{equation} \frac{u^2}{2} + g\eta = \text{constant}. \end{equation}

Equation (2.3) lacks mechanisms for wave breaking (e.g. wave overturning and air-entertainment), and for transferring organized energy to turbulence and heat. With breaking waves, the NLSWE-based energy balance does not close and the balance residual is interpreted as the dissipation rate of organized energy.

2.2. Frequency-resolved energy balance

A frequency-resolved energy balance details energy exchanges between different frequencies and identifies the frequency bands with the largest residual (dissipation). Following Henderson & Bowen (2002), a variable $X(t)$ is Fourier decomposed over frequencies $f$ spaced $\Delta f$ apart,

(2.9)\begin{equation} X(t) = \sum_f \bar{X}^f {\rm e}^{\mathrm{i} \omega t}, \end{equation}

where $\omega =2{\rm \pi} f$, $\bar {X}^f(x,y,t)$ denotes the complex amplitude of the frequency decomposition. To ensure real valued functions $(\bar {X}^f)^* = \bar {X}^{-f}$, where ${*}$ indicates the complex conjugate. The temporal variation of spectral amplitudes is associated with slow scale changes in mean energy (on scale $T$), whereas the complex exponential term accounts for oscillatory behaviour (with scale $f^{-1}$). The lowest considered frequency $f_0$ in principal separates wave-like dynamics from those associated with ‘very low frequency’ (VLF) vortical flows that are not included. A meaningful decomposition requires that means change on a much slower time scale $T$ than the longest wave frequencies considered (i.e. $f_0^{-1} T \ll 1$). The assumptions of scale separation and negligible vortical flow in the IG band might be problematic in the strongly nonlinear inner surfzone, but relaxing these assumptions is beyond the present scope.

Frequency balances follow from substitution of (2.9) for all time-dependent variables into (2.1)–(2.3), collecting terms at like frequencies, and multiplication by ${\rm e}^{-\mathrm {i} \omega t}$:

(2.10)\begin{gather} \frac{\partial\bar{\eta}^f}{\partial t} + \mathrm{i}\omega \bar{\eta}^f +\frac{\partial \overline{{D u_m}}^f }{\partial x_m} = 0, \end{gather}

(2.11)\begin{gather}\frac{\partial \overline{u_m}^f}{\partial t} + \mathrm{i}\overline{u_m}^f + \overline{ u_n \frac{\partial u_m}{\partial x_n}}^f + g\frac{\partial \bar{\eta}^f }{\partial x_m} ={-}\overline{\left(\frac{\tau_m}{D}\right)}^f\!, \end{gather}

(2.12)\begin{gather}\frac{\partial \overline{D u_m}^f}{\partial t} + \mathrm{i}\omega \overline{D u_m}^f + \frac{\partial \overline{D u_n u_m}^f}{\partial x_n} + g \overline{D \frac{\partial \eta}{\partial x_m} }^f ={-}\overline{\tau_m}^f. \end{gather}

A balance for the temporal mean of the kinetic energy density follows from combining (2.11) and (2.12),

(2.13)\begin{equation} \left( \frac{1}{4}\overline{D u_m}^f \times \mathrm{Eq.} (2.11)^* +\mathrm{C.C.} \right) + \left( \frac{1}{4}(\overline{u_m}^f)^* \times \mathrm{Eq.}(2.12) + \mathrm{C.C.} \right).\end{equation}

Dividing by $\Delta f$, considering the limit $\Delta f \rightarrow 0$, taking expected values, and using the definition of the cospectrum,

(2.14)\begin{equation} {C}_f \left(X;Y\right) = \lim_{\Delta f \rightarrow 0} {Re}\left\{\left\langle\frac{ \bar{X}^f \left(\bar{Y}^f\right)^*}{\Delta f}\right\rangle\right\}, \end{equation}

results in the following kinetic energy balance:

(2.15)\begin{align} &\frac{1}{2} \frac{\partial}{\partial t} {C}_f \left( Du_m;u_m \right) + g \; {C}_f \left( Du_m;\frac{\partial \eta}{\partial x_m} \right) + \frac{1}{2} \frac{\partial}{\partial x_m} \left[ {C}_f\left( D u_n u_m;u_m \right) \right]\nonumber\\ &\quad = \frac{1}{2} {C}_f \left( D u_m u_n; \frac{\partial u_m}{\partial x_n} \right) - \frac{1}{2} {C}_f \left( D u_m; u_n \frac{\partial u_m}{\partial x_n} \right)\nonumber\\ &\qquad + \frac{1}{2} g {C}_f \left( D u_m; \frac{\partial \eta }{\partial x_m} \right) \nonumber\\ &\qquad - \frac{1}{2} g {C}_f \left( u_m; D \frac{\partial\eta}{\partial x_m} \right) - \frac{1}{2} {C}_f \left( D u_m; \frac{\tau_m}{D} \right) - \frac{1}{2} {C}_f \left(u_m;\tau_m\right). \end{align}

A potential energy balance is derived by multiplying (2.10) with $\frac {1}{2}g \bar {\eta }^f$, followed by steps similar to the kinetic balance,

(2.16)\begin{equation} \frac{1}{2} g \frac{\partial}{\partial t} {C}_f \left( \eta; \eta \right) + g {C}_f \left( \eta; \frac{\partial Du_m}{\partial x_m} \right) = 0. \end{equation}

Combining the potential and kinetic energy balances yields a frequency balance for the total organized energy,

(2.17)\begin{equation} \frac{\partial E_f}{\partial t} + \frac{\partial }{\partial x_n}\left( F^{L}_{n,f}+F^{{NL}}_{n,f}\right) = S^{{NL}}_f + S^{\tau}_f, \end{equation}

where the flux $F_{n,f}$ is decomposed into a (quasi-)linear (superscript $L$) and nonlinear (superscript $NL$) contribution (i.e. $F_{n,f} = F_n^{{L}} + F_n^{{NL}}$). In a linear approximation the correlation between $u$ and $\eta$ is the only contribution to the flux. For this reason we will consider $F_n^{{L}}$ the ‘linear’ part of the flux. However, both $u$ and $\eta$ contain nonlinear corrections, and the correlation does not exactly equal the $E c_{g}$ fully linear approximation. The other frequency-dependent terms are given by

(2.18)\begin{align} E_f &= \frac{1}{2}g {C}_f \left( \eta; \eta \right) + \frac{1}{2} {C}_f \left( D u_m; u_m \right), \end{align}

(2.19)\begin{align} F^{{L}}_{n,f} &= g {C}_f \left( du_n ; \eta \right), \end{align}

(2.20)\begin{align} F^{{NL}}_{n,f} &= g {C}_f\! \left( \eta u_n ; \eta \right) + \frac{1}{2} {C}_f \left( D u_n u_m; u_m \right), \end{align}

(2.21)\begin{align} S_f^{{NL}} &= \frac{1}{2} \left[ {C}_f\left( D u_m u_n; \frac{\partial u_m}{\partial x_n} \right) - {C}_f\left( D u_m; u_n \frac{\partial u_m}{\partial x_n} \right) \right]\nonumber\\ &\quad + \frac{1}{2} g \left[ {C}_f \left( \eta u_m; \frac{\partial\eta}{\partial x_m} \right) - {C}_f \left( u_m; \eta \frac{\partial\eta}{\partial x_m} \right) \right],\end{align}

(2.22)\begin{align} S_f^{\tau} &={-} \frac{1}{2} {C}_f\! \left( D u_m; \frac{\tau_m}{D} \right) - \frac{1}{2} {C}_f \left(u_m;\tau_m\right). \end{align}

Note that cross-shore variation of the still water depth is included at lowest order in (2.17) through gradients of $F^{{L}}_{n,f}$ (2.19).

This frequency-resolved energy balance is the central result of this work. Apart from the change to a spectral formulation, the balance derivation mirrors that of the bulk balance. The bulk energy, energy flux and frictional dissipation terms ((2.4)–(2.7)) all have clear frequency-domain counterparts ((2.18)–(2.21)). The bulk and frequency-resolved formulations are internally consistent because

(2.23)\begin{equation} \int C_f \left( X; Y \right) \,\mathrm{d}f = \left\langle XY\right\rangle \end{equation}

and

(2.24)\begin{equation} \int {C}_f \left( X; YZ \right)\, \mathrm{d}f = \int {C}_f \left( XY;Z \right) \,\mathrm{d}f = \left\langle XYZ\right\rangle. \end{equation}

Bulk expressions for $E$, $F_n$ and $S^\tau$ are recovered when $E_f$, $F_{n,f}$ and $S^\tau$ are integrated over all frequencies. Further, integrated contributions from the nonlinear interaction term $S^{{NL}}_f$ vanish, consistent with the expectation that nonlinear interactions redistribute energy across frequencies but conserve total energy.

This spectral energy balance accounts for all mechanisms potentially transferring energy between SS and IG waves, including the nonlinear shoaling of bound waves (e.g. List 1992; Janssen et al. 2003; Battjes et al. 2004), excitation of free IG waves over a sloping bed (e.g. Mei & Benmoussa 1984; Nielsen & Baldock 2010; Contardo et al. 2021; Liao et al. 2021) and breakpoint generation (Symonds et al. 1982). However, within the spectral framework all these mechanisms are represented as either contributions to nonlinear flux gradients or nonlinear interactions, so that identifying the dominant mechanism that drives the interaction is in general not possible. Here, we therefore quantify the spectral energy flow, but otherwise will not attempt to interpret these in the context of mechanisms proposed in the literature.

3. Methodology

3.1. Numerical model

The highly detailed flow variables required to evaluate the frequency-resolved energy balance (2.17) in the surfzone are not readily available from field or laboratory experiments. Instead, we use the nonlinear and fully dispersive wave model SWASH to obtain the required variables. The wave model SWASH is a multilayer non-hydrostatic wave-flow model, and essentially a direct numerical implementation of the Reynolds-averaged Navier–Stokes equations (Zijlema, Stelling & Smit 2011). Previous studies have shown that SWASH accurately describes the nonlinear transformation and breaking of SS waves (e.g. Smit, Zijlema & Stelling 2013; Smit et al. 2014), the IG wave dynamics (e.g. Rijnsdorp, Smit & Zijlema 2014; de Bakker et al. 2016; Fiedler et al. 2019) and run-up oscillations at the beach (e.g. Ruju, Lara & Losada 2014; Nicolae Lerma et al. 2017).

3.2. Model set-up

The energy balance is studied for SS waves normally incident on mild (1/100) and moderately (1/30) sloping beaches. Irregular waves were generated in 15 m depth using a weakly reflective wavemaker to avoid rereflection at IG frequencies. The irregular wave field had a JONSWAP (Joint North Sea Wave Project) spectral shape with significant wave height $H_{m0}=3$ m and peak period $T_p=10$ s. The wavemaker signal was based on weakly nonlinear theory (Hasselmann 1962) to suppress generation of free IG waves (see Rijnsdorp et al. (2014) for further details).

A horizontal resolution of 1 m was used, which is approximately $1/70$ of the peak wavelength at the edge of the surfzone, in combination with a time-step of 0.025 s (corresponding to a Courant number $CFL\approx 0.3$). A fine vertical resolution (10 vertical layers) was used to ensure that the kinematic condition for the onset of wave breaking (i.e. the orbital velocity exceeds the wave celerity) is captured accurately and wave breaking occurs at the correct location (Smit et al. 2013). The tangential bottom stress $\tau _b$ was estimated using the law of the wall with roughness height $4\times 10^{-4}$ m, a representative value for smooth concrete (e.g. Chow 1959). The $k-\epsilon$ turbulence model accounts for vertical mixing due to shear in the water column but does not describe overturning waves, air entrainment, breaking generated turbulence or the transfer of organized energy into turbulence and heat. In NLSWE-based models like SWASH, waves steepen until a jump-discontinuity (or shock) develops that represents a broken wave. Now SWASH uses the momentum-conserving (shock-capturing) numerical scheme of Stelling & Duinmeijer (2003) to solve the momentum equations (2.2) in their conservative form (2.3). As a result, jumps-discontinuities dissipate total energy in accordance with a hydraulic jump. Such shock-capturing numerical schemes allows NLSWE-based wave models to simulate the bulk dissipation of a breaking wave without accounting for wave-breaking generated turbulence (e.g. Tissier et al. 2012; Smit et al. 2013).

The model depth-averaged horizontal velocities and surface elevation were sampled with a horizontal resolution of 1 m and a temporal resolution of 20 Hz, sufficient to estimate accurately the terms of the energy balance at IG frequencies (Appendix A). The run-up signal was defined as the location where the instantaneous water depth ($d+\eta$) was smaller than 5 cm. Variables were output for 240 min (${>}1000$ peak wave periods) after a spin-up time of 10 min.

3.3. Data analysis

Time series of velocity or sea surface elevation $X$ are separated with

(3.1)\begin{equation} X=X_{vlf}+X_{ig}+X_{ss}, \end{equation}

where $X_{ig}$ represents the IG frequency band (defined as $\frac {1}{20}f_p< f\leq 0.5f_p$, with $f_p$ the peak period at the numerical wavemaker) and $X_{ss}$ corresponds to the SS frequency band (defined as $f > 0.5 f_p$). Frequencies below the IG band are considered part of the VLF band ($f\leq \frac {1}{20}f_p$). In the following, we use labels with a normal font (in subscripts, for readability) and capital font (in-line text) to distinguish between variables belonging to the VLF, IG and SS frequency band and from combinations thereof (separated by a comma in subscripts and by a spaced en-dash for in-line text).

Spectral moments were computed by integrating surface elevation and run-up spectra over their respective frequency bands ($m_n=\int E f^n \,{\rm d}f$). The significant wave height $H_{m0}$ and run-up height $R_{m0}$ were computed as $4 \sqrt {m_0}$ with $m_0$ the zeroth-order moment of the surface elevation or run-up spectra. The mean wave period was computed as $T_{m01}={m0}/{m1}$. To quantify the bulk wave nonlinearity, the Ursell parameter $\text {Ur}=\delta /\mu ^2$ was computed based on the significant SS wave height ($\delta ={a}/({d+\bar {\eta }})$, with the wave amplitude $a=\frac {1}{2}H_{m0,{SS}}$) and the mean wave period ($\mu = k d$, with wavenumber $k$ computed based on $T_{m01,{SS}}$). Similarly, we also estimated the nonlinearity of the IG waves by estimating the wave amplitude as $a=\frac {1}{2}H_{m0,{IG}}$ and the wavenumber based on the mean IG wave period $T_{m01,{IG}}$.

The frequency dependent NLSWE energy balance (2.17) was evaluated based on the modelled (depth-averaged) flow and free surface variables. To compare the NLSWE energy balance with existing theories, we evaluated the new theory and the energy balance terms of Henderson & Bowen (2002), Henderson et al. (2006) and Herbers & Burton (1997) in a consistent manner. All spectra and cross-spectra in this work were computed with $50\,\%$ overlapping Hanning windows and a segment length of 8000 samples, resulting in a frequency resolution $\Delta f=0.0025$ Hz.

To quantify the contribution from IG and SS frequencies to the NLSWE energy balance, terms of the frequency dependent energy balance were integrated over both the SS and IG frequency bands. Furthermore, we quantified how different nonlinear contributions affected the nonlinear flux (2.20) and nonlinear interaction term (2.21) (see Appendix B). The time-domain analysis (Fiedler et al. 2019) avoids the cumbersome bookkeeping of traditional estimates based on bispectral integrations (e.g. de Bakker et al. 2015). With this approach, the frequency-dependent nonlinear flux term $F^{NL}_f$ and nonlinear interaction term $S_f^{NL}$ were decomposed as

(3.2)\begin{gather} F^{NL}_f=F^{NL}_{f,{ig,ig,ig}}+F^{NL}_{f,{ig,ig,ss}}+F^{NL}_{f,{ig,ss,ss}}+F^{NL}_{f,{ss,ss,ss}}+F^{NL}_{f,{vlf}}, \end{gather}

(3.3)\begin{gather}S_f^{NL}=S^{NL}_{f,{ig,ig,ig}}+S^{NL}_{f,{ig,ig,ss}}+S^{NL}_{f,{ig,ss,ss}}+S^{NL}_{f,{ss,ss,ss}}+S^{NL}_{f,{vlf}}. \end{gather}

In these equations, terms with subscript $[\cdots ]_{ig,ig,ig}$ and $[\cdots ]_{{ss,ss,ss}}$ represent correlations between three IG and three SS wave components, respectively. Terms with $[\cdots ]_{{ig,ig,ss}}$ represent correlations between two IG and one SS wave component, terms with $[\cdots ]_{{ig,ss,ss}}$ represent correlations between a single IG and two SS wave components, and the terms with subscript $[\cdots ]_{vlf}$ account for correlations between three wave components of which at least one is a VLF component.

Integrating the above frequency dependent $F^{NL}_f$ and $S_f^{NL}$ over the IG frequencies results in the decomposed bulk IG flux and bulk IG nonlinear interaction terms:

(3.4)\begin{gather} F^{NL}_{IG}=F^{NL}_{IG,ig,ig,ig}+F^{NL}_{IG,ig,ig,ss}+F^{NL}_{IG,ig,ss,ss}+F^{NL}_{IG,vlf}, \end{gather}

(3.5)\begin{gather}S_{IG}^{NL}=S^{NL}_{IG^\pm,ig,ig,ig}+S^{NL}_{IG,ig,ig,ss}+S^{NL}_{IG,ig,ss,ss}+S^{NL}_{IG,vlf}. \end{gather}

In these decomposed terms, correlations between three SS components are not included as they do not contribute to IG fluxes. Furthermore, $S^{NL}_{IG,ig,ig,ig}$ integrates to zero over the IG band as interactions between three IG components redistribute energy but do not account for a net energy transfer to or away from IG frequencies. In order to quantify the energy flow within the IG band by $S^{NL}_{IG,ig,ig,ig}$, we integrated $S^{NL}_{f,{ig,ig,ig}}$ over the IG band by considering either positive or negative interactions (denoted as $S^{NL}_{IG^\pm,ig,ig,ig}=S^{NL}_{IG^+,ig,ig,ig}+S^{NL}_{IG^-,ig,ig,ig}$).

To evaluate the biphase $\beta$ for the individual IG and SS band and the correlations between both bands, we computed the asymmetry ($As$) and skewness ($Sk$) in the time-domain for all possible triads (excluding the VLF motions) following the approach of Fiedler et al. (2019). Subsequently, we computed the biphase $\beta$ and bicoherence for each combination of triads. For example, the biphase of triads between three SS components was computed as

(3.6)\begin{equation} \beta_{ss,ss,ss}=\tan^{{-}1}\left(\frac{{As}_{ss,ss,ss}}{{Sk}_{ss,ss,ss}}\right), \end{equation}

with a corresponding bicoherence of

(3.7)\begin{equation} \sqrt{{As}_{ss,ss,ss}^2+{Sk}_{ss,ss,ss}^2}. \end{equation}

4. Bulk wave evolution

The cross-shore evolution of the wave height ($H_{IG}$, $H_{SS}$) and wave spectrum illustrate the well known transition from SS dominance offshore (depth $d=10$ m) to IG dominance near the shoreline and in run-up (figure 1). On both slopes, the wave height $H_{m0}$ increased slightly during shoaling, followed by breaking at $d\approx 6$ m (with the breakpoint $x_b$ located at $x{_b}\approx 200$ and $700$ m for the 1/30 and 1/100 slope, respectively) (figure 1a,b). During shoaling, energy at the higher harmonics increased gradually, associated with nonlinear wave-interactions between SS components (e.g. $2f_p$ in figure 1g,h). This spectral signature is consistent with gradual steepening of SS waves as they shoal and break (figure 2).

Figure 1. Overview of bulk wave evolution. (a–f) Total, SS, IG and VLF wave height $H$ and run-up $R$, set-up $\bar {\eta }$ and still water depth $d$ (see legends) versus cross-shore location $x$. The still water depth is zero at $x=0$, and $x>0$ offshore. Filled circles in (c,d) indicate significant run-up (with values relative to the right axis). In (a–f), cells that are always or intermittently wet are shown with solid and dashed curves, respectively. Dashed vertical line indicates the outer edge of the surfzone, where the total wave height is largest). (g–j) Power spectra of surface elevation and run-up on 1/30 and 1/100 slopes (blue and orange curves, respectively) at (g) $d= 10$ m, (h) $d=6$ m (approximate SS break point), (i) $x^\prime =0$ (the last cell that is always wet) and (j) of the run-up signal. In (j), black lines indicate $f^{-4}$ slopes separated by distance $\beta ^{3}$ (dash–dot) and $\beta ^4$ (dashed).

Figure 2. Time series of the total surface elevation (thin grey) and the band-passed IG surface elevation signal (blue) from $d=10$ m depth (top) to $d=0$ m depth (bottom) on slopes of (a) 1/30 and (b) 1/100. At the bottom of both (a) and (b), the total (thin grey), band-passed IG (blue) and VLF (green) run-up signal is shown. For improved visualization, the IG signal of the surface elevation and run-up is translated vertically to oscillate around the mean of the respective total signal. The surfzone, where SS wave energy decreases, starts around $d=6\unicode{x2013}7$ mm on both slopes (indicated by the green dashed line). The IG energy begins to decrease (flux gradient $\partial _x F_{IG}<0$) in approximately 1.5 and 4 m depth on the 1/30 and 1/100 slope, respectively (indicated by the red dashed line).

On both slopes, the IG wave height $H_{m0,{IG}}$ was smaller than the SS wave height $H_{m0,{SS}}$ seaward of and within most of the surfzone. Seaward of the surfzone, IG waves were out of phase with the forcing SS wave groups (figure 2), consistent with bound IG waves (e.g. Longuet-Higgins & Stewart 1962). On the 1/100 slope, $H_{m0,{IG}}> H_{m0,{SS}}$ for $x<100$ m, whereas on the 1/30 slope, IG motions were largest only near the shoreline ($x<20$ m) (figure 1c,d,i). On both slopes, motions at IG frequencies dominated the run-up signal (figures 1 and 2). The filtered time signals indicate that IG waves gradually pitch forward on the 1/100 slope to a bore-like front with small SS waves riding on their crest (figure 2).

On both slopes, IG energy levels of the run-up spectra were much larger compared with the surface elevation spectra at $x^\prime =0$ (the shallowest location where a cell is always wet), whereas run-up energy levels at SS frequencies were lower compared with SS spectral energy levels at $x^\prime =0$ (figure 1i,j). The run-up energy levels roll-off as $f^{-4}$ at SS frequencies, and the roll-off (often associated with saturation) extends well into the IG band (figure 1j). For the same incident waves, the run-up energy levels depend on the beach slope $\beta$ and are separated by a distance proportional to $\beta ^{-l}$ where $l=3$ and $4$ for the present slopes of 1/30 and 1/100, respectively.

During shoaling ($x>x_{b}$), the bulk SS Ursell parameter $\text {Ur}_{SS}\approx 0.5$. Here, dispersive and nonlinear effects are equally important in the SS band, indicating that the classical Boussinesq scaling holds. The present energy balance based on the NLSWE is non-dispersive and does not fully explain the physics of shoaling SS waves. However, the NLSWE energy balance is potentially accurate in the surfzone and especially the inner surfzone, where nonlinearity dominates over dispersion ($\text {Ur}_{SS}>1$). The IG Ursell parameter $\text {Ur}_{IG}>1$ for $d<10$ m (not shown), and the NLSWE energy balance is expected to explain the IG wave dynamics throughout the domain.

5. Energy balance

5.1. Sea-swell balance

The cross-shore variation of the (frequency dependent) NLSWE energy balance terms shows that seaward of the surfzone the energy flux increases ($\partial _x F_f>0$) at harmonics $2f_p$ and $3f_p$ and decreases at $f_p$ (blue and red shading in figure 3c,d for $x>x_{b}$, respectively). These flux gradients are largely balanced by nonlinear interactions $S^{NL}_f$ with negligible bottom stress $S^\tau _f$ (figure 3e–h). Integrating the energy balance terms over the primary SS ($0.5 f_p < f \leq 3 f_p$) and over the superharmonic frequencies ($f>3f_p$) shows that the energy balance approximately closes seaward of the surfzone ($x>x_b$) on both slopes (grey curves in figure 4c,d).

Figure 3. Cross-shore variation of the surface elevation and run-up spectra (a,b) and frequency-dependent NLSWE energy balance terms (c–h) on a 1/30 (a,c,e,g) and 1/100 slope (b,d, f,h). The horizontal grey dotted line indicates the offshore peak frequency $f_p$ and black dashed lines indicate the limits of the IG ($\frac {1}{20}f_p< f\leq f_p/2$), primary SS ($\frac {1}{2}f_p< f\leq 3f_p$) and superharmonic ($f>3f_p$) frequency bands. The dashed vertical line indicates the breakpoint $x_b$ ($d\approx 6$ m).

Figure 4. Cross-shore variation of the bulk fluxes (a,b) and non-zero bulk NLSWE energy balance terms (c,d) of the SS band on slopes of 1/30  (a,c) and 1/100 (b,d). The balance terms in (c,d) are further decomposed into primary ($\frac {1}{2}f_p< f\leq 3f_p$) and superharmonic ($f>3f_p$) SS frequencies. The dashed vertical line indicates the breakpoint $x_b$ ($d\approx 6$ m).

Within the surfzone ($x< x_b$), nonlinear interactions strengthen and transfer energy towards increasingly high superharmonic frequencies ($f>3f_p$, figure 3e, f). The energy flux $F_f$, however, did not increase at the superharmonic frequencies (figure 3c,d). Lacking significant dissipation from bottom stress (figure 3g,h), the superharmonic residual ($f>3f_p$) was negative and relatively large in the surfzone (grey dashed lines in figure 4c,d). The residual was smaller in the more energetic primary band (grey solid lines in figure 4c,d). The transfer of SS energy to superharmonic frequencies with breaking waves is generally consistent with previous studies based on weakly nonlinear but dispersive balances (e.g. Herbers et al. 2000; Smit et al. 2014).

5.2. Infragravity balance

On both slopes nonlinear interactions increase the energy flux at the IG frequencies in the shoaling region and much of the surfzone (pale blue, $S_f^{NL}>0$ and $\partial _x F_f>0$ in figure 3c–f). Deeper inside the surfzone, nonlinear interactions changed sign, resulting in $\partial _x F_f<0$ at higher IG frequencies. On the 1/30 slope, lower IG frequencies continued to receive energy up to close to the shoreline. In contrast, on the 1/100 slope, $\partial _x F_f$ gradually changes sign over the whole IG band and all IG frequencies lose energy for $x<300$ m. The IG dissipation from bottom friction $S_f^{\tau }$ was small except for the lowest frequencies in very shallow water (figure 3g,h).

On the 1/100 slope, bulk IG fluxes $F_{IG}$ were simplified by the relatively low amplitude of seaward going IG waves ($F_{IG}^{{L}^+} \gt \gt F_{IG}^{{L}^-}$, figure 5b). At $d=10$ m, the relative radiation coefficient $R^2_{IG}$ = $F_{IG}^{{L}^-}/F_{IG}^{{L}^+} \approx 0.05$. With this negligible seaward IG radiation, $F_{IG}^{{L}} \approx F_{IG}^{{L}^+} \approx E c_g$ as assumed by de Bakker et al. (2016) and others. On the 1/30 slope, $R^2_{IG} \approx 0.6$ at $d=10$ m, and $E c_g$ is inaccurate (figure 5a). The accuracy of the present linear estimates of seaward and shoreward IG fluxes for these nonlinear waves is unknown.

Figure 5. Cross-shore variation of the (decomposed) bulk IG fluxes (a–d), the bulk IG energy balance terms (e, f) and the (decomposed) bulk nonlinear interaction terms (g,h) on a 1/30 (a,c,e,g) and 1/100 slope (b,d, f,h). The shoreline is at $x=0.$ In (a,b), the decomposition of the linear IG flux $F^{L}_{IG}$ into shoreward and seaward components ($F^{L^+}_{IG}$ and $F^{L^-}_{IG}$, respectively) is not shown when the WKB (Wentzel–Kramers–Brillouin) assumption is violated (${\omega ^2 d}/{g}<10\,h_x^2$, with $\omega$ from $T_{m01,{IG}}$). Vertical lines indicate the location of (dashed) the seaward surfzone edge and (dash–dotted) where the IG flux gradient $\partial _x F_{IG}$ changes sign.

On both slopes, the nonlinear contribution $F_{IG}^{{NL}}$ to the total IG flux ($F_{IG}$) is significant (figure 5a,b). In the shoaling region and a large part of the surfzone, $F_{IG}^{NL}<0$ due to a negative correlation between IG waves and SS components, as occurs with bound IG waves that are $180^\circ$ out of phase with wave group forcing. In the surfzone, the phase drifts away from $180^\circ$ (e.g. Janssen et al. 2003). The cross-shore location of the maximum $F_{IG}$ (where $\partial _x F_{IG}$ changed sign) is around $x=450$ m ($\approx 0.65 x_b$) on the 1/100 and $x=50$ m ($\approx 0.25x_{b}$) on the 1/30 slope (figure 5a,b). Farther shoreward in the inner surfzone, IG waves are losing energy on both slopes. These IG energy losses ($\partial _x F_{IG}<0$) were largely explained by $S^{NL}_{IG}$ and a small bottom stress $S^{\tau }_{IG}$. The IG residual is relatively small everywhere (grey lines in figure 5e, f).

5.3. Contributors to the flux and nonlinear interactions at IG frequencies

To quantify the contribution from different triads to the nonlinear interactions, we decomposed $S^{NL}_{IG}$ into correlations between different combinations of the SS wave signal, the IG wave signal and the remaining VLF signal (see § 3.3 and Appendix B). Some previous studies estimated the total flux from the linear flux alone (e.g. Thomson et al. 2006; de Bakker et al. 2015), which can be problematic in the nearshore (e.g. Henderson et al. 2006). In the present simulations, the linear flux $F_{IG}^{L}$ under-predicted the total flux $F_{IG}$ in the inner surfzone shoreward of the location where $F_{IG}$ peaked, and over-predicted $F_{IG}$ seaward of the location where $F_{IG}$ peaked (compare black with blue curves, figure 5a,b). Seaward of the breakpoint, the dominant contributor to $F_{IG}^{NL}$ is due to interactions between one IG and two SS wave components $F^{NL}_{IG,ig,ss,ss}$. This is consistent with numerous previous studies (e.g. Longuet-Higgins & Stewart (1962), Henderson et al. (2006), Ruju et al. (2012), Mendes et al. (2018) and many others). The interactions between one IG and two SS waves, associated with the forcing of IG waves by the SS waves (e.g. Hasselmann 1962; Longuet-Higgins & Stewart 1962), also (nearly) completely explained $S_{IG}^{NL}$ (figure 5g,h). Within the surfzone, $S^{NL}_{IG,ig,ss,ss}$ did not explain all the negative work by $S_{IG}^{NL}$, and additional interactions are required to close the IG energy balance. Nonlinear interactions with at least one VLF component ($S^{NL}_{IG,vlf}$) did not provide a significant contribution (figure 5g,h). Instead, the inclusion of interactions between two IG and one SS component $S^{NL}_{IG,ig,ig,ss}$ was required to accurately estimate $S_{IG}^{NL}$ inside the surfzone. Separating $S^{IG}_{ig,ig,ig}$ (which integrated to zero over the IG band) into positive and negative contributions ($S^{IG^+}_{ig,ig,ig}$ and $S^{IG^-}_{ig,ig,ig}$, respectively) quantifies the energy flow within the IG band due to triad interactions among IG waves. For both slopes, $S^{IG^+}_{ig,ig,ig}$ and $S^{IG^-}_{ig,ig,ig}$ became non-zero shoreward of the location where $F_{IG}$ peaked. Their contribution was largest for the mild-slope (figure 5h), where they resulted in a flow of energy from lower IG to higher IG frequencies (not shown). This energy gain at higher IG frequencies is subsequently balanced by $S^{NL}_{ig,ig,ss}$, which transports the IG flux to SS frequencies.

5.4. Infragravity wave dissipation

For both slopes, a strong increase followed by an intense reduction of the energy flux occurred at the IG frequencies. To quantify the contribution of the different processes that contributed to the net gain and loss of IG flux inside the surfzone, we cross-shore integrated the IG energy balance terms (only considering cells that were always wet) over the region where the individual terms were either positive or negative.

The net nearshore gain of IG flux was approximately two times larger for the 1/100 slope than the 1/30 slope ($1.71\,\mathrm {m}^4\,\mathrm {s}^{-3}$ versus $0.99\,\mathrm {m}^4\,\mathrm {s}^{-3}$). This gain was nearly completely explained by nonlinear interactions (figure 6a), although $S_{IG}^{NL}$ overestimated the net gain by approximately $25\,\%$ on the mild slope (i.e. a net residual of $\approx -25\,\%$). For both slopes, the nonlinear interactions were dominated by interactions between one IG and two SS components, associated with the well known forcing of IG waves by SS wave groups.

Figure 6. Relative contribution of different energy balance terms to the IG energy flux (a) gain and (b) loss on a 1/30 (blue) and 1/100 (orange) slope. Triad interactions among three IG waves (IG–IG–IG) integrate to zero. Instead, loss terms were cross-shore integrated over cells with negative values, and gain terms were integrated over positive cells (subscript ${IG}^-$ and ${IG}^+$, respectively).

The net nearshore loss of IG flux was approximately six times larger for the 1/100 slope than the 1/30 slope ($-1.93$ versus $-0.32 \,\mathrm {m}^4\,\mathrm {s}^{-3}$). For the mild slope, the net loss exceeded the net gain ($-1.93$ versus $1.71\,\mathrm {m}^4\,\mathrm {s}^{-3}$). This mismatch can be explained by the incoming flux contribution from IG waves generated at the wavemaker (estimated as $F^{L^+}_{IG}+F^{NL}_{IG}=0.28\, \mathrm {m}^4\,\mathrm {s}^{-3}$), indicating near complete dissipation of shoreward propagating IG motions on the mild slope that is consistent with negligible reflections at IG frequencies (figure 5b). For the 1/30 slope, the net loss was only a fraction of the net gain ($-0.32$ versus $0.99\,\mathrm {m}^4\,\mathrm {s}^{-3}$), associated with non-negligible reflection of IG waves. However, the cross-shore integration was limited to cells that were always wet and additional dissipation that occurred in shallower water (including the swash zone) is missing from this analysis and will contribute to this mismatch.

On both slopes, the nearshore loss of the IG flux was nearly completely explained by nonlinear interactions $S_{IG}^{NL}$ and bottom friction $S_{IG}^\tau$ (figure 6b), although the combined effect of $S^{NL}$ and $S_{IG}^\tau$ overestimated the net IG flux loss resulting in a residual of ${\approx }10\,\%$–$20\,\%$ (figure 6b). For both slopes, bottom friction accounted for at most ${\approx }25\,\%$ of the nearshore IG flux losses and nonlinear interactions were responsible for the majority of the flux loss. Separating the nonlinear interactions into correlations between different combinations of VLF, IG and SS wave components (§ 5.3) indicates that interactions between one IG and two SS components ($S^{NL}_{IG,ig,ss,ss}$) resulted in the largest transfer of energy away from IG frequencies. The relative contribution of these interactions was comparable for both slopes, and they explained up to ${\approx }65\,\%$ of the nearshore loss in IG flux. Interactions between one IG and two SS components ($S^{NL}_{IG,ig,ig,ss}$) were typically weaker but still significant on both slopes (${\approx }20\,\%$). Energy transfers to frequencies below the IG band ($S^{NL}_{IG,vlf}$) were small on both slopes. Consistent with interpretation of laboratory experiments (e.g. Baldock et al. 2000; Baldock 2006, 2012) and results from weakly nonlinear energy balances (e.g. Thomson et al. 2006; de Bakker et al. 2015, 2016), these results confirm that interactions with two IG components became significant inside the inner surfzone and result in substantial energy transfer from IG to SS frequencies, where the wave energy is ultimately dissipated.

6. Discussion

6.1. Comparison NLSWE energy balance with existing theories

Theories for nonlinear interactions range from Boussinesq models for resonant waves in shallow water (e.g. Herbers & Burton 1997) up to more generalized theories that account for full linear dispersive effects (e.g. Bredmose et al. 2005; Janssen 2006). Several previous studies used the nonlinear interaction term based on the Boussinesq scaling from Herbers & Burton (1997), denoted as $S^{NL}_{HB97}$, to evaluate the nearshore IG energy balance (e.g. Thomson et al. 2006; de Bakker et al. 2015, 2016). Assuming weakly nonlinear waves by retaining terms up to $O(\delta ^3)$, and relating $u$ to $\eta$ using $u=({c}/{d})\eta$ and ${\partial \eta }/{\partial x}=-({1}/{d})({\partial \eta }/{\partial t})$ with $c=\sqrt {g d}$ (which combine to the linearized flat-bed assumption, ${\partial u}/{\partial x}=-({1}/{c})({\partial \eta }/{\partial t})$), the nonlinear interaction $S^{NL}$ term from the NLSWE balance reduces to $S^{NL}_{HB97}$ (Appendix C).

The energy balance with Boussinesq scaling of Herbers & Burton (1997) is

(6.1)\begin{equation} \frac{\partial F_{f,{HB97}}}{\partial x}+S^{NL}_{f,{HB97}}=S^\tau_f, \end{equation}

in which $F_{f,{HB97}}=E c_g$ and $S^{NL}_{f,{HB97}}$ depends on a coupling coefficient and $Q$, the imaginary part of the bispectrum. We computed $S^{NL}_{HB97}$ following Henderson et al. (2006),

(6.2)\begin{equation} S^{NL}_{f,{HB97}}=\frac{3 g \omega}{2 h} Q_f ( \eta^2, \eta), \end{equation}

where $h$ is the mean water depth ($h=d+\bar {\eta }$). The flux term $F$ of the Boussinesq scaled energy balance breaks down in the case of significant shoreline reflections (as is the case on the 1/30 slope, figure 5a). As an alternative, we used the linear flux estimates to evaluate this term, which yields similar results in the case of weak shoreline reflections (figure 5b).

The Boussinesq scaled energy flux and nonlinear interaction term are in balance in the shoaling region ($x>x_b$, figure 7a,b), which confirms that in this region the weakly nonlinear and dispersive balance closes to the order considered. However, this balance deteriorated significantly in the surfzone on both slopes, where wave nonlinearity (amplitude over depth) is large. Here, $S^{NL}_{IG,HB97}$ overestimated the linear flux gradient (figure 7a,b). In contrast, the fully nonlinear but non-dispersive energy balance based on the NSLWE closed with much smaller residual compared with the Boussinesq scaled energy balance throughout the surfzone (compare figures 7e and 7f).

Figure 7. Nearshore IG energy balance for different theories on (a,c) 1/30 and (b,d) 1/100 slope. Panels (a–d) show the cross-shore evolution of the cross-shore energy flux gradient minus the bottom stress ($\partial _x F_{IG}-S^{\tau }_{{IG}}$) and the nonlinear interaction term ($S^{{NL}_{IG}}$) of (a,b) NLSWE balance (this paper) and Boussinesq scaled balance Herbers & Burton (1997), and (c,d) balance of Henderson & Bowen (2002) and Henderson et al. (2006). Panels (e–g) are scatter plots of $S^{NL}_{IG}$ versus $\partial _x F_{IG}-S^{\tau }_{{IG}}$ for (e) the NLSWE balance ( f) Herbers & Burton (1997) and (g) Henderson & Bowen (2002) and Henderson et al. (2006). In each panel, the relative bias and skill scores are given for the respective theory. Solid lines indicate closure of the IG energy balance. For the Boussinesq scaled balance, results are not plotted when the underlying WKB assumption is violated (${\omega ^2 d}/{g}<10 h_x^2$, with $\omega$ from $T_{m01,{IG}}$).

Henderson & Bowen (2002) derived a nonlinear energy balance for IG wave frequencies that considers the mean potential and total pseudokinetic energy derived from the Lagrangian mean energy in the water column (see Appendix D). Henderson et al. (2006) subsequently adapted this theory by assuming weak nonlinearity. The weakly nonlinear theory of Henderson et al. (2006) did not provide a closing balance between the energy flux gradient $\partial _x F_{IG,H06}$ and nonlinear interaction term $S^{NL}_{IG,H06}$ (figures 7c,d and 7g). The residual was largest in the inner surfzone, where the assumption of weak nonlinearity breaks down (e.g. $S^{{NL}}_{ig,ig,ss}$ becomes non-negligible, figure 5g,h), and $S^{NL}_{H06}$ under estimated the loss of IG energy flux. The nonlinear version of this balance (Henderson & Bowen 2002) describes a closer balance between the flux gradient $\partial _x F_{H02}$ and nonlinear interaction term $S^{{NL}}_{H02}$ throughout most of the surfzone (figure 7c,d). However, $S^{NL}_{H02}$ blows up close to the shoreline (where $d\rightarrow 0$ m) due to its dependence on the Lagrangian velocity ($\tilde {u}=u+({1}/{d})\eta u$, Appendix D).

6.2. Infragravity wave breaking

Past studies on mild slopes have linked substantial IG energy dissipation with IG breaking (e.g. Van Dongeren et al. 2007; de Bakker, Tissier & Ruessink 2014; de Bakker et al. 2015, 2016). Combining energy balances with bispectral analysis, several previous studies found that significant transfers of IG energy to higher harmonics coincided with skewed and asymmetric shapes at the IG frequencies (as was observed in the time signals on the 1/100 slope, figure 2) and linked this to the possibility that IG waves were breaking (e.g. de Bakker et al. 2015, 2016). These studies relied on weakly nonlinear energy balances that do not close at IG frequencies, complicating a quantitative understanding of the nearshore IG wave dissipation. The fully nonlinear non-dispersive energy balance derived in this work did (nearly) completely close at IG wave frequencies, and allowed for a meticulous study into the nearshore dissipation of IG waves (§ 5.4).

The IG wave dynamics on the 1/30 and 1/100 slope showed remarkably similar patterns up to the inner surfzone. In the shoaling region, energy was continuously transferred from SS wave groups to the IG band by nonlinear interactions (figure 5e, f). As the waves approach the surfzone, the biphase between these IG–SS–SS triads gradually moved away from the equilibrium response $\beta _{ig,ss,ss}=180^\circ$ (figure 8c). The IG band received energy through the first part of the outer surfzone (figure 5e, f), up to the location where $\beta _{ig,ss,ss}\approx 0^\circ$ (figure 8c), indicating that here the largest SS waves propagate on the crest of the IG waves (see the time signals near the dashed red line in figure 2). Farther shoreward, the IG band lost energy to the SS band due to $S^{NL}_{ig,ss,ss}$ interactions. These interactions persisted until deep inside the surfzone, and were responsible for ${\approx }65\,\%$ of the net IG flux loss (figure 6b). In the inner surfzone, nonlinear interactions between two IG and one SS harmonic became significant on both slopes (figure 5g,h), accounting for approximately $20\,\%$ of the net IG flux loss (figure 6b).

Figure 8. Cross-shore variation of the (a–e) biphase $\beta$ and (e–j) bicoherence $C$ versus total depth for (a,b) the full signal (all frequencies), and (b–j) triads between different combinations of IG and SS components based on the total signal (full line) and signal corresponding to shoreward propagating waves (dotted line) for the 1/30 (blue lines) and 1/100 (orange lines) slope. Vertical lines indicate the location of (dashed) the seaward surfzone edge and (dash–dotted) where the IG flux gradient $\partial _x F$ changes sign (see figure 5). To suppress noise, biphases are not shown when the corresponding bicoherence $<0.05$.

Only on the 1/100 slope, triad interactions involving at least two IG components became increasingly significant in the inner surfzone ($x<250$ m), with increasing bicoherence (figure 8i,j). The biphase of these IG triads and of the other triads involving at least one IG harmonic tended to $\beta =-90^\circ$, indicative of a forward-pitching (saw-tooth shaped) wave ($0 < d < 1$ m, figure 2). Surprisingly, the IG–IG–SS triads only accounted for approximately $25\,\%$ of the net IG flux loss, as $S^{NL}_{ig,ss,ss}$ drained the majority of the IG flux (figure 6b). Seaward of the inner surfzone, $S^{NL}_{ig,ss,ss}$ resembled interactions between IG waves and SS groups that caused an energy transfer from the IG waves back to the SS waves. In the inner surfzone IG waves have a saw-tooth shape (figure 2) and such IG–SS–SS triads can also represent IG self-interactions. Although we cannot establish which of these two alternatives are responsible for the energy transfer, we expect that IG self-interactions may dominate $S^{NL}_{ig,ss,ss}$ when the water depth is sufficiently small ($d<1$ m, figure 1d). Seaward of this location, $S^{NL}_{ig,ss,ss}$ accounted for $\approx 46\,\%$ of the total IG flux loss (equivalent to ${\approx }70\,\%$ of the loss induced by IG–SS–SS triads). This indicates that even though IG waves were breaking on the mild slope, nonlinear interactions between IG waves and SS wave groups were responsible for at least half of the net IG flux loss.

6.3. Validity of NLSWE balance at SS frequencies

The energy balance based on the NLSWE provided a good description of the surfzone IG dynamics simulated with the fully nonlinear and dispersive SWASH model. Despite neglecting dispersive effects in NLSWE (not in SWASH), the IG energy balance approximately closes across the surfzone. Furthermore, IG energy losses do occur, but are mostly explained through nonlinear transfers to the SS band. Residuals at SS frequencies are harder to directly interpret. Outside the surfzone, the small residual (figure 4c,d) is likely due to the exclusion of work performed by the non-hydrostatic pressure in the flux terms. Inside the surfzone, a combination of neglected physics in the analysis (i.e. turbulent stresses) and dissipative losses due to breaking are responsible for the residual.

In SWASH, breaking is modelled as a discontinuity, and dissipation occurs at the scale of the spatial resolution $\Delta x$. The associated frequency scale of O($\sqrt {gd}/\Delta x$) falls within the superharmonic range ($f>3f_{p}$) throughout the domain. As a consequence, the residual in the primary SS band is likely primarily the result of inaccuracies in the NLSWE derived balance (figure 4c,d). Although we suspect that the large residual in the superharmonic range is a dissipation signature, the required extension to three-dimensional non-hydrostatic flow is outside the current scope of work. We further found that inaccuracies in the NLSWE balance become large when spatial and temporal resolution are mismatched in the analysis. For example, when a high spatial resolution captures the dissipation at superharmonics but superharmonics are not resolved due to a low temporal sampling resolution (Appendix A).

7. Conclusions

A fully nonlinear non-dispersive energy balance for surfzone waves based on the NLSWE was derived and applied to unidirectional spectral waves on a moderately and a mildly sloping beach. The energy balance was evaluated based on output from a fully nonlinear and dispersive wave model (the non-hydrostatic wave-flow model SWASH).

The new theory predicts a closed energy balance at the IG frequencies. Throughout the shoaling region and most of the surfzone, the gain and loss of the IG energy flux was primarily balanced by nonlinear wave interactions. The nonlinear interactions explained the energy growth through the shoaling region up to the outer surfzone, and also caused most of the energy loss deeper inside the surfzone. Bottom friction caused only small dissipation in very shallow water. The nonlinear energy transfer was primarily explained by interactions between a single IG and two SS wave components. Only deep inside the surfzone on the mildly sloping beach, interactions between two IG and a single SS wave component (suggesting IG wave breaking) became significant and contributed to the IG energy flux loss.

Existing theories to estimate the nearshore IG energy balance that assume weak nonlinearity do not account for all nonlinear contributions to the energy flux and/or nonlinear interaction terms. This degrades the energy balance closure and complicates the interpretation of the IG energy balance. The fully nonlinear non-dispersive balance was shown to provide an improved description of the IG energy balance inside the surfzone.