Dissipation of wave energy ($S_{\rm ds}$)

Whitecapping: Komen et al. (1984) formulation


The processes of whitecapping in the SWAN model is represented by the pulse-based model of Hasselmann (1974). Reformulated in terms of wave number (rather than frequency) so as to be applicable in finite water depth (cf. the WAMDI group, 1988), this expression is:

  $$S_{\rm ds,w} (\sigma,\theta) = -\Gamma \tilde{\sigma} \frac{k}{\tilde{k}} E (\sigma,\theta)$$ (2.44)

where ${\tilde{\sigma}}$ and ${\tilde{k}}$ denote the mean frequency and the mean wave number, respectively, and the coefficient $\Gamma$ depends on the overall wave steepness. This steepness dependent coefficient, as given by the WAMDI group (1988), has been adapted by Günther et al. (1992) based on Janssen (1991a) (see also (Janssen, 1991b)):

  $$\Gamma = \Gamma_{\rm KJ} = C_{\rm ds} ((1-\delta) + \delta \frac{k}{\tilde{k}}) \left(\frac{\tilde{s}}{\tilde{s}_{\rm PM}} \right)^p$$ (2.45)

For $\delta=0$ the expression of $\Gamma$ reduces to the expression as used by the WAMDI group (1988). The coefficients $C_{\rm ds}$, $\delta$ and $p$ are tunable coefficients, ${\tilde{s}}$ is the overall wave steepness, ${\tilde{s}}_{\rm PM}$ is the value of ${\tilde{s}}$ for the Pierson-Moskowitz spectrum (1964): ${\tilde{s}}_{PM} = \sqrt{3.02 \times 10^{-3}}$. The overall wave steepness ${\tilde{s}}$ is defined as

  $$\tilde{s} = \tilde{k} \sqrt{E_{\rm tot}}$$ (2.46)

The mean frequency ${\tilde{\sigma}}$, the mean wave number ${\tilde{k}}$ and the total wave energy $E_{\rm tot}$ are defined as (cf. the WAMDI group, 1988):

  $$\tilde{\sigma} = \left( E^{-1}_{\rm tot} \int_{0}^{2\pi} \int_{0}^{\infty} \frac{1}{\sigma}E(\sigma,\theta) d\sigma d\theta \right)^{-1}$$ (2.47)

  $$\tilde{k} = \left(E^{-1}_{\rm tot} \int_{0}^{2\pi} \int_{0}^{\infty} \frac{1}{\sqrt{k}}E(\sigma,\theta) d\sigma d\theta\right)^{-2}$$ (2.48)

  $$E_{\rm tot} = \int_{0}^{2\pi} \int_{0}^{\infty} E(\sigma,\theta)d\sigma d\theta$$ (2.49)

The values of the tunable coefficients $C_{\rm ds}$ and $\delta$ and exponent $p$ in this model have been obtained by Komen et al. (1984) and Janssen (1992) by closing the energy balance of the waves in idealized wave growth conditions (both for growing and fully developed wind seas) for deep water. This implies that coefficients in the steepness dependent coefficient $\Gamma$ depend on the wind input formulation that is used. Since two different wind input formulations are used in the SWAN model, two sets of coefficients are used. For the wind input of Komen et al. (1984; corresponding to WAM Cycle 3; the WAMDI group, 1988): $C_{\rm ds} = 2.36 \times 10^{-5}$, $\delta=0$ and $p=4$. Janssen (1992) and also Günther et al. (1992) obtained (assuming $p=4$) $C_{\rm ds} = 4.10 \times 10^{-5}$ and $\delta=0.5$ (as used in the WAM Cycle 4; Komen et al., 1994).


It is well-known that SWAN underestimates structurally the mean (or peak) wave periods by 10 to 20%. This has also been observed in the SWAN hindcasts as described by Rogers et al. (2003). Investigations of Rogers et al. (2003) showed that adjusting the parameter $\delta$ from 0 to 1 leads to an improved prediction of the wave energy at lower frequencies. Because of this, $\delta$ is set to 1 as default since version 40.91A. However, it should be mentioned that adapting $\delta$ without retuning $C_{\rm ds}$ may lead to exceedence of the theoretical limits on wave height proposed by Pierson and Moskowitz (1964).


Whitecapping: saturation-based model; and wind: Yan model


