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 Reference Lamb1932)). Nonlinear processes become increasingly important in shallow water, where waves pitch forward and break (Stoker Reference Stoker1957). During shoaling, nonlinear interactions cause a significant transfer of energy from the spectral peak frequency towards superharmonic and subharmonic frequencies (Phillips Reference Phillips1960; Hasselmann Reference Hasselmann1962; Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962; Hasselmann et al. Reference Hasselmann, Sell, Ross and Müller1976). 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 Reference Elgar and Guza1985; Elgar et al. Reference Elgar, Guza, Raubenheimer, Herbers and Gallagher1997), 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 Reference Webb, Zhang and Crawford1991; Aucan & Ardhuin Reference Aucan and Ardhuin2013) but can reach 
$O$(m) in shallow water during storms (e.g. Sheremet et al. Reference Sheremet, Staples, Ardhuin, Suanez and Fichaut2014; Matsuba, Shimozono & Sato Reference Matsuba, Shimozono and Sato2020).
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. (Reference Bertin2018) 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 Reference Symonds, Huntley and Bowen1982; Schäffer Reference Schäffer1993; Janssen, Battjes & Van Dongeren Reference Janssen, Battjes and Van Dongeren2003; Nielsen & Baldock Reference Nielsen and Baldock2010; Contardo et al. Reference Contardo, Lowe, Hansen, Rijnsdorp, Dufois and Symonds2021; Liao et al. Reference Liao, Li, Liu and Zou2021, and many others). Such theoretical background combined with numerical modelling (e.g. Reniers et al. Reference Reniers, Van Dongeren, Battjes and Thornton2002; Van Dongeren et al. Reference Van Dongeren, Reniers, Battjes and Svendsen2003; Lara, Ruju & Losada Reference Lara, Ruju and Losada2011), laboratory experiments (e.g. Boers Reference Boers1996; Baldock et al. Reference Baldock, Huntley, Bird, O'Hare and Bullock2000; Baldock & Huntley Reference Baldock and Huntley2002) and field campaigns (e.g. Herbers, Elgar & Guza Reference Herbers, Elgar and Guza1994; Herbers et al. Reference Herbers, Elgar, Guza and O'Reilly1995b; Herbers, Elgar & Guza Reference Herbers, Elgar and Guza1995a; Okihiro, Guza & Seymour Reference Okihiro, Guza and Seymour1992), 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. Reference Battjes, Bakkenes, Janssen and Van Dongeren2004; Van Dongeren et al. Reference Van Dongeren, Battjes, Janssen, Van Noorloos, Steenhauer, Steenbergen and Reniers2007; Baldock Reference Baldock2012; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015). 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 Reference Henderson and Bowen2002; Henderson et al. Reference Henderson, Guza, Elgar, Herbers and Bowen2006; Ruju, Lara & Losada Reference Ruju, Lara and Losada2012; de Bakker, Tissier & Ruessink Reference de Bakker, Tissier and Ruessink2016). 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 Reference Phillips1970),
\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 Reference Freilich and Guza1984; Agnon & Sheremet Reference Agnon and Sheremet1997; Herbers & Burton Reference Herbers and Burton1997). For example, Herbers, Russnogle & Elgar (Reference Herbers, Russnogle and Elgar2000) 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 Reference Smit and Janssen2016).
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. Reference Thomson, Elgar, Raubenheimer, Herbers and Guza2006; Van Dongeren et al. Reference Van Dongeren, Battjes, Janssen, Van Noorloos, Steenhauer, Steenbergen and Reniers2007; Torres-Freyermuth, Lara & Losada Reference Torres-Freyermuth, Lara and Losada2010; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015; Liao et al. Reference Liao, Li, Liu and Zou2021). A subset of nonlinear contributions to the energy flux can be included by using a fully nonlinear low-frequency energy balance (e.g. Phillips Reference Phillips1970; Schäffer Reference Schäffer1993), 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. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006) 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 Reference Longuet-Higgins and Stewart1962). 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. Reference Henderson, Guza, Elgar, Herbers and Bowen2006; Ruju et al. Reference Ruju, Lara and Losada2012; Guedes, Bryan & Coco Reference Guedes, Bryan and Coco2013; Rijnsdorp, Ruessink & Zijlema Reference Rijnsdorp, Ruessink and Zijlema2015; Mendes et al. Reference Mendes, Pinto, Pires-Silva and Fortunato2018). de Bakker et al. (Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015, Reference de Bakker, Tissier and Ruessink2016) 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. Reference Van Dongeren, Battjes, Janssen, Van Noorloos, Steenhauer, Steenbergen and Reniers2007; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015, Reference de Bakker, Tissier and Ruessink2016).
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. Reference Ruju, Lara and Losada2012; Guedes et al. Reference Guedes, Bryan and Coco2013; Rijnsdorp et al. Reference Rijnsdorp, Ruessink and Zijlema2015; Mendes et al. Reference Mendes, Pinto, Pires-Silva and Fortunato2018). 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 (Reference Henderson and Bowen2002) 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,
\begin{gather} \frac{\partial \eta}{\partial t} + \frac{\partial D u_n}{\partial x_n} = 0, \end{gather}
\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$,
\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,
\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
\begin{gather} E = \tfrac{1}{2}g \langle\eta^2\rangle + \tfrac{1}{2}\left\langle Du_m u_m \right\rangle, \end{gather}
\begin{gather}F_n = \left\langle Q_n \left[ \frac{1}{2}u_m u_m + g\eta \right]\right\rangle, \end{gather}
\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,
\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 (Reference Henderson and Bowen2002), a variable 
$X(t)$ is Fourier decomposed over frequencies 
$f$ spaced 
$\Delta f$ apart,
\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}$:
\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}
\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}
\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),
\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,
\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:
\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,
\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,
\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
\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}
\begin{align} F^{{L}}_{n,f} &= g {C}_f \left( du_n ; \eta \right), \end{align}
\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}
\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}
\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
\begin{equation} \int C_f \left( X; Y \right) \,\mathrm{d}f = \left\langle XY\right\rangle \end{equation}and
\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 Reference List1992; Janssen et al. Reference Janssen, Battjes and Van Dongeren2003; Battjes et al. Reference Battjes, Bakkenes, Janssen and Van Dongeren2004), excitation of free IG waves over a sloping bed (e.g. Mei & Benmoussa Reference Mei and Benmoussa1984; Nielsen & Baldock Reference Nielsen and Baldock2010; Contardo et al. Reference Contardo, Lowe, Hansen, Rijnsdorp, Dufois and Symonds2021; Liao et al. Reference Liao, Li, Liu and Zou2021) and breakpoint generation (Symonds et al. Reference Symonds, Huntley and Bowen1982). 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 Reference Zijlema, Stelling and Smit2011). Previous studies have shown that SWASH accurately describes the nonlinear transformation and breaking of SS waves (e.g. Smit, Zijlema & Stelling Reference Smit, Zijlema and Stelling2013; Smit et al. Reference Smit, Janssen, Holthuijsen and Smith2014), the IG wave dynamics (e.g. Rijnsdorp, Smit & Zijlema Reference Rijnsdorp, Smit and Zijlema2014; de Bakker et al. Reference de Bakker, Tissier and Ruessink2016; Fiedler et al. Reference Fiedler, Smit, Brodie, McNinch and Guza2019) and run-up oscillations at the beach (e.g. Ruju, Lara & Losada Reference Ruju, Lara and Losada2014; Nicolae Lerma et al. Reference Nicolae Lerma, Pedreros, Robinet and Sénéchal2017).
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 Reference Hasselmann1962) to suppress generation of free IG waves (see Rijnsdorp et al. (Reference Rijnsdorp, Smit and Zijlema2014) 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. Reference Smit, Zijlema and Stelling2013). 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 Reference Chow1959). 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 (Reference Stelling and Duinmeijer2003) 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. Reference Tissier, Bonneton, Marche, Chazel and Lannes2012; Smit et al. Reference Smit, Zijlema and Stelling2013).
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
\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 (Reference Henderson and Bowen2002), Henderson et al. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006) and Herbers & Burton (Reference Herbers and Burton1997) 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. Reference Fiedler, Smit, Brodie, McNinch and Guza2019) avoids the cumbersome bookkeeping of traditional estimates based on bispectral integrations (e.g. de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015). With this approach, the frequency-dependent nonlinear flux term 
$F^{NL}_f$ and nonlinear interaction term 
$S_f^{NL}$ were decomposed as
\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}
\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:
\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}
\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. (Reference Fiedler, Smit, Brodie, McNinch and Guza2019). 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
\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
\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 Reference Longuet-Higgins and Stewart1962). 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. Reference Herbers, Russnogle and Elgar2000; Smit et al. Reference Smit, Janssen, Holthuijsen and Smith2014).
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. (Reference de Bakker, Tissier and Ruessink2016) 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. Reference Janssen, Battjes and Van Dongeren2003). 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. Reference Thomson, Elgar, Raubenheimer, Herbers and Guza2006; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015), which can be problematic in the nearshore (e.g. Henderson et al. Reference Henderson, Guza, Elgar, Herbers and Bowen2006). 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 (Reference Longuet-Higgins and Stewart1962), Henderson et al. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006), Ruju et al. (Reference Ruju, Lara and Losada2012), Mendes et al. (Reference Mendes, Pinto, Pires-Silva and Fortunato2018) 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 Reference Hasselmann1962; Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962), 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. Reference Baldock, Huntley, Bird, O'Hare and Bullock2000; Baldock Reference Baldock2006, Reference Baldock2012) and results from weakly nonlinear energy balances (e.g. Thomson et al. Reference Thomson, Elgar, Raubenheimer, Herbers and Guza2006; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015, Reference de Bakker, Tissier and Ruessink2016), 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 Reference Herbers and Burton1997) up to more generalized theories that account for full linear dispersive effects (e.g. Bredmose et al. Reference Bredmose, Agnon, Madsen and Schäffer2005; Janssen Reference Janssen2006). Several previous studies used the nonlinear interaction term based on the Boussinesq scaling from Herbers & Burton (Reference Herbers and Burton1997), denoted as 
$S^{NL}_{HB97}$, to evaluate the nearshore IG energy balance (e.g. Thomson et al. Reference Thomson, Elgar, Raubenheimer, Herbers and Guza2006; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015, Reference de Bakker, Tissier and Ruessink2016). 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 (Reference Herbers and Burton1997) is
\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. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006),
\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 (Reference Herbers and Burton1997), and (c,d) balance of Henderson & Bowen (Reference Henderson and Bowen2002) and Henderson et al. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006). 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 (Reference Herbers and Burton1997) and (g) Henderson & Bowen (Reference Henderson and Bowen2002) and Henderson et al. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006). 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 (Reference Henderson and Bowen2002) 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. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006) subsequently adapted this theory by assuming weak nonlinearity. The weakly nonlinear theory of Henderson et al. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006) 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 Reference Henderson and Bowen2002) 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. Reference Van Dongeren, Battjes, Janssen, Van Noorloos, Steenhauer, Steenbergen and Reniers2007; de Bakker, Tissier & Ruessink Reference de Bakker, Tissier and Ruessink2014; de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015, Reference de Bakker, Tissier and Ruessink2016). 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. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015, Reference de Bakker, Tissier and Ruessink2016). 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.
Acknowledgements
We are very grateful to J. Fiedler for her contributions during the initial phase of this work.
Funding
R.T.G. was supported by USACE and California State Parks. P.B.S. was supported by the Office of Naval Research through grants N00014-20-1-2439.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Sensitivity of energy balance to sampling in space and time
The spatial and temporal sampling required to estimate accurately the nearshore IG energy balance was explored by evaluating the energy balance terms with variable spatial (
$\Delta x=1\unicode{x2013}8$ m) and temporal (
$\Delta t = 2\unicode{x2013}10$ Hz) resolution (the SWASH numerical model output is 
$\Delta x=1$ m and 
$\Delta t = 40$ Hz). For a constant sampling resolution of 
$\Delta t=10$ Hz, the bulk IG energy flux gradient 
$\partial _x F_{IG}$ and nonlinear interaction term 
$S^{NL}_{IG}$ were insensitive to the spatial resolution (red lines in figure 9a,c). In contrast, estimates of both terms were sensitive to the temporal resolution for a fixed grid resolution of 
$\Delta x=1$ m (blue lines). Sampling errors were quantified with the nRMSE,
\begin{equation} \mathrm{nRMSE}=\frac{ \sqrt{ \left\langle \left( Q - Q_r \right)^2 \right\rangle }}{ \left|Q_r\right|_{max}}, \end{equation}
where angular brackets indicates averaging over all cross-shore output locations, 
$|\cdots |_{max}$ indicates the maximum absolute value and 
$Q$ is a term of the bulk IG energy balance (with subscript r indicating the reference results corresponding to the finest spatial and temporal sampling). Errors are generally small when space and time resolution are relatively high 
$L_p/\Delta x, T_p/\Delta t > 30$ (figure 9). Errors in IG flux gradients and nonlinear terms are largest when a high spatial resolution (
$L_p/\Delta x >20, \Delta x<5$ m) is mismatched with low temporal resolution (
$T_p/\Delta t<20, \Delta t>0.5\,{\rm s}$).
Figure 9. Sensitivity of the bulk IG energy balance terms to sampling resolution in space (
$\Delta x$) and time (
$\Delta t$) on 1/30 slope. Cross-shore variation of (a) energy flux gradient and (b) nonlinear interaction for a varying grid size given a fixed time sampling (
$\Delta t=10$ Hz, 
$T_p/\Delta t=100$, red lines), and a varying time sampling given a fixed grid size (
$\Delta x=1$ m, 
$L_p/ \Delta x= 100$, blue lines). Light to dark colours indicate coarse to fine sampling resolutions, and the result for a varying time step are mirrored around the 
$y$-axis for visualization. Panels (b) and (d) show the normalized root mean square error (nRMSE) for the full parameter space, with the spatial and temporal sampling related to the peak off-shore wavelength 
$L_p(\approx 100$ m) and wave period 
$T_p(\approx 10$ s), respectively.
Appendix B. Contributors to the nonlinear energy flux and nonlinear interaction terms
Bispectral analysis can be used to estimate the contributions of correlations between particular wave components to third-order wave statistics An alternative to the bispectrum (de Bakker et al. Reference de Bakker, Herbers, Smit, Tissier and Ruessink2015) is time-domain analysis (band-passed correlations) between filtered signals (Fiedler et al. Reference Fiedler, Smit, Brodie, McNinch and Guza2019). Using time domain band-passing, a wave or velocity signal 
$X(t)$ is decomposed into VLF, IG and SS motions:
\begin{equation} X=X_{vlf}+X_{ig}+X_{ss}. \end{equation}
Substitution of (B1) in the general cospectrum between three signal 
$X$, 
$Y$ and 
$Z$ yields
\begin{align} {C}_f \left( X Y; Z \right) &= {C}_f \left( \left(X_{vlf}+X_{ig}+X_{ss}\right) \left(Y_{vlf}+Y_{ig}+Y_{ss}\right); \left( Z_{vlf}+Z_{ig}+Z_{ss} \right)\right) \nonumber\\ &= {C}_f \left(X_{ss} \left(Y_{vlf}+Y_{ig}+Y_{ss}\right); \left( Z_{vlf}+Z_{ig}+Z_{ss} \right)\right) \nonumber\\ &\quad +{C}_f \left( X_{ig} \left(Y_{vlf}+Y_{ig}+Y_{ss}\right); \left( Z_{vlf}+Z_{ig}+Z_{ss} \right)\right) \nonumber\\ &\quad +{C}_f \left( X_{vlf} \left(Y_{vlf}+Y_{ig}+Y_{ss}\right); \left( Z_{vlf}+Z_{ig}+Z_{ss} \right)\right). \end{align}
Decomposing the nonlinear interaction term 
$S^{NL}_f$ (2.21), and collecting terms that are permutations of the same signals yields
\begin{equation} S^{{NL}}_f = S^{{NL}}_{f,{ss,ss,ss}} +S^{{NL}}_{f,{ig,ig,ig}} + S^{{NL}}_{f,{ig,ig,ss}} + S^{{NL}}_{f,{ig,ss,ss}} + S^{{NL}}_{f,{vlf}}, \end{equation}
in which 
$S^{{NL}}_{f,{ig,ig,ig}}$ and 
$S^{{NL}}_{f,{ss,ss,ss}}$ represent triad interactions between three IG and three SS wave components, respectively, 
$S^{{NL}}_{f,{ig,ig,ss}}$ represents interactions between two IG and one SS wave component, and 
$S^{{NL}}_{f,{ig,ss,ss}}$ represents interactions between a single IG and two SS wave components. Finally, 
$S^{{NL}}_{f,{vlf}}$ is short for all interactions that include at least one VLF component (the last term in (B2)). The nonlinear flux 
$F^{NL}_{f}$ can be decomposed in a similar fashion:
\begin{equation} 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{equation}Appendix C. Comparison with weakly nonlinear theory
We consider unidirectional waves over constant depth with a typical ratio 
$\delta$ between amplitude and depth that is small but finite. Neglecting friction and retaining terms up to 
$O(\delta ^3)$ in nonlinearity, the NLSWE energy balance (2.17) can be written as
\begin{equation} \frac{\partial E_f}{\partial t} + \frac{\partial F_f}{\partial x} = S^{nl}_f\end{equation}with
\begin{align} E_f &= \frac{g}{2}{C}_f(\eta,\eta) + \frac{1}{2} {C}_f (Du,u), \end{align}
\begin{align} F_{f} &= g {C}_f \left( D u ; \eta \right) + \frac{1}{2} {C}_f ( d u^2; u ), \end{align}
\begin{align} S_f^{{NL}} &= \frac{1}{2} \left[ {C}_f\left( {\rm d} u^2; \frac{\partial u}{\partial x} \right) - {C}_f\left( {\rm d} u; u \frac{\partial u}{\partial x} \right) \right] \nonumber\\ &\quad + \frac{1}{2} g \left[ {C}_f \left( \eta u; \frac{\partial\eta}{\partial x_m} \right) - {C}_f \left( u; \eta \frac{\partial\eta}{\partial x} \right) \right]. \end{align}
With waves propagating in the positive 
$x$ direction, the negative Riemann characteristic of the one-dimensional NLSWE has the constant value of 
$u-2\sqrt {gD}=-2c$ (with 
$c=\sqrt {gd}$). The nonlinear depth averaged velocity and surface elevation are related by
\begin{equation} u = 2\sqrt{gD} - 2c \approx \frac{c}{d}\left[\eta - \frac{\eta^2}{2d} + O(\delta^2) \right]. \end{equation}
Using (C5) and retaining terms through 
$O(\delta ^3)$ the energy and flux ((C2) and (C3), respectively) reduce to
\begin{gather} E_f = g {C}_f \left( \eta ; \eta \right), \end{gather}
\begin{gather} F_{f}= gc\left[ {C}_f \left( \eta ; \eta \right)+ \frac{1}{d}{C}_f \left( \eta ; \eta^2 \right) \right]= c\left[ E_f + \frac{g}{d} \int \mathcal{B}^{{Re}}_{f-f^\prime,f^\prime}\,\mathrm{d}f^\prime \right]. \end{gather}
The real part of the surface elevation bispectrum 
$\mathcal {B}^{{Re}}_{f-f^\prime,f^\prime }$ is typically neglected in weakly nonlinear balances, under the assumption mean quantities change on a slow spatial scale. Consequently gradients (
${\partial F_f}/{\partial x}$) are of order 
$O(\delta$), so the flux correction enters dynamically at 
$O(\delta ^4$) in (C1), and can be neglected (Herbers & Burton Reference Herbers and Burton1997). Lastly, to the leading order in nonlinearity and neglecting bottom slope 
$\partial _x \approx -({1}/{c})\partial _t$ and using that in the frequency domain 
$\partial _t \rightarrow \mathrm {i}2{\rm \pi} f$, yields
\begin{equation} S_f^{{NL}} = \frac{gc}{d} \left[ {C}_f \left( \eta^2; \frac{\partial \eta}{\partial x} \right) - {C}_f \left( \eta; \eta \frac{\partial \eta}{\partial x} \right) \right] ={-}\frac{ 3\omega gc}{2d} \int \mathcal{B}^{{Im}}_{f-f^\prime,f^\prime} \mathrm{d}f^\prime, \end{equation}
with 
$\mathcal {B}^{{Im}}_{f-f^\prime,f^\prime }$ the imaginary bispectrum, consistent with previous weakly nonlinear WKB balances for shoreward propagating waves (e.g. Herbers & Burton Reference Herbers and Burton1997).
Appendix D. Comparison with Henderson & Bowen (Reference Henderson and Bowen2002)
The balance considered in Henderson & Bowen (Reference Henderson and Bowen2002) can be found from combination of the frequency domain conservative momentum and mass equations alone as
\begin{equation} \lim_{\Delta f \rightarrow 0} \frac{1}{\Delta f}\langle(\overline{\tilde{u}}_m^f/2 \times \mathrm{Eq.}\;(2.12)^* + \mathrm{C.C.}) + (g\bar{\eta}^f/2 \times \mathrm{Eq.}\;(2.10)^* +\mathrm{C.C.})\rangle\!,\end{equation}
where 
$\tilde {u}_m = (({d+\eta })/{d})u_m$ is the Lagrangian (or mass transport) velocity. The result can be written as a balance of the form
\begin{equation} \frac{\partial \tilde{E}_f}{\partial t} + \frac{\partial }{\partial x_n}(\tilde{F}^{L}_{n,f}+\tilde{F}^{{NL}}_{n,f}) = \tilde{S}^{{NL}}_f + \tilde{S}^{\tau}_f, \end{equation}in which
\begin{gather}\tilde{E}_f = \frac{1}{2}g {C}_f \left( \eta; \eta \right) + \frac{d}{2} {C}_f \left(\tilde{u}_m; \tilde{u}_m \right), \end{gather}
\begin{gather}\tilde{F}^{{L}}_{n,f} = g {C}_f \left( du_n ; \eta \right), \end{gather}
\begin{gather}\tilde{F}^{{NL}}_{n,f} = g {C}_f \left( \eta u_n ; \eta \right) + {C}_f \left( D u_n u_m + \frac{\delta_{m,n}}{2}\eta^2; \tilde{u}_n \right), \end{gather}
\begin{gather}\tilde{S}_f^{{NL}} = {C}_f\left( D u_n u_m + \frac{\delta_{m,n}}{2}\eta^2; \frac{\partial \tilde{u}_n}{\partial x_m} \right), \end{gather}
\begin{gather}\tilde{S}_f^{\tau} ={-}{C}_f \left( \tilde{u}_m; \tau_m \right). \end{gather}Though qualitatively similar to the NLSWE energy balance (2.17), this is not an energy balance. Integrating (D3)–(D7) over all frequencies, and manipulating the result to resemble closely the bulk energy balance ((2.4)–(2.7)), results in
\begin{equation}
\frac{\partial}{\partial t}\left[ E +
\left\langle\frac{\eta}{d}D u_m u_m\right\rangle\right] +
\frac{\partial F_n}{\partial x_n} = S^\tau - \left\langle
\frac{\eta u_m \tau_m}{d} - \frac{\eta
u_m}{d}\frac{\partial D u_m u_n}{\partial x_n} - g
\frac{\tilde{u}_m}{2} \frac{\partial \eta^2}{\partial
x_m}\right\rangle. \end{equation}
In effect, this balance considers the time evolution of a quantity that is the sum of mean potential and a total pseudokinetic energy derived from the Lagrangian mean energy in the water column. This differs from the standard energy measure by the term in brackets on the left-hand side. The above form (spectral or bulk) is non-unique, and the right-hand side terms can be manipulated into flux-like contributions and source terms through the chain rule. This is reflected in the different forms of the balance considered in, for example, Henderson & Bowen (Reference Henderson and Bowen2002) and Henderson et al. (Reference Henderson, Guza, Elgar, Herbers and Bowen2006). However, in all cases the bulk balance takes a non-conserved form, and is generally not strictly interpretable as an energy balance. That said, for mild slopes and mild nonlinearity the balance is consistent with weakly nonlinear theory (Henderson et al. Reference Henderson, Guza, Elgar, Herbers and Bowen2006), and significant differences are only expected at high Ursell numbers.
Open access
