The Performance of a Spectral Wave Model at Predicting Wave Farm Impacts

For renewable ocean wave energy to support global energy demands, wave energy converters (WECs) will likely be deployed in large numbers (farms), which will necessarily change the nearshore environment. Wave farm induced changes can be both helpful (e.g., beneficial habitat and coastal protection) and potentially harmful (e.g., degraded habitat, recreational, and commercial use) to existing users of the coastal environment. It is essential to estimate this impact through modeling prior to the development of a farm, and to that end, many researchers have used spectral wave models, such as Simulating WAves Nearshore (SWAN), to assess wave farm impacts. However, the validity of the approaches used within SWAN have not been thoroughly verified or validated. Herein, a version of SWAN, called Sandia National Laboratories (SNL)-SWAN, which has a specialized WEC implementation, is verified by comparing its wave field outputs to those of linear wave interaction theory (LWIT), where LWIT is theoretically more appropriate for modeling wave-body interactions and wave field effects. The focus is on medium-sized arrays of 27 WECs, wave periods, and directional spreading representative of likely conditions, as well as the impact on the nearshore. A quantitative metric, the Mean Squared Skill Score, is used. Results show that the performance of SNL-SWAN as compared to LWIT is “Good” to “Excellent”.


Introduction
Renewable energy from ocean waves has the potential to play a large role in the world's future energy supply; the theoretically available wave power resource is around 2.1 TW [1]. For wave energy to be extracted economically at a large scale, wave energy converters (WECs) need to be deployed in arrays, or farms, consisting of many devices clustered together, facilitating operations and maintenance, and maximizing the use of infrastructure, like subsea cables.
By the conservation of energy, extracting energy from the wave field means that there must be a net reduction of wave heights in the region of the array. It is well understood within the field of coastal engineering that variations in wave heights can result in wave-driven currents and sediment transport. As such, wave farm driven environmental changes have the potential to both positively (e.g., coastal resilience) and negatively (e.g., adverse habitat change) impact adjacent beaches and coastal ecosystems. Evaluation tools must be able to optimize array designs to balance energy generation with environmental impacts.
Researchers have investigated the impact of wave farms on nearshore processes, often using phase-averaged, spectral numerical models, such as SWAN (Simulating WAves Nearshore) [2]. seas. In McNatt et al. [20], this was extended to small WEC arrays. Qualitative and quantitative comparisons were made of the wave field in what is typically referred to as the "far field", that is, at some distance from the array. Findings showed that WEC reflection, which is modeled in WAMIT, but not in SNL-SWAN resulted in higher SNL-SWAN error, and that the error also was higher for regular waves and unidirectional waves.
The short-comings of previous comparative studies of WECs in spectral wave models are that they considered only small WEC arrays, and the verification cases consisted of only the near-array region modeled over a flat bathymetry.
Here, SNL-SWAN is further verified by comparison with linear wave-body theory through additional considerations aimed at making the study more representative of real conditions: 1. Medium-sized arrays of 27 WECs are modeled, representing a medium-sized wave farm.
Smaller arrays of 5 WECs were compared in McNatt et al. [20], and, while larger arrays could be considered in the future, an array of 27 WECs represents a medium-sized wave farm and is a stepping stone in both the development of commercial wave farms and the verification of SNL-SWAN. To model these arrays, an approach referred to as Linear Wave Interaction Theory (LWIT) was used. 2. Only directionally spread sea conditions are modeled, based on literature on the range of realistic directional spreading [21]. The range of directional spreading includes highly focused storm waves, however it does not include idealized purely unidirectional seas. In previous studies, SNL-SWAN did not show good agreement with linear wave theory for unidirectional seas because without the SWAN diffraction option on, it is not capably of wave diffraction [19]. However, directional spreading is a form of wave diffraction and unidirectional seas are not representative of real sea conditions. 3. A realistic beach bathymetry is used. In previous studies, the authors considered only a flat sea bottom. However, the aim of SNL-SWAN is to assess the impact of wave farms on nearshore processes. As such, the wave fields produced by each method were fed in as boundary conditions for a SWAN model that propagates the wave field to shore.
The verification of SNL-SWAN by comparison with LWIT is appropriate because LWIT solves the linear-wave boundary value problem and is a type of 'wave-structure interaction solver' referred to by Stratigaki et al. [14]. Linear-wave-body theory is widely accepted as a good method for computing wave-body interactions [22], and it has been validated experimentally (e.g., Reference [23]); LWIT in particular has been verified by code-to-code comparison (e.g., Reference [24]).
The novelty in this approach includes the verification of a spectral wave model, again, a medium-sized array (previous work considered at most 5 WECs), the focus on realistic wave conditions, and the verification over beach (non-flat) bathymetry.
The results produced by SNL-SWAN and LWIT are compared as plots of the wave field, as well as a quantitative comparative metric: the Mean Squared Skill Score (MSSS). Results show that the main relative errors between SNL-SWAN and LWIT due to wave reflection in shorter waves, which are modeled by LWIT but not in SNL-SWAN. Despite this, the performance of SNL-SWAN when compared to LWIT is "Good" to "Excellent".