An alternative description for whitecapping in SWAN is given by Van der Westhuysen et al. (2007) and Van der Westhuysen (2007), which is an adapted form of the expression of Alves and Banner (2003). The latter is based on the apparent relationship between wave groups and whitecapping dissipation. This adaption is due to the fact that it can also be applied to mixed sea-swell conditions and in shallow water. This was done by removing the dependencies on mean spectral steepness and wavenumber in the original expression, and by applying source term scaling arguments for its calibration (see below). This led to the following expression for whitecapping dissipation

  $$S_{\rm ds,break} (\sigma,\theta) = -C'_{\rm ds} \left( \frac{B(k)}{B_{\rm r}} \right)^{p/2} (\tanh(kh))^{(2-p_0)/4} \sqrt{gk} E(\sigma,\theta)$$ (2.50)

in which the density function $B(k)$ is the azimuthal-integrated spectral saturation, which is positively correlated with the probability of wave group-induced breaking. It is calculated from frequency space variables as follows

  $$B(k) = \int_{0}^{2\pi} c_g k^3 E(\sigma,\theta) d\theta$$ (2.51)

and $B_{\rm r} = 1.75 \times 10^{-3}$ is a threshold saturation level. The proportionality coefficient is set to $C'_{\rm ds} = 5.0 \times 10^{-5}$. When $B(k) > B_{\rm r}$, waves break and the exponent $p$ is set equal to a calibration parameter $p_0$. For $B(k) \leq B_{\rm r}$ there is no breaking, but some residual dissipation proved necessary. This is obtained by setting $p = 0$. A smooth transition between these two situations is achieved by (Alves and Banner, 2003)

  $$p = \frac{p_0}{2} + \frac{p_0}{2} \tanh \left[ 10 \left( \sqrt{\frac{B(k)}{B_{\rm r}}} - 1 \right) \right]$$ (2.52)

In SWAN, however, $p$ is simply set to $p_0$ (see below).


In Van der Westhuysen (2007) the dissipation modes of breaking and non-breaking waves are separated, so that they are active over different parts of the spectrum:

  $$S_{\rm ds,w}(\sigma,\theta) = f_{\rm br}(\sigma)S_{\rm ds,break} + \left[ 1 - f_{\rm br}(\sigma) \right] S_{\rm ds,non-break} \ ,$$ (2.53)

where $S_{\rm ds,break}$ is the contribution by breaking waves (2.50), and $S_{\rm ds,non-break}$ dissipation by means other than breaking (e.g. turbulence). The changeover between the two modes is made with a smooth transition function $f_{\rm br}$ similar to (2.52):

  $$f_{\rm br}(\sigma) = \frac{1}{2} + \frac{1}{2} \tanh \left[ 10 \left( \sqrt{\frac{B(k)}{B_{\rm r}}} - 1 \right) \right]$$ (2.54)

Since relatively little is known about the dissipation mechanisms of the non-breaking low-frequency waves, their dissipation is not modelled in detail. Instead, the expression (2.44) is used for $S_{\rm ds,non-break}$, to provide general background dissipation of non-breaking waves. For this component, the parameter settings of Komen et al. (1984) are applied.


