<figure><img src="image_26.png" /><figcaption>**Figure 1.** Schematic representation of a deep convolutional neural network. The input layer consists of an active wave breaking candidate and has shape 256 × 256 × 3 (image height, image width and number of RGB channels). In the case of a grey-scale image, the single grey channel was triplicated. The red dashed rectangle in the analyzed image shows the region of interest which, in this particular case, was equal to the stereo video reconstruction area. The intermediary convolutional layers (or backbones) are interchangeable with different options available (see text for details). The function of the convolutional layers is to extract features from the input image by using convolutions (3 × 3 in this example) and max pooling⁴⁹ (that is, selecting the brightest pixel in a given window). The last convolutional layer is flattened (that is, turned into a one-dimensional vector) and is connected to two fully-connected (that is, multi-layer perceptron-like) layers and one final classification layer. A 50% dropout (that is, random selection of neurons in a layer) is applied after each fully connected layer. The final classification layer has one unit and uses sigmoid activation with a threshold of 0.5.</figcaption></figure>

<table><caption>**Table 1.** Data characterization summary table. *f* is the sampling frequency in Hertz, *D* is the duration of the experiment in seconds, *H*<sub>s</sub> is the significant wave height computed from the wave spectrum, *T*<sub>p1</sub> and *T*<sub>p2</sub> are, respectively, the peak wave period of the first and second spectral partitions computed following Hanson and Jensen<sup>50</sup>, *D*<sub>p</sub> is the peak wave direction, and *U*<sub>10</sub> is the wind speed measured or converted (denoted by the *) to a height of 10 m above the sea surface using Large and Pond's<sup>51</sup> formula. Tr. S. and Ts. S. are the sample size of the train and test datasets, respectively.</caption><thead><tr><th>Location</th><th>Date and time</th><th>f [Hz]</th><th>D. [min]</th><th>H<sub>s</sub> [m]</th><th>T<sub>p1</sub> [s]</th><th>T<sub>p2</sub> [s]</th><th>U<sub>10</sub>[ms<sup>-1</sup>]</th><th>D<sub>p</sub></th><th>Tr. S.</th><th>Ts. S.</th></tr></thead><tbody><tr><td>Black Sea</td><td>2011/10/01 14:18</td><td>12</td><td>07</td><td>0.3</td><td>6.20</td><td>3.10</td><td>10.7*</td><td>WSW</td><td>1000</td><td>100</td></tr><tr><td>Black Sea</td><td>2011/10/04 09:38</td><td>12</td><td>20</td><td>0.36</td><td>6.10</td><td>2.63</td><td>10.1*</td><td>WSW</td><td>1000</td><td>100</td></tr><tr><td>Black Sea</td><td>2011/10/04 11:07</td><td>12</td><td>30</td><td>0.45</td><td>6.10</td><td>3.16</td><td>12.2*</td><td>WSW</td><td>1000</td><td>100</td></tr><tr><td>Black Sea</td><td>2011/10/04 13:30</td><td>12</td><td>30</td><td>0.55</td><td>6.60</td><td>3.71</td><td>12.9*</td><td>WSW</td><td>1000</td><td>100</td></tr><tr><td>Black Sea</td><td>2013/09/22 13:00</td><td>10</td><td>15</td><td>0.66</td><td>4.30</td><td>30</td><td>8.7*</td><td>E</td><td>1000</td><td>100</td></tr><tr><td>Black Sea</td><td>2013/09/25 12:15</td><td>12</td><td>15</td><td>0.41</td><td>4.10</td><td>1.20</td><td>6.1*</td><td>N</td><td>1000</td><td>100</td></tr><tr><td>Black Sea</td><td>2013/09/30 10:20</td><td>12</td><td>15</td><td>0.65</td><td>5.70</td><td>1.80</td><td>15.2*</td><td>N</td><td>1000</td><td>100</td></tr><tr><td>Adriatic Sea</td><td>2014/03/27 09:10</td><td>12</td><td>60</td><td>1.36</td><td>5.02</td><td>-</td><td>9.9</td><td>ENE</td><td>2000</td><td>100</td></tr><tr><td>Adriatic Sea</td><td>2015/03/05 10:35</td><td>12</td><td>60</td><td>1.33</td><td>6.10</td><td>5.02</td><td>9.0</td><td>ENE</td><td>2000</td><td>100</td></tr><tr><td>Yellow Sea</td><td>2017/05/13 05:00</td><td>10</td><td>05</td><td>1.93</td><td>5.20</td><td>4.02</td><td>13.4</td><td>NW</td><td>2000</td><td>100</td></tr><tr><td>La Jument</td><td>2017/12/15 14:20</td><td>10</td><td>30</td><td>5.88</td><td>10.00</td><td>6.80</td><td>13.3</td><td>W</td><td>2000</td><td>100</td></tr><tr><td>La Jument</td><td>2018/01/03 09:39</td><td>10</td><td>30</td><td>10.03</td><td>12.80</td><td>9.30</td><td>17.9</td><td>W</td><td>2000</td><td>100</td></tr><tr><td>La Jument</td><td>2018/01/04 11:43</td><td>10</td><td>30</td><td>7.52</td><td>11.10</td><td>-</td><td>14.7</td><td>W</td><td>2000</td><td>100</td></tr><tr><td colspan="8">Total:</td><td>19,000</td><td>1300</td></tr></tbody></table>

sample (label 1, or positive) there are nine instances that were not active wave breaking (label 0, or negative). Data augmentation<sup>38</sup> (rotation, vertical and horizontal flips, and zoom) was employed during training to increase the variety of samples of the positive label.