SNL-SWAN
SWAN is a spectral (phase-averaged) wave transformation model acting over the two horizontal spatial dimensions that can be used to generate wind waves and transform wave conditions to a nearshore area. SWAN simulates wind-wave generation, transformation on the current, refraction, shoaling, nonlinear wave-wave interactions, and energy dissipation through white-capping and breaking. Obstacles in SWAN can absorb or reflect energy, but the program does not model wave interactions between bodies.
Implementing WECs within SNL-SWAN is done by introducing obstacles within the model domain. The obstacles are set to be approximately equal to the width of the device in the direction of wave travel. The computational grid resolution is set so that the width can be closely approximated. To effectively capture the width of the device, the minimum recommended computational grid spacing is partially dependent on the ratio of the device size to the array spacing. For the arrays used in this analysis, the computational grid spacing was selected as 1/2 the device diameter in order to resolve the devices in each row of the array. Thus, a dense computational grid is typically required near the array; this could be nested inside a larger environmental grid, but for the purposes of this study, it was not.
Transmission of wave energy through the obstacle is set to be frequency dependent. The amount of energy removed from the incident wave spectra at the seaward side of the obstacle is set based on a WEC's power capture curve which is calculated from linear wave theory and WAMIT. The amount of energy absorbed by the WEC at each frequency is set in the SWAN obstacle so that both linear wave theory calculation and SNL-SWAN simulations remove equal amounts of energy, at each frequency, from the wave spectra.

Linear Wave Interaction Theory
Although not presented explicitly as a conclusion, a key outcome of the authors' previous work [20] was that, due to the computation time required, WAMIT was not practical for computing the wave fields of large arrays. Instead, herein, the so-called Linear Wave Interaction Theory (LWIT) is used to compute the linear-wave array wave fields.
LWIT was developed by Kagemoto and Yue [25], who stated that the wave interactions between multiple, separated, floating bodies could be computed by knowing a set of characteristics or wave-amplitude coefficients, which describe (1) the waves that each body radiates while moving and (2) the waves that each body scatters due to incident waves. The wave coefficients of each body can be found in a variety of ways: analytically for cylindrical shapes (e.g., Reference [26]), numerically [27], and from numerical or physical measurements of the wave field [28].
As an example, using an array of two bodies: the waves radiated and scattered by Body A become incident waves on Body B, which in turn scatters those waves and moves radiating waves; both the scattered and radiated waves of body B then become incident waves on Body A, and so on.
Unknown amplitudes of the scattered and radiation waves, as well as of the device forces and motions, can be found in a matrix computation. This means that knowing each body's wave coefficients allows for the computation of hydrodynamic array interactions without precise knowledge of each body's geometry. This is opposed to the direct method of array computation by BEM codes (e.g., WAMIT) in which the exact surface of each body in the array is represented by panels. There are fewer unknown values in the Interaction Theory computation than in the direct BEM computation, which is important because matrix solution computational time is proportional to the square of the number of unknowns. Typically, the number of unknown wave amplitudes (LWIT) is two orders of magnitude less than the number of unknown panels amplitudes (BEM), which means that LWIT computation is on the order of 10,000 times faster than a direct computation [28].
An important aspect of the LWIT is that it is equally as accurate as a direct BEM computation as long as spacing requirements are met, and a sufficient number of wave coefficients are included. Figure 1 shows a comparison of LWIT computation with direct BEM computation using WAMIT from McNatt [28]; except very close to each WEC, there is excellent agreement-the errors in the Interaction Theory compared to WAMIT are less than 0.1% for most of the wave field.