The wind input expression used in saturation-based model is based on that by Yan (1987). This expression embodies experimental findings that for strong wind forcing, $u_*/c > 0.1$ say, the wind-induced growth rate of waves depends quadratically on $u_*/c$ (e.g. Plant 1982), whereas for weaker forcing, $u_*/c < 0.1$ say, the growth rate depends linearly on $u_*/c$ (Snyder et al., 1981). Yan (1987) proposes an analytical fit through these two ranges of the form:

  $$\beta_{\rm fit} = D \left( \frac{u_*}{c} \right)^2 \cos(\theta-\a... ... \left( \frac{u_*}{c} \right) \cos(\theta-\alpha) + F \cos(\theta-\alpha) + H$$ (2.55)

where $D$,$E$,$F$ and $H$ are coefficients of the fit. Yan imposed two constraints:

  $$\beta_{\rm fit} \approx \beta_{\rm Snyder}\, \quad \mbox{for} \quad \frac{U_5}{c} \approx 1\,\,\, (\mbox{or} \,\, \frac{u_*}{c} \approx 0.036)$$ (2.56)

and

  $$\lim \limits_{u_*/c \rightarrow \infty} \beta_{\rm fit} = \beta_{\rm Plant}$$ (2.57)

in which $\beta_{\rm Snyder}$ and $\beta_{\rm Plant}$ are the growth rates proposed by Snyder et al. (1981) and Plant (1982), respectively. Application of Eqs. (2.56) and (2.57) led us to parameter values of $D=4.0 \times 10^{-2}$, $E=5.52 \times 10^{-3}$, $F=5.2 \times 10^{-5}$ and $H=-3.02 \times 10^{-4}$, which are somewhat different from those proposed by Yan (1987). We found that our parameter values produce better fetch-limited simulation results in the Pierson and Moskowitz (1964) fetch range thant the original values of Yan (1987).


Finally, the choice of the exponent $p_0$ in Eqs. (2.50) and (2.52) is made by requiring that the source terms of whitecapping (Eq. 2.50) and wind input (Eq. 2.55) have equal scaling in frequency, after Resio et al. (2004). This leads to a value of $p_0 = 4$ for strong wind forcing ($u_*/c > 0.1$) and $p_0 = 2$ for weaker forcing ($u_*/c < 0.1$). A smooth transition between these two limits, centred around $u_*/c = 0.1$, is achieved by the expression

  $$p_0(\sigma) = 3 + \tanh \left[ w \left( \frac{u_*}{c} - 0.1 \right) \right]$$ (2.58)

where $w$ is a scaling parameter for which a value of $w = 26$ is used in SWAN. In shallow water, under strong wind forcing ($p_0 = 4$), this scaling condition requires the additional dimensionless factor ${\tanh(kh)}^{-1/2}$ in Eq. (2.50), where $h$ is the water depth.


Whitecapping: dissipation on opposing current


When a wave field meets an adverse current with a velocity that approaches the wave group velocity, waves are blocked, which may cause steepness-induced breaking and reflection. Ris and Holthuijsen (1996) show that SWAN underestimates wave dissipation in such situations, leading to a strong overestimation in the significant wave height. Van der Westhuysen (2012) proposes a saturation-based whitecapping expression for the required enhanced dissipation. The dissipation due to current influence is taken to be proportional to the relative increase in steepness due to the opposing current, expressed in terms of the relative Doppler shifting rate $c_\sigma/\sigma$ (2.13):

  $$S_{\rm wc,curr} (\sigma,\theta) = -C''_{\rm ds} \max \left [ \fra... ...ma}, 0 \right ] \,\left( \frac{B(k)}{B_{\rm r}} \right)^{p/2} E(\sigma,\theta)$$ (2.59)

The proportionality coefficient is set to $C''_{\rm ds}$ = 0.8. The remaining parameters are given by Van der Westhuysen (2007), namely $B_{\rm r} = 1.75 \times 10^{-3}$ and $p = p_0$ according to (2.58).


Wind input, whitecapping, and non-breaking dissipation by ST6


The “ST6” source term package was implemented in an unofficial (NRL) version of SWAN starting in 2008 and initial development was documented in Rogers et al. (2012). (At the time, it was referred to as “Babanin et al. physics” rather than “ST6”.) ST6 was implemented in the official version of WAVEWATCH III(WW3) starting in 2010, and this implementation was documented in Zieger et al. (2015). Since 2010, developments in the two models have largely paralleled each other, insofar as most notable improvements are implemented in both models. As such, the documentation for WW3 (public release versions 4 or 5) is largely adequate documentation of significant changes to the source terms in SWAN since the publication of Rogers et al. (2012), and do not need to be repeated here. We point out three notable exceptions to this.


The first notable difference relates to the SSWELL ZIEGER option in SWAN, for representation of non-breaking dissipation. The steepness-dependent coefficient for this term, introduced to WW3 in version 5, has not yet been implemented in SWAN. The non-breaking dissipation instead follows ST6 in WW3 version 4. These two methods are contrasted by Zieger et al. (2015) in their equations 23 and 28.


The second notable difference is that ST6 in SWAN permits the use of the non-breaking dissipation of Ardhuin et al. (2010), the SSWELL ARDHUIN option in SWAN. This option is not available in WW3/ST6, but is instead the non-breaking dissipation used in WW3/ST4.


The third notable difference is that the wind speed scaling which was $U=28 u_\star$, following Komen et al. (1984), has been replaced with $U=S_{ws} u_\star$, where $S_{ws}$ is a free parameter. Use of $S_{ws}>28$ (we use $S_{ws}=32$) yields significant improvements to the tail level, correcting overprediction mean square slope. This necessitates tuning of the $a_1$ and $a_2$ coefficients. Settings are suggested in the SWAN User Manual. At time of writing, this feature has not yet been ported to WW3.