Five state-of-the-art neural network architectures, or backbones (VGG16<sup>39</sup>, ResNet50V2<sup>40</sup>, InceptionResNetV2<sup>41</sup>, MobileNetV2<sup>42</sup> and EfficientNetB5<sup>43</sup>), were implemented and can be chosen by end users
(see our Github code repository https://github.com/caiostringari/deepwaves for usage guidance). All of these backbones make use of convolutional layers to gradually extract information from the input images. The main differences between them are how many layers they have, the size of the convolution windows, how data is normalized in each layer and how each layer connects to each other (or to previous layers in the stack). For example, VGG16 is 16 layers deep, uses 3 × 3 convolution windows and has no normalization. InceptionResNetV2 is 164 layers deep, uses a mix of 5 × 5, 3 × 3 and 1 × 1 convolution windows and uses batch normalization⁴⁴ to help avoiding overfitting. Residual Networks (namely, ResNet50V2) not only connect adjacent layers but also take into account errors (residuals) from previous layers. In general, more modern architectures (namely, Efficient-Net) are wider (mainly by having parallel convolution windows) as well as much deeper than older architectures (namely, VGG16). The final top layers of the network were fixed for all backbones and consisted of flattening the last convolutional layer, two fully-connected layers with 50% dropout⁴⁵, and a final classification layer with sigmoid activation.

The optimization step (training) was done using the *Adam*⁴⁶ implementation of the stochastic gradient descent method⁴⁷ and binary cross-entropy was used as the loss function. Note that this step must be computed using a graphics processing unit (GPU) in order to achieve feasible computation times. The models took from three (VGG16) to twelve (EfficientNetB5) hours to train using a NVIDIA GTX 1080 GPU with a batch size of sixty-four images. After the neural networks were sufficiently trained, the best performing network (VGG16, see below) was used to classify all the naïvely identified wave breaking candidates and only the events classified as active wave breaking were kept for further analyses. Although VGG16 was chosen for presenting the results, the performance of the other architectures is nearly identical to VGG16 on real-world applications. Finally, note that user-tunable parameters for the neural networks (for example, learning rate and neuron activation thresholds) were kept unchanged from their default values in the TensorFlow library⁴⁸. This implicates that more aggressive parameter optimization could improve the results presented here even further. For sake of brevity, we refer the reader to the project's code repository for guidance on how to select hyper-parameters.

The last step of the method consisted of grouping the active wave breaking events in space and time (note that at this stage bright pixels are only grouped in space). Time-discrete wave breaking events had their bounding ellipse region filled in pixel space (i, j) and were stacked in time (t) resulting in a point-cloud-like three-dimensional (i, j, t) structure. The DBSCAN³⁴ algorithm was then used to cluster these data and obtain uniquely labelled clusters. The two parameters for the DBSCAN algorithm, that is, the minimum number of samples in each cluster ($n_{min}$) and the minimum distance allowed between two samples (eps), were set to be equals to the minimum number of points inside a fitted ellipse among all ellipses ($n_{min}$) and equals to the sampling frequency in Hertz (eps). These values for $n_{min}$ and eps were constant among the analyzed datasets. Note that this final step can be replaced by any other density-based clustering algorithm or other more sophisticated algorithms such as SORT⁵² (which is also available on the project's code repository but was not used here).

Note that up to the clustering step, all calculations were done in pixel domain and “time” was obtained from sequential frame numbers. To convert between pixel and metric coordinates, the grids constructed using the stereo-video dataset available from Guimarães et al.³⁶ were used in the present study. If stereo video data are not available, such conversions could be done knowing the camera position in real-world (that is, metric) coordinates, which is usually done using surveyed ground control points³³,³⁴. Conversion from sequential frame numbers to time (in seconds) can be done by knowing the sample rate (in frames per second, that is, Hertz) for the data. The total amount of computational time required to process 20 min of raw video recorded at 10Hz with an image size of 5 megapixels is approximately two hours on a modern six-core computer (assuming that a pre-trained neural network is available). Much shorter processing times are achievable for smaller image sizes and higher number of computation threads.

**Evaluation metrics.** Due to the imbalanced characteristic of the dataset, the classification accuracy score in isolation is not an appropriate metric to evaluate the present test dataset. For instance, a classifier that guesses that all labels are negative (that is, 0) would automatically obtain a high score (≈ 90%). To properly assess the performance of the classifier, more robust metrics were defined. These metrics are defined as follows:

* True Positives (TF) and True Negatives (TN) are samples that were correctly classified. False Positives (FP) or Type I errors are false samples that were incorrectly classified as true, and False Negatives (FN) or Type II errors are true samples that were classified as false.

* Accuracy is the percentage of examples correctly classified considering both classes, that is, $Accuracy = \frac{TP+TN}{T+N}$, where T and N are the total number of positive and negative samples, respectively.

* Precision is the percentage of predicted positives that were correctly classified, that is, $Precision = \frac{TP}{TP+FP}$.

* Recall is the percentage of actual positives that were correctly classified, that is, $Recall = \frac{TP}{TP+FN}$.

* Area under the curve (AUC) is the area defined by plotting the FP rate against the TP rate (also referred to as receiver operating characteristic curve). This metric indicates the probability of a classifier ranking a random positive sample higher than a random negative sample.

All metrics described above were monitored during the training process and were individually assessed to rank model performance. The training curves shown in Fig. 2a were used to assess when overfitting, that is, decreases in the loss function value for the training dataset that were not reflected in the validation dataset, started to occur. Training epochs after overfitting started were discarded. Here we favor the AUC curves as shown in Fig. 2b to indicate better performing models because AUC is a more robust metric than the classification score and at the same time presents a smooth evolution with training epochs. Finally, a confusion matrix (or table of confusion) as shown in Fig. 2c was plotted for each model to assess Type I and Type II errors.
<figure><img src="image_27.png" /><figcaption>**Figure 2.** Examples of training curves and confusion matrix for the best overall performing model (VGG16). (a) Loss function value for training and validation data. (b) AUC value for training and validation data. In (a) and (b) the hatched area indicates the epochs after which the model started to overfit, the thick colored lines show smoothed loss or AUC values (average at every 10 epochs), and the transparent lines show raw loss or AUC values. (c) confusion matrix calculated using test data only.</figcaption></figure>

<table><caption>**Table 2.** Evaluation metrics for all tested backbone architectures. Refer to the main manuscript text for definition of evaluation metrics. The best performing models are shown in boldface. Results are sorted by AUC.</caption><thead><tr><th>Model</th><th>Accuracy</th><th>TP</th><th>FP</th><th>TN</th><th>FN</th><th>Precision</th><th>Recall</th><th>AUC</th></tr></thead><tbody><tr><td><strong>Train</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ResNetV250</td><td>0.97</td><td>1414</td><td>198</td><td>13978</td><td>280</td><td>0.877</td><td>0.835</td><td>0.989</td></tr><tr><td>VGG16</td><td>0.93</td><td>855</td><td>273</td><td>13911</td><td>831</td><td>0.758</td><td>0.507</td><td>0.943</td></tr><tr><td>InceptionResnetV2</td><td>0.927</td><td>886</td><td>359</td><td>13823</td><td>802</td><td>0.712</td><td>0.525</td><td>0.932</td></tr><tr><td>EfficientNetB5</td><td>0.772</td><td>1403</td><td>3346</td><td>10920</td><td>297</td><td>0.295</td><td>0.825</td><td>0.874</td></tr><tr><td>MobileNet</td><td>0.904</td><td>436</td><td>268</td><td>13916</td><td>1250</td><td>0.619</td><td>0.259</td><td>0.848</td></tr><tr><td><strong>Validation</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>VGG16</td><td>0.932</td><td>221</td><td>65</td><td>3478</td><td>204</td><td>0.773</td><td>0.52</td><td>0.946</td></tr><tr><td>InceptionResnetV2</td><td>0.921</td><td>190</td><td>81</td><td>3466</td><td>231</td><td>0.701</td><td>0.451</td><td>0.93</td></tr><tr><td>EfficientNetB5</td><td>0.809</td><td>353</td><td>687</td><td>2856</td><td>72</td><td>0.339</td><td>0.831</td><td>0.897</td></tr><tr><td>MobileNet</td><td>0.908</td><td>123</td><td>64</td><td>3479</td><td>302</td><td>0.658</td><td>0.289</td><td>0.878</td></tr><tr><td>ResNetV250</td><td>0.919</td><td>197</td><td>97</td><td>3450</td><td>224</td><td>0.67</td><td>0.468</td><td>0.873</td></tr><tr><td><strong>Test</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>VGG16</td><td>0.876</td><td>106</td><td>80</td><td>945</td><td>69</td><td>0.57</td><td>0.606</td><td>0.855</td></tr><tr><td>ResNetV250</td><td>0.881</td><td>95</td><td>63</td><td>962</td><td>80</td><td>0.601</td><td>0.543</td><td>0.843</td></tr><tr><td>InceptionResnetV2</td><td>0.882</td><td>91</td><td>57</td><td>968</td><td>84</td><td>0.615</td><td>0.52</td><td>0.839</td></tr><tr><td>EfficientNetB5</td><td>0.873</td><td>88</td><td>65</td><td>960</td><td>87</td><td>0.575</td><td>0.503</td><td>0.827</td></tr><tr><td>MobileNet</td><td>0.875</td><td>30</td><td>5</td><td>1020</td><td>145</td><td>0.857</td><td>0.171</td><td>0.768</td></tr></tbody></table>

## Results

**Classifier performance.** From the analysis of all training curves, confusion matrices, and from Table 2, the best performing backbone architecture during training was ResNet50V2 by a considerable margin (AUC = 0.989). These results however, did not translate to the validation dataset (AUC = 0.873). Considering only the validation data, VGG16 was the best performing backbone with AUC = 0.946. Considering only the test dataset, the best performing model was also VGG16 (AUC = 0.855). Overall, VGG16 was selected as the best performing model and the results presented in the next sections will use this backbone. Other evaluation metrics such as the accuracy score, precision, and recall closely followed the evolution of AUC with training epochs (compare Fig. 2a and b, for example). In general, as the loss value decreased, the number of false positives decreased, which made the precision and recall to increase. This behavior was consistent for all models. Given that different architectures may perform better for other datasets and further optimization could change the model ranking presented here, all pre-trained models and their training metrics evolution are made available in the project's code repository. Finally, it is worth mentioning that larger models (for example, VGG19, ResNET152 and EfficientNetB7) could achieve better results but this was not attempted there due to hardware limitations (that is, these models required more memory than what was available on the NVIDIA GTX 1080 GPU used in this study).
<figure><img src="image_28.png" /><figcaption>**Figure 3.** Example of the application of the method. (a) Results of the naïve wave breaking detection (thresholding + DBSCAN) for La Jument data (03/01/2018 09:39). Note the great amount of passive foam being detected as active wave breaking. (b) Results of active wave breaking detection using VGG16 as backbone for La Jument data (03/01/2018 09:39). Note the significant reduction in the amount of passive foam being detected. In both plots, number of clusters refers to the number of clusters identified by DBSCAN. The red dashed rectangle indicates the region of interest which, in these examples, was the same as the stereo-video reconstruction area. (c) Results of the naïve wave breaking detection (thresholding + DBSCAN) for Black sea data (04/10/2011 09:38). (d) Results of active wave breaking detection using VGG16 as backbone for Black Sea data (04/10/2011 09:38). Animations of these results are available at https://github.com/caiostringari/deepwaves. In this particular example, the image was subdivided into blocks of 256 × 256 pixels for processing. Note that identical results were seen using other architectures other than VGG16 to classify these data.</figcaption></figure>

**Test case with real-world data.** Figure 3 shows the results of the application of the best performing model architecture (VGG16) on La Jument (2018/01/03 09:39) and Black Sea (2011/10/04 09:38) data. Visual inspection of these results confirmed the ability of the neural network to correctly classify the naïvely identified wave breaking candidates and only keep instances that represented active wave breaking. Moreover, the same neural network was able to correctly classify active wave breaking events for the rogue waves seen at La Jument²⁹ and for the small wind-generated breakers seen in the Black Sea data. This result highlights the ability of the neural network to generalize well on the dataset, which is a difficult result to achieve. From the analysis of the training curves, the averaged classification error (accounting for the imbalance in the data) should be of the order of ≈ 10%, which to the authors' knowledge it was not assessed by other wave breaking detection methods. We strongly recommend that future research should report active wave breaking detection errors when used to assess or further develop models that depend on wave breaking data.
<figure><img src="image_29.png" /><figcaption>**Figure 4.** Comparison between averaged classification error for active wave breaking detection between Mironov and Dulov's (2008)³⁰ method and all the neural network architectures implemented here. The error bars represent one standard deviation from the mean. The data used for this comparison are from the 2013 Black Sea experiments described in Table 1.</figcaption></figure>

**Comparison with Mironov and Dulov (2008).** To highlight the improvements obtained by the present method, this section presents a comparison between Mironov and Dulov's (2008)³⁰ automatic active wave breaking detection method and the present method. To perform this task, data from the 2013 Black Sea experiments in Table 1 that had previously been classified and investigated in detail by Guimarães⁵⁵ were used. All the active breaking events (that is, considering all 15 min of data) detected using Mironov and Dulov's (2008)³⁰ were manually classified as true or false and compared to the labels predicted by our method. To the authors' knowledge, these data are the only currently available data that has been classified by both methods as well as manually. Figure 4 shows the result of the compassion between models. On average, our method has relatively ≈ 50% less error than Mironov and Dulov's (2008) method with an averaged absolute reduction in error of ≈ 15%. The results in Fig. 4 are also consistent with the results seen in Table 2 which showed that our model had errors in the order of ≈ 15% when considering the validation and test datasets. Note that all tested neural network architectures performed very similarly with only a slight advantage for InceptionResnetV2. It is also worth mentioning that Mironov and Dulov's (2008) method was designed and optimized to work specifically with data from the Black Sea and it is not guaranteed that it will generalize to other datasets. Our method, on the contrary, has been shown to generalize well for very distinct datasets (see Fig. 3, for example) and with further optimization could achieve even better performance.

**Wave breaking statistics.** In this section we briefly present examples of wave breaking statistics that can be directly derived from the classified wave breaking data. For brevity, the analysis is limited to data from the Black Sea because our classifier performed best at this location (classification errors < 10% considering both 2011 and 2013 experiments). Five quantities will be analysed: the wave breaking duration ($T_{br}$), the wave breaking area ($A_{br}$), the major (a) and minor (b) axis of the fitted ellipses (representative of the wave breaking lengths at their maximum during the active part of the wave breaking), and Phillips' distribution $\Lambda(c)dc$. These quantities were obtained directly from space-time clustered wave breaking events with the exception of the cumulative wave breaking area ($A_{br}$) which was calculated from the projections of pixels clustered in the first step of the method to metric coordinates. The results of this analyses are shown in Fig. 5.

The wave breaking duration ($T_{br}$) normalized by wave peak period ($T_p$), Fig. 5a) roughly followed a shifted Gamma probability density function (PDF) and had a mean value of 0.12 and most frequent value (mode) of 0.13. This result shows that the active part of the wave breaking process happens very quickly. The wave breaking area ($A_{br}$, Fig. 5b) closely followed a Pareto distribution which indicates that large wave breaking events are relatively rare in the data. The ratio between the major and minor axis of the fitted ellipses (a/b, Fig. 5c) followed a Beta PDF and had a mean of 2.5 and mode of 1.9, which indicates that the ellipses' major axis is approximately double the size of the minor axis. Assuming a negligible wave front angle, the wave breaking area scaling relation from Duncan⁵⁶ ($A_{br}/b^2$, Fig. 5d) also followed a Beta PDF and had mean of 0.1 and mode of 0.8, which is consistent with the previously reported value⁵⁶ (0.1 ± 0.01). The wave breaking area (Fig. 5e) showed a quadratic increase with wave breaking event duration, which is trivial but seems not to have been previously directly shown before. Finally, Fig. 5f shows the $\Lambda(c)dc$ distributions obtained using our method and considering that the ellipse major axis is representative of the wave breaking length. The observed distributions greatly deviated from the theoretical $c^{-6}$ relation¹². Note that all PDF fits shown here were statistically significant at the 95% confidence level using the two-tailed Kolmogorov-Smirnov test. See Table 3 for the description of the parameters of the PDFs presented in this section and the Discussion section for the contextualization of these results and the possible implications that they have for future research.
<figure><img src="image_30.png" /><figcaption>**Figure 5.** Examples of statistical properties of breaking waves that can be directly obtained from the proposed method. (a) Probability distribution of the wave breaking duration ($T_{br}$) normalized by wave peak period ($T_p$). The blue line shows the Gamma fit to the data, the purple dashed line shows the mean value (0.13) and the orange dashed line shows the mode value (0.12). (b) Probability distribution of the wave breaking area ($A_{br}$) normalized by wavelength ($\frac{g}{2\pi}T_p^2$). The blue line shows the Pareto fit to the data. (c) Probability distribution of the ratio between major and minor axis (a/b). The blue line shows the Beta fit to the data, the purple dashed line shows the mean value (2.4) and the orange dashed line shows the mode value (1.8). (d) Probability distribution of the wave breaking area scaling parameter ($A_{br}/b^2$). The blue line shows the Beta fit to the data, the purple dashed line shows the mean value (0.1), the orange dashed line shows the mode value (0.08), and the red line shows Duncan's constant value (0.11). (e) Evolution of the wave breaking area over time. The black markers show the observed values, the error bars show values binned at 0.15 s intervals, the blue line shows the quadratic regression to the data, the blue swath shows the 95% confidence interval for the regression, and $r^2$ is the correlation coefficient. (f) Phillips$^{12}$ $\Lambda(c)dc$ distributions assuming that the wave breaking speed is the same as the phase speed of the carrier wave (that is, $c = c_{br}$) grouped by wind speed. The black dashed line shows the $c^{-6}$ theoretical decay, the colored markers show the average crest length binned at $0.1 \text{ms}^{-1}$ wave speed intervals, and the colored lines show a running average with a window size of 10 bins.</figcaption></figure>

<table><caption>**Table 3.** Fitted parameters for the distributions seen in Fig. 5. All fits were done using *Scipy*$^{57}$ and, consequently, the values reported here follow *Scipy's* conventions.</caption><thead><tr><td></td><td>Distribution</td><td>Parameter 1</td><td>Parameter 2</td><td>Location</td><td>Scale</td></tr></thead><tbody><tr><td>Figure 5a</td><td>Gamma</td><td>3.001</td><td>-</td><td>0.053</td><td>0.0260</td></tr><tr><td>Figure 5b</td><td>Pareto</td><td>4.339</td><td>-</td><td>-0.007</td><td>0.007</td></tr><tr><td>Figure 5c</td><td>Beta</td><td>6.117</td><td>73.948</td><td>0.207</td><td>30.440</td></tr><tr><td>Figure 5d</td><td>Beta</td><td>6.389</td><td>40.470</td><td>-0.760</td><td>24.438</td></tr></tbody></table>

Discussion

We have presented a new method to detect active wave breaking in video imagery data that is robust and easily transferable to future studies. Overall, VGG16 was the best performing architecture which is a surprising result given that VGG is a considerably older architecture that has been superseded by more recent architectures such as ResNets and EfficientNets$^{43}$. Also surprisingly, EfficientNet was one of the worst performers considering the test dataset despite the state-of-the-art results that this architecture achieved in recent years$^{43}$. One explanation could be that given VGG16 is the model with the highest number of parameters (despite being the shallowest network), it adapted better to the characteristics of the present dataset. Another explanation could be that EfficientNet is currently overfit on the ImageNet$^{58}$ dataset that was used to train this network hence its high reported classification scores. Another aspect to take into consideration is the speed in which the models can predict on new data (inference). While VGG16 was the best performing model, its size slows down inference time. In this regard, MobileNetV2 offers real-time inference speed at 10Hz for small image sizes (for example 512 x 512px). Consequently, a model based on MobileNetV2 could be used to detect active wave breaking in an operational
<figure><img src="image_31.png" /><figcaption>Passive foam (class 0)</figcaption></figure>

<figure><img src="image_32.png" /></figure>

**Figure 6.** Results of Grad-CAM<sup>59</sup> for all unique experiment locations described in Table 1 applied to actively breaking waves (top row) and to passive foam (bottom row). (a) Actively breaking wave example recorded at Adriatic Sea (2015/03/05 10:35). (b) Actively breaking wave example recorded at the Black Sea (2011/10/04 11:07). (c) Actively breaking wave example recorded at La Jument (2018/01/03 09:39). (d) Actively breaking wave example recorded at the Yellow Sea (2017/05/13 05:00). (e) Passive foam example recorded at Adriatic Sea (2015/03/05 10:35). (f) Passive foam example recorded at the Black Sea (2011/10/04 11:07). (g) Passive foam example recorded at La Jument (2018/01/03 09:39). (h) Passive foam example recorded at the (2017/05/13 05:00). In all panels, the color scale indicates the weights of the class activation map with brighter colors showing regions of the image which the neural network used to classify each particular sample.

fashion in, for example, offshore floating wind turbines that are susceptible to wave breaking impacts and used to adjust anchor and cable settings in real-time.

From the analysis of the application of the method on real-world data, it was visually observed that small active wave breaking events were not always detected, particularly for the Black Sea data. There are two possible explanations for this error. The simplest explanation could be that this type of events is under-represented in the training dataset. The solution for this issue consists of adding more data representative of these instances to the training dataset. The second possibility is that the image size becomes too small in the deeper layers of the network which makes it impossible for the network to learn these events (see below for further discussion). A solution for this issue could be to increase the input image size but this was not attempted here due to hardware constraints (that is, memory limitation, as discussed above).

Neural networks have been historically seen as black boxes in which only the final classification outputs are relevant. There has been, however, an increase in interest to understanding how neural networks are learning. One technique that is of particular interest for the present paper is the concept of Gradient-weighted Class Activation Mapping (Grad-CAM)⁵⁹. Briefly, this technique shows which regions of a convolutional layer are more important for the results obtained by the classifier. Figure 6 shows the results of Grad-CAM applied to examples for all unique locations from Table 1 considering VGG16’s last convolutional layer. When considering only actively breaking waves (Fig. 6a to d) it is evident that VGG16 closely mimicked how a human would classify these data, that is, it directly searched for the regions of the image that corresponded to active wave breaking. In the case of passive foam (Fig. 6e to h), VGG16 seemed to use larger portions of the image but at the same time focused on the flocculent foam seen in the images as a human classifier would do. In general, these results show that our model truly learned how to classify the images and is not merely guessing the labels.

The promising results presented here indicate that the current method should be extended to an object detection framework. Such a framework would eliminate the need for the image thresholding and DBSCAN steps. This implementation could be done by using a strongly-supervised architecture such as UNet⁶⁰, or by using a weakly-supervised method derived from Grad-CAM, for example. As a final recommendation regarding the machine learning aspect of this paper, we strongly encourage future researchers to add samples to the training dataset that matches their specific research needs instead of blindly applying the provided pre-trained models.

We have also presented examples of wave breaking statistics that can be obtained using the proposed method. In general, the patterns observed here agreed with previously reported scaling factors that support the idea of
wave breaking self-similarity. For example, Fig. 5f directly showed that the scaling parameter $A_{br}/b^2$ approaches the constant 0.1 value from Duncan's⁵⁶ laboratory experiments. Another variable that showed signs of a self-similar behavior was the wave breaking area ($A_{br}$) which was very well described by a Pareto distribution and presented a steady quadratic increase with wave breaking duration ($T_{br}$). Extensive research¹¹,¹⁵ has been grounded on the assumption that wave breaking is self-similar, but inconsistencies with this approach have been reported before¹⁴. Contrarily to other studies¹¹,³¹, however, the $\Lambda(c)$ distributions obtained here did not closely match the theoretical $c^{-6}$ for a sea-state in equilibrium. As reported before¹³, these differences may be due to the fact that here we only considered actively breaking waves for our analysis whereas other studies seem to track the speeds of both actively breaking waves and passive foam³¹,³². Another possibility is that other phenomena not accounted by Phillips'¹² theory (for example, wave modulation⁶¹) play important role on wave breaking. A future publication that investigates the mechanisms related to the wave breaking kinematics using the method and data obtained here will soon follow.

## Conclusion

We described a novel method to detect and classify actively breaking waves in video imagery data. Our method achieved promising results when assessed using several different metrics. Further, by analyzing the deeper layers of our neural network, we showed that the model mimicked how a human classifier would perform a similar classification task. As an application of the method, we presented wave breaking statistics and breaking wave crest length distributions. Our method can thus be useful for the investigation of several ocean-atmosphere interaction processes and should, in the future, lead to advancements in operational wave forecast models, gas-exchange models, and general safety at seas. Finally, we strongly recommend that future research should focus on standardized methods to study wave breaking so that a consistent dataset can be generated and used freely and unambiguously by the community.

## Data availability

Data, pre-trained networks and auxiliary programs necessary to reproduce the results in this paper are available at Zenodo (https://doi.org/10.5281/zenodo.4495599) and at our Github repository (https://github.com/caios-tringari/deepwaves).

Received: 18 September 2020; Accepted: 28 January 2021

Published online: 11 February 2021

## References

1. Battjes, J. A. & Janssen, J. Energy loss and set-up due to breaking of random waves. *Coast. Eng.* **32**, 569–587 (1978).

2. Thornton, E. B. & Guza, R. T. Transformation of wave height distribution. *J. Geophys. Res.* **88**, 5925–5938 (1983).

3. Banner, M. L., Babanin, A. V. & Young, I. R. Breaking probability for dominant waves on the sea surface. *J. Phys. Oceanogr.* **30**, 3145–3160. https://doi.org/10.1175/1520-0485(2000)030<3145:BPFDWO>2.0.CO;2 (2000).

4. Banner, M. L., Gemmrich, J. R. & Farmer, D. M. Multiscale measurements of ocean wave breaking probability. *J. Phys. Oceanogr.* **32**, 3364–3375. https://doi.org/10.1175/1520-0485(2002)032<3364:MMOOWB>2.0.CO;2 (2002).

5. Cavaleri, B. Y. L. Wave modeling: where to go in the future. *Bull. Am. Meteorol. Soc.* **87**, 207–204. https://doi.org/10.1175/BAMS-87-2-207 (2006).

6. The Wavewatch Development Group (WW3DG). User manual and system documentation of WAVEWATCH III R version 6.07. Technical Report, NOAA/NWS/NCEP/MMAB, College Park, MD, USA (2019).

7. Boooij, N., Ris, R. C. & Holthuijsen, L. H. A third-generation wave model for coastal regions 1. Model description and validation. *J. Geophys. Res.* **104**, 7649–7666 (1999).

8. The WANDI Group. The WAN model: a third generation ocean wave prediction model. *J. Phys. Oceanogr.* **18**, 1775–1810 (1988).

9. Filipot, J. F., Ardhuin, F. & Babanin, A. V. A unified deep-to-shallow water wave-breaking probability parameterization. *J. Geophys. Res. Oceans* **115**, 1–15. https://doi.org/10.1029/2009JC005448 (2010).

10. Stringari, C. E. & Power, H. E. The fraction of broken waves in natural surf zones. *J. Geophys. Res. Oceans* **124**, 1–27. https://doi.org/10.1029/2019JC015213 (2019) (*arXiv:1904.06821v1*).

11. Melville, W. K. & Matusov, P. Distribution of breaking waves at the ocean surface. *Nature* **417**, 58–63. https://doi.org/10.1038/41705-8a (2002).

12. Phillips, O. M. Spectral and statistical properties of the equilibrium range in wind-generated gravity waves. *J. Fluid Mech.* **156**, 505–531. https://doi.org/10.1017/S0022112085002221 (1985).

13. Banner, M. L., Zappa, C. J. & Gemmrich, J. R. A note on the Phillips spectral framework for ocean whitecaps. *J. Phys. Oceanogr.* **44**, 1727–1734. https://doi.org/10.1175/JPO-D-13-0126.1 (2014).

14. Gemmrich, J. R., Banner, M. L. & Garrett, C. Spectrally resolved energy dissipation rate and momentum flux of breaking waves. *J. Phys. Oceanogr.* **38**, 1296–1312. https://doi.org/10.1175/2007JPO3762.1 (2008).

15. Romero, L. Distribution of surface wave breaking fronts. *Geophys. Res. Lett.* **46**, 10463–10474. https://doi.org/10.1029/2019GL0834-08 (2019).

16. Yurovsky, Y. Y., Kudryavtsev, V. N., Chapron, B. & Grodsky, S. A. Modulation of Ka-band doppler radar signals backscattered from the sea surface. *IEEE Trans. Geosci. Remote Sens.* **56**, 2931–2948. https://doi.org/10.1109/TGRS.2017.2787459 (2018).

17. Huang, L., Liu, B., Li, X., Zhang, Z. & Yu, W. Technical evaluation of Sentinel-1 IW mode cross-pol radar backscattering from the ocean surface in moderate wind condition. *Remote Sens.* **9**, 1–21. https://doi.org/10.3390/rs9080854 (2017).

18. Hwang, P. A., Zhang, B., Toporkov, J. V. & Perrie, W. Comparison of composite Bragg theory and quad-polarization radar backscatter from RADARSAT-2: with applications to wave breaking and high wind retrieval. *J. Geophys. Res. Oceans* **115**, 1–12. https://doi.org/10.1029/2009JC005995 (2010).

19. Monahan, E. C. Oceanic whitecaps: sea surface features detectable via satellite that are indicators of the magnitude of the air-sea gas transfer coefficient. *J. Earth Syst. Sci.* https://doi.org/10.1007/BF02701977 (2002).

20. Reul, N. & Chapron, B. A model of sea-foam thickness distribution for passive microwave remote sensing applications. *J. Geophys. Res. C Oceans* **108**, 1–19. https://doi.org/10.1029/2003jc001887 (2003).

21. Carini, R. J., Chickadel, C. C., Jessup, A. T. & Thompson, J. Estimating wave energy dissipation in the surf zone using thermal infrared imagery. *J. Geophys. Res. Oceans* **120**, 3937–3957. https://doi.org/10.1002/2014JC010561.Received (2015).
22. Wang, S. et al. Improving the upper-ocean temperature in an ocean climate model (FESOM 1.4): shortwave penetration versus mixing induced by nonbreaking surface waves. *J. Adv. Model. Earth Syst.* **11**, 545–557. https://doi.org/10.1029/2018MS001494 (2019).

23. Komori, S. et al. Laboratory measurements of heat transfer and drag coefficients at extremely high wind speeds. *J. Phys. Oceanogr.* **48**, 959–974. https://doi.org/10.1175/JPO-D-17-0243.1 (2018).

24. Buscombe, D. & Carini, R. J. A data-driven approach to classifying wave breaking in infrared imagery. *Remote Sens.* **11**, 1–10. https://doi.org/10.3390/RS11070859 (2019).

25. Buscombe, D., Carini, R. J., Harrison, S. R., Chickadel, C. C. & Warrick, J. A. Optical wave gauging using deep neural networks. *Coast. Eng.* **155**, 103593. https://doi.org/10.1016/j.coastaleng.2019.103593 (2020).

26. Kim, J., Kim, J., Kim, T., Huh, D. & Caires, S. Wave-tracking in the surf zone using coastal video imagery with deep neural networks. *Atmosphere* **11**, 1–13. https://doi.org/10.3390/atmos11030304 (2020).

27. Bieman, J. P. D., Ridder, M. P. D. & Gent, M. R. A. V. Deep learning video analysis as measurement technique in physical models. *Coast. Eng.* **158**, 103689. https://doi.org/10.1016/j.coastaleng.2020.103689 (2020).

28. Zheng, C. W., Chen, Y. G., Zhan, C. & Wang, Q. Source tracing of the swell energy: a case study of the Pacific Ocean. *IEEE Access* **7**, 139264–139275. https://doi.org/10.1109/ACCESS.2019.2943903 (2019).

29. Filipot, J.-F. et al. La Jument lighthouse: a real-scale laboratory for the study of giant waves and their loading on marine structures. *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **377**, 20190008. https://doi.org/10.1098/rsta.2019.0008 (2019).

30. Mironov, A. S. & Dulov, V. A. Detection of wave breaking using sea surface video records. *Meas. Sci. Technol.* https://doi.org/10.1088/0957-0233/19/1/015405 (2008).

31. Sutherland, P. & Melville, W. K. Field measurements and scaling of ocean surface wave-breaking statistics. *Geophys. Res. Lett.* **40**, 3074–3079. https://doi.org/10.1002/grl.50584 (2013).

32. Kleiss, J. M. & Melville, W. K. Observations of wave breaking kinematics in fetch-limited seas. *J. Phys. Oceanogr.* **40**, 2575–2604. https://doi.org/10.1175/2010JPO4383.1 (2010).

33. Bradski, G. The OpenCV library. Dobb J. Softw. Tools. https://github.com/opencv/opencv/wiki/CiteOpenCV (2000).

34. Ester, M., Kriegel, H.-P., Sander, J. & Xu, X. Density-based clustering methods. *Compr. Chemometr.* **2**, 635–654. https://doi.org/10.1016/B978-044452701-1.00067-3 (1996).

35. Mostagh, N. Minimum volume enclosing ellipsoids. *Convex Optim.* **111**, 112–118 (2005).

36. Guimarães, P. V. et al. A data set of sea surface stereo images to resolve space-time wave fields. *Sci. Data* **7**, 1–12. https://doi.org/10.1277/af599f42-2770-4d6d-8209-13f40e2c292f (2020).

37. Hastie, T., Tibshirani, R. & Friedman, J. *The Elements of Statistical Learning* 2nd edn. (Springer International Publishing, Berlin, 2009).

38. Shorten, C. & Khoshgoftaar, T. M. A survey on image data augmentation for deep learning. *J. Big Data* https://doi.org/10.1186/s40537-019-0197-0 (2019).

39. Simonyan, K. & Zisserman, A. Very deep convolutional networks for large-scale image recognition. In *3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings* 1–14 (2015). arXiv:1409.1556v6.

40. He, K., Zhang, X., Ren, S. & Sun, J. Identity mappings in deep residual networks. *Lecture Notes in Computer Science (including sub-series Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)* **9908 LNCS**, 630–645, https://doi.org/10.1007/978-3-319-46493-0_38 (2016). arXiv:1603.05027.

41. Szegedy, C., Ioffe, S., Vanhoucke, V. & Alemi, A. A. Inception-v4, inception-ResNet and the impact of residual connections on learning. In *31st AAAI Conference on Artificial Intelligence*, AAAI 2017 4278–4284 (2017). arXiv:1602.07261.

42. Howard, A. G. et al. MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications (2017). arXiv:1704.04861.

43. Tan, M. & Le, Q. V. EfficientNet: Rethinking model scaling for convolutional neural networks. In *36th International Conference on Machine Learning, ICML 2019 2019-June*, 10691–10700 (2019). arXiv:1905.11946.

44. Ioffe, S. & Szegedy, C. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In *32nd International Conference on Machine Learning, ICML 2015* **1**, 448–456 (2015) arXiv:1502.03167.

45. Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. & Salakhutdinov, R. Dropout: a simple way to prevent neural networks from overfitting. *J. Mach. Learn. Res.* **15**, 1929–1958. https://doi.org/10.1109/ICAEES.2016.7888100 (2014).

46. Kingma, D. P. & Ba, J. Adam: A Method for Stochastic Optimization. In *3rd International Conference for Learning Representations*, 1–15, http://doi.acm.org.ezproxy.lib.ucf.edu/10.1145/1830483.1830503 (San Diego, California, 2014). arXiv:1412.6980.

47. Ruder, S. An overview of gradient descent optimization algorithms. arXiv:1609.04747 (2016).

48. Abadi, M. et al. TensorFlow: Large-scale machine learning on heterogeneous systems (2015). Software available from tensorflow.org.

49. Nagi, J. et al. Max-Pooling Convolutional Neural Networks for Vision-based Hand Gesture Recognition. In *2011 IEEE International Conference on Signal and Image Processing Applications, ICSIPA 2011*, 342–347 (2011).

50. Hanson, J. L. & Jensen, R. Wave system diagnostics for numerical wave models. In *8th International Workshop on Wave Hindcasting and Forecasting, Oahu, Hawaii, November* (2004).

51. Large, W. G. & Pond, S. Open ocean momentum flux measurements in moderate to strong winds. *J. Phys. Oceanogr.* https://doi.org/10.1175/1520-0485(1981)0112.0.CO;2 (1981).

52. Bewley, A., Ge, Z., Ott, L., Ramos, F. & Upcroft, B. Simple online and realtime tracking. In *Proceedings - International Conference on Image Processing, ICIP 2016-August*, 3464–3468, https://doi.org/10.1109/ICIP.2016.7533003(2016). arXiv:1602.00763.

53. Holman, R. A. & Stanley, J. The history and technical capabilities of Argus. *Coast. Eng.* **54**, 477–491. https://doi.org/10.1016/j.coastaleng.2007.01.003 (2007).

54. Schwendeman, M., Thomson, J. & Gemmrich, J. R. Wave breaking dissipation in a Young Wind Sea. *J. Phys. Oceanogr.* **44**, 104–127. https://doi.org/10.1175/JPO-D-12-0237.1 (2014).

55. Guimarães, P. V. *Sea surface and energy dissipation*. Ph.D. thesis, Université de Bretagne Loire (2018).

56. Duncan, J. H. An experimental investigation of breaking waves produced by a Towed Hydrofoil. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **377**, 331–348. https://doi.org/10.1098/rspa.1981.0127 (1981).

57. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in python. *Nat. Methods* https://doi.org/10.1038/s41592-019-0686-2 (2020).

58. Russakovsky, O. et al. ImageNet large scale visual recognition challenge. *Int. J. Comput. Vis.* **115**, 211–252. https://doi.org/10.1007/s11263-015-0816-y (2015) (arXiv:1409.0575).

59. Selvaraju, R. R. et al. Grad-CAM: visual explanations from deep networks via gradient-based localization. *Int. J. Comput. Vis.* **128**, 336–359. https://doi.org/10.1007/s11263-019-01228-7 (2020) (arXiv:1610.02391).

60. Ronneberger, O., P. Fischer & Brox, T. U-net: Convolutional networks for biomedical image segmentation. In *Medical Image Computing and Computer-Assisted Intervention (MICCAI)*, vol. 9351 of LNCS, 234–241 (Springer, 2015). (available on arXiv:1505.04597 [cs.CV]).

61. Donelan, M. A., Haus, B. K., Plant, W. J. & Troianowski, O. Modulation of short wind waves by long waves. *J. Geophys. Res. Oceans* **115**, 1–12. https://doi.org/10.1029/2009JC005794 (2010).
**Acknowledgements**

This work benefited from France Energies Marines and State financing managed by the National Research Agency under the Investments for the Future program bearing the reference numbers ANR-10-IED-0006-14 and ANR-10-IEED-0006-26 for the projects DiME and CARAVELE.

**Author contributions**

C.E.S. idealized and implemented the method, C.E.S. drafted the manuscript, P.V.G., F.L., J.F.F. and R. D. conducted the field experiments, C.E.S. and P.V.G. analyzed the data, J.F.F. provided funding. All authors reviewed the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

Correspondence and requests for materials should be addressed to C.E.S.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_33.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021
# Chapter 4

## Improvement of spectral wave models: wave energy dissipation parameterizations and validation in coastal water

### 4.1 Introduction

#### 4.1.1 The spectral wave models

Spectral wave models are the only tool able to simulate the wave conditions at a given offshore wind site in real time and over a period long enough (typically 10 to 20 years) to be useful for the definition of the current, near-future (10-day horizon), average and extreme wave conditions, required for the operation and design of offshore wind turbines.

The main reason for this is that they solve the evolution in time and space of the wave energy spectrum. This statistical variable carries the immense advantage to evolve much slower in time (and space in deep water) than the individual waves themselves. This allows large numerical time (and space in deep water) steps, resulting in a relatively low computational cost.

More specifically, spectral wave models solve the evolution equation of the wave energy spectrum [Gelci et al., 1957]:

$$ \frac{dE(f, \theta)}{dt} = S_{atm}(f, \theta) + S_{nl}(f, \theta) + S_{bk}(f, \theta) + S_{bot}(f, \theta) + \dots \quad (4.1) $$

This reflects that following a wave packet in physical and spectral spaces, the spectral wave energy density, $E$, evolves due to wind input ($S_{atm}$), energy exchanges between spectral components ($S_{nl}$), wave breaking dissipation ($S_{bk}$) and bottom friction dissipation ($S_{bot}$) and other possible source terms out of the scope of the present discussion.
One of the main difficulty of the practical implementation of equation 4.1 re-
sides in the lack of knowledge regarding the physics underlying the source terms.
Ardhuin et al. [2010] gives a detailed overview of the source term parameteriza-
tions and of the empirical adjustments required to optimize the model accuracy.
Their effort to develop the WAVEWATCHIII (hereinafter WW3) model still rep-
resents the state of the art of this discipline.

In storm conditions and in shallow water, wave breaking generally dominates
the dissipation of wave energy. Therefore, in theory, improving its parameteri-
zation may improve the accuracy of the spectral models. Being able to extract
meaningful wave breaking statistics from spectral models is also critical for the
definition of the design sea-state (if slamming loads are the design loads), for
air-sea coupled models and for wind retrieval for space-borne sensors.

In shallow water, pioneering efforts of Battjes and Janssen [1978] and Thorn-
ton and Guza [1983] have inspired the development of breaking dissipation source
terms coupling the parameterization of two physical variables : the breaking wave
height distribution and the dissipation rate per unit crest length.

This sound approach was applied only decades later to the "whitecapping"
(deep water wave breaking) dissipation that is traditionally parameterized through
an additional source term. Indeed, in the first spectral wave models, whitecapping
dissipation was initially treated as a tuning knob loosely related to the average
sea states steepness [WAMDI Group, 1988]. One of the reasons for this, was the
absence of reliable observations and hence parameterizations for the statistics of
breaking waves. In the early 2000's, several research efforts on the development of
new breaking statistics models [Banner et al., 2000; Babanin et al., 2001; Filipot
et al., 2010], motivated the development of more physically-related source terms
applicable to deep water only [Alves and Banner, 2003; Van der Westhuysen,
2009; Banner and Morison, 2010] or tackling wave breaking as a unique process
from deep-to-shallow water [Filipot and Ardhuin, 2012].

The improvements of the wave breaking source terms behaviour in deep water
has motivated many studies in the recent years, by Rogers et al. [2012], Zieger
et al. [2015], Stopa et al. [2016], Liu et al. [2019] and Romero [2019]. In these
papers, the wave breaking dissipation parameterization again relies on saturation-
type approaches, derived from the empirical results from Banner et al. [2000]
and Banner et al. [2002], relating wave breaking statistics to the local spectral
steepness (see Chapter 3).

In shallow water, most of the recent efforts were produced by the SWAN
(Simulating WAves Nearshore) development team through a series of papers led
by van der Westhuysen [van der Westhuysen, 2010, 2012] and Salmon [Salmon
and Holthuijsen, 2015; Salmon et al., 2015]. They introduced dependencies on
the spectral biphase, bottom slope and non dimensional water depth. Their wave
breaking source terms were however only indirectly validated through comparison
with observations of standard integrated wave parameters.

The next section presents efforts conducted from my PhD thesis up to now
to improve the parameterization of the wave breaking-induced dissipation. The
results may be useful to :

1. Improve wave forecasting and hindcasting at offshore wind farms

2. Develop appropriate spectral models to be coupled to hydrodynamic and sediment transport models at shallow water offshore wind farm (and power cable landing zone)

3. Provide accurate breaking statistics that can be used in the design procedure of offshore wind turbines or to feed Ocean-Wave-Atmospheric Coupled models.

# 4.2 Development of spectral wave models dissipation source terms

## 4.2.1 Development of the WW3 model

Here we shall detail research efforts carried out to improve the WW3 spectral wave model [Tolman et al., 2009] developed at NOAA/NCEP, through implementation of more accurate and physical wave breaking dissipation source terms.

As mentioned earlier, the shallow water approach from Thornton and Guza [1983], although almost 40-year old, appears as a solid physical basis to develop the wave breaking source terms in wave models. In deep water, a major difficulty arises because wave fields are generally made of different wave scales that travel at different phase speeds and with possibly different nonlinearity levels, hence breaking occurrences and dissipation rates. It seems therefore wise to attempt to separate these wave scales to compute, for each of them, their breaking statistics and their related energy dissipation rate. The definition of the wave scales can however be quite arbitrary. Longuet-Higgins [1985] mentioned that band-pass filtering the wave spectrum between [0.5, 1.5]f_p allows to conserve the group structure of the wave field, governed by the so-called "dominant wave" scale. Later, Banner et al. [2000] and Babanin et al. [2001] revised the dominant wave bandwidth to [0.7, 1.3]f_p. Applying this concept to the entire spectrum, Filipot et al. [2010] developed a spectral wave breaking probability parameterization inspired from the definition of Thornton and Guza [1983] extended to deep water and multiple wave scales.

This breaking probability parameterization of Filipot et al. [2010] recalled in Chapter 3 was used for deriving a new wave breaking dissipation source term, introduced by Filipot and Ardhuin [2012] and implemented in WW3. The source term is based on a spectral dissipation rate per unit area, for each wave scale $f_s$ as:
$$D(f_s) = Q(f_s)\epsilon(f_s)\Pi(f_s) \qquad (4.2)$$

With $Q(f_s)$ the wave breaking probability model from Filipot et al. [2010], $\epsilon(f_s)$ the dissipation rate per unit crest length and $\Pi(f_s)$ the crest length density per unit area.

Parameterizing the term $\epsilon(f_s)$ represents an important challenge, most of all because measuring the wave breaking energy dissipation rate of real waves is nearly impossible. In Filipot and Ardhuin [2012] (as in Chawla and Kirby [2002]) an analogy with the bore model extended to deep water was applied. Recent advances regarding the observations of nonlinear waves and dissipation in the field, brought by Grasso et al. [2012]; Martins et al. [2018, 2020] provide guidance for improving further the parameterization of $\epsilon$ as will be discussed later in this document.

This effort by Filipot and Ardhuin [2012], implemented in the source term package named TEST500 of the WW3 model, has been used by the National Oceanic and Atmospheric Administration (NOAA) for forecasting waves over the Great Lakes in North America. Its capabilities to capture storm sea states, that is critical for the survivability of offshore wind turbines was examined, in the frame of the DiMe project by Ruju et al. [2020] and will be detailed in next section.

### 4.2.2 Evaluation of WW3 source terms accuracy in extreme conditions

A recent paper by Ruju et al. [2020] (see Appendix 4.A) investigates the effect of the different source term packages and forcing wind fields on the wave hindcasting in storm conditions using the WW3 model. This work focuses on the winter of 2013-2014, known for the frequency and intensity of storms that struck the North East Atlantic region. It uses satellite and wave data for estimating model performance offshore and in coastal area, respectively. The model is based on a grid covering the globe with a regular mesh of 0.5 deg. horizontal resolution, to capture waves in the offshore domain. It forces an unstructured grid covering the continental shelf and allowing to represent more finely the coastal processes affecting the transformation of waves (e.g.: bottom friction, bathymetric and current-induced refraction).

The forcing wind fields come from a blending of different satellite-derived wind fields from Bentamy et al. [2019] and from the ERA interim reanalysis [Hersbach et al., 2020]. This work concludes that, among the different combinations of source term packages and forcing wind sources, T500 with ERA-interim represents the best alternative to capture storms sea states, confirming the results from Filipot and Ardhuin [2012] (see Figure 4.2.1).

Further improvements of the wave breaking source terms based on work carried out at FEM are introduced in Chapter 6.
<figure><img src="image_34.png" /><figcaption>Figure 4.2.1: Left panel: normalized root mean square error for $H_s$, computed through a comparison with altimeter observations over the entire 2007-year. The top panel shows the results using the source terms packages TEST441 [Ardhuin et al., 2010], the middle one with TEST500 [Filipot and Ardhuin, 2012] and the bottom one using Bidlot et al. [2005] source terms. Right panel : $H_s$ NRMSE of the three source term packages. Illustrations reproduced from Filipot and Ardhuin [2012].</figcaption></figure>

### 4.2.3 Development of the SWAN model

The SWAN model, developed by TU Delft, is the most widely used spectral model for coastal and nearshore applications. Most of the research effort regarding its capacities in representing nearshore wave fields has been devoted to simulate waves over beaches with mild or moderate slopes and with rather low bottom roughness [e.g. Booij et al., 1997; Castelle et al., 2007]. However, there are different environments where the bottom slope and roughness departs significantly from the situations generally investigated with the SWAN model. Tropical islands, as those of the Hawaiian Archipelago are part of them. They are generally surrounded by steep reefs covered by a rough substrate, either made of coral or lava rocks. When the reef barrier is connected to the shore by a lagoon, as in the Hawaiian islands, the reef is referred to as a fringing-reef. In such environment, the wave transformation is largely affected by the strong bottom slope that affects the nonlinear wave evolution and breaking processes. Furthermore, the strong bottom roughness results in high bottom friction dissipation rates that can commensurate with those induced by breaking, or even exceeds them over the down-wave face of the coral barrier (i.e. over the lagoon). This statement has motivated the development of SWAN, described in [Filipot and Cheung, 2012] (see Appendix 4.B), to address the challenging fringing reef environment. The corresponding effort focuses on the adaptation and development of wave breaking and bottom friction dissipation source terms.
One of the main challenge faced by *Filipot and Cheung* [2012] is that it is impossible to isolate the effect of wave breaking and bottom friction on the wave energy dissipation rates observed in the field. Bringing improvements on the bottom friction and wave energy dissipation, using field data for calibration is therefore highly speculative. To circumvent this problem, we relied on wave tank experiments, carried out at U. Michigan by *Demirbilek et al.* [2007] to study the evolution of waves above a steep bathymetry representative of the Guam tropical Island. Because the wave tank bottom was made of a smooth plastic material, this dataset was used to examine the behaviour of the *SWAN* model to capture wave breaking dissipation over steep bottom slopes, assuming that the bottom friction was negligible. Since no observation of breaking statistics was available from the wave tank experiment, the work had merely consisted in adjusting the existing wave breaking parameterization by *Battjes and Janssen* [1978] already implemented in *SWAN*. The analysis shows that the source term used with the default parameters generally employed for sandy beaches with mild slopes [*Battjes and Janssen, 1978*] fails in reproducing the nonlinear wave evolution and breaking over such steep slope. An optimization procedure, led to increase the breaker index parameter, and dissipation coefficient by 20 and 9%, respectively. These empirical adjustments are consistent with the observations reported in the literature [e.g. *Svendsen, 1987*; *Tsai et al., 2005*], suggesting that steeper bottom slopes tend to cause breaking onset to occur in shallower water and generate steeper (plunging) breakers with higher dissipation rate [*Roeber et al., 2010*; *Derakhti et al., 2018*].

On top of wave breaking (as already introduced), it is well acknowledged that coral or rocky reef can be an important source of wave energy dissipation in shallow water. *Lowe et al.* [2005] studied the bottom-friction induced wave energy dissipation above the Kaneohe Bay reef, in Oahu, Hawaii. They further used their data to develop a new bottom friction dissipation source term, extending the work of *Madsen* [1995]. They were able to show that the adjustment of the bottom roughness scale in their model led to a value consistent with observations at their study site. These findings, motivated *Filipot and Cheung* [2012] to implement their parameterization into *SWAN*. The corresponding parameterization version, was further calibrated using field observations collected at Mokuleia beach fringing reef, located on the North Shore of Oahu island, Hawaii.

Since then, this *SWAN* version has been used by research teams of the University of Hawaii to study the wave processes above fringing reefs [e.g. *Li et al., 2016*].
## 4.3 Methodology to validate the ability of spectral models to capture areas of high and potentially breaking waves

In coastal areas, uneven bathymetries cause waves to focus through the well-known refraction effects. The areas of energy focusing can be seen as hot spots for wave energy harvesting or as dangerous regions by the offshore wind industry. The latter will try to avoid these locations because of the potential presence of high and/or breaking waves and the subsequent increased wave loads on the structures. Capturing this wave energy variability at the scales of interest; $O(10-100m)$, over a domain large enough, $O(1-10km)$ was beyond any observational systems so far. The only available tool was the numerical wave models, but because of the lack of appropriate spatial observations, their ability to capture the wave variability was uncertain. The paper from Varing et al. [2021a] (Appendix 4.C) presents a methodology to use Synthetic Aperture Radar (SAR) satellite observations to capture the spatial variability of the wave energy induced by bathymetric refraction. One key finding of the paper lies in the approach developed to enable a consistent comparison between numerical models outputs and SAR observations. Indeed, SAR images give snapshot of the normalized radar cross section from which a map of the instantaneous wave amplitude can be inferred, after application of a spectral filter [described in Varing et al., 2021a]. However, the wave field resulting from this procedure involves the effect of wave groups which cannot be accounted for, neither by the Boussinesq model BOSZ since the phases of the wave components inputed at the boundary are randomly drawn nor by the spectral model SWAN, that propagates directly the wave energy spectrum, losing hence the wave phase information. To go around these difficulties, the wave groupiness effect was filtered out by applying an ensemble-like averaged using a set of SAR images corresponding to wave conditions assumed to be sufficiently alike. This study provides two major results:

1. For the studied site, the spectral and phase-resolving models disagree on the spatial energy gradients location and magnitude

2. SAR images may bring insightful results to identify regions of potential dangerous breaking waves and help qualifying the limitations of the numerical wave models

This suggests that the classical procedure (based on the implementation of spectral models) to investigate the wave conditions at future bottom-fixed offshore wind farms should be reconsidered if the up-wave bathymetry is strongly uneven. SAR images in complement to numerical models information should be exploited to optimize the siting of the wind turbines inside the farms. Further thoughts and potential research axes are developed in Chapter 6.
## 4.4 Discussions & conclusions

Thanks to their computational efficiency, spectral wave models are still valuable tools for characterizing the wave climate at offshore renewable energy sites. Their main drawbacks reside in their often-empirical parameterizations of the different physical processes that control wave growth and dissipation. We further note that this limitation is also shared with phase-resolving models for most of the physical processes of interest, including wave breaking dissipation that remains a common blocking point. The research efforts reported here to improve the wave breaking dissipation source terms, mostly relied on better constraining the wave breaking statistics, should be complemented with research addressing the energy dissipation by unit crest length, which is a major unknown of the problem at hand. For this, the inspiring research of Derakhti et al. [2018] and Martins et al. [2020] should be further exploited in the aim of improving the wave breaking source terms in spectral models. Revisiting the work of Thornton and Guza [1983] at the light of these recent results, should offer a significant potential of innovation for shallow water spectral wave modeling. Again, here complementary observations of breaking waves in shallow water, combining breaking wave height distribution together with dissipation rate estimation, following Martins et al. [2020] could be a first important step in this direction (see additional discussion in Chapter 6). Further efforts should also be directed to develop new observational means to capture the coastal and nearshore wave field variability for spectral models validation.
# Appendix

## 4.A Paper from Ruju & al. (2020)
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

# Spectral wave modelling of the extreme 2013/2014 winter storms in the North-East Atlantic

Andrea Ruju<sup>a,d,*</sup>, Jean-François Filipot<sup>a</sup>, Abderrahim Bentamy<sup>b</sup>, Fabien Leckler<sup>c,a</sup>

<figure><img src="image_37.png" /></figure>

## ARTICLE INFO

**Keywords:**

Spectral wave modelling

Wind forcing

Wave energy dissipation

Wave breaking

Extreme storms

Storm tracking

## ABSTRACT

This work aims to investigate the impact of wind forcing datasets and wave breaking parameterizations on spectral wave model performance under extremely energetic conditions. For this purpose we used the wave model WaveWatch III to simulate the evolution of the highly energetic storms that occurred in winter 2013/2014 in the North-East Atlantic. We forced the wave model with two different wind datasets: one proceeding from the ECMWF ERA5 reanalysis dataset and the other from satellite observations. Moreover, two wave energy dissipation parameterizations were tested: Test471 and Test500. The model accuracy was assessed by comparing the output datasets with buoy data both in deep and coastal water. Moreover, wave height measurements from satellite were used to assess the model accuracy along storm tracks across the ocean. The accuracy of simulated results shows a significant dependence on the wind forcing and wave dissipation parameterization used. Error metrics computed under storm conditions at wave buoys are consistent with those computed along storm tracks. At the wave buoy locations, all datasets tend to underestimate wave parameters at the peaks of the storms.

## 1. Introduction

Recent work has reported that extreme sea state conditions have increased in terms of frequency and intensity in the last decades (Young and Ribal, 2019; Reguero et al., 2018). This trend, related to climate change and possibly involved in a long-term tendency, has significance for engineering applications: among them we can mention coastal hazard assessment, offshore ship operations and the design of marine structures. Marine engineers and scientists often combine datasets proceeding from different sources in an effort to achieve an accurate and exhaustive description of extreme events and their impacts (O'Reilly et al., 2016; Castelle et al., 2015; Masselink et al., 2016). In this context, by integrating in situ and remote measurements, third-generation spectral wave models and their output make a fundamental contribution towards a better understanding and prediction of extreme wave events.

Third-generation spectral wave models are widely used nowadays for wave hindcast and forecast at global and regional scales (Bernier et al., 2016; Besio et al., 2016; Perez et al., 2017; Sandhya et al., 2018; Ruju et al., 2019). These models solve the wave action balance equation with a set of source terms encompassing the effects of physical processes from wave generation to dissipation (Tolman et al., 2013). Although the recent implementation of physical-based parameterizations

has led to an increase of model output accuracy, simulating extreme wave events remains a challenge (van Vledder et al., 2016; Holthuijsen et al., 2012; Zieger et al., 2015; Campos et al., 2019). This is mainly due the paucity of observations available during the evolution and at the peak of extreme events with respect to moderate and more frequent conditions. As a result of the data used during the parameterization development and model calibration processes, model uncertainties are generally higher for rare wave conditions. For instance, Filipot and Ardhuin (2012) reported a deterioration of error statistics associated with different parameterizations for significant wave heights above 8 m.

Under energetic and storm conditions characterized by large wave steepness values, the wave energy dissipation parameterization takes a key role in spectral evolution and wave growth limitation. Despite the significant attention received, it is likely to represent the least understood source term (Ardhuin et al., 2010). In addition to parameterizations, is well acknowledged that wave model accuracy strongly depends on the accuracy of the wind forcing dataset (Stopa et al., 2016). This work aims to investigate the impact of two different wave breaking parameterizations and two wind forcing datasets under extremely energetic wave conditions. We use the third generation wave

<footnote><sup>a</sup> France Energies Marines, Batiment Cap Ocean, Technopole Brest Iroise, 525 Avenue Alexis de Rochon, 29280 Plouzane, France</footnote>

<footnote><sup>b</sup> Laboratoire d'Oceanographie Spatiale, Institut Francais pour la Recherche et l'Exploitation de la Mer, 29280 Plouzane, France</footnote>

<footnote><sup>c</sup> Service Hydrographique et Oceanographique de la Marine, 29200 Brest, France</footnote>

<footnote><sup>d</sup> Department of Chemical and Geological Sciences, University of Cagliari, 09042 Monserrato, Italy</footnote>

<footnote>* Corresponding author at: Department of Chemical and Geological Sciences, University of Cagliari, 09042 Monserrato, Italy.
E-mail address: rujua@unica.it (A. Ruju).</footnote>
<figure><img src="image_38.png" /><figcaption>**Fig. 1.** Global (blue dots) and local (green dots) grid configuration over the Eastern Atlantic region. Red dots indicate the buoy locations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

model WaveWatch III (WWIII), version 5.16, to simulate the sequence of severe storms occurred that in the North-East Atlantic during the winter 2013/2014.

Previous work has recognized the winter of 2013/2014 as one of the most exceptional in terms of storm sequence and intensity in the North-East Atlantic Ocean (Wadey et al., 2014; Masselink et al., 2016). Due to the relatively south paths of these extra-tropical cyclones, extreme energetic wave conditions were recorded by coastal monitoring systems of Western European countries, from Portugal to Ireland. On coastal areas, these storms drove extreme surge, runup and overtopping causing large morphological changes and strong damage to infrastructures (Castelle et al., 2015; Scott et al., 2016; Autret et al., 2016).

The two wave energy dissipation parameterizations tested in this work are Test471 and 500. They are both included in the parameterization group ST4 available in WWIII version 5.16. Moreover, we assess the impact of two wind forcing datasets. One of them is constituted by the wind analysis obtained through the use of various remotely sensed wind observations (Bentamy et al., 2019; Desbiolles et al., 2017). The other is the ERA5 reanalysis dataset (Hersbach et al., 2019). Model accuracy is assessed by comparing simulated results with the measurements from buoys located in the North-East Atlantic as well as satellite observations along storm tracks.

## 2. Methods

### 2.1. Data collection and storm identification at wave buoys

We collected in situ wave parameters from eight North-East Atlantic wave buoys belonging to different observational networks. Two of them (62163 and 62001) are offshore buoys, located in water depths exceeding 2500 m. The other six (62069, 62103, 62064, 4403, 5602, DW5) are coastal buoys deployed in mean water depths ranging from 30 to 68 m (see Fig. 1 showing the geographical setting). Table 1 lists the wave buoys with mean water depth and main wave parameters.

These buoys are exposed to a combination of long-period Atlantic swells and locally-generated wind waves. Due to the shelter offered

<table><caption>**Table 1**
Wave buoys with mean water depth and significant wave height $H_s$ statistics, for the period considered.</caption><thead><tr><th>Buoy</th><th>depth [m]</th><th>mean $H_s$ [m]</th><th>$H_{s,99}$ [m]</th><th>$H_{s,70}$ [m]</th></tr></thead><tbody><tr><td>62163</td><td>2526</td><td>5.0</td><td>11.9</td><td>5.9</td></tr><tr><td>62001</td><td>4554</td><td>4.5</td><td>10.9</td><td>5.4</td></tr><tr><td>62069</td><td>66</td><td>3.9</td><td>9.6</td><td>4.7</td></tr><tr><td>62103</td><td>68</td><td>2.6</td><td>7.1</td><td>3.1</td></tr><tr><td>62064</td><td>54</td><td>3.4</td><td>8.2</td><td>4.1</td></tr><tr><td>4403</td><td>30</td><td>2.4</td><td>5.8</td><td>3.0</td></tr><tr><td>5602</td><td>45</td><td>3.6</td><td>8.7</td><td>4.3</td></tr><tr><td>DW5</td><td>42</td><td>2.5</td><td>6.6</td><td>3.0</td></tr></tbody></table>

by the surrounding coastline and small islands and the dissipation on the continental shelf, wave height at buoys 62103, 4403 and DW5 is significantly smaller than that at other locations. Buoy 62103 lies in the British channel and it is thus partially sheltered by the Brittany and Cornwall peninsulas. On the other hand, the presence of the islands of Ushant and Belle-Ile dampens the incoming wave energy hitting buoys DW5 and 4403, respectively. Wave propagation at coastal buoy locations is affected not only by topographic features but also by tidal dynamics (currents and water levels) that can be particularly intense in proximity of Brittany shores. All the buoys chosen in this work provide a high time coverage of nearly 100% for the winter 2013/2014 on which this work focuses on.

We used the peak-over-threshold (POT) (Mathiesen et al., 1994) method to identify the 24-h independent storms occurred during the 2013/2014 winter at buoy 62163. The $H_s$ threshold was chosen equal to the 30% exceedance $H_s$ ($H_{s,70}$) calculated over the 2013/2014 winter period. We retained only the storms with a duration larger than 12 h that met the independence criterion with more than 24 h between the end of a storm and the beginning of the following one. Although the threshold value of 30% may seem low for extreme event analysis, due to the highly-energetic period considered, this method allowed the identification of 13 storms in the winter period comprised between the 21st of December and the 21st of March (dates usually taken as of the meteorological start and end of winter). However, we extended the winter period up to the 31 of March to include the 14th storm occurred on the 24th of March; see upper panel of Fig. 2. The extreme wave parameters representative of each storm of the sample were selected as the values occurring at the time in which the maximum wave height was observed during the storm duration.

The adoption of the same method, used for the event identification at buoy 62163, would have led to a different number of storms at each buoy location. For consistency, we recognized at buoys locations the same storms first identified at the offshore buoy 62163. Since this buoy lies at the westernmost location and North-East Atlantic storms are mainly moving eastward (Dodet et al., 2010), they are likely to hit first buoy 62163 and then continue propagating until they reach the other buoys. For this reason, at the other buoy location we expect that both the beginning and the end of a storm happen later than at buoy 62163. Therefore, we identified the beginning of a storm at each buoy location as the time at which $H_s$ firstly increases over $H_{s,70}$ after the beginning of the storm at 62163. Analogously, the end of the storm was set at the time at which $H_s$ falls below $H_{s,70}$ after the end of the storm at 62163. Note that $H_{s,70}$ is different at each location.

The criteria of storm independence and minimal duration prescribed at buoy 62163 are not always met at the other locations. This is particularly evident at buoy 4403 where, just before the February 10th, storms S7 and S8 are contiguous since $H_s$ remains above the threshold for a considerable amount of time from the start of storm S7 to the end of the storm S8. Nevertheless, this procedure has the main benefit of allowing the identification of the same storms (14 in number) at each buoy location, each of them being related to the same synoptic system (see Fig. 2 that highlights the storms over the time series of $H_s$ measured by buoys).
<figure><img src="image_39.png" /><figcaption>Fig. 2. Time series of significant wave height $H_s$ at the buoy locations. Yellow rectangles extend over the storm duration. Grey dashed lines indicate the $H_s$ thresholds $H_{s,70}$ used for storm identification. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 2.2. Data collection and storm tracking from atmospheric pressure

We used the fifth generation ECMWF atmospheric reanalysis ERA5 (Hersbach et al., 2019) as a database to track the low-pressure systems propagation across the Atlantic Ocean during the 2013/2014 winter. We identified the low-pressure systems from the atmospheric pressure at the sea level. First, at each ERA5 output time instant the active low-pressure systems (that we can classify as extra-tropical cyclones) were identified as those systems that have a pressure value lower than 980 hPa and imposing a minimum distance of 5° between different systems. Moreover, the evolution of system propagation in time and space was made assuming a maximum velocity of 120 km/h (33.3 m/s) of the low-pressure system.

This method led to the identification of a large number of low-pressure systems whose life duration spanned from few hours up to several days for the most persistent. To focus on the same events recognized at the buoy locations, we looked for the active systems at the time arrival of the 14 energetic storms, in terms of $H_s$, recorded at the buoy 62163. We named these low-pressure systems with the same name of the storms they drove at the buoy locations (S1, S2, etc.). Since an energetic low-pressure system drove the $H_s$ peak occurred in the last stage of storm S4 at buoy 62163, we added this one (calling it S4B) to the sample constituted by the 14 systems that were active at the beginning of the 14 storms. Once identified, the propagation path of these 15 systems, responsible for the largest wave heights in the winter 2013/2014, was tracked back from its generation in the Western Atlantic and forward to its dissolution in the Eastern Atlantic. Fig. 3

shows the paths of the low-pressure systems highlighting the intense extra-tropical cyclones driving the five highest $H_s$ at the storm peak at buoy 62163.

## 2.3. Modelling techniques

We used the numerical model WAVEWATCH III (WWIII) (Tolman, 2016) version 5.16 to simulate the energetic wave dynamics that occurred in the winter 2013/2014. WWIII is a spectral wave model able to reproduce the physical processes governing wave motion over a wide range of water depths. Its physical and numerical configurations make it suitable to perform hindcast and forecast at global and regional scales. The governing equation of the model is the wave action balance equation in which the source and sink of wave energy is taken into account by means of a set of source terms:

$$ \frac{DN}{Dt} = \frac{S}{\sigma}, \qquad (1) $$

where $D/Dt$ represents the total derivative (moving with a wave component) and $S$ represents the net effect of sources and sinks for the wave action spectrum $N = E/\sigma$ (where $E$ is the energy spectrum and $\sigma$ is the intrinsic frequency of the wave). Parameterizations are usually divided into four main source terms: atmospheric $S_{atm}$, nonlinear $S_{nl}$, ocean $S_{oc}$ and bottom $S_{bt}$. This work focuses on the parameterization of wave energy dissipation by breaking included in the ocean term $S_{oc}$, that is assumed to be the most important sink of wave energy in storm seas (van Vledder et al., 2016).
<figure><img src="image_40.png" /><figcaption>Fig. 3. Paths of the low-pressure systems recorded in the winter 2013/2014 in the North Atlantic. The paths of the 5 most intense systems are highlighted by coloured lines. Markers are 6-h spaced. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

In this study, a multigrid approach allowed the optimization of computational cost given the wide range of physical process scales we are focusing on: from the long scales of wave and swell generation in the deep ocean to the small scales of wave-current interaction and depth-induced processes in coastal water. The model ran over a rectangular grid with a constant spatial resolution of 0.5° covering the entire globe. Wave spectra computed over this grid represented the boundary conditions for the coastal simulations performed over an unstructured grid extending over coastal water from Northern Spain in the South to the British channel in the North. This unstructured mesh was developed in the scope of the HOMONIM project, funded by the French government, in order to improve the operational wave surge forecasting system along French Atlantic coast (Michaud et al., 2015). It is made up of 92 757 nodes with a decreasing resolution from 10 km at offshore boundaries to about 200 m at the coastline and is supported by an accurate and recent 100 m resolution bathymetry also developed in the HOMONIM project (see Biscara et al., 2014). The triangle-based grid is used in WAVEWATCH III with the explicit N-scheme based on contour residual distribution (see Roland, 2009 for a review). Initially implemented in the Wind Wave Model-II (WWM-II), this numerical scheme have then been successfully validated in WAVEWATCH III on an unstructured mesh closely similar to ours (Ardhuin et al., 2009; Boudiere et al., 2013). See Fig. 1 that shows the model grid configuration together with the location of the buoys. It is worth mentioning that, the mean water depths of 4 coastal buoys given in Table 1 are between the minimum and maximum water depths extracted from the model simulations, as expected. The exceptions are buoys 5602 and 4403 whose mean water depth provided by the responsible entity (CEREMA in this case) is smaller than the depth ranges from the model: water depth ranges 56–61 m at 5602 and 35–41 m at 4403.

We discretized the WWIII wave spectra into 32 frequencies and 24 directions. The frequency range extended from 0.0373 to 0.716 Hz, with a frequency increment factor of 1.1. Wave directions were linearly spaced resulting in an angular resolution of 15°. Resonant nonlinear wave-wave interactions occurring between four wave components (quadruplets) were computed with the Discrete Interaction Approximation (DIA) method. Triad wave interactions, accounting for nonlinear energy transfer in the nearshore, were also included through the LTA model. Wave dissipation was simulated with the parameterizations of Ardhuin et al. (2010) and Filipot and Ardhuin (2012) in Test471 and

<table><caption>**Table 2**
Source term treatment in WWIII.</caption><thead><tr><th></th><th>S<sub>in</sub> + S<sub>ds</sub></th><th>S<sub>nl</sub></th><th>S<sub>ir</sub></th><th>S<sub>bot</sub></th><th>S<sub>db</sub></th><th>S<sub>bs</sub></th></tr></thead><tbody><tr><td>Parameterization</td><td>ST4</td><td>NL1</td><td>TR1</td><td>BT4</td><td>DB1</td><td>BS1</td></tr></tbody></table>

Test500, respectively (see Section 2.4). Table 2 lists the parameterizations used for the main source terms, see the WWIII manual (Tolman, 2016) for an exhaustive description of these terms. Note that both Test471 and Test500 are included in parameterization group ST4.

## 2.4. Parameterization of the dissipation induced by wave breaking

The two wave dissipation parameterizations assessed in this work are those by Ardhuin et al. (2010) and Filipot and Ardhuin (2012). Consistent with previous literature (Filipot and Ardhuin, 2012; Leckler et al., 2013), they are referred to as Test471 and Test500, respectively. These formulations recognize that wave energy can be dissipated by the breaking process in two ways: a spontaneous breaking in which the energy of a wave packet is dissipated by the breaking of that very wave packet and a cumulative breaking dissipation in which energy dissipation is the result of the breaking of longer waves wiping out shorter waves. Test471 and Test500 differ in the way the spontaneous breaking source term $S_{bk}$ is computed. For this reason we briefly outline here the computation process for $S_{bk}$.

### 2.4.1. Test471

Following the work of Phillips (1984), several wave parameterizations related breaking probability to spectral saturation. Ardhuin et al. (2010) (therein after ARD10) introduced a saturation-based semiempirical wave breaking parameterization with a larger dissipation rate in the mean wave direction, consistent with the observations of Mironov and Dulov (2007). The directional saturation spectrum $B'(k, \theta)$ is defined as:

$$B'(k, \theta) = \int_{\theta-\Delta}^{\theta+\Delta} k^3 \cos^2(\theta - \theta') E(k, \theta') d\theta', \quad (2)$$

with $\Delta = 80°$, $E(k, \theta')$ is the frequency spectrum and $k$ is the wave number. ARD10 extrapolated the theory of Banner et al. (2000), originally
formulated for dominant waves, over the entire directional spectrum to obtain the breaking probability parameterization $Q_b$:

$$Q_b(k, \theta) = 28.16 \cdot \max[\sqrt{B'(k, \theta)} - \sqrt{B'_r, 0}]^2, \quad (3)$$

where $B'_r$ is the breaking threshold with a correction providing a constant ratio of the root-mean-square orbital velocity and phase speed at different water depths d (Filipot et al., 2010):

$$B'_r = B_r Y(M_4 Y^3 + M_3 Y^2 + M_2 Y + M_1), \quad (4)$$

where $Y = \tanh(kd)$. The deep water threshold $B_r$ and the other constants in the polynomial fit can be found in ARD10. The factor 28.16 comes from the original factor of 22 of Banner et al. (2000), modified by taking into account that wave steepness is on the order of $1.6\sqrt{B'}$ and that the wave counting analysis for a given wave scale from Banner et al. (2000) tends to give a number of waves twice less than that expected for monochromatic waves. The dissipation term of spontaneous breaking $S_{bk}$ is:

$$S_{bk}(k, \theta) = \sigma \frac{C_{ds}}{B_r^2} \{\delta_d \max[B(k) - B'_r, 0]^2 + (1 + \delta_d) \max[B'(k, \theta) - B'_r, 0]^2\} E(k, \theta), \quad (5)$$

in which $C_{ds}$ is a dissipation constant, $\delta_d$ is a coefficient that controls the directionality of breaking and $B(k)$ is the maximum value of $B'(k, \theta)$ for $\theta$ in the range $[0, 2\pi]$. Although this formulation is able to address both deep water and depth-induced breaking, ARD10 warned about the uncertainties involved in its application in shallow water environments.

### 2.4.2. Test500

With the main aim of overcoming the limitations of previous wave breaking parameterizations, Filipot et al. (2010) and Filipot and Ardhuin (2012) (therein after FAB12) made a significant effort towards a unified breaking parameterization valid from the deep ocean to the surf zone. Filipot et al. (2010) divided the frequency spectrum into wave scales with finite bandwidth centred at frequency $f_i$. Then, following Thornton and Guza (1983), they assumed that the breaking wave height distribution for each scale is given by the product of a Rayleigh distribution $P_R(H, f_i)$ and a weight function $W_b(H, f_i)$. In order to extend the formulation outside shallow water, they replaced the breaking criterion of Thornton and Guza (1983), based on the relative water depth $H/d$, with the breaking parameter defined by Miche (1944):

$$\beta_r = \frac{\bar{k}_r H_r}{\tanh(\bar{k}_r d)}, \quad (6)$$

where $\bar{k}_r$ and $H_r$ are the representative wave number and wave height for each wave scale $f_i$. The breaking wave height function $W_b$ is:

$$W_b(H, f_i) = 1.5 \left[ \frac{\beta_r}{\beta_{t,lin}} \right]^2 \left\{ 1 - \exp\left[ - \left( \frac{\beta}{\beta_{t,lin}} \right)^4 \right] \right\}, \quad (7)$$

where $\beta_{t,lin}$ is the breaking threshold defined by Miche (1944) but that takes into account the wave linearization (Filipot et al., 2010), inherent to the wave scale decomposition. The breaking probability for the wave scale with central frequency $f_i$ is:

$$Q_b(f_i) = \int_0^\infty P_R(H, f_i) \cdot W_b(H, f_i) dH. \quad (8)$$

The dissipation source term $S_{bk,i}$ for the component involved in the wave scale i is then given by:

$$S_{bk,i}(f) = \frac{D(f_i)E(f)}{\int_0^\infty E(f)df}, \quad (9)$$

where $D(f_i)$ is the dissipation rate per unit area

$$D(f_i) = Q_b(f_i) \Pi(f_i) \epsilon(f_i), \quad (10)$$

being $\Pi(f_i)$ and $\epsilon(f_i)$ the crest length density per unit area and the dissipation rate per unit length of breaking crest, respectively (Filipot and Ardhuin, 2012).

Since the frequency windows overlap, each spectral component is associated with several wave scales. The dissipation source term is expressed as:

$$S_{bk}(f) = \frac{1}{N} \sum_{i=1}^{N} S_{bk,i}(f), \quad (11)$$

in which N is the number of wave scales involving the frequency f.

## 2.5. Atmospheric forcing

In this work we tested the impact of two atmospheric forcing datasets on WWIII model performance, under a winter characterized by a sequence of exceptionally energetic storm conditions in the North-East Atlantic (Masselink et al., 2016). The first forcing dataset proceeds from the wind analysis obtained through the use of various remotely sensed wind observations. On the other hand, ERA5 reanalysis represents the second forcing dataset. These forcing datasets are briefly described in the following sections.

Besides wind forcing, water levels and flow velocities computed by the Model for Applications at Regional Scales (MARS) (Lazure and Dumas, 2008) were included in the WWIII simulations. MARS simulations were carried out over three nested grids with spatial resolution ranging from 2 km to 250 m in the shallower areas. MARS output was included only in the simulations over the unstructured grid, allowing the computation of wave-current interactions in the coastal environment.

### 2.5.1. Satellite winds

The remotely sensed data, also referred to as satellite wind analyses, used in this study are mostly derived from scatterometer wind retrievals in combination with radiometer observations (Bentamy et al., 2019; Desbiolles et al., 2017). The main sources of remotely sensed wind data are from scatterometers onboard Metop-A (2007-present) and Metop-B (2012-present), and named ASCAT-A and ASCAT-B. Ancillary remotely sensed data are derived from radiometers Special Sensor Microwave Imager Sounder (SSMI/S) onboard the Defence Meteorological Satellite Program (DMSP) F16 (2003-present) and F17 (2006-present), and from WindSat onboard Coriolis satellite (2003-present).

The scatterometer retrieval in combination with radiometer wind observations, and with the European Center of Medium Weather Forecasts (ECMWF) re-analysis model ERA Interim (Simmons et al., 2007), are used for determining regular in space and time surface wind analyses (Desbiolles et al., 2017). These are available at synoptic times (00h:00, 06h:00, 12h:00, and 18h:00 UTC), over the global oceans with a spatial resolution of 0.25°0.25°. Their accuracy, determined through comprehensive comparisons with 6-hourly averaged buoy winds, is of same order of scatterometer retrieval accuracy.

Regarding the study topic, it is of interest to determine some statistics aiming at the characterization of the remotely wind speed and direction analyses at regional scale. To achieve such purpose, satellite wind analyses are compared to collocated (in space and time) 6-hourly averaged wind speed and direction measured by buoys 62103, 62163 and 62001 (see their location in Fig. 1). Scatter plots (not shown) indicate that satellite wind analyses agree well with buoy estimates for all wind speed and direction ranges, including high wind conditions. The correlation between buoy and satellite is almost 1, while symmetrical regression slope and intercept parameter are 1 and of 0.1 m/s, respectively. Furthermore, the low Root Mean Square Error RMSE values (lower than 1 m/s and 20° for wind speed and wind direction, respectively) attest the quality of satellite wind analyses. Satellite data improve the comparisons with in situ wind measurements with respect to the ERA Interim model. In fact, the satellite wind dataset contributes
<figure><img src="image_41.png" /><figcaption>Fig. 4. Observed (grey line) and predicted time series of significant wave height Hs at the deep water buoy 62163 (a) and at the coastal buoy 5602 (b). The four computed datasets are shown. Storms are coloured in yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

to reduce the bias and RMSE values, improving the comparisons for high wind conditions and confirming the results of Bentamy et al. (2017) and Desbiolles et al. (2017).

### 2.5.2. ECMWF winds

In 2018 ECMWF released the ERA5 reanalysis with spatial resolution of 0.25° and 1-h intervals (Hersbach et al., 2019). This dataset combines worldwide observations with model data collected from the 1979 until present. Atmospheric variables are given at the surface and on model levels. The variables used as forcing for the WWIII simulations in this work are the horizontal components of the wind speed at 10 m above the sea level (Tolman, 2016).

Although the ERA5 dataset was originally released with an hourly output resolution, we reduced the time resolution to 6-h intervals in the WWIII forcing. This has been done with the purpose of having the two forcing datasets assessed in this study with same spatial and time resolution, thus ensuring an insightful result comparison. The influence of time resolution of forcing winds on model results will be further addressed in Section 4.

## 3. Results

The impact of forcing winds and energy dissipation parameterizations on model performance is assessed by comparing the simulated datasets with the datasets produced by buoy and altimeter observations. Here, we adopt normalized statistics with the main aim of comparing a large range of wave conditions. The normalized root-mean-square-error NRMSE and normalized bias NBIAS are defined as follows:

$$ \text{NRMSE} = \sqrt{\frac{\sum (O_i - M_i)^2}{\sum O_i^2}}. \qquad (12) $$

$$ \text{NBIAS} = \frac{\sum (O_i - M_i)}{\sum O_i}. \qquad (13) $$

where Oi and Mi are the observed and modelled variables.

For the sake of clarity we assign the names of each model output dataset by specifying the dissipation parameterization used followed by the wind forcing. For instance, the output including Test500 and the ECMWF forcing is called Test500ECMWF.

### 3.1. Wave buoys

Fig. 4 shows the comparison between observed and modelled significant wave height Hs time series at the offshore buoy 62163 and at the coastal buoy 5602. Despite some differences among the four model outputs, overall the WWIII datasets are able to capture the main evolution of the observed dataset. However, evident discrepancies can be found at storm peaks where the modelled Hs underestimate the observations, especially in coastal water. An exception to this trend is represented by Test500Satellite at buoy 62163 that seems to better represent the Hs evolution at the peak of the main storms.

To achieve a quantitative assessment of model performance, modelled time series are linearly interpolated over the observed time series. The scatterplots of the total number of samples N, divided into deep water and coastal water, are shown in Fig. 5, for the wave parameter Hs. NRMSE is on the order 0.1 at deep water buoy locations, ranging between 0.098 for Test500ECMWF and 0.134 for Test500Satellite. Test500Satellite slightly overestimates Hs (NBIAS is 0.031), whereas a small underestimation is given by Test471ECMWF (NBIAS is -0.025). The other two datasets are practically unbiased (|NBIAS|<0.003).

In general terms, for all datasets the NRMSE of Hs increases by few points percentage at coastal buoy locations. This is an expected result given the additional modelling challenges represented by the coastal environment (van Vledder et al., 2016), such as complex bathymetries and tidal currents, with respect to the deep ocean. A remarkable result comes from the observation of coastal water scatterplots and the associated NBIAS of Hs. In fact, whereas all datasets overestimate Hs with positive NBIAS between 0.01 and 0.042, the points associated to higher Hs (Hs>10 m) fall below the line of perfect agreement, meaning that those large Hs are underestimated. This is in agreement with the Hs underestimation at storm peaks already observed in Fig. 4.

The analysis of model accuracy is integrated by the Taylor diagrams in Fig. 6. All datasets have Correlation Coefficients larger than 0.95. In deep water, Test471Satellite and Test500ECMWF have the larger agreement with the observations in terms of Standard Deviation and RMSD. In coastal water, the differences between datasets are less marked, with Test500Satellite slightly improving the prediction in terms of Standard Deviation.
<figure><img src="image_42.png" /><figcaption>Fig. 5. Scatter plot of observed versus modelled significant wave height $H_s$ in deep water (a–d) and coastal water (e–h). The four datasets are shown. The grey line represents the mean of the model values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Fig. 6. Taylor diagrams of Significant wave height in deep (a) and coastal (b) water. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Fig. 7. Scatter plot of observed versus modelled extreme significant wave height $H_s$ in deep water (a–d) and coastal water (e–h). The four datasets are shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 3**
$H_s$ statistics.</caption><thead><tr><th rowspan="3">Dataset</th><th colspan="6">H<sub>s</sub> NRMSE</th></tr><tr><th colspan="3">Deep water</th><th colspan="3">Coastal water</th></tr><tr><th>All</th><th>Extreme</th><th>Storm</th><th>All</th><th>Extreme</th><th>Storm</th></tr></thead><tbody><tr><td>T471Satellite</td><td>0.11</td><td>0.12</td><td>0.11</td><td>0.13</td><td>0.19</td><td>0.13</td></tr><tr><td>T471ECMWF</td><td>0.10</td><td>0.16</td><td>0.11</td><td>0.12</td><td>0.19</td><td>0.13</td></tr><tr><td>T500Satellite</td><td>0.13</td><td>0.09</td><td>0.14</td><td>0.14</td><td>0.16</td><td>0.14</td></tr><tr><td>T500ECMWF</td><td>0.10</td><td>0.11</td><td>0.10</td><td>0.13</td><td>0.16</td><td>0.13</td></tr></tbody><thead><tr><th rowspan="3">Dataset</th><th colspan="6">H<sub>s</sub> NBIAS</th></tr><tr><th colspan="3">Deep water</th><th colspan="3">Coastal water</th></tr><tr><th>All</th><th>Extreme</th><th>Storm</th><th>All</th><th>Extreme</th><th>Storm</th></tr></thead><tbody><tr><td>T471Satellite</td><td>-0.00</td><td>-0.09</td><td>-0.01</td><td>0.01</td><td>-0.12</td><td>-0.00</td></tr><tr><td>T471ECMWF</td><td>-0.03</td><td>-0.15</td><td>-0.05</td><td>0.03</td><td>-0.13</td><td>0.00</td></tr><tr><td>T500Satellite</td><td>0.03</td><td>-0.01</td><td>0.05</td><td>0.03</td><td>-0.08</td><td>0.02</td></tr><tr><td>T500ECMWF</td><td>-0.00</td><td>-0.09</td><td>-0.01</td><td>0.04</td><td>-0.09</td><td>0.03</td></tr></tbody></table>

### 3.1.1. *Storm evolution and peaks*

To explore in more detail the model performance in addressing extreme $H_s$, Fig. 7 provides the scatterplots of observed and modelled $H_s$ at the storm peak. In this case the number of samples $N$ is simply given by the product of the number of storms times the number of locations. For all datasets the NRMSE increases in coastal water. Moreover, NBIAS is always negative confirming the $H_s$ underestimation at storm peaks. Both in deep and coastal water, Test500Satellite gives the lowest

<table><caption>**Table 4**
$T_2$ and $F_e$ statistics in coastal water.</caption><thead><tr><th rowspan="3">Dataset</th><th colspan="6">T<sub>02</sub> NRMSE</th></tr><tr><th colspan="3">Deep water</th><th colspan="3">Coastal water</th></tr><tr><th>All</th><th>Extreme</th><th>Storm</th><th>All</th><th>Extreme</th><th>Storm</th></tr></thead><tbody><tr><td>T471Satellite</td><td>0.12</td><td>0.11</td><td>0.08</td><td>0.33</td><td>0.51</td><td>0.34</td></tr><tr><td>T471ECMWF</td><td>0.09</td><td>0.14</td><td>0.08</td><td>0.34</td><td>0.55</td><td>0.36</td></tr><tr><td>T500Satellite</td><td>0.14</td><td>0.11</td><td>0.09</td><td>0.32</td><td>0.44</td><td>0.32</td></tr><tr><td>T500ECMWF</td><td>0.10</td><td>0.13</td><td>0.07</td><td>0.31</td><td>0.49</td><td>0.32</td></tr></tbody><thead><tr><th rowspan="3">Dataset</th><th colspan="6">T<sub>02</sub> NBIAS</th></tr><tr><th colspan="3">Deep water</th><th colspan="3">Coastal water</th></tr><tr><th>All</th><th>Extreme</th><th>Storm</th><th>All</th><th>Extreme</th><th>Storm</th></tr></thead><tbody><tr><td>T471Satellite</td><td>0.04</td><td>-0.05</td><td>0.01</td><td>-0.01</td><td>-0.35</td><td>-0.04</td></tr><tr><td>T471ECMWF</td><td>-0.01</td><td>-0.10</td><td>-0.03</td><td>-0.05</td><td>-0.42</td><td>-0.11</td></tr><tr><td>T500Satellite</td><td>0.06</td><td>-0.03</td><td>0.03</td><td>0.05</td><td>-0.27</td><td>0.05</td></tr><tr><td>T500ECMWF</td><td>0.00</td><td>-0.08</td><td>-0.02</td><td>-0.01</td><td>-0.36</td><td>-0.05</td></tr></tbody></table>

NBIAS in absolute value (the underestimation is less pronounced). Whereas the largest $H_s$ underestimation (minimum NBIAS) is provided by Test471ECMWF.

In contrast with Fig. 7, we do not observe a systematic negative NBIAS for $H_s$ in Fig. 8. This Figure shows the scatterplots of observed and modelled $H_s$ collected during storms in deep water and at coastal buoy locations. In coastal water, the positive NBIAS values of datasets involving the parameterization Test500 are likely to be driven by the large number of $H_s$ data below 6 m. However, underestimation is still noticeable for more energetic conditions ($H_s>6$ m). Fig. 10 aims at
<figure><img src="image_45.png" /><figcaption>Fig. 8. Scatter plot of observed versus modelled significant wave height $H_s$ in deep water (a–d) and in coastal water (e–h) during storm duration. The four datasets are shown. The grey line represents the mean of the model values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

pointing out the difference between the identification of the extreme $H_s$ at the peak of a storm (the circle, in this case storm S10) and the identification of the $H_s$ values collected during a storm (the thick line). Note that the time instants at which extreme $H_s$ occurs for the observed and modelled dataset do not necessarily coincide. For instance, with the ECMWF and satellite datasets, the modelled maximum is slightly ahead and delayed, respectively.

The Taylor diagrams of Fig. 9 integrate the information provided by the scatter plots of Fig. 8. Correlation Coefficients are on the order of 0.9 for all datasets. In deep, Test500Satellite gives the largest Standard Deviations. In coastal water, the low values of the modelled Standard Deviations seem to confirm the peak underestimation already suggested by the scatterplots of Fig. 8.

Fig. 11 shows the extreme wave energy flux $F_e$ computed at the peak of the storm and during the storm occurrence in coastal water. The wave energy flux $F_e$ is calculated from linear theory, assuming a Rayleigh distribution of wave heights (Longuet-Higgins, 1952), as the product between the wave energy density $E$ and the groupcelerity $c_g$:

$$F_e = E \cdot c_g, \qquad (14)$$

in which

$$E = \frac{1}{8} \rho g H_{\text{RMS}}^2. \qquad (15)$$

$$c_g = \frac{1}{2} c \left( 1 + \frac{2kh}{\sinh(2kh)} \right), \qquad (16)$$

where $\rho$ is the water density, $g$ is the acceleration of gravity, $H_{\text{RMS}}$ is the root mean square wave height ($H_{\text{RMS}} = H_s/1.4$), $c$ is the wavecelerity, $k$ is the wave number and $h$ is the water depth. Both $c$ and $k$ are computed from linear wave theory using the mean period $T_{02}$. $F_e$ plots tend to be more scattered than those of $H_s$, with Test471ECMWF providing the largest NRMSE for $F_e$ both at the peak and during the storm. Test471ECMWF also gives the largest underestimation of $F_e$. Test500Satellite is the only dataset that overestimates $F_e$ (NBIAS = 0.05) during storms although, analogously the $H_s$ trend commented in Fig. 8, large values ($F_e > 0.5 \frac{\text{Mj}}{\text{ms}}$) are clearly underestimated. This point will be discussed in Section 4. Tables 3 and 4 list the error statistics of the four datasets at the wave buoys.

### 3.1.2. Spectral wave analysis and sea/swell decomposition

The availability of the two-dimensional spectra at the four coastal buoys managed by CEREMA (62069, 62064, 4403, 5602) allows the spectral analysis of storms occurred in the winter 2013/2014 at those locations. An example of spectral evolution during a storm is displayed in Figs. 12 and 13 for storm S10 at the coastal buoy 62069. At the beginning of the storm, the measured spectrum shows a variegated shape with multiple peaks (see for instance the secondary peak at 0.17 Hz and
<figure><img src="image_46.png" /><figcaption>Fig. 9. Taylor diagrams of significant wave height during storms in deep (a) and coastal (b) water. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>Fig. 10. Observed (grey) and modelled (black) significant wave height $H_s$ during storm S10 at buoy 5602. Thick lines highlight the $H_s$ evolution during the storm. The circles indicate the extreme $H_s$ values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

200°) that are less marked in the modelled spectrum (T471ECMWF). At the end of the storm, the computed spectrum reproduces the secondary peak at frequencies lower than 0.1 Hz and direction nearly opposite with respect to mean storm direction, probably due to wave reflection at the shoreline. Fig. 13 highlights as the wide 1-D frequency spectrum observed at the beginning of the storm tends to a more narrow shape

as the storm attenuates towards the end. At the peak of the storm, the energy gap between the modelled spectra and the measured one is particularly evident. At the storm beginning, both the atmospheric $S_{in}$ and the energy dissipation $S_{ds}$ source terms proceeding from datasets using satellite forcing are larger that those of the ECMWF dataset. Moreover, Test500 seems to give a smoother dissipation spectrum.
<figure><img src="image_48.png" /><figcaption>Fig. 11. Scatter plot of observed versus modelled extreme (a–d) and collected during storms (e–h) energy flux $F_e$ (a–d) in coastal water. The four datasets are shown. The grey line represents the mean of the model values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

This is consistent with Leckler et al. (2013) and is likely due to the averaging over wave scales of Eq. (11). On the contrary, as a result of its dissipation rate that is local in frequency, Test471 gives a higher dissipation rate at the peak frequency at the beginning and at the peak of the storm. At the end of the storm, Test500 gives the higher dissipation rates. However, their impact on storm evolution seems to be limited since these dissipation rates at the end of the storm are three orders of magnitudes smaller than those at the storm peak.

The identification of the wind and swell components of the spectrum is carried out by means of the wave age criterion first introduced by Hanson and Phillips (2001). The wind sea component W is defined as:

$$W = E^{-1}E|_{U_p>c}, \qquad (17)$$

where E is the total spectral energy and $E|_{U_p>c}$ is the energy of the region of the spectrum under the direct influence of the wind. $U_p$ is the projection of the wind speed, with direction $\theta$, along the mean wave direction $\theta_w$:

$$U_p = CU_{10} \cos(\theta - \theta_w), \qquad (18)$$

where C has been set equal to 1.7. Fig. 14 shows the scatter plots of the wind and swell components of the significant wave height for

<table><caption>Table 5
Sea and swell waves statistics in coastal water.</caption><thead><tr><th rowspan="3">Dataset</th><th colspan="6">H<sub>p</sub> NRMSE</th></tr><tr><th colspan="3">Wind waves</th><th colspan="3">Swell waves</th></tr><tr><th>All</th><th>Extreme</th><th>Storm</th><th>All</th><th>Extreme</th><th>Storm</th></tr></thead><tbody><tr><td>T471Satellite</td><td>0.13</td><td>0.13</td><td>0.11</td><td>0.15</td><td>0.14</td><td>0.11</td></tr><tr><td>T471ECMWF</td><td>0.12</td><td>0.14</td><td>0.10</td><td>0.15</td><td>0.15</td><td>0.11</td></tr><tr><td>T500Satellite</td><td>0.13</td><td>0.10</td><td>0.12</td><td>0.15</td><td>0.12</td><td>0.12</td></tr><tr><td>T500ECMWF</td><td>0.12</td><td>0.10</td><td>0.12</td><td>0.15</td><td>0.13</td><td>0.13</td></tr></tbody><thead><tr><th rowspan="3">Dataset</th><th colspan="6">H<sub>p</sub> NBIAS</th></tr><tr><th colspan="3">Wind waves</th><th colspan="3">Swell waves</th></tr><tr><th>All</th><th>Extreme</th><th>Storm</th><th>All</th><th>Extreme</th><th>Storm</th></tr></thead><tbody><tr><td>T471Satellite</td><td>-0.05</td><td>-0.07</td><td>0.04</td><td>-0.01</td><td>-0.08</td><td>-0.08</td></tr><tr><td>T471ECMWF</td><td>0.03</td><td>-0.10</td><td>0.01</td><td>-0.02</td><td>-0.10</td><td>-0.07</td></tr><tr><td>T500Satellite</td><td>-0.03</td><td>-0.01</td><td>0.06</td><td>-0.00</td><td>-0.05</td><td>-0.04</td></tr><tr><td>T500ECMWF</td><td>0.05</td><td>-0.06</td><td>0.03</td><td>-0.02</td><td>-0.08</td><td>-0.04</td></tr></tbody></table>

the four datasets. Swell waves are characterized by a larger NRMSE than wind waves, with all the datasets that tend to underestimate extreme values larger than 6 m. The extreme values of wind waves
<figure><img src="image_49.png" /><figcaption>**Fig. 12.** Computed (T471ECMWF) and measured frequency-directional wave spectra at the beginning (a and d), peak (b and e) and end (c and f) of the storm S10 at buoy 62069. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>**Fig. 13.** Frequency wave spectra at the beginning (a), peak (b) and end (c) of the storm S10 at buoy 62069 (T471ECMWF). The values of $H_s$ and $T_{02}$ are those measured by the buoy; the wind speed at 10 m $U_{10}$ is from the ECMWF dataset. Atmospheric source terms (d–f). Dissipation source term (g–i). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

($>10$ m) are larger than those of swell waves. These extreme wind wave conditions are slightly underestimated, although the ECMWF forcing leads to a positive NBIAS due to overestimation of moderate values. Table 5 reports the error statistics of the sea and swell components of

the four datasets at the coastal wave buoys. It is worth noting that all datasets have, respectively, positive and negative NBIAS for wind and swell waves under storm duration. Therefore, the virtually unbiased total $H_s$ values during storms may be the result of a balance between
<figure><img src="image_51.png" /><figcaption>Fig. 14. Scatter plot of observed versus modelled significant wave height $H_s$ of the sea component (a–d) and the swell component (e–h). The four datasets are shown. The grey line represents the mean of the model values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

a small overestimation of the wind component combined with a small underestimation of the swell component.

## 3.2. Storm tracking

Fig. 15 shows the result of the storm tracking process for storm S10. The storm path from its generation to its dissipation is superposed to the atmospheric pressure field at the time of storm arrival at buoy 62163 (see Fig. 15b). The centre of the low pressure system lies few hundreds km west of the buoy. The large pressure gradients in the southern part of the system are capable of driving strong westerly winds, as displayed in Fig. 16a.

The intense wind forcing in the southern part of the low pressure system controls the wave storm propagation across the Atlantic. Fig. 16 highlights the spatial relationship between the low pressure system and the wind and wave height fields as the storm hits the offshore buoy 62163. This Figure displays the results obtained by the Test471ECMWF. The largest winds (in excess of 25 m/s) and significant wave heights (in excess of 10 m) are predicted to occur inside the half circle of 800 km radius depicted in Fig. 16.

To assess the model performance along the storm propagation, we retain the $H_s$ altimeter measurements falling inside a half-circle of

800 km of radius south of the centre of the recognized low-pressure systems. See Figs. 3 and 16. Fig. 17 shows the comparison between the $H_s$ from altimeter measurements and from model computations along the storm propagation paths. The model results have been interpolated from the regular grid over the altimeter path. The observation of the scatterplots of Fig. 17 draws the attention to the combined role played by model forcing and parameterization used. Whereas Test471Satellite and Test500ECMWF are characterized by minimal NBIAS values, Test500Satellite and Test471ECMWF give positive and negative NBIAS values, respectively. Test500ECMWF is the one showing the lowest NRMSE.

Fig. 18 compares the wind from the ECMWF and satellite forcing datasets. The wind speed values are extracted along the low-pressure systems path propagation, thus corresponding to the time and location of $H_s$ values of Fig. 17. The NBIAS is slightly negative meaning that, along the extra-tropical cyclone paths, the winds of ECMWF forcing are smaller than those of the satellite forcing. This underestimation becomes more evident for strong winds above 20 m/s. This is consistent with Fig. 17, in which for a given parameterization, NBIAS of $H_s$ is lower with the adoption of ECMWF forcing. In fact, the smaller ECMWF storm winds are likely to yield lower energy transfer rates from
<figure><img src="image_52.png" /><figcaption>**Fig. 15.** (a) Atmospheric pressure field at the moment of the S10 storm arrival at the deep-water buoy 62163 (grey point). The dashed line indicates the low-pressure system path. (b) Time series of significant wave height $H_s$ recorded at the buoy 62163. The red dashed line marks the time instant of panel (a). The storms at the buoy 62163 are highlighted in yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the atmosphere to the wave motion, eventually reducing the sea state growth along storm tracks.

## 4. Discussion

In this work, wave spectral numerical simulations under storm wave conditions show a substantial dependence on the wind forcing and wave dissipation parameterization used. Roland and Arduin (2014) suggested that the quality of wind data and source term parameterizations are the main factors defining the accuracy of spectral wave results. Here, we address this subject under extreme storm conditions. This section highlights the main outcomes of the present study and discusses its results in the light of previous work.

A first analysis assesses the model performance separately at deep water buoys and coastal buoys. Model performance decreases when computed data are compared with buoy measurements in coastal water. This result is consistent with previous studies (Ravdas et al., 2018). In contrast with deep water waves, coastal waves are controlled by the combined effect of irregular shorelines, uneven bathymetry and mean water level oscillations. In addition, the interaction with strong tidal currents is not negligible. Besides the challenges in modelling the complex physics of coastal wave processes, the quality of bathymetric and ocean circulation data play a significant role. In fact, tidal currents, mean water levels and bathymetry data are inevitably affected by errors that may propagate into the wave model and therefore decrease its accuracy in the nearshore.

The NRMSE and NBIAS values of $H_s$ along storm tracks across the Atlantic are consistent with those at deep water buoys during storm duration. For a given parameterization, the use of the wind forcing from satellite data tends to increase the NBIAS with respect to the use of the ECMWF wind forcing. This is likely to be related to an underestimation of extreme winds by the ECMWF reanalysis dataset (Rascle and Arduin, 2013), see Fig. 18. Analogously, for a given wind forcing, NBIAS rises when Test500 is used. This result is consistent with Filipot and

Arduin (2012). In general terms, Test500ECMWF seems to be the most robust choice for the simulation of storm evolution both along storm tracks (NRMSE = 0.1, NBIAS = 0.02) and at buoy locations (NRMSE = 0.1, NBIAS = -0.01).

The spectral analysis shows that all datasets tend to underestimate the swell component oh $H_s$ at four coastal buoys during storm conditions. Fig. 12 reveals that during an energetic storm a considerable amount of low-frequency energy is propagating with an opposite direction with respect to the main storm direction. This highlights the importance of shallow water processes, such as wave reflection at the shoreline, in determining the total energy budget in coastal water. To have a first assessment of the importance of wave reflection during storms, we isolate the seaward component of spectrum (retaining the seaward wave components that are more than 90° apart from the main wave direction). The comparison of the $H_s$ with the measurements shows that in the model this component is largely underestimated (with NBIAS ranging between -0.19 and -0.25 depending on the dataset). Although the great part of wave energy is propagating shoreward, this strong underestimation of the seaward component of the spectrum is likely to contribute to the underestimation of the swell component under storms in coastal water.

In contrast with its values during storm events, $H_s$ at the storm peak is systematically underestimated as revealed by its negative NBIAS values. This result is valid at both offshore and coastal buoys. Although all datasets share this trend, this underestimation is particularly marked using the combination of ECMWF forcing and Test471 parameterization: NBIAS = -0.15 and -0.12 in deep and coastal water, respectively. The dataset Test500Satellite provides the smallest NBIAS, in absolute values, thus reducing the underestimation. Although Test500Satellite provides the smallest NRMSE and absolute NBIAS, it seems that it has the main drawback of overestimating $H_s$ larger than 15 m (see Figs. 7 and 17).

When comparing the results coming from the different model settings described in this work, it is worth mentioning that the calibration
<figure><img src="image_53.png" /><figcaption>**Fig. 16.** (a) Wind speed field at the moment of the S10 storm arrival at the deep-water buoy 62163. (b) Significant wave height $H_s$ field at same time of panel (a). The dashed line indicates the low-pressure system path. The solid black line shows the half circle with radius of 800 km. Altimeter measurement locations are shown by the grey dots. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

of the two wave breaking parameterizations used here has been carried out with the ECMWF wind forcing (Tolman, 2016). Nevertheless, both parameterizations perform well with satellite forcing showing minimal bias for the entire winter timeseries. Since satellite data are expected to improve the characterization of high wind conditions with respect to the ECMWF products (Bentamy et al., 2017), we argue that a new calibration of parameterizations T471 and T500 with the satellite wind forcing may lead to an improvement of the prediction of extreme sea states. However, the calibration task is beyond the scope of this work, focusing on the winter 2013/2014, as longer simulations are required.

Another aspect to be taken into account is that, as can be seen in Fig. 10, $H_s$ measurements are more noisy than simulated results. These spikes are likely to enhance the $H_s$ peak underestimation by the model that could then be mitigated by applying a despiking filter to the measured timeseries. However, we prefer to avoid this alteration due to its involved subjectivity. It can be worth mentioning that Castelle et al. (2015) simply applied a linear coefficient equal to 1.07 to adjust the $H_s$ values from WWIII simulations to the measurements from buoy 62064 under the same period.

An underestimation is found also for extreme mean wave period $T_{02}$ and energy flux $F_e$ in coastal water. NBIAS for $T_{02}$ is slightly smaller, in absolute value, than that for $H_s$. On the other hand, the stronger underestimation of peak values of $F_e$ (NBIAS between -0.27 and -0.42)

is due to the $F_e$ parameter definition, resulting from the product of E which is function of $H_s^2$ and $c_g$ which is function of wave period.

The observed underestimation of extreme wave parameters highlights the importance of the choice of an accurate hindcast product for extreme wave analysis purposes. In fact, long-term hindcasts affected by errors in extreme sea state conditions can strongly impact the probabilistic moments and the tail of the distributions used for extreme event analysis (Campos et al., 2019). This may have a crucial importance especially in calculation design.

### 4.1. Impact of time resolution of forcing winds

Due to the observed rapid evolution of sea states under extreme weather conditions, it appears plausible that the time resolution of the wind forcing might have an impact on wave model result accuracy. In this paper we have tested two wind datasets with the same time resolution of 6 h. Here, we assess a possible negative impact of relatively low time resolutions of the wind forcing dataset. The accuracy of the wave output proceeding from simulations with the Test500 parameterization and forced by ECMWF ERA5 winds at one hour resolution is discussed.

In our case the high-resolution wind forcing does not lead to a reduction of the NRMSE of $H_s$ for the entire dataset in deep water (NRSME = 0.10). The $H_s$ data stays unbiased (NBIAS < 0.005). The main benefit of using a high-resolution wind forcing seems to be related
<figure><img src="image_54.png" /><figcaption>Fig. 17. Scatterplots of the observed versus modelled significant wave height $H_s$ along the storm tracks. Panels (a), (b), (c) and (d) show the comparison for different forcing and parameterizations. The grey line represents the mean of the model values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

to the ability of catching extreme $H_s$ at the storm peak. In fact, the NRMSE for extreme $H_s$ experiences a small decrease (from 0.111 to 0.106). Moreover, the one hour time resolution forcing leads to a smaller underestimation of extreme $H_s$: NBIAS passes from -0.09 to -0.08. Although detectable, the impact of the time resolution of the forcing winds is undoubtedly limited. The limited magnitude of this improvement seems to be related to the lengthy evolution of Atlantic swells that progressively gain energy along their tracks across the ocean. These swells characterized by a large wave age are likely to dominate the sea state at the considered wave buoys, thus reducing the impact of fast wind oscillations included in a high time resolution wind forcing.

What we have observed here in terms of the impact of the increased time resolution of winds on wave modelling is consistent with Men-taschi et al. (2015) who suggested that a resolution increase of the forcing wind field (in their case it was a spatial resolution increase) may not lead to an improvement of single point statistics. According to previous studies (Cavaleri, 2009; Ardhuin et al., 2007; Bertotti and Cavaleri, 2009) they attributed this result to the so-called double penalty effect: some features and patterns may be missed or reproduced in a wrong place in space and time by the model.

## 4.2. Parameterization ST6

Version 5.16 of WWIII includes the new package ST6 which is designed for the parameterization of wind input, wave breaking and swell dissipation. We comment here the results obtained by activating ST6. For our dataset, setting the FAC parameter equal to 1.09, that means increasing the value of the wind drag by 9%, yields a reduction

of absolute NBIAS with respect to the default value of FAC = 1. This is in agreement with Zieger et al. (2015) who used the same value in combination with CFSR wind reanalysis. In fact, our results shows a clear under- and overestimation of $H_s$ with the other two values proposed by Zieger et al. (2015): FAC = 1 in combination with CFS winds and FAC = 1.23 in combination with NOGAPS winds. A more detailed sensitivity analysis of the FAC or other parameters included in the parameterization ST6 is beyond the scope of this work.

Although the use of ST6 with FAC = 1.09 leads to a small NBIAS for the entire dataset (<0.04), $H_s$ values at storm peaks remain underes-timated. With both Satellite and ECMWF wind forcing dataset, the use ST6 leads to NBIAS values of extreme (at the storm peak) $H_s$ comprised between the values associated with Test471 and Test500. For instance, at deep water buoys with ECMWF forcing, the NBIAS for extreme $H_s$ is equal to -0.13, a value that lies between the those associated with Test471 (-0.15) and Test500 (-0.09), see Table 3. This results suggest that, despite ineluctable differences, the general behaviour of ST6 in predicting $H_s$ under moderate and extreme conditions is analogous to what we have already seen and commented for Test471 and Test500 of the parameterization group ST4.

## 5. Conclusions

The aim of this work was to assess the impact of wave breaking parameterizations and wind forcing datasets on the accuracy of spec- tral wave model results under storm wave conditions. We used the WWIII model to simulate the storm sequence occurred in the winter 2013/2014 on the North-East Atlantic. This work focused on two wave breaking parameterizations included in the parameterization group
<figure><img src="image_55.png" /><figcaption>**Fig. 18.** Scatterplot of satellite wind speed versus ECMWF wind speed along the storm tracks. The grey line represents the mean of the model values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

ST4: Test471 and Test500. Moreover, we tested two forcing datasets with six-hour time resolution winds: one based on satellite observations and another based on the ECMWF ERA5 reanalysis. The analysis was carried out firstly by identifying the individual storms at North-East Atlantic buoy locations and then following the storm tracks across the ocean. The main findings are summarized here:

1. The choice of the combination of the wave breaking parameterization and the wind forcing dataset significantly affects the model results in terms of NBIAS and NRMSE of wave parameters. This is valid for wave parameters computed over the entire time series, during storm evolution as well as at the storm peaks. The change of a given breaking parameterization or wind forcing dataset leads to changes in NBIAS and NRMSE of $H_s$ and $T_{02}$ that are on the order of 5%. Due to its definition involving the product between group parameters function of $H_s^2$ and $T_{02}$, $F_e$ suffers more variability.

2. For a given wave dissipation parameterization and wind forcing, the NBIAS and NRMSE values of $H_s$ computed under storm conditions at wave buoys are consistent with those computed along storm tracks across the Atlantic. Test500 together with the satellite wind forcing gives higher $H_s$ values, thus increasing the NBIAS with respect to, respectively, Test471 and ECMWF wind forcing. By improving the error metrics, Test500ECMWF seems to represent the most robust choice for simulating the storm evolution.

3. Negative NBIAS values of $H_s$ at the storm peaks reveal a significant underestimation of extreme wave conditions that is

particularly marked at the coastal buoy locations. This underestimation, common to all the tested datasets, is reduced by using the Test500Satellite dataset.

4. The spectral analysis shows that at the coastal buoys a considerable amount of energy is propagating seaward during storms, possibly as a result of wave reflection at the shoreline. This seaward component is strongly underestimated by the model (NBIAS of the order of -0.2), thus contributing to the underestimation of the swell component at the coastal buoy locations.

5. The use of the high-resolution wind forcing (one-hour resolution) ERA5 does not significantly improve the error statistics computed over the entire time series at the wave buoys. The main benefit of using a high-resolution forcing resides in the (limited, on the order of 1%) improvement of NRMSE and NBIAS values of extreme wave conditions at storm peaks.

### CRediT authorship contribution statement

**Andrea Ruju:** Conceptualization, Methodology, Data curation, Writing - original draft, Visualization, Investigation, Writing - review & editing. **Jean-Francois Filipot:** Conceptualization, Methodology, Writing - review & editing. **Abderrahim Bentamy:** Resources, Analysis, Writing - review & editing. **Fabien Leckler:** Resources, Analysis, Writing - review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements

This work was financially supported by the ARCWIND project "Adaptation and implementation of floating wind energy conversion technology for the Atlantic region" (EAPA 34/2016), which is co-financed by the European Regional Development Fund through the Interreg Atlantic Area Programme. Moreover, it benefited from government support managed by the Agence Nationale de la Recherche, France under the program Investissement d'Avenir with the reference ANR-10-IEED-0006-14 and ANR-10-IEED-0006-26 related to the projects DiMe, France and CARAVELE, France. We thank two anonymous reviewers for their insightful suggestions that improved this paper.

References

Arduin, F., Bertotti, L., Bidlot, J.-R., Cavaleri, L., Filipetto, V., Lefevre, J.-M., Wittmann, P., 2007. Comparison of wind and wave measurements and models in the Western Mediterranean Sea. Ocean Eng. 34 (3), 526-541. http://dx.doi.org/10.1016/j.oceaneng.2006.02.008, URL: http://www.sciencedirect.com/science/article/pii/S0029801806001193.

Arduin, F., Marié, L., Rascle, N., Forget, P., Roland, A., 2009. Observation and estimation of Lagrangian, Stokes, and eulerian currents induced by wind and waves at the sea surface. J. Phys. Oceanogr. 39 (11), 2820-2838. http://dx.doi.org/10.1175/2009JPO4169.1, arXiv:https://doi.org/10.1175/2009JPO4169.1.

Arduin, F., Rogers, E., Babanin, A.V., Filipot, J.-F., Magne, R., Roland, A., van der Westhuysen, A., Queffelou, P., Lefevre, J.-M., Aouf, L., Collard, F., 2010. Semipirical dissipation source functions for ocean waves. Part I: Definition, calibration, and validation. J. Phys. Oceanogr. 40 (9), 1917-1941. http://dx.doi.org/10.1175/2010JPO4324.1, arXiv:https://doi.org/10.1175/2010JPO4324.1.

Autret, R., Dodet, G., Fichaut, B., Suanez, S., David, L., Leckler, F., Arduin, F., Ammann, J., Grandjean, P., Allemand, P., Filipot, J.-F., 2016. A comprehensive hydro-geomorphic study of cliff-top storm deposits on banng island during winter 2013-2014. Marine Geol. 382, 37-55. http://dx.doi.org/10.1016/j.margeo.2016.09.014, URL: http://www.sciencedirect.com/science/article/pii/S0025322716302201.

Banner, M.L., Babanin, A.V., Young, I.R., 2000. Breaking probability for dominant waves on the sea surface. J. Phys. Oceanogr. 30 (12), 3145-3160. http://dx.doi.org/10.1175/1520-0485(2000)030<3145:BFPDWO>2.0.CO;2, arXiv:https://doi.org/10.1175/1520-0485(2000)030<3145:BFPDWO>2.0.CO;2.

Bentamy, A., Grodsky, S.A., Elyouncha, A., Chapron, B., Desbiolles, F., 2017. Homogenization of scatterometer wind retrievals. Int. J. Climatol. 37 (2), 870-889. http://dx.doi.org/10.1002/joc.4746, URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.4746, arXiv:https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/joc.4746.

Bentamy, A., Mouche, A., Grouazel, A., Moujane, A., Mohamed, A.A., 2019. Using sentinel-1A SAR wind retrievals for enhancing scatterometer and radiometer regional wind analyses. Int. J. Remote Sens. 40 (3), 1120-1147. http://dx.doi.org/10.1080/01431161.2018.1524174, arXiv:https://doi.org/10.1080/01431161.2018.1524174.

Bernier, N.B., Alves, J.-H.G.M., Tolman, H., Chawla, A., Peel, S., Pouliot, B., Bélanger, J.-M., Pellerin, P., Lépine, M., Roch, M., 2016. Operational wave prediction system at environment Canada: Going global to improve regional forecast skill. Weather Forecast. 31 (2), 353-370. http://dx.doi.org/10.1175/WAF-D-15-0087.1, arXiv:https://doi.org/10.1175/WAF-D-15-0087.1.

Bertotti, L., Cavaleri, L., 2009. Wind and wave predictions in the Adriatic Sea. J. Mar. Syst. 78, S227 - S234. http://dx.doi.org/10.1016/j.jmarsys.2009.01.018, URL: http://www.sciencedirect.com/science/article/pii/S0924796309001511.

Besio, G., Mentaschi, L., Mazzino, A., 2016. Wave energy resource assessment in the Mediterranean Sea on the basis of a 35-year hindcast. Energy 94, 50-63. http://dx.doi.org/10.1016/j.energy.2015.01.044, URL: http://www.sciencedirect.com/science/article/pii/S0360544215014127.

Biscara, L., Schmitt, T., Correard, S., Créach, R., 2014. Modèles numériques de bathymétrie pour la prévision hydrodynamique du dispositif vigilance vaguessubmersions. In: Journées Nationales Génie Côtier - Génie Civil. pp. 547-556. http://dx.doi.org/10.5150/jngcgc.2014.060.

Boudiere, E., Maisondieu, C., Arduin, F., Accensi, M., Pineau-Guillou, L., Lepesqueur, J., 2013. A suitable metocean hindcast database for the design of Marine energy converters. Int. J. Marine Energy 3-4, e40 - e52. http://dx.doi.org/10.1016/j.ijomae.2013.11.010, URL: http://www.sciencedirect.com/science/article/pii/S2214166913000362. Special Issue Selected Papers - EWTEC2013.

Campos, R., Soares, C.G., Alves, J., Parente, C., Guimaraes, L., 2019. Regional long-term extreme wave analysis using hindcast data from the South Atlantic Ocean. Ocean Eng. 179, 202-212. http://dx.doi.org/10.1016/j.oceaneng.2019.03.023, URL: http://www.sciencedirect.com/science/article/pii/S0029801818318110.

Castelle, B., Marieu, V., Bujan, S., Splinter, K.D., Robinet, A., Senechal, N., Ferreira, S., 2015. Impact of the winter 2013-2014 series of severe Western Europe storms on a double-barred sandy coast: Beach and dune erosion and megacusp embayments. Geomorphology 238, 135-148. http://dx.doi.org/10.1016/j geomorph.2015.03.006, URL: http://www.sciencedirect.com/science/article/pii/S0169555X15001385.

Cavaleri, L., 2009. Wave modeling—Missing the peaks. J. Phys. Oceanogr. 39 (11), 2757-2778. http://dx.doi.org/10.1175/2009JPO4067.1, arXiv:https://doi.org/10.1175/2009JPO4067.1.

Desbiolles, F., Bentamy, A., Blanke, B., Roy, C., Mestas-Nuñez, A.M., Grodsky, S.A.,HBette, S., Cambon, G., Maes, C., 2017. Two decades [1992-2012] of surface wind analyses based on satellite scatterometer observations. J. Mar. Syst. 168, 38-56. http://dx.doi.org/10.1016/j.jmarsys.2017.01.003, URL: http://www.sciencedirect.com/science/article/pii/S0924796316302068.

Dodet, G., Bertin, X., Taborda, R., 2010. Wave climate variability in the Northeast Atlantic Ocean over the last six decades. Ocean Model. 31 (3), 120-131. http://dx.doi.org/10.1016/j.ocemod.2009.10.010, URL: http://www.sciencedirect.com/science/article/pii/S1463500309002066.

Filipot, J.-F., Arduin, F., 2012. A unified spectral parameterization for wave breaking: From the deep ocean to the surf zone. J. Geophys. Res.: Oceans 117 (C11), http://dx.doi.org/10.1029/2011JC007784, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JC007784, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2011JC007784.

Filipot, J.-F., Arduin, F., Babanin, A.V., 2010. A unified deep-to-shallow water wave-breaking probability parameterization. J. Geophys. Res.: Oceans 115 (C4), http://dx.doi.org/10.1029/2009JC005448, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2009JC005448, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2009JC005448.

Hanson, J.L., Phillips, O.M., 2001. Automated analysis of ocean surface directional wave spectra. J. Atmos. Ocean. Technol. 18 (2), 277-293. http://dx.doi.org/10.1175/1520-0426(2001)018<0277:AAOOSD>2.0.CO;2, arXiv:https://doi.org/10.1175/1520-0426(2001)018<0277:AAOOSD>2.0.CO;2.

Hersbach, H., Bell, W., Berrisford, P., Horányi, A., J., M.-S., Nicolas, J., Radu, R., Schepers, D., Simmons, A., Soci, C., Dee, D., 2019. Global reanalysis: goodbye ERA-Interim, hello ERAS. pp. 17-24. http://dx.doi.org/10.21957/vf291hehd7, URL: https://www.ecmwf.int/node/19027.

Holthuijsen, L.H., Powell, M.D., Pietrzak, J.D., 2012. Wind and waves in extreme hurricanes. J. Geophys. Res.: Oceans 117 (C9), http://dx.doi.org/10.1029/2012JC007983, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2012JC007983, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2012JC007983.

Lazure, P., Dumas, F., 2008. An external-internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS). Adv. Water Resour. 31 (2), 233-250. http://dx.doi.org/10.1016/j.adwatres.2007.06.010, URL: http://www.sciencedirect.com/science/article/pii/S0309170807001121.

Leckler, F., Arduin, F., Filipot, J.-F., Robinov, A., 2013. Dissipation source terms and whitecap statistics. Ocean Model. 70, 62-74. http://dx.doi.org/10.1016/j.ocemod.2013.03.007, URL: http://www.sciencedirect.com/science/article/pii/S1463500313000474. Ocean Surface Waves.

Longuet-Higgins, M.S., 1952. On the statistical distribution of the heights of sea waves. J. Mar. Res. 11 (3), 245-266.

Masselink, G., Scott, T., Poate, T., Russell, P., Davidson, M., Conley, D., 2016. The extreme 2013/2014 winter storms: hydrodynamic forcing and coastal response along the southwest coast of England. Earth Surf. Process. Landf. 41 (3), 378-391. http://dx.doi.org/10.1002/esp.3836, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/esp.3836, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/esp.3836.

Mathiesen, M., Goda, Y., Hawkes, P.J., Mansard, E., Martin, M.J., Peltier, E., Thompson, E.F., Vledder, G.V., 1994. Recommended practice for extreme wave analysis. J. Hydraul. Res. 32 (6), 803-814. http://dx.doi.org/10.1080/00221689409498691, arXiv:https://doi.org/10.1080/00221689409498691.

Mentaschi, L., Besio, G., Cassola, F., Mazzino, A., 2015. Performance evaluation of Wavewatch III in the Mediterranean Sea. Ocean Model. 90, 82-94. http://dx.doi.org/10.1016/j.ocemod.2015.04.003, URL: http://www.sciencedirect.com/science/article/pii/S1463500315000578.

Michaud, H., Dalphinet, A., Huchet, M., Pasquet, A., Baraille, R., Leckler, F., Aouf, L., Roland, A., Sikiric, M., Arduin, F., Filipot, J.-F., 2015. Implementation of the next French operational coastal wave forecasting system and application to a wave-current interaction study. In: 14th International Workshop on Wave Hindcasting and Forecasting and 5th Coastal Hazard Symposium.

Miche, A., 1944. Mouvements ondulatoires de la mer en profondeur croissante ou décroissante. Forme limite de la houle lors de son déferlement. Application aux digues maritimes. Troisième partie. Forme et propriétés des houles limites lors du déferlement. Croissance des vitesces vers la rive. Ann. Ponts Chaussées Tome 114, 369-406.

Mironov, A.S., Dulov, V.A., 2007. Detection of wave breaking using sea surface video records. Meas. Sci. Technol. 19 (1), 015405. http://dx.doi.org/10.1088/0957-0233/19/1/015405.

O'Reilly, W., Olfe, C.B., Thomas, J., Seymour, R., Guza, R., 2016. The California coastal wave monitoring and prediction system. Coast. Eng. 116, 118-132. http://dx.doi.org/10.1016/j.coastaleng.2016.06.005, URL: http://www.sciencedirect.com/science/article/pii/S0378383916301120.
Perez, J., Menendez, M., Losada, I.J., 2017. GOW2: A global wave hindcast for coastal applications. Coast. Eng. 124, 1–11. http://dx.doi.org/10.1016/j.coastaleng.2017.03.005, URL: http://www.sciencedirect.com/science/article/pii/S0378383917300443.

Phillips, O.M., 1984. On the response of short ocean wave components at a fixed wavenumber to ocean current variations. J. Phys. Oceanogr. 14 (9), 1425–1433. http://dx.doi.org/10.1175/1520-0485(1984)014<1425:OTROSO>2.0.CO;2, arXiv: https://doi.org/10.1175/1520-0485(1984)014<1425:OTROSO>2.0.CO;2.

Rascle, N., Ardhuin, F., 2013. A global wave parameter database for geophysical applications. Part 2: Model validation with improved source term parameterization. Ocean Model. 70, 174–188. http://dx.doi.org/10.1016/j.ocemod.2012.12.001, URL: http://www.sciencedirect.com/science/article/pii/S1463500312001709. Ocean Surface Waves.

Ravdas, M., Zacharioudaki, A., Korres, G., 2018. Implementation and validation of a new operational wave forecasting system of the Mediterranean Monitoring and Forecasting Centre in the framework of the Copernicus Marine Environment Monitoring Service. Nat. Hazards Earth Syst. Sci. 18 (10), 2675–2695. http://dx.doi.org/10.5194/nhess-18-2675-2018, URL: https://www.nat-hazards-earth-syst-sci.net/18/2675/2018/.

Reguero, B.G., Losada, I.J., Méndez, F.J., 2018. A recent increase in global wave power as a consequence of oceanic warming. In: Nature Communications.

Roland, A., 2009. Development of WWM II: Spectral Wave Modelling on Unstructured Meshes (Ph.D. thesis). Technische Universität Darmstadt, p. 212.

Roland, A., Ardhuin, F., 2014. On the developments of spectral wave models: numerics and parameterizations for the coastal ocean. Ocean Dyn. 64 (6), 833–846. http://dx.doi.org/10.1007/s10236-014-0711-z.

Ruju, A., Passarella, M., Trogu, D., Buosi, C., Ibba, A., De Muro, S., 2019. An operational wave system within the monitoring program of a mediterranean beach. J. Marine Sci. Eng. 7 (2), http://dx.doi.org/10.3390/jmse7020032, URL: https://www.mdpi.com/2077-1312/7/2/32.

Sandhya, K., Murty, P., Deshmukh, A.N., Nair, T.B., Shenoi, S., 2018. An operational wave forecasting system for the east coast of India. Estuar. Coast. Shelf Sci. 202, 114–124. http://dx.doi.org/10.1016/j.ecss.2017.12.010, URL: http://www.sciencedirect.com/science/article/pii/S0272771417310417.

Scott, T., Masselink, G., O'Hare, T., Saulter, A., Poate, T., Russell, P., Davidson, M., Conley, D., 2016. The extreme 2013/2014 winter storms: Beach recovery along the southwest coast of England. Marine Geol. 382, 224–241. http://dx.doi.org/10.1016/j.margeo.2016.10.011, URL: http://www.sciencedirect.com/science/article/pii/S0025322716302766.

Simmons, A., Uppala, S., Dee, D., Kobayashi, S., 2007. ERA-Interim: New ECMWF Reanalysis Products from 1989 Onwards. ECMWF, pp. 25–35. http://dx.doi.org/10.21957/pocnex23c6, URL: https://www.ecmwf.int/node/17713.

Stopa, J.E., Ardhuin, F., Babanin, A., Zieger, S., 2016. Comparison and validation of physical wave parameterizations in spectral wave models. Ocean Model. 103, 2–17. http://dx.doi.org/10.1016/j.ocemod.2015.09.003, URL: http://www.sciencedirect.com/science/article/pii/S1463500315001614. Waves and coastal, regional and global processes.

Thornton, E.B., Guza, R.T., 1983. Transformation of wave height distribution. J. Geophys. Res.: Oceans 88 (C10), 5925–5938. http://dx.doi.org/10.1029/JC088iC10p05925, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JC088iC10p05925, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JC088iC10p05925.

Tolman, H.L., 2016. User Manual and System Documentation of WAVEWATCH III Version 5.16. Tech. Rep., NOAA/NWS/NCEP: 5830 University Research Court, College Park, MD 20740, USA.

Tolman, H.L., Banner, M.L., Kaihatu, J.M., 2013. The NOPP operational wave model improvement project. Ocean Model. 70, 2–10. http://dx.doi.org/10.1016/j.ocemod.2012.11.011, URL: http://www.sciencedirect.com/science/article/pii/S1463500312001722. Ocean Surface Waves.

van Vledder, G.P., Hulst, S.T.C., McConochie, J.D., 2016. Source term balance in a severe storm in the Southern North Sea. Ocean Dyn. 66 (12), 1681–1697. http://dx.doi.org/10.1007/s10236-016-0998-z.

Wadey, M.P., Haigh, I.D., Brown, J.M., 2014. A century of sea level data and the UK's 2013/14 storm surges: an assessment of extremes and clustering using the Newlyn tide gauge record. Ocean Sci. 10 (6), 1031–1045. http://dx.doi.org/10.5194/os-10-1031-2014, URL: https://www.ocean-sci.net/10/1031/2014/.

Young, I.R., Ribal, A., 2019. Multiplatform evaluation of global trends in wind speed and wave height. Science 364 (6440), 548–552. http://dx.doi.org/10.1126/science.aav9527, URL: https://science.sciencemag.org/content/364/6440/548, arXiv:https://science.sciencemag.org/content/364/6440/548.full.pdf.

Zieger, S., Babanin, A.V., Rogers, W.E., Young, I.R., 2015. Observation-based source terms in the third-generation wave model WAVEWATCH. Ocean Model. 96, 2–25. http://dx.doi.org/10.1016/j.ocemod.2015.07.014, URL: http://www.sciencedirect.com/science/article/pii/S1463500315001237. Waves and coastal, regional and global processes.
4.B Paper from Filipot & Cheung (2012)
<figure><img src="image_56.png" /></figure>

# Spectral wave modeling in fringing reef environments

Jean-François Filipot *, Kwok Fai Cheung

Department of Ocean and Resources Engineering, University of Hawaii at Manoa, Honolulu, HI 96822, USA

## ARTICLE INFO

**Article history:**

Received 26 August 2011

Received in revised form 5 April 2012

Accepted 10 April 2012

Available online 22 May 2012

**Keywords:**

Bottom friction

Fringing reefs

Phase-averaged model

Spectral model

SWAN

Wave breaking

## ABSTRACT

This study examines the implementation and validity of phase-averaged spectral models in describing wave transformation over fringing reefs. The wave breaking mechanism over the abrupt reef crest, dissipation due to bottom friction, and nonlinear energy transfer over the reef flat set the problem apart from general applications in non-tropical settings. The Simulating Wave Nearshore (SWAN) model, which is widely used in coastal engineering applications, has embedded parameterizations of these surf-zone processes for examination. In addition, we extend the model to include a recently implemented bottom-friction source term for tropical reefs. Comparisons of model predictions with laboratory and field observations over fringing reefs identify promising wave-breaking and bottom-friction parameterizations for calibration. The calibrated parameterizations enable SWAN to reproduce the observed wave height and setup and extend the applicability of spectral models to fringing reef environments. SWAN, however, lacks the parameterization to describe the nonlinear energy transfer toward the infra-gravity band that was observed in the laboratory and field studies.

© 2012 Elsevier B.V. All rights reserved.

## 1. Introduction

Spectral wave models are widely used around the world for numerous applications from operational wave predictions at global and regional scales to site-specific coastal engineering applications. These models provide a phase-averaged description of the generation and propagation of ocean waves subject to wind forcing. Let $f$ denote the frequency and $\theta$ the wave direction. The evolution of the wave energy density spectrum $E(f, \theta)$ in time and space is defined by the energy balance equation

$$ \frac{dE(f, \theta)}{dt} = S \qquad (1) $$

where $d/dt$ denotes material derivative and $S$ the source term (Gelci et al., 1957). In third generation spectral wave models, the spectrum can describe multi-modal sea states and the source term contains parameterizations of the physical processes that include energy input by winds, nonlinear energy transfer among frequency components, and dissipation through wave breaking and bottom friction (e.g., Tolman, 2008; Komen et al., 1994; Booj et al., 1999; Zijlema, 2010). These models owe their success to their low computational costs and proven accuracy for typical ocean conditions due in part to the extensive and continuing efforts on parameterization of the physics involved (e.g., Ardhuin et al., 2010). The early parameterizations have been mostly designed for the open ocean, but efforts have been made to extend their capabilities

and implementations to shallow water (Eldeberky, 1996; Ris et al., 1999; Filipot and Ardhuin, in press). Although the accuracy of these parameterizations has been extensively examined for continental beach environments, very little is known on their ability to model wave transformation over fringing reefs.

Commonly found in the tropics and subtropics, fringing reefs are characterized by a steep slope with an abrupt transition at the reef crest to a shallow platform attached to the shoreline. Intense wave breaking typically occurs around the reef crest enhancing the oxygen content and circulation that support the coral ecosystem (Abelson et al., 1993). Wave breaking at fringing reefs protects the coastline from severe wave conditions, but increases the water level due to wave setup (e.g. Vetter et al., 2010). Over the flat, the bottom friction associated with the reef formation and corals may dissipate as much energy as wave breaking (Lugo-Fernández et al., 1998; Lowe et al., 2005). Although the nonlinear triad interactions of the frequency components are energy conservative, they cause dramatic alteration of the spectral shape across the reef. Energy is transferred from the spectral peak to the super-harmonics during the shoaling phase (Eldeberky, 1996), making the wave crest sharper and the trough flatter, and then toward the infra-gravity band in the surf zone (Sheremet et al., 2011). Seiching occurs if the natural period of the reef flat matches the infra-gravity band (Nwogu and Demirbilek, 2010; Roeber et al., 2010a, 2010b; Péquignet et al., 2011). These processes, which have important implications for coastal protection as well as for the fragile coral ecosystem, have drawn a growing interest in the community.

Because of the low computational cost, the capability to account for wind effects, and ready coupling with a circulation model for storm surge computations, researchers and engineers have applied spectral wave models to tropical and subtropical coastal regions without

<footnote>* Corresponding author at: Service Hydrographique et Océanographique de la Marine, Brest, France. Tel.: +33 298 149 936.</footnote>

<footnote>E-mail addresses: filipot@shom.fr (J-F. Filipot), cheung@hawaii.edu (KF. Cheung).</footnote>
knowing the full impact of their limitations. The motivation of this study is to assess the validity of spectral wave models for application in fringing reef environments. We examine the parameterizations for wave breaking and bottom-friction in the commonly-used Simulating Wave Nearshore (SWAN) model (Booij et al., 1999; Ris et al., 1999). In addition, we implement the bottom friction parameterization of Lowe et al. (2005) for a barrier reef in SWAN. A comparison with the laboratory data from Demirbilek et al. (2007) for wave transformation and setup over a smooth reef configuration provides an initial assessment of the wave breaking parameterizations and identification of a promising scheme for calibration. A general assessment of the bottom-friction parameterizations in SWAN is made with 10 days of observational data over a fringing reef on North Shore, Oahu, Hawaii. The dataset allows calibration of the parameterization of Lowe et al. (2005) for the fringing reef at the site. In both case studies, the ability of SWAN in reproducing both the total energy and the spectral shape is examined and the advantages and limitations of spectral wave modeling are discussed.

## 2. Modeling of coastal processes

The governing Eq. (1) for phase-averaged spectral models already accounts for wave shoaling and refraction in the evolution of the energy density spectrum over time and space. Other coastal processes such as depth-induced wave breaking, dissipation due to bottom friction, and triad nonlinear energy transfer among frequency components are included in the source term through parameterization of the physical processes. This section provides a brief review of these parameterizations in SWAN, the formulation of Lowe et al. (2005) for bottom friction in reef environments, and the computation of wave setup.

### 2.1. Depth-induced wave breaking

The physical processes of wave breaking are complex and its spectral signature is still an active field of research (e.g., Vink, 2001; Filipot et al., 2010). The parameterization of the processes includes the wave breaking probability and the subsequent dissipation rate. In shallow water, the breaker index limits the wave height $H$ to a fraction of the water depth $d$ as

$$ \gamma = \frac{H}{d} \qquad (2) $$

Longuet-Higgins (1974) derived $\gamma=0.78$ from solitary wave theory. For irregular waves on a sloping bottom, the breaker index varies over a wide range. For instance, Raubenheimer et al. (1996) observed values ranging from 0.3 to 1.2 based on the significant wave height $H_s$. Several hypotheses were advanced to explain the variability of the breaker index in terms of the offshore wave steepness (Battjes and Stive, 1985), the bottom slope (Raubenheimer et al., 1996; Nelson, 1987), the water depth parameter $kd$, where $k$ is the wave number (Ruessink et al., 2003), and the phase angle of the complex bi-spectrum (van der Westhuysen, 2010). The dissipation rate is typically formulated by approximation of the broken waves as bores with heights inferred from probability density functions based on the Rayleigh distribution (Goda, 1975; Battjes and Janssen, 1978; Thornton and Guza, 1983).

SWAN includes five parameterizations to describe depth-induced wave breaking processes. The formulation by Battjes and Janssen (1978), hereafter BJ78, provides the framework for three of the parameterizations, in which the total dissipation rate is given by

$$ D_{tot} = \frac{1}{4} B Q_B \bar{f} H_m^2 \qquad (3) $$

where $B$ is a constant controlling the level of energy dissipation, $\bar{f}$ is the mean frequency, $H_m$ is the maximum wave height in terms of the breaker index (2), and $Q_B$ is the breaking probability of the random sea state. The breaking probability $Q_B$ can be estimated by means of

$$ \frac{1-Q_B}{-\ln Q_B} = \left( \frac{H_{rms}}{H_m} \right)^2 \qquad (4) $$

where $H_{rms}$ is the root-mean-square wave height. The breaker index is the key parameter in the energy dissipation formulation. Battjes and Janssen (1978) provided $\gamma=0.73$ from a wide range of wave breaking observations on beaches. From laboratory and field data, Nelson (1987, 1994, 1994), hereafter NE87, introduced a dependence of the breaker index on the bottom slope $\alpha$ as

$$ \gamma = 0.55 + 0.88 \exp[-0.012 \cot(\alpha)] \qquad (5) $$

in which $\gamma$ varies from 0.55 for horizontal bottoms to 1.3 for very steep slopes. From analysis of two datasets gathered in barred beach environments, Ruessink et al. (2003), hereafter RU03, reported a dependence of the breaker index on the water depth parameter as

$$ \gamma = 0.76kd + 0.29 \qquad (6) $$

In BJ78, NE87 and RU03, the breaker index defines a single maximum wave height after the onset of breaking and thus also controls the dissipation rate in the surf zone.

Thornton and Guza (1983), hereafter TG83, parameterized the dissipation rate with a distribution of breaking wave heights in contrast to the approaches based on BJ78. The total dissipation rate in the random wave field is given by

$$ D_{tot} = - \frac{B^3 \bar{f}}{4 d} \int_{0}^{\infty} H^3 P(H) W(H) dH \qquad (7) $$

in which $P$ is the Rayleigh distribution and $W$ is a weight function given by

$$ W(H) = \left( \frac{H_{rms}}{\gamma d} \right)^2 \qquad (8) $$

where a breaker index of $\gamma=0.42$ is predefined in the parameterization. Van der Westhuysen (2010), hereafter WE10, assumes a constant distribution of the weight function and parameterizes the biphase, which accounts for wave asymmetry and hence the nonlinearity and breaking susceptibility. The resulting dissipation rate is

$$ D_{tot} = - \frac{3\sqrt{\pi} B^3 \bar{f}}{16 d} \left( \frac{\beta}{\beta_{ref}} \right)^2 H_{rms}^3 \qquad (9) $$

in which the biphase is expressed in terms of the Ursell number $Ur$ as

$$ \beta = \frac{-\pi}{2} + \frac{-\pi}{2} \tanh\left(\frac{0.2}{Ur}\right) \qquad (10) $$

and $\beta_{ref}=-4\pi/9$. Note that, in relating the breaking probability to wave nonlinearity, this approach shares some similarities with the work of Filipot et al. (2010).

The breaking probability and dissipation rate produce reasonable results for continental shelf beaches with gentle slopes. In SWAN, a default value of $B=1$ is typically used. The total dissipation is distributed over the frequency range and direction bins in proportion to the spectral density. The source term for depth-induced breaking dissipation is given by

$$ S_{bk}(f, \theta) = D_{tot} \frac{E(f, \theta)}{E_{tot}} \qquad (11) $$
where $D_{tot}$ is the dissipation rate from Eqs. (3), (7) or (9) depending on the parameterization and $E_{tot}$ is the total spectral energy from integration of $E(f, \theta)$. This approach, which does not alter the spectral shape during the wave breaking process, is consistent with the findings of Eldeberky and Battjes (1996) from flume experiments involving single-peaked spectra.

## 2.2. Bottom friction

When the water depth becomes less than one half of the wavelength, the wave orbital motions extend to the bottom resulting in energy dissipation via various linear or nonlinear processes. Shemdin et al. (1978) identified and examined four dissipation processes: wave orbital velocity friction on the bed, percolation of the water particles through the bottom substrate, wave-induced motion of a soft muddy bottom, and scattering of surface waves by bottom irregularities. Friction dominates in sandy and hard substrate conditions commonly found along the world's coastlines. The shear stress and the dissipation rate are typically parameterized in terms of a drag law or an eddy viscosity model (Mirfenderesk and Young, 2003). In spectral wave models, the dissipation source term for bottom friction can be expressed as

$$S_{bf}(f, \theta) = -C_b \frac{(2\pi f)^2}{g^2 \sinh^2 kd} E(f, \theta) \quad (12)$$

where g is acceleration due to gravity and $C_b$ is a bottom friction coefficient derived from parameterization of the physical processes (Bertotti and Cavaleri, 1994).

SWAN includes three bottom-friction parameterizations based on the empirical JONSWAP model of Hasselmann et al. (1973), the drag law model of Collins (1972), and the eddy viscosity model of Madsen (1988). The empirical JONSWAP model assumes a constant $C_b$ of 0.038 m$^2$ s$^{-3}$ from observations of swell energy decay in the North Sea even though the observed value varies over two orders of magnitude during the field study. The drag law model, which is primarily a simplification of Hasselmann and Collins (1968) for regular waves, can be written as

$$C_b = C_f g U_{rms} \quad (13)$$

where $C_f=0.015$ and $U_{rms}$ is the root-mean-square wave orbital velocity at the bottom. Madsen (1988) proposed a formulation based on the eddy viscosity model yielding:

$$C_b = f_w \frac{g}{\sqrt{2}} U_{rms} \quad (14)$$

in which the friction coefficient $f_w$ is given by Jonsson (1967) as

$$\frac{1}{4\sqrt{f_w}} + \log_{10}\left(\frac{1}{4\sqrt{f_w}}\right) = -0.08 + \log_{10}\left(\frac{a_{b,r}}{K_w}\right) \quad (15)$$

where $a_{b,r}$ is a representative bottom orbital amplitude and $K_w$ denotes the hydraulic roughness scale. SWAN does not include the effect of a mean current on the bottom friction (Booij et al., 1999). According to Tolman (1992), the impact of the current is negligible compared to the error in the estimation of a correct hydraulic roughness scale.

In addition, we implemented the bottom-friction parameterization of Lowe et al. (2005), hereafter LW05, into SWAN. Their work is based on Madsen (1994), in which the bottom-friction coefficient is expressed in terms of a frequency-dependent dissipation factor $f_e(\sigma)$ as

$$C_b = \frac{gU_{rms}}{4} f_e(\sigma) \quad (16)$$

The energy dissipation factor is in turn given by

$$f_e(\sigma) = \sqrt{f_w(\sigma_r)} \sqrt{f_w(\sigma)} \cos\varphi(\sigma) \quad (17)$$

in which the wave friction factor, the phase relation between shear stress and velocity, and the representative radian-frequency are given by

$$f_w(\sigma) = \exp\left[a_1 \left(\frac{U_{rms}}{\sigma K_w}\right)^{a_2} + a_3\right] \quad (18)$$

$$\varphi(\sigma) = 33 - 6.0 \log_{10} \left( \frac{U_{rms}}{K_w \sigma} \right) \quad (19)$$

$$\sigma_r = \int_0^\infty \sigma u_b^2(\sigma) d\sigma / \int_0^\infty u_b^2(\sigma) d\sigma \quad (20)$$

This formulation involves a number of free parameters for calibration. LW05 used the empirical coefficients $a_1=5.5$, $a_2=-0.2$, and $a_3=-6.3$ from Nielsen (1992) to facilitate calculation of the hydraulic roughness scale from measurements. They determined an optimal value of $K_w=0.16$ m through comparison of an energy flux model with observed dissipation rates, in the absence of breaking waves, along a cross-shore transect on a barrier reef at Kaneohe, Hawaii. It is remarkable that physical measurements with a profiler at 14 locations on the reef flat revealed very similar values of $K_w=0.14 \pm 0.03$ m, suggesting nearly homogenous surface roughness over the study area and the consistence of the method itself.

## 2.3. Triad interactions

The triad interactions dominate the nonlinear energy transfer among frequency components in coastal waters with $kd<1$. During shoaling, wave energy cascades down from the spectral peak to the super-harmonics resulting in wave profiles with steeper crests and flatter troughs. The crest loses its vertical symmetry and continues to steepen prior to wave breaking. Eldeberky (1996) provided a description of these effects resulting from the triad interactions through an energy density formulation of the Boussinesq model of Madsen and Sørensen (1993). Although the formulation only includes self interactions of col-linear waves, it provides a reasonable description of the energy transfer to the super-harmonics for general application.

The parameterization adopted in SWAN, known as the lumped triad approximation (LTA) in Eldeberky (1996), relies on an empirical expression for the biphase, derived from the laboratory experiments of Battjes and Beji (1992) and Arcilla et al. (1994). The source term for the triad interactions is written as

$$S_{nl3}(f, \theta) = S_{nl3}^{+}(f, \theta) + S_{nl3}^{-}(f, \theta) \quad (21)$$

in which

$$S_{nl3}^{+}(f, \theta) = \max\left\{0, \alpha_{EB} 2\pi c_c g J^2 \left|\sin(\beta)\right| \left[E^2(f/2, \theta) - 2E(f/2, \theta)E(f, \theta)\right]\right\} \quad (22)$$

$$S_{nl3}^{-}(f, \theta) = -2S_{nl3}^{+}(2f, \theta) \quad (23)$$

where the biphase $\beta$ is given by the parameterization (10), $c$ and $c_g$ are the celerity and group velocity, respectively, $\alpha_{EB}$ is a predefined constant, and $J$ is a coupling coefficient of the triad interactions derived from the energy density formulation of Madsen and Sørensen (1993).

## 2.4. Wave setup

The depth-induced breaking and bottom-friction parameterizations play an important role in defining the coastal wave height, which in turn modifies the water level through the radiation stresses. The water-level increase along the shore, known as wave setup, is important to flood hazard assessment and beach safety. SWAN may provide the radiation stresses to a circulation model for computation of wave-induced currents and water-level variations in the surf zone. This is commonly
done in modeling of storm surge due to hurricane landfalls (e.g., Cheung et al., 2007; Dietrich et al., 2010). Alternatively, SWAN may directly provide an estimate of the wave setup by assuming stationary wave conditions.

Let (x, y) define the Cartesian coordinate system in the horizontal plane. In one dimension, SWAN solves the momentum balance equation, in which the hydrostatic pressure gradient balances the phase-averaged force induced by the waves. The variation of the mean surface elevation $\bar{\eta}$ is expressed in terms of the radiation stress $S_{xx}$ as

$$ \frac{d\bar{\eta}}{dx} = - \frac{1}{\rho g(d + \bar{\eta})} \frac{dS_{xx}}{dx} \qquad (24) $$

For two-dimensional cases, SWAN estimates the variation of the mean surface elevation from the Poisson equation

$$ \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial}{\partial x} \left[ \frac{\rho g(d + \bar{\eta}) \partial \bar{\eta}}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{\rho g(d + \bar{\eta}) \partial \bar{\eta}}{\partial y} \right] = 0 \qquad (25) $$

in which

$$ F_x = - \frac{\partial S_{xx}}{\partial x} - \frac{\partial S_{xy}}{\partial y} \qquad (26) $$

$$ F_y = - \frac{\partial S_{yy}}{\partial y} - \frac{\partial S_{yx}}{\partial x} \qquad (27) $$

are the rotation-free part of the wave-induced forces determined from the radiation stresses in the x-y plane (Dingemans, 1987). The reader may refer to Holthuijsen (2007) for the full expressions of the radiation stresses and for more details on the wave setup calculations in SWAN.

### 3. Comparison with wave flume data

Demirbilek et al. (2007) reported a flume experiment on wave transformation over a fringing reef that allows comparison and calibration of the parameterizations for wave breaking dissipation. The flume, which is 35 m long, 0.7 m wide, and 1.6 m high, features a 1:64 scaled model of a reef typical of the southeast shore of Guam. Fig. 1 shows the locations of 9 wave gages and the profile of the model that includes a composite slope for the reef face, a reef flat, and a beach. The measurements at G4 are not used here as they were found to be inaccurate. The toe of the reef is located at 15.5 m from a plunger-type wave maker, which generates irregular sea states based on the JONSWAP spectrum with random phases. As the reef model was built of polyvinyl chloride (PVC), bottom friction is negligible for the smooth and impervious material. We may assume, at least in the vicinity of the reef crest, that depth-induced breaking is the only dissipation mechanism. Table 1 lists the 19 test cases, which cover the input significant wave height $H_s = 3.2–8.5$ cm, peak period $T_p = 1–2.5$ s, and water depth $h_r = 1.6–5.1$ cm over the reef flat. These correspond to offshore waves with heights and periods up to 5.4 m and 20 s respectively and at least 1 m of water over the reef flat in the prototype.

The five parameterizations for depth-induced wave breaking in SWAN are examined through the 19 test cases. The computational domain extends from G1 to the end of the flume at 1 cm resolution. The incident wave spectrum recorded at G1 provides the boundary condition for the calculation with 42 frequency bins from 0.05 to 3 Hz. SWAN is run in the stationary and one-dimensional Cartesian modes with the triad source term and the setup calculation activated. We first focus on the ability of the five parameterizations in their default settings to predict the significant wave height across the reef. Fig. 2 shows the results from 6 cases representative of the incident wave heights, periods, and water levels considered. The incident waves typically break between G5 and G7 and evolve into turbulent

<figure><img src="image_57.png" /><figcaption>Fig. 1. Schematic and instrument locations for the reef experiment of Demirbilek et al. (2007). The reference elevation at zero corresponds to a water depth of 5.1 cm on the reef flat for illustration. Circles: wave gages.</figcaption></figure>

bores on the reef flat. The model with the various parameterizations performs reasonably well in describing the wave height at the onset of the breaking process and in the surf zone for the low-energy test 20. The discrepancy between the measurements and the computed data increases with the incident wave height. For the energetic wave conditions in tests 19, 48, and 32, all the parameterizations with default settings result in excessive energy dissipation prior to wave breaking leading to underestimation of the wave height at the onset near the reef crest. The model over-predicts the extent of the surf zone resulting in excessive dissipation and underestimation of the energy on the reef flat. The increase of the measured wave height from G8 to G9 in tests 48 and 32 might be attributed to reflection from the beach and formation of partial standing waves.

The overall performance of the five parameterizations in the 19 tests is assessed with four metrics. The Root Mean Square Error (RMSE), the Normalized Root Mean Square Error (NRMSE), the bias

<caption>Table 1</caption>

Water depth over the reef flat and input wave conditions for the laboratory tests of Demirbilek et al. (2007).

<table><thead><tr><td>Test</td><td>h<sub>r</sub> (cm)</td><td>H<sub>s</sub> (cm)</td><td>T<sub>p</sub> (s)</td></tr></thead><tbody><tr><td>15</td><td>5.1</td><td>6.2</td><td>1.00</td></tr><tr><td>20</td><td>5.1</td><td>6.1</td><td>1.25</td></tr><tr><td>16</td><td>5.1</td><td>5.2</td><td>1.50</td></tr><tr><td>17</td><td>5.1</td><td>7.8</td><td>1.50</td></tr><tr><td>21</td><td>5.1</td><td>8.2</td><td>1.75</td></tr><tr><td>18</td><td>5.1</td><td>8.5</td><td>2.00</td></tr><tr><td>19</td><td>5.1</td><td>8.3</td><td>2.50</td></tr><tr><td>26</td><td>1.6</td><td>5.8</td><td>1.00</td></tr><tr><td>27</td><td>1.6</td><td>5.5</td><td>1.25</td></tr><tr><td>28</td><td>1.6</td><td>4.7</td><td>1.50</td></tr><tr><td>29</td><td>1.6</td><td>7.1</td><td>1.50</td></tr><tr><td>30</td><td>1.6</td><td>7.6</td><td>1.75</td></tr><tr><td>31</td><td>1.6</td><td>8.5</td><td>2.00</td></tr><tr><td>32</td><td>1.6</td><td>7.9</td><td>2.50</td></tr><tr><td>44</td><td>3.1</td><td>3.2</td><td>1.00</td></tr><tr><td>45</td><td>3.1</td><td>6.1</td><td>1.00</td></tr><tr><td>46</td><td>3.1</td><td>5.9</td><td>1.25</td></tr><tr><td>47</td><td>3.1</td><td>5.0</td><td>1.50</td></tr><tr><td>48</td><td>3.1</td><td>7.5</td><td>1.50</td></tr><tr><td>57</td><td>3.1</td><td>7.7</td><td>1.75</td></tr><tr><td>58</td><td>3.1</td><td>8.5</td><td>2.00</td></tr><tr><td>59</td><td>3.1</td><td>8.2</td><td>2.50</td></tr></tbody></table>
<figure><img src="image_58.png" /><figcaption>**Fig. 2.** Comparison of observed and computed significant wave heights for 6 representative test cases. Black squares: measurements; black line: BJ78 with $\gamma=0.73$ and $B=1$ (default); green line: TG83; magenta line: NE87; blue line: RU03; orange line: WE10.</figcaption></figure>

($B$), and the normalized bias ($NB$) between the modeled and observed significant wave heights ($H_m, H_o$) at each gage are defined as

$$RMSE(m) = \sqrt{\frac{\sum_{i=1}^{n} (H_{m,i} - H_{o,i})^2}{n}} \qquad (28)$$

$$NRMSE(\%) = 100\% \times \sqrt{\frac{\sum_{i=1}^{n} (H_{m,i} - H_{o,i})^2}{\sum_{i=1}^{n} (H_{o,i})^2}} \qquad (29)$$

$$B(m) = \frac{\sum_{i=1}^{n} (H_{m,i} - H_{o,i})}{n} \qquad (30)$$

$$NB(\%) = 100\% \times \frac{\sum_{i=1}^{n} (H_{m,i} - H_{o,i})}{\sum_{i=1}^{n} (H_{o,i})} \qquad (31)$$

where $n=19$ denotes the number of tests. Table 2 summarizes the results from the error analysis for the five parameterizations. BJ78 yields the lowest errors overall. This outcome is unexpected as, for instance, NE87 seems to be more suitable because it is a refined version of BJ78 featuring the effect of the bottom slope. The reason is that the breaking index (5) in NE87 is capped at 0.81 in SWAN corresponding to a bottom slope of 0.01 (Booij et al., 1999), which is much smaller compared to that on the reef face in the experiment. However, uncapping the breaker index would lead to a non-physical value of $\gamma=1.35$ for the reef slope of 0.106 in the experiment. WE10 appears to be a good candidate because the parameterization of the biphase is expected to relate the wave nonlinearity and the onset of wave breaking more precisely. However, WE10 as well as TG83 and RU03 initiate the breaking process earlier than BJ78 resulting in a strong bias of the predictions across the reef.

The results from BJ78, which already have the lowest errors, might be further improved by increasing the breaker index $\gamma$. This delays the breaking onset and reduces the width of the surf zone. However, waves tend to break as a plunger over the reef crest and propagate as

turbulent bores over the flat (Roeber et al., 2010a, 2010b). The breaking process is more dissipative and the breaking intensity $B$ needs adjustment as well. We conducted a sensitivity analysis by varying $\gamma$ and $B$.

<table><caption>**Table 2**

Error metrics for depth-induced wave breaking parameterizations.</caption><thead><tr><th>Parameterization</th><th>Gage</th><th>RMSE (cm)</th><th>NRMSE (%)</th><th>BIAS (cm)</th><th>NBIAS (%)</th></tr></thead><tbody><tr><td rowspan="7">BJ78<br>(γ=0.73, B=1)</td><td>G2</td><td>0.10</td><td>1.43</td><td>0.07</td><td>1.00</td></tr><tr><td>G3</td><td>0.16</td><td>2.37</td><td>0.09</td><td>1.43</td></tr><tr><td>G5</td><td>0.43</td><td>5.95</td><td>-0.30</td><td>-4.39</td></tr><tr><td>G6</td><td>0.63</td><td>8.97</td><td>-0.57</td><td>-8.28</td></tr><tr><td>G7</td><td>0.65</td><td>13.41</td><td>-0.56</td><td>-12.00</td></tr><tr><td>G8</td><td>0.92</td><td>32.67</td><td>-0.79</td><td>-30.07</td></tr><tr><td>G9</td><td>1.24</td><td>42.95</td><td>-1.03</td><td>-38.56</td></tr><tr><td rowspan="7">TG83</td><td>G2</td><td>0.10</td><td>1.43</td><td>0.07</td><td>1.00</td></tr><tr><td>G3</td><td>0.16</td><td>2.36</td><td>0.09</td><td>1.41</td></tr><tr><td>G5</td><td>0.97</td><td>13.41</td><td>-0.77</td><td>-11.15</td></tr><tr><td>G6</td><td>1.46</td><td>20.68</td><td>-1.36</td><td>-19.81</td></tr><tr><td>G7</td><td>1.80</td><td>37.02</td><td>-1.69</td><td>-36.40</td></tr><tr><td>G8</td><td>1.25</td><td>44.56</td><td>-1.16</td><td>-44.07</td></tr><tr><td>G9</td><td>1.62</td><td>56.04</td><td>-1.44</td><td>-54.06</td></tr><tr><td rowspan="7">NE87</td><td>G2</td><td>0.10</td><td>1.43</td><td>0.07</td><td>1.00</td></tr><tr><td>G3</td><td>0.16</td><td>2.37</td><td>0.09</td><td>1.43</td></tr><tr><td>G5</td><td>0.38</td><td>5.30</td><td>-0.26</td><td>-3.77</td></tr><tr><td>G6</td><td>0.64</td><td>9.00</td><td>-0.57</td><td>-8.35</td></tr><tr><td>G7</td><td>1.10</td><td>22.68</td><td>-0.98</td><td>-20.98</td></tr><tr><td>G8</td><td>1.06</td><td>37.80</td><td>-0.94</td><td>-35.75</td></tr><tr><td>G9</td><td>1.34</td><td>46.37</td><td>-1.13</td><td>-42.42</td></tr><tr><td rowspan="7">RU03</td><td>G2</td><td>0.10</td><td>1.43</td><td>0.07</td><td>1.00</td></tr><tr><td>G3</td><td>0.16</td><td>2.37</td><td>0.09</td><td>1.42</td></tr><tr><td>G5</td><td>0.80</td><td>11.02</td><td>-0.51</td><td>-7.35</td></tr><tr><td>G6</td><td>1.20</td><td>16.95</td><td>-1.05</td><td>-15.24</td></tr><tr><td>G7</td><td>1.70</td><td>35.01</td><td>-1.58</td><td>-33.92</td></tr><tr><td>G8</td><td>1.53</td><td>54.49</td><td>-1.39</td><td>-52.71</td></tr><tr><td>G9</td><td>1.80</td><td>62.49</td><td>-1.58</td><td>-59.06</td></tr><tr><td rowspan="7">WE10</td><td>G2</td><td>0.10</td><td>1.43</td><td>0.07</td><td>1.00</td></tr><tr><td>G3</td><td>0.16</td><td>2.37</td><td>0.09</td><td>1.43</td></tr><tr><td>G5</td><td>0.80</td><td>11.14</td><td>-0.48</td><td>-6.88</td></tr><tr><td>G6</td><td>0.88</td><td>12.48</td><td>-0.78</td><td>-11.36</td></tr><tr><td>G7</td><td>0.67</td><td>13.80</td><td>-0.60</td><td>-12.98</td></tr><tr><td>G8</td><td>1.44</td><td>51.32</td><td>-1.31</td><td>-49.47</td></tr><tr><td>G9</td><td>1.92</td><td>66.51</td><td>-1.67</td><td>-62.56</td></tr></tbody></table>
<table><caption>**Table 3**
Error metrics for the depth-induced wave breaking parameterization BJ78 with $\gamma = 0.94$ and $B = 1.09$.</caption><thead><tr><th>Gage</th><th>RMSE (cm)</th><th>NRMSE (%)</th><th>BIAS (cm)</th><th>NBIAS (%)</th></tr></thead><tbody><tr><td>G2</td><td>0.10</td><td>1.44</td><td>0.07</td><td>1.01</td></tr><tr><td>G3</td><td>0.16</td><td>2.38</td><td>0.10</td><td>1.45</td></tr><tr><td>G5</td><td>0.24</td><td>3.30</td><td>-0.12</td><td>-1.85</td></tr><tr><td>G6</td><td>0.37</td><td>5.25</td><td>-0.08</td><td>-1.12</td></tr><tr><td>G7</td><td>0.37</td><td>7.62</td><td>0.06</td><td>1.39</td></tr><tr><td>G8</td><td>0.57</td><td>20.17</td><td>-0.30</td><td>-11.19</td></tr><tr><td>G9</td><td>0.85</td><td>29.62</td><td>-0.55</td><td>-20.72</td></tr></tbody></table>

from 0.7 to 1.3 at 0.01 increments for each test case. The optimal values of $\gamma = 0.94$ and $B = 1.09$ minimize the average NRMSE at G5, G6, G7 and G8 around the reef crest, where intensive wave breaking occurs. This confirms, *a posteriori*, the logic of the optimization. Table 3 shows an appreciable improvement of the error metrics in comparison to the default BJ78 parameterization in SWAN. Fig. 3 shows that the optimal BJ78 parameterization gives better overall agreement with the observed wave height in comparison to the default setting. Most importantly, the present parameterization captures the onset of the breaking process at the reef crest as observed in the experiment. The breaking process and the subsequent wave height attenuation defines the wave setup in the surf zone. Fig. 4 shows the computed and observed wave setup at G8 over the reef flat and Table 4 summarizes the error metrics. The optimal BJ78 parameterization represents an improvement over the default setting to give the best agreement with the observed data. The other parameterizations underestimate the setup comparing to the BJ78 default setting because of the wider predicted surf zones. The optimal values of $\gamma$ and $B$ provide a convenient setting in BJ78 for practical applications with general wave conditions in fringing reef environments, even though their values as well as other parameterizations can be optimized on a case by case basis (Massel and Gourlay, 2000; Apotsos et al., 2008; Sheremet et al., 2011).

The spectral transformation provides further insights into the wave processes. Fig. 5 compares the observed spectra across the reef with the computed results from the optimal BJ78 parameterization

for representative conditions in the experiment. Test 44 with input $H_s = 3.4$ cm and $T_p = 1$ s represents low energy wave conditions with the breaking onset at the reef crest between G6 and G7, while the more energetic test 32 with $H_s = 8.4$ cm and $T_p = 2.5$ s produces breaking waves between G5 and G6, before reaching the reef crest. The wave breaking parameterization captures the energy reduction across the reef in both cases. Nonlinear energy transfer is evident as the waves shoal over the reef slope with development of super and sub-harmonics. The triad source term reproduces the growth of the super-harmonics at twice the peak frequency, but does not describe the sub-harmonics due to the lack of proper parameterization in SWAN. The computed spectrum at G8 shows the same peaks as it does before breaking, while the

<figure><img src="image_59.png" /><figcaption>**Fig. 4.** Scatter plot for wave setup at G8. Red crosses: optimal BJ78 with $\gamma = 0.94$ and $B = 1.09$; Black squares: BJ78 with $\gamma = 0.73$ and $B = 1$ (default); green asterisks: TG83; magenta circles: NE87; blue diamonds: RU03; orange triangles: WE10.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>**Fig. 3.** Comparison of observed and computed significant wave heights for 6 representative test cases. Black squares: measurements, black line: BJ78 with $\gamma = 0.73$ and $B = 1$ (default), red line: optimal BJ78 with $\gamma = 0.94$ and $B = 1.09$.</figcaption></figure>
<table><caption>**Table 4**
Error metrics of computed wave setup at G8 on the reef flat.</caption><thead><tr><th></th><th>RMSE (cm)</th><th>NRMSE (%)</th><th>BIAS (cm)</th><th>NBIAS (%)</th></tr></thead><tbody><tr><td>TG83</td><td>0.44</td><td>10.57</td><td>-0.38</td><td>-42.82</td></tr><tr><td>NE87</td><td>0.26</td><td>6.18</td><td>-0.23</td><td>-25.48</td></tr><tr><td>RU03</td><td>0.39</td><td>9.43</td><td>-0.33</td><td>-36.34</td></tr><tr><td>WE10</td><td>0.24</td><td>5.70</td><td>-0.21</td><td>-22.86</td></tr><tr><td>BJ78 (γ=0.73, B=1)</td><td>0.23</td><td>5.56</td><td>-0.21</td><td>-23.48</td></tr><tr><td>BJ78 (γ=0.94, B=1.09)</td><td>0.14</td><td>3.42</td><td>-0.13</td><td>-13.99</td></tr></tbody></table>

observed spectrum broadens with a significant part of the energy shifted toward the infra-gravity band. For the same set of laboratory data, Sheremet et al. (2011) managed to reproduce both the sub and super-harmonics with fully nonlinear, triad interaction models.

4. Comparison with field data

With the depth-induced breaking processes calibrated, we focus our attention on dissipation due to bottom friction in natural reef environments. The US Army Corps of Engineers PILOT project provided wave measurements across a fringing reef off Mokuleia beach, Oahu, Hawaii. Fig. 6 shows the site location, bathymetry, and instrument transect profile. The site is open to north swells during the winter months. The open ocean wave conditions are recorded by Buoy 51201 at 200 m water depth. The bathymetry of the Hawaiian Islands is well defined by numerous LiDAR and multi-beam surveys. The slope of the insular shelf on Oahu's north shore is quite gentle from the 200-m depth contour to the shore. The reef, which starts at approximately 20 m water depth, rapidly rises to 10 m and the slope becomes gentler until the reef crest. Six sensors were deployed along a transect across the reef to 20 m depth.

<figure><img src="image_61.png" /><figcaption>**Fig. 5.** Comparison of observed and computed wave energy spectra across the reef for Test 44 (low energy) and Test 32 (high energy). Black line: measurements; red line: optimal BJ78 with $\gamma=0.94$ and $B=1.09$.</figcaption></figure>
<figure><img src="image_62.png" /><figcaption>**Fig. 6.** Site information for the field study. (a) Study site on North Shore, Oahu. (b) Study site at Mokuleia. (c) Instrument transect. Square: buoy 51201; circles: coastal wave sensors.</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>**Fig. 7.** Recorded offshore wave conditions at buoy 51201 during the field study.</figcaption></figure>

Two of the sensors were located in a shallow lagoon of approximately 1.5 m deep between the reef crest and the shore. Although the profile differs from the laboratory model of Demirbilek et al. (2007), the reef slopes in front of the crest are comparable to produce similar depth-induced breaking processes.

The six sensors provided continuous measurements of the coastal wave conditions from December 28, 2007 to January 8, 2008, during which Buoy 51201 recorded a series of wave events. The recorded significant wave height, peak period, and direction from the National Data Buoy Center are shown in Fig. 7. A north swell event with $H_s=2–2.5$ m and $T_p=12–16$ s occurred during the first three days. This is followed by waves with $H_s=2.5–3.5$ m and $T_p=10–12$ s due to a system passing through from the north. Strong trade winds were developed on January 3 and 4 augmenting the significant wave height with wind seas. The significant wave height then slowly decreased to about 1 m toward the end of the measurements, while the wave period remained relatively constant at $T_p=11–13$ s. The site is open to these wave events, which had peak directions varying

<figure><img src="image_64.png" /><figcaption>**Fig. 8.** Comparison of recorded and computed significant wave heights at Mokuleia. Blue dots: observations; black line: JONSWAP; green line: Madsen (1988); magenta line: Collins (1972); and red line: LW05 with $K_w = 0.16$ m.</figcaption></figure>
slightly about the north. The instruments at the site recorded a set of 270 semi-hourly spectra and wave parameters for examination of the model parameterizations. The computation utilizes two levels of nested grids covering the regions shown in Fig. 6a and b. The level-1 grid extends from buoy 51201 to the shore at 1 second resolution, while the level-2 grid covers the instrument transect from 30 m water depth to the shore at 3 m resolution. The spectral data at the Buoy 51201 provides the input conditions to the level-1 grid, which in turn defines the boundary conditions for the level-2 grid. The global tidal model TPXO.6 of Egbert and Erofeeva (2002) defines the water level for the computation. SWAN was run in the non-stationary and Cartesian modes with 36 direction bins and 72 frequency bins from 0.005 to 1 Hz. The optimal BJ78 parameterization for depth-induced wave breaking, the LTA parameterization for the triad interactions, and the wave setup calculation are activated in the level-2 grid.

We first examine the three parameterizations for bottom friction in SWAN as well as Lowe et al. (2005) in their default settings. The significant wave height is estimated from integration of the computed energy density up to 0.2 Hz consistent with the cut-off frequency of the observed spectra. Fig. 8 compares the observed and computed significant wave heights at buoy 51201 and the six sensors at the site. The good agreement between the computed and recorded data at buoy 51201 indicates that the offshore wave conditions are properly input into SWAN. The wave height decreases slightly from the buoy to A7 corresponding to the initial stage of the shoaling process. An average peak period of 12 s during the case study gives a depth parameter of $kd \approx 0.8$ at A7. The wave height should increase shoreward of A7 given the shore-normal, peak wave direction, but the observations and computed results from A7 to A4 show little variation, suggesting a balance between shoaling and dissipation by bottom friction. The model trends to give slightly larger wave heights on the reef slope during the trade wind event on January 3 and 4. Although the north swell is dominant as indicated by the peak wave direction, the sea state is bimodal with a component from the east to northeast associated with the trade winds.

While buoy 51201 recorded the full event from an offshore location for model input, the northeast headland of Oahu might have partially shielded the site from the wind waves. This is evident in the reduced signature of the wind waves at the sensors in relation to that at buoy 51201. In contrast, other features recorded at the buoy extend across the reef slope with little attenuation.

The waves break and their height decreases dramatically across the reef crest between A4 and A2. The depth-limited wave heights computed at A2 and A1, which are independent of the offshore conditions, reflect the tide signals quite clearly. The four parameterizations give rise to very similar wave heights across the reef and provide good agreement with the observed data up to A4, but overestimate the significant wave height at A2 and A1 in the shallow lagoon. This is expected for the three default bottom friction parameterizations in SWAN as they were calibrated for sand beds, which are much smoother than the reef at Mokuleia. The LW05 parameterization also over-predicts the significant wave heights at A2 and A1, indicating that the bottom roughness at Mokuleia is likely larger than that at Kaneohe. A closer look at the results suggests that a hydraulic roughness scale of 0.16 m yields acceptable dissipation rates over the reef slope from A7 to A4, but appears to be too small shoreward of A4. Pawlak and MacCready (2002) and Bandet (2009) identified different turbulent flow regimes and dissipation rates over large bedform roughness depending on the wave orbital amplitude and the roughness scale. As stressed by Tolman (1992), most of the errors in the bottom friction parameterization are related to estimation of the appropriate bottom roughness. A single roughness scale thus cannot fully describe the dissipation due to bottom friction at this site.

An aerial photo of the site in Fig. 9 shows an abrupt change of the surface texture and color between A4 and A2 at approximately 5 m depth. The substrate transitions from bare reefs covered by thin layers of light color sand to rough reefs with branching and fossil corals (Jaramillo and Pawlak, 2011). This corresponds to an increase of the bottom roughness from the reef slope to the crest and the inner region.

<figure><img src="image_65.png" /><figcaption>Fig. 9. Aerial image of the reef and sensor locations at Mokuleia from Google Earth.</figcaption></figure>
Despite being qualitative, the observed bottom roughness suggests that the model results might be improved by using a higher roughness scale for the nearshore region. Since our knowledge on the spatial variability is limited, we assume an abrupt transition at the 5-m contour from $K_w = 0.16$ m to a higher value. Through a sensitivity analysis, we obtained $K_w = 0.5$ m in the nearshore region to match the recorded data at A2 and A1. Fig. 10 compares the recorded and computed significant wave heights with and without the abrupt increase of the nearshore roughness scale. The energy dissipation from A4 to A2 is from a combination of wave breaking and bottom friction, but as the offshore significant wave height decreases to less than 1 m toward the end the case study, bottom friction becomes the dominant process. The use of $K_w = 0.5$ m appears to work well over the range of significant wave heights during the event, rendering support for both the wave-breaking and bottom-friction parameterizations. The use of the large $K_w$ here is not to offset underestimations of wave breaking dissipation as verified toward the end of the case study when wave breaking becomes unimportant and bottom friction dominates the dissipation from A4 to A2. In addition, the selected roughness scale correctly accounts for the dissipation from A2 to A1 with no wave

breaking involved. This bottom roughness scale is higher than that at Kaneohe, but corroborates visual observations made at the site and measurements on the south shore of Oahu by Nunes and Pawlak (2008).

The comparison with the flume experiment in Section 3 has already shown that an accurate description of wave height attenuation across the reef is important for estimation of the wave setup. We thus utilize the composite roughness scale of $K_w = 0.16$ and 0.5 m in the LW05 parameterization to examine the wave setup in natural reef environments. Fig. 11 compares the computed and recorded wave setup at A1 for the entire case study and along the instrument transect during the rise, peak, and decay of the swell event from December 31, 2009 to January 5, 2008. The observed wave setup is estimated from the surface elevation time series using the global tidal model TPXO.6 of Egbert and Erofeeva (2002). This might introduce errors into the deduced data as the global tidal model does not resolve island-scale processes. Nevertheless, the model reproduces the observed wave setup reasonably well especially for the energetic events in which the signal to noise level is high. The swell reaches 1.95 m significant wave height at buoy 51201 on December 31. The waves break at the reef crest causing an increase of the mean water level in the lagoon. The peak arrives in the morning of January 2 with an offshore wave height of 3 m and produces breaking waves on the upper reef slope resulting in a wider surf zone. The setup retreats

<figure><img src="image_66.png" /><figcaption>**Fig. 10.** Comparison of recorded and computed significant wave heights at Mokuleia. Blue dots: observations; red line: LW05 with $K_w = 0.16$ m; green line: LW05 with $K_w = 0.16$ and 0.5 m.</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>**Fig. 11.** Comparison of recorded and computed wave setup at Mokuleia. (a) Wave setup at A1 during the entire case study. (b) Wave setup on 31/12/2007 10:34 am. (c) Wave setup on 1/2/2008 6:34 am. (d) Wave setup on 1/5/2008 1:34 am. (e) Reef profile along the instrument transect. Blue dots and circles: recorded data; green line: LW05 with $K_w = 0.16$ and 0.5 m.</figcaption></figure>