WEC Arrays
A simple canonical WEC was modeled for the purposes of the study. It is a heaving point absorber consisting of a spherical buoy of diameter, d = 10 m, connected via a taught mooring line to the sea bed, Figure 2. Power is extracted through the heave motion of the buoy. The performance of the WEC was modeled with linear wave theory [22], where the hydrodynamic coefficients were computed with the boundary element method (BEM) code WAMIT [16], in which only the submerged hemisphere of the WEC is modeled. WEC power absorption is characterized by the relative capture width (RCW), which is defined as the ratio of the power absorbed by the WEC to the incident power that passes across the WEC's width.
where ω is the radial wave frequency, d is the device width or diameter, D pto is the Please define if appropriate. damping, χ(ω) is the WEC's motion response, ρ is the water density, g is gravitational constant, and c g (ω, h) is the wave group velocity, which is a function of ω and water depth, h. The WEC's RCW curve is shown in Figure 3. The WEC's estimated annual energy production at the European Marine Energy Centre (EMEC)'s Billia Croo test site is 520 MWh, which is equivalent to a 150 kW device with a 40% capacity factor.
Medium-sized wave farms were created with arrays of 27 WECs arranged in two array configurations: (a) two rows and (b) a cluster; see Figure 4. The spacing between the WECs in each row and between the rows is set to either 5d or 10d where d is the device diameter. For a farm of 27 WECs, each with a rated power of 150 kW, the farm size is 4 MW.

Environment
The environment is defined by the bathymetry and wave conditions. The bathymetry consists of a flat bottom at a water depth of 60 m in the region of the WEC arrays, extending to 1 km in the lee of the array. The flat bottom leads into an idealized beach profile that is uniform in the longshore direction and follows a equilibrium beach profile over 3 km to the shoreline [29], as in Figure 5. The WEC array is located at x = 0. Idealized bathymetry was selected in order to simplify the problem and view and isolate differences between the the nearshore results. Figure 5. The bathymetry modeled in the study consisted of a flat bottom in the region of the WEC array (WEC array centered at x = 0), and an idealized beach profile [29] in the nearshore.
Wave conditions were modeled as directional spectra at three periods and three spreading parameters. The spectral shape considered was a Pierson-Moskovitz spectrum, where the relationship between the peak period, T p , and the significant wave height, H s , is given by: H s = 0.04T 2 p . Three peak periods were modeled, T p = 6, 8, 12 s representing short, medium, and long waves, respectively, at potential WEC deployment sites in the U.S. [30].
The directional spreading was modeled using a cosine-squared spreading function, parameterized by the spreading coefficient, s [31]. Three spreading values were modeled, s = 2, 10, 40, representing high directional spreading, moderate spreading, and highly unidirectional waves. The selection of these values was based on correlating the cosine-squared spreading value, s, against ranges of the spreading parameter, φ = 0.74, 0.87, 0.95, representing high, moderate, and low directional spreading from Forristall and Ewans [21].

Computational Domains
The regions under consideration were divided into two computational domains (see Figure 6):  The Array Domain has two boundary conditions (BCs), labeled in Figure 6 as (1) and (2). BC (1) is set to be the incident wave spectra as specified in Section 2.3. However, it should be noted that this BC is only relevant in SNL-SWAN. BCs are treated differently in linear wave theory [22]. BC (2) forms the connection between the Array Domain and the Nearshore Domain. Significant wave heights extracted from the Array domain in both SNL-SWAN and LWIT are set as the input BC for SWAN in the Nearshore Domain.
In the Nearshore Domain, BC (3) was treated to represent the impact of nearshore processes of shoaling and wave breaking as the offshore conditions moved shoreward. This was done by carrying out a simulation on the Nearshore Domain where BC (2) was simply the offshore incident wave conditions (i.e., no WEC array present) and then extracting the wave heights along the centre of the domain and setting these to BC (3).
BC (4) is the zero depth shoreline which has no incident waves.