Other less notable differences include: changes to linear wind input (Cavaleri and Malanotte-Rizzoli 1981) in SWAN/ST6, changes to calculation of viscous stress in SWAN/ST6, and dissipation by viscosity in the water, added to SWAN/ST6.


Bottom friction


The bottom friction models that have been selected for SWAN are the empirical model of JONSWAP (Hasselmann et al., 1973), the drag law model of Collins (1972) and the eddy-viscosity model of Madsen et al. (1988). The formulations for these bottom friction models can all be expressed in the following form:

  $$S_{\rm ds,b} = -C_{\rm b} \frac{\sigma^2}{g^2 \sinh^2 kd} E (\sigma,\theta)$$ (2.60)

in which $C_{\rm b}$ is a bottom friction coefficient that generally depends on the bottom orbital motion represented by $U_{\rm rms}$:

  $$U^2_{\rm rms} = \int_{0}^{2\pi} \int_{0}^{\infty} \frac{\sigma^2}{\sinh^2 kd} E(\sigma,\theta) d\sigma d\theta$$ (2.61)

Hasselmann et al. (1973) found $C_{\rm b}=C_{\rm JON}=0.038$m$^{2}$s$^{-3}$ which is in agreement with the JONSWAP result for swell dissipation. However, Bouws and Komen (1983) suggest a value of $C_{\rm JON}~=~0.067$m$^{2}$s$^{-3}$ for depth-limited wind-sea conditions in the North Sea. This value is derived from revisiting the energy balance equation employing an alternative deep water dissipation. Recently, in Zijlema et al. (2012) it was found that a unified value of $0.038$m$^{2}$s$^{-3}$ can be used if the second order polyomial fit for wind drag of Eq. (2.36) is employed. So, in SWAN 41.01 this is default irrespective of swell and wind-sea conditions.


The expression of Collins (1972) is based on a conventional formulation for periodic waves with the appropriate parameters adapted to suit a random wave field. The dissipation rate is calculated with the conventional bottom friction formulation of Eq. (2.60) in which the bottom friction coefficient is $C_{\rm b} = C_f g U_{\rm rms}$ with $C_f = 0.015$ (Collins, 1972)^2.2.


Madsen et al. (1988) derived a formulation similar to that of Hasselmann and Collins (1968) but in their model the bottom friction factor is a function of the bottom roughness height and the actual wave conditions. Their bottom friction coefficient is given by

  $$C_{\rm b} = f_w \frac{g}{\sqrt{2}} U_{\rm rms}$$ (2.62)

in which $f_w$ is a non-dimensional friction factor estimated by using the formulation of Jonsson (1966) cf. Madsen et al. (1988):

  $$\frac{1}{4\sqrt{f_w}} + \log_{10} (\frac{1}{4\sqrt{f_w}}) = m_f + \log_{10} (\frac{a_b}{K_{\rm N}})$$ (2.63)

in which $m_f = -0.08$ (Jonsson and Carlsen, 1976) and $a_b$ is a representative near-bottom excursion amplitude:

  $$a^2_{b} = 2\int_{0}^{2\pi} \int_{0}^{\infty} \frac{1}{\sinh^2 kd} E(\sigma,\theta) d\sigma d\theta$$ (2.64)

and $K_{\rm N}$ is the bottom roughness length scale. For values of $a_b/K_{\rm N}$ smaller than 1.57 the friction factor $f_w$ is 0.30 (Jonsson, 1980).


Depth-induced wave breaking