Means of Comparison
Wave fields produced by SNL-SWAN are compared to those produced by LWIT using the significant wave height, H s , because H s is the driver of nearshore physical processes. Several quantitative metrics are used. Measures referring to SNL-SWAN wave fields are given the subscript SNL; those referring to LWIT wave fields are given the subscript LW IT; incident significant wave height is specified as H s,i .

Normalized H s
The normalized H s is the ratio of H s at a field point and the incident H s,i : Values ofĤ s = 1 mean that the field H s is equal to the incident H s ; values greater than 1 indicate wave heights are higher than the incident, and values less that 1 indicate lower wave heights.

Normalized Error
The normalized Error is the measure of the difference between the SNL-SWAN wave field and the LWIT wave field normalized by the incident H s,i : Values of equal to 0 indicate an exact match between SNL-SWAN and LWIT; values greater than 0 indicate that SNL-SWAN over-predicts the wave height compared to LWIT, and vice versa for values less than 0.

Mean-Squared Error
The Mean-Squared Error (MSE) is a measure of the error over a region of N points: Values of MSE have the same mean as that of but over a region.

Mean-Squared Skill Score
A measure of the "skill" of a model is described in Sutherland et al. [32] as "non-dimensional measure of the accuracy of a prediction relative to the accuracy of a baseline prediction, which could, for example, be an initial condition, an average value, a random choice (from the possible range of outcomes) or a simple empirical predictor".
Herein, a Mean-Squared Skill Score (MSSS) was defined as a comparison of SNL-SWAN against the LWIT relative to the baseline case, in which SNL-SWAN results were compared to an undisturbed wave field (i.e., no WEC array): where H s,0 is the significant wave height of the undisturbed wave field, that is, the wave field in the absence of WECs. This was found by running SWAN for the Nearshore Domain where BC (2) was the offshore BC. MSE 0 gives the measure of the impact of (or error caused by) the WEC array as assessed by SNL-SWAN. Our definition of MSSS is sometimes referred to as the Brier Skill Score (BSS) in coastal engineering and meteorology [32].
The value of MSSS gives us a quantitative measure of the performance of SNL-SWAN at predicting the wave field as compared to the LWIT baseline: • MSSS = 0 would be produced by not modeling the WEC array at all. • MSSS < 0: SNL-SWAN produces worse results compared to not modeling the WEC array at all. • MSSS = 1: SNL-SWAN exactly matches LWIT results (no error). • 0 < MSSS < 1: SNL-SWAN makes some improvement to the wave field (compared to not modeling the array at all).
As pointed out by Sutherland et al. [32]: "However, as skill scores are used more and more, it will become apparent what value of skill score is needed before a model result can be considered "good"".
To that end, we related the quantitative MSSS to a qualitative interpretation from Sutherland et al. [32], as shown in Table 1. Throughout the text, the terms "Good" or "Excellent" will be used, as capitalized, to refer to the defined ranges given in Table 1.

Results
This section presents the results of the study. First, an example is given to illustrate how results are presented and to address general considerations. Following that, the performance of SNL-SWAN is considered in terms of the MSSS, and, finally, the performance of SNL-SWAN is considered in terms of the dimensional errors in wave height.  ). The predominant incident wave direction is left to right, i.e., from −x to +x. Here, one can see the WEC array, given as small black dots centered at (x, y) = (0, 0). In the LWIT wave field, there are also small white circles near the WECs where the wave field could not be resolved. The flat-bottom domain extends from x = −2 km to x = 1 km. Thereafter, the sloped beach extends from a depth of 60 m (z = −60 m) to a depth of 0 at x = 4 km.

Example Results
In the LWIT wave field, subplot (a), there is some reflection,Ĥ s > 1 in the region directly in front of the array. This is not present in the SNL-SWAN wave field, subplot (b). Both wave fields show a wave shadow,Ĥ s < 1, in the lee of the array; it is narrow near the array but appears to broaden as it approaches the shoreline. Subplot (c) shows the normalized error ( ) between the two wave fields, highlighting regions where SNL-SWAN under-predicts the wave heights (− ) and over-predicts wave heights, + . In addition, indicated are the locations of the transects shown in the next row.
The middle row, subplots (d)-(g), shows longshore transects at different locations in the cross-shore ofĤ s ) for the undisturbed (no WEC array) wave field (black dashed), the LWIT wave field (blue), and the SNL-SWAN wave field (red). They are selected at cross-shore intervals of 1 km, except for the transect closest to shore; at the shoreline, due to wave breaking, there is no wave field, and the cross-shore location of x = 3.8 km was selected as the nearest 100 m interval where waves were present for all results.
Finally, the bottom row, subplot (h), shows the MSSS for longshore transects extending across the domain as a function of the cross-shore location, x. The MSSS calculation begins just behind the location of the array and extends to a distance of x = 3.8 km.

Wave Breaking
As the propagation of waves is modeled over the nearshore bathymetry in SWAN, three main processes impact it: shoaling, refraction, and wave breaking. Wave breaking, in particular, has a dominant impact on the nearshore wave field. One can see, in Figure 7, that wave heights along the longshore direction reduce as the the waves approach the shoreline. In the wave field plots, this looks like a sort of widening of the shadow as the waves approach the shore.

Boundary Condition Effects
As discussed in Section 2.4, only H s (not the entire spectrum) was used as an input to the Nearshore Domain BCs. One can see that this creates a small kink as a boundary condition effect at the boundary between the Array Domain and the Nearshore Domain, which is located at x = 1 km; this is particularly noticeable in the LWIT results.
Furthermore, there is a similar effect on the cross-shore boundaries in Nearshore Domain (BC (3) in Figure 6). In the plot of the 3.8 km transect in Figure 7, the undisturbed wave field line (black dashed) should be flat all the way across, but it shows a slight up-turn near the boundaries. While worth noting, neither of these boundary condition effects made a significant difference to the conclusions drawn from the results.

Garden Sprinkler Effect
Another result is that the SNL-SWAN wave fields show a ripple in the lee of the array, which can be seen most clearly in the transects. The ripple worsens as the waves approach shore. This is due to the so-called Garden Sprinkler Effect, which is often found in spectral models and is caused by discretization of the incident wave directions and interaction with obstacles.
The ripple could likely be reduced by increasing the discretization in wave direction in SWAN. However, this was not practical due to memory allocation limits for these simulations. And, like the BC issue discussed above, it did not have an impact on the conclusions, which is notable as it frequently arises in practical coastal engineering work.

Mean-Squared Skill Score Results
Thirty-six cases were considered: two array configurations, two array spacings, three wave periods, and three directional spreadings.
For each case, the MSSS for the longshore transect across the entire domain at 3.8 km was used as the basis of comparison. Various other regions of assessment for the MSSS were examined, including transects not extending the width of the domain and two-dimensional areas. However, it was found that, while the region did have an impact on the value of the MSSS, this did not change the value categories as given in Table 1 nor the trends in the results.
Of the 36 cases, the range of MSSS scores was from MSSS = 0.30-0.97, equating to performance of SNL-SWAN as "Good" to "Excellent". None of the results were worse than "Good". A simple visual scan revealed that only the cases for T p = 6 s resulted in "Good" rather than "Excellent" fits, as shown in Table 2. Table 2. Number of cases scored as "Excellent" and "Good" in terms of the MSSS from Sutherland et al. [32] for each wave period.