To model the energy dissipation in random waves due to depth-induced breaking, the bore-based model of Battjes and Janssen (1978) is used in SWAN. The mean rate of energy dissipation per unit horizontal area due to wave breaking $D_{\rm tot}$ is expressed as

  $$D_{\rm tot} = - \frac{1}{4} \alpha_{\rm BJ} Q_b (\frac{\tilde{\si... ...2_{\rm max} = - \alpha_{\rm BJ} Q_b \tilde{\sigma} \frac{H^2_{\rm max}}{8\pi}$$ (2.65)

in which $\alpha_{\rm BJ} = 1$ in SWAN, $Q_b$ is the fraction of breaking waves determined by

  $$\frac{1 - Q_b}{\ln Q_b} = -8 \frac{E_{\rm tot}}{H_{\rm max}^2}$$ (2.66)

in which $H_{\rm max}$ is the maximum wave height that can exist at the given depth and ${\tilde{\sigma}}$ is a mean frequency defined as

  $$\tilde{\sigma} = E^{-1}_{\rm tot} \int_{0}^{2\pi} \int_{0}^{\infty} \sigma E(\sigma,\theta) d\sigma d\theta$$ (2.67)

The fraction of depth-induced breakers ($Q_b$) is determined in SWAN with

  $$Q_b = \left\{ \begin{array}{ll} 0 \, , & \mbox{for } \beta \le... ....2 < \beta < 1\\ \\ 1 \, , & \mbox{for } \beta \geq 1 \end{array} \right.$$ (2.68)

where $\beta = H_{\rm rms}/H_{\rm max}$. Furthermore, for $\beta \leq 0.5$, $Q_0 = 0$ and for $0.5 < \beta \leq 1$, $Q_0 = (2\beta-1)^2$.


Extending the expression of Eldeberky and Battjes (1995) to include the spectral directions, the dissipation for a spectral component per unit time is calculated in SWAN with:

  $$S_{\rm ds,br} (\sigma,\theta) = \frac{D_{\rm tot}}{E_{\rm tot}} E... ...a) = -\frac{\alpha_{\rm BJ} Q_b \tilde{\sigma}}{\beta^2 \pi} E(\sigma,\theta)$$ (2.69)

The maximum wave height $H_{\rm max}$ is determined in SWAN with $H_{\rm max} = \gamma d$, in which $\gamma$ is the breaker parameter and $d$ is the total water depth (including the wave-induced set-up if computed by SWAN). In the literature, this breaker parameter $\gamma$ is often a constant or it is expressed as a function of bottom slope or incident wave steepness (see e.g., Galvin, 1972; Battjes and Janssen, 1978; Battjes and Stive, 1985; Arcilla and Lemos, 1990; Kaminsky and Kraus, 1993; Nelson, 1987, 1994). In the publication of Battjes and Janssen (1978) in which the dissipation model is described, a constant breaker parameter, based on Miche's criterion, of $\gamma=0.8$ was used. Battjes and Stive (1985) re-analyzed wave data of a number of laboratory and field experiments and found values for the breaker parameter varying between 0.6 and 0.83 for different types of bathymetry (plane, bar-trough and bar) with an average of 0.73. From a compilation of a large number of experiments Kaminsky and Kraus (1993) have found breaker parameters in the range of 0.6 to 1.59 with an average of 0.79.


An alternative to the bore-based model of Battjes and Janssen (1978) is proposed by Thornton and Guza (1983). This model can be regarded as an alteration of Battjes and Janssen with respect to the description of the wave height probability density function. The total dissipation due to depth-induced breaking is formulated as

  $$D_{\rm tot} = - \frac{B^3 \tilde{\sigma}}{8 \pi d} \int_{0}^{\infty} H^3 \, p_b(H)\,dH$$ (2.70)

in which $B$ is a proportionality coefficient and $p_b(H)$ is the probability density function of breaking waves times the fraction of breakers, $Q_b$. Based on field observations, the wave heights in the surf zone are assumed to remain Rayleigh distributed, even after breaking. This implies that all waves will break, not only the highest as assumed by Battjes and Janssen (1978). The function $p_b(H)$ is obtained by multiplying the Rayleigh wave height probability density function $p(H)$, given by

  $$p(H) = \frac{2H}{H^2_{\rm rms}}\exp \left( - \left( \frac{H}{H_{\rm rms}} \right)^2 \right)$$ (2.71)

by a weighting function $W(H)$ defined so that $0 \leq W(H) \leq 1$, to yield

  $$p_b(H) = W(H)\,p(H)$$ (2.72)

Thornton and Guza (1983) proposed the following weighting function in which the fraction of breaking waves is independent of the wave height:

  $$W(H) = Q_b = \left( \frac{H_{\rm rms}}{\gamma d} \right)^n$$ (2.73)

with a calibration parameter $n$(=4) and a breaker index $\gamma$ (not to be confused with the Battjes and Janssen breaker index!). The integral in expression (2.70) can then be simplified, as follows

  $$\int_{0}^{\infty} H^3 \, p_b(H)\,dH = Q_b\,\int_{0}^{\infty} H^3 \, p(H)\,dH = \frac{3}{4}\sqrt{\pi}\,Q_b\,H^3_{\rm rms}$$ (2.74)

Hence,

  $$D_{\rm tot} = - \frac{3 B^3 \tilde{\sigma}}{32 \sqrt{\pi} d} \, Q_b\, H^3_{\rm rms}$$ (2.75)