Visualization of Results
The performance of SNL-SWAN for all cases was at least "Good", and, in most cases, "Excellent". This can be further understood by considering plots of the results. Figure 8 shows wave fields for a 2-row array, of spacing = 10d with moderate spreading for the three periods: T p = 6, 8, 12 s; the MSSS at the 3.8 km transect are: 0.43, 0.95, and 0.97, respectively. One can see that LWIT and SNL-SWAN wave fields have visually good matches for each of these runs, particularly for T p = 8, 12 s, where the patterns in each pair of wave fields are of comparable shape and magnitude; this is corroborated by the Error plot, in which the relative error is less than 1% ( < 0.01). At T p = 6 s, one can see that the LWIT wave field shows some wave reflection in front of the array (higher wave heights), which is not captured in SNL-SWAN. Correspondingly, LWIT shows a deeper wave shadow than SNL-SWAN; that is, some of the LWIT wave shadow is the results of wave reflection, not just power absorption, and SNL-SWAN only captures a wave shadow due to power absorption. This is again corroborated by the Error plot, which shows errors of around +/− 1-2%. Figure 9 shows the corresponding longshore transect for the above set of conditions. Again, one can see that SNL-SWAN over-predicts the wave height (i.e., under-predicts the depth of the shadow) as compared to LWIT for T p = 6 s for all transects, while the match between LWIT and SNL-SWAN is very good for T p = 8, 12 s. However, as before, the error between LWIT and SNL-SWAN is only around = 0.01-0.02, or 1-2% of the incident significant wave height. In Figure 9, particularly for T p = 12 s at the 3.0 km and 3.8 km transects, one can also see the impact of the lateral BCs-the undisturbed wave field line becomes curved up at the boundaries; this impact is also manifested in the LWIT and SNL-SWAN wave fields. Figure 10 shows the MSSS for longshore transects as a function of the cross-shore location for the above set of conditions. For T p = 6, 8 s, the MSSS starts very high, nearly 1, and remains high throughout the domain; SNL-SWAN performs very well compared to LWIT throughout the domain.
Interestingly, the MSSS for T p = 6 s starts near 1 but decays to less than 0.5 in the nearshore, resulting in a "Good" rather than "Excellent" score. The reason for this is illustrated in Figure 11: the MSE (i.e., a measure of the error between SNL-SWAN and LWIT) remains constant throughout the cross-shore, but, for MSE 0 , the impact of the array on the wave field decreases towards the shore. Therefore, the relative error increases, i.e., the MSSS worsens.  Mean-Squared Error (MSE) and the mean-squared error of Sandia National Laboratories-Simulating WAves Nearshore (SNL-SWAN) against the no-WEC condition (MSE 0 ) for T p = 6 s and the conditions given in Figure 8.
Considering that, for T p = 6 s, the MSSS got worse towards the shore simply because the impact of the wave farm decreased suggests that one should ask: • How significant are the impacts of the array? • Are the errors between SNL-SWAN and LWIT meaningful in the context of nearshore physical processes?
To answer this, consider simply the difference in meters between the SNL-SWAN and LWIT wave fields: ∆H s = |H s,SNL − H s,LW IT |. ∆H s is considered for the entire Nearshore Domain for each period. Table 3 gives the maximum and mean values of ∆H s for each period. For T p = 6 s, where the MSSS values were lowest, the maximum ∆H s = 0.089 m, or 9 cm. As a dimensional value, this is very low and well within the range of error of met-ocean models typically used to predict nearshore wave conditions (e.g., Reference [33] or Reference [34]).

Conclusions
The performance of SNL-SWAN for modeling WEC array impacts on the nearshore was assessed by comparing it to analogous results from LWIT, LWIT being more accurate at modeling wave-body interactions. The analysis considered medium-sized arrays in four configurations for a range of wave conditions representative of real ocean conditions in terms of T p , H s , and directional spreading.
Quantitative results showed that the agreement between SNL-SWAN and LWIT was "Good" to "Excellent"; MSSS > 0.3 for all conditions, and MSSS > 0.5 for most conditions. Shorter waves, T p = 6 s, performed worse, because wave reflection is not captured in SNL-SWAN, where wave reflection accounts for a significant proportion of the wave shadow for shorter waves. However, this could be addressed by including the impact of reflection in the transmission coefficient curve rather than just wave energy absorption as was done in this study.
Nevertheless, errors between SNL-SWAN and LWIT were small, resulting in only maximum differences in H s between the models of 9 cm for short and medium waves and 20 cm for long waves, which is well within the error of typical met-ocean models used to predict nearshore wave conditions. Further work to investigate other types of WECs, particularly non-axis-symmetric ones, even larger arrays, and the inclusion reflection in SNL-SWAN could be carried out.
Users of SNL-SWAN should be aware that, if large amounts of wave reflection are present, SNL-SWAN, using only the power capture, will over-predict down-wave wave heights, that is, under-predict the shadow. However, this difference may still be relatively small in absolute terms. Where the wave shadow is predominantly due to power absorption, users should have high confidence in the outputs of SNL-SWAN for realistic sea conditions.