Modeling Sea Ice Effects for Wave Energy Resource Assessments

: Wave-generated power has potential as a valuable coastal resource, but the wave climate needs to be mapped for feasibility before wave energy converters are installed. Numerical models are used for wave resource assessments to quantify the amount of available power and its seasonality. Alaska is the U.S. state with the longest coastline and has extensive wave resources, but it is affected by seasonal sea ice that dampens the wave energy and the full extent of this dampening is unknown. To accurately characterize the wave resource in regions that experience seasonal sea ice, coastal wave models must account for these effects. The aim of this study is to determine how the dampening effects of sea ice change wave energy resource assessments in the nearshore. Here, we show that by combining high-resolution sea ice imagery with a sea ice/wave dampening parameterization in an unstructured grid, the Simulating Waves Nearshore (SWAN) model improves wave height predictions and demonstrates the extent to which wave power decreases when sea ice is present. The sea ice parametrization decreases the bias and root mean square errors of wave height comparisons with two wave buoys and predicts a decrease in the wave power of up to 100 kW/m in areas around Prince William Sound, Alaska. The magnitude of the improvement of the model/buoy comparison depends on the coefﬁcients used to parameterize the wave–ice interaction.


Introduction
Wave energy conversion is an emerging source of renewable energy that has the potential to supply one quarter of the total U.S. annual electricity demand [1]. Alaska has the longest coastline of the 50 U.S. states and sizable waves whose energy could be captured with wave energy converters. The southern Alaskan coastline has been identified as a coastline with potential for wave energy power generation due to the size of the waves and their spectral and directional characteristics [2,3]. A 32-year numerical modeling resource characterization study confirmed the wave energy potential of the southern Alaskan coastline [3] but did not assess the resource in areas with known seasonal sea ice. To accurately assess these regions, numerical models must include the dampening effects of sea ice on the waves that may decrease the available power.
The dampening effects of sea ice on waves have been observed in field campaigns in both the Arctic and Antarctic and the observations have been used to parameterize the ice/wave interactions [4][5][6][7][8][9][10]. The parameterization of the dampening is an active area of research with the dependence on ice type and distance from the ice edge being investigated [11]. Large-scale numerical models have begun to add ice/wave parameterizations to examine the relationship between ocean waves and sea ice [7,9,10,12,13]. This study applies an ice/wave parameterization to waves in nearshore ice along the Alaskan coastline to determine how wave energy resource assessments are affected by the presence of sea ice.
Sea ice occurs seasonally along the Alaskan coastline. In early spring, sea ice is found along most of the coastline (Figure 1a), but by early summer, significant amounts only remain in the north (Figure 1b). In fall, there are no areas with ice coverage (Figure 1c) and by mid-December, ice has been returning to most of the coastline (Figure 1d). The coastline along the northern Gulf of Alaska is a region that has both a very energetic sea state [2] and seasonal sea ice. Prince William Sound (PWS), along the northern Gulf of Alaska, has seasonal sea ice and two wave buoys maintained by the National Data Buoy Center (NDBC), which make it a good site to study the effects of sea ice on the wave energy resource (Figure 2a). Data from a previously published wave model at this location [14] show that the modeled wave heights are overestimated when sea ice is present (Figure 2b).  Wave resource assessments that conform to international standards, such as International Electrotechnical Commission (IEC) TS 62600-101, require consideration of the effect of sea ice. One of the most widely used numerical models for wave resource assessment is the Simulating Waves Nearshore (SWAN) model [3,[14][15][16][17][18]; yet, until recently, SWAN has not included the effects of sea ice on waves. With the release of version 41.31, a wave-ice interaction source term has been added to the model. This prompts a reevaluation of model-based wave resource assessments in areas where sea ice and large wave energy resources are present. This paper aims to examine the impacts of sea ice on the local wave field in PWS using the latest SWAN model with the sea ice effect and understand how the wave resource characterization changes when sea ice is considered.

Materials and Methods
The sea ice effects on the wave energy resource were examined with a high-resolution, unstructured-grid numerical model in Alaska, where seasonal sea ice occurs and two wave buoys are available for model validation. Sea ice data were derived from high-resolution satellite images acquired by a passive microwave sensor that measures sea ice regardless of cloud conditions.

Study Area
The study area of PWS is located along the southern coastline of Alaska (red box in Figure 1c). The depths in PWS are less than 1000 m, but outside the sound, they increase to over 4000 m (Figure 3a). Prince William Sound is subject to seasonal sea ice from late fall through spring. The ice forms around the edges of the sound and occasionally forms across the mouth between Hinchinbrook Island and Montague Island (Figure 2a). Two NDBC buoys recorded wave height data that were used for model validation. Buoy 46060 is located inside PWS and buoy 46061 is outside the sound between Hinchinbrook Island and Montague Island (Figures 2a and 3a). Buoy data were available almost continuously from 2016-2020, except from October through mid-December of 2018 at buoy 46061 and February through July of 2016 at buoy 46060.

Wave Model
The SWAN model is a third-generation model that solves the spectral action balance equation [19,20]. In the spectral action balance equation, the wave action, N, balances sources and sinks of energy as:

∂N ∂t
where c i is the group velocity in the ith dimension, θ is the wave direction, and σ is the radian frequency of the waves. The wave action is defined as the ratio of the wave energy, E, and the frequency, σ. The sources and sinks of energy are due to wind, S in , whitecapping, S wc , non-linear quadruplet interactions, S nl4 , dissipation due to bottom friction, S bot , depth-induced wave breaking, S brk , and dissipation due to sea ice, S ice . The sea ice dissipation term is shown in bold because it is a new addition to SWAN in version 41.31, which is used in this study. Version 41.31 includes the effects of sea ice as direct wave dissipation and a reduction of wave growth due to wind when ice is present [21]. It is an empirical/parametric method by which the energy decay is implemented as coefficients affecting different wave frequencies. S ice is defined as: where D ice is the energy decay rate, c g is the wave group velocity, and k i is a 6th order polynomial describing the dissipation due to sea ice: The coefficients C 1 , C 2 , . . . C 6 are constants that depend on the type of ice at the model grid location. Table 1 lists four sets of coefficients that have been derived from field studies in the Arctic and Antarctic. All field-derived sets of coefficients discussed in this paper are from experiments that reported values for C 2 and C 4 and zeros for the other coefficients. The Rogers et al., 2018 [22] coefficients are from a study conducted in the Arctic, where the ice was in the form of frazil or pancake ice. The Meylan et al., 2014 [6] coefficients were derived with data from Antarctica, where the ice was in the form of 10-25 m floes. The Hosekova et al., 2020 [11] coefficients were derived using data from the Arctic, where the waves were measured in nearshore frazil or pancake ice at least 500 m from the ice edge. The determination of these coefficients is an area of active research; thus, evaluating the model sensitivity to the selection of coefficients is desired. To this end, the Hosekova et al., 2020 coefficients were doubled to evaluate the effects of larger coefficients on the wave dampening. These coefficients are discussed in more detail in Section 3.2. Table 1. Coefficients used in SWAN to model sea ice/wave interactions. Hereafter, the sets of coefficients will be referred to by the bold initials. SWAN has been widely used for wave resource assessments because it can be used with an unstructured grid to capture wave dynamics in nearshore areas with complicated coastlines and bathymetry [23][24][25][26][27]. The model mesh developed in this study implements a variable resolution with lower resolution offshore (~2.5 km) and higher resolution near the coast (~300 m) ( Figure 3b). It had 274,929 nodes and 538,146 elements with open boundary nodes spaced 2.5 km apart. The bathymetry used in the grid was obtained from the Southern Alaska Coastal Relief Model [28] and the National Oceanic and Atmospheric Administration (NOAA) Fisheries digital elevation models [29][30][31]. In this SWAN implementation, dissipation due to bottom friction is modeled based on JONSWAP [32] and depth-induced wave breaking is modeled based on Battjes and Janssen [33]. Energy transfers due to triad interactions were also included [34].
The model was run from 2016-2020 with and without the sea ice effect except during the months of July, August, and September, when it was only run without ice because there was no ice present. Seven days preceding each period were run as model spin up and were disregarded in the analysis. The month of December was run for sensitivity tests from 2016-2019, with four sets of coefficients in the sea ice/wave implementation to determine the impact of the coefficients on the wave resource due to ice. December 2020 was not run because the wind data were incorrect during that month.

Boundary Conditions and Sea Ice Input
The SWAN model is forced with wind, ice, and wave data. The wind data are from the Climate Forecast System Reanalysis (CFSv2) [35]. Wave hourly spectral boundary conditions were obtained from a reimplementation of the multi-grid Wavewatch III™ model run by the NOAA National Centers for Environmental Prediction [36]. Sea ice concentration data were downloaded from the University of Hamburg's Integrated Climate Data Center (ftp://ftp-projects.cen.uni-hamburg.de/seaice/AMSR2/3.125km/ accessed on 30 April 2020). The daily sea ice concentration maps have a resolution of 3.125 km and are calculated from the Advanced Microwave Scanning Radiometer 2 (AMSR2) Level-1R brightness temperature measurements [37,38]. The coverage is not continuous along the southern coastline of Alaska but does include PWS (Figures 1 and 2a).

Results
The model was validated with significant wave height data from NDBC buoy 46060 inside PWS and buoy 46061 near the mouth of PWS. The model parameterization of the interaction between sea ice and waves was investigated by testing four sets of coefficients in the sea ice implementation. Maps of wave power with and without the sea ice implementation show the areas around PWS that were most affected by sea ice.

Model Validation
The modeled significant wave height, H s , values agree dwell with the buoy values inside and outside PWS and the agreement improved when sea ice was accounted for in SWAN. Time series of the month of December for 2016-2019 show the agreement during times of large and small waves at both the buoy near the mouth of PWS, which encountered large waves, and the buoy inside PWS, where the waves were smaller (Figure 4). When no sea ice was present, the SWAN without ice output agreed with the SWAN with ice implementation and both models agreed well with the buoy data (14 December 2016, Figure 4a). During storms with sea ice, H s was sometimes overestimated when SWAN was run without ice and H s agreed better with the buoy data when SWAN was run with ice (8 December 2019, Figure 4g,h). Occasionally, the SWAN with ice implementation overly dampened wave heights due to the ice effect, and therefore, underestimated H s (2 December 2017, Figure 4c).
The waves at buoy 46061 near the mouth of PWS were almost twice as large as the waves at buoy 46060 inside PWS, but at both buoys, the larger waves were overestimated when sea ice was not included in the model. This overestimation is shown in Figure 4g Figure 5 and Table 2. The RMSE was improved at buoy 46061 from 0.64 m for the model without ice to 0.57 m for the model with ice (Figure 5a). At buoy 46060, the RMSE did not change when sea ice was added to the model (Figure 5b). The ice implementation lowered the bias at buoy 46061 from 0.44 m to 0.35 m and at buoy 46060 from 0.12 to 0.02 m ( Figure 5). The linear correlation coefficient, R, decreased slightly when ice was added to the model at both locations. The percentage error decreased at both buoys with the ice implementation, and the scatter index decreased at 46061 but increased at 46060. The scatter index improved at buoy 46060 when sea ice was added to the model, with different coefficients describing the wave/ice interaction, as discussed in Section 3.2. Overall, the comparison statistics improved when ice was added to the model.   The effect of sea ice on the waves depends on the amount of sea ice present. Figure 6 shows how the modeled and measured H s values changed with the sea ice concentration. The sea ice concentration was calculated at the location of buoy 46061 near the mouth of PWS. At that location, ice directly dampened waves measured by buoy 46061 and dampened waves that propagated into PWS and were measured by buoy 46060. Modeled versus measured H s scatterplots that were separated by sea ice concentration (SIC) showed changes in the comparisons as the SIC increased from below 10% to above 40%. At buoy 46061, the linear regression was closest to the one-to-one line when the SIC was above 20% and below 30% (Figure 6c). At buoy 46061, the wave heights were overestimated by both the model without ice and with ice until the SIC was above 30%. When the SIC was above 30% the SWAN model with ice underestimated H s at buoy 46061, whereas the SWAN model without ice still overestimated H s (Figure 6d,e). At buoy 46060, the linear regression was above the one-to-one line for the model without ice until the SIC was greater than 40%. The model with ice lowered the linear regression close to the one-to-one line for all sea ice conditions when the SIC was less than 40%. When the SIC was greater than 40%, there were very few measurements at buoy 46060, and adding sea ice to the model did not change the comparison.

Parameterization of Sea Ice/Wave Interaction in SWAN
The sea ice/wave interaction parameterization in SWAN 41.31 uses sets of coefficients to model the dampening of waves by ice. The dependence of the wave heights on these coefficients was explored by running the model with the four sets of coefficients listed in Table 1 Table 3 and Figure 7. The RMSE, bias, R, PE, and SI values all decreased as the coefficients increased. A decrease in all parameters corresponded to an improvement in the model/buoy comparison except for the R value. The R value represented the correlation of the wave height signals in time. The sea ice dampening effect decreased the wave heights, and therefore, decreased the wave height changes in time, which decreased the R value of the comparison. The parameter that improved the most was the bias at buoy 46060. It decreased to zero when the 2H coefficients were used. The RMSE decreased the most at both buoys when using the 2H coefficients.
A longer sensitivity test was carried out by running the model with the four sets of coefficients for the month of December from 2016-2019. The results of this sensitivity test are shown in Figure 8 and Table 4. The improvements of the parameters at buoy 46061 are similar in pattern and magnitude to those from the single-month test. The parameters do not show as much change or improvement at buoy 46060 as they did in the single-month test. The bias and PE both are negative at buoy 46060 when the largest (2H) coefficients are used. This indicates that the model is dampening the waves too much.     The time series of the comparisons with buoy data are shown in Figure 9. When no ice was present, all five model implementations agreed (1 December 2017, Figure 9c). When ice was present, the wave heights decreased in size as the coefficients increased from R to M, H, and 2H coefficients (Figure 9). The agreement with buoy H s values varied when the 2H coefficient implementation agreed best with the buoy and the other implementations overestimated the waves (7 December 2017, Figure 9d) when the 2H coefficient implementation underestimated the waves and the R coefficient implementation agreed well with the buoy (3 December 2017, Figure 9c). Statistical comparisons of H s predicted by the model implementations and measured by the buoy show the model implementation with the 2H ice coefficients predicted waves closest in size to those measured by buoy 46061 near the mouth of PWS. At buoy 46060, inside of PWS, the one-month sensitivity test showed the most improvement in the wave height comparison when the 2H coefficients were used, but the four-month sensitivity test found very little change in the RMSE with any of the coefficients and a negative bias when the 2H coefficients were used. The improvement in the comparison statistics when using the 2H statistics at buoy 46061 and at buoy 46060 in the one-month test demonstrated that the 2H coefficients best describe the wave/ice interaction at these two buoy locations. The 2H coefficients were used for the long model runs whose statistical comparisons to buoy data were reported in Table 2. Future research could improve the understanding of the coefficients by examining the ice type and thickness around PWS.

Effect of Sea Ice on the Wave Resource
Sea ice dissipates waves, and therefore, decreases the available power that could be extracted by a wave energy converter. Figure 10 shows examples of this decrease on 7 December 2017 when the AMSR2 satellite instrument measured ice around the edges and across the mouth of PWS (Figure 2a). The model that was run without sea ice showed wave heights above 1 m inside PWS (Figure 10a). When ice was added to the model simulation, the wave heights decreased significantly at the mouth of and inside PWS (Figure 10b). The decrease is shown as a difference in Figure 10c, where the difference is largest at the mouth of PWS. The decrease in H s due to ice corresponds to a decrease in power. The power from waves is defined in terms of the density of water, ρ, the acceleration due to gravity, g, the wavenumber, k, the water depth, h, and the significant wave height, H s . It is calculated as the mean energy density, E, times the group wave velocity, c g , P = Ec g . The mean energy density is given by linear wave theory as: and the group velocity as: where the phase velocity is: The wave power inside PWS was estimated to be larger when the model was run without ice than when it was run with ice (Figure 10d,e). The power difference is largest at the mouth of PWS and along the coasts of Montague and Hinchinbrook Islands (Figure 10f). At the mouth of PWS and along the coasts of Montague and Hinchinbrook Islands, the power was reduced from 100 kW/m to almost zero by the presence of sea ice. This decrease of wave power by up to 100 kW/m is a significant decrease that should be accounted for when assessing the wave power resource at that location. The wave power was also reduced landward of the mouth of PWS. This is associated with swells being dissipated by sea ice near the mouth because sea ice was only present near the shore and not in the middle of PWS during the study period. This will not only affect the magnitude of the wave power but also the frequency distribution at PWS.
The seasonal changes in available power were investigated by comparing the average monthly power computed from 2016-2020 ( Figure 11). The maximum power difference due to sea ice occurred in December. Significant power differences occurred in January and February, minimal differences occurred in March, April, and May, and no power differences were detected in June, October, and November. The model was not run with the ice implementation for July, August, and September because no sea ice was present around PWS during those months.

Discussion and Conclusions
The sea ice/wave interaction addition to a high-resolution unstructured-grid SWAN model was used to investigate the dampening effects of sea ice on wave heights and wave power in PWS in Alaska. Wave heights and power decreased throughout PWS when sea ice was present, but the largest effects were seen at the mouth of the sound and along the coast of Montague Island. Sea ice was not observed to cover the sound entirely but the waves in the sound decreased when ice across the mouth impeded wave propagation from the ocean. This reduction in wave power was found to occur even on the northern coast of PWS.
Wave power was significantly reduced when sea ice was present. The strongest effect occurred at the mouth of PWS due to the large waves and high concentrations of sea ice. The power decrease changed seasonally with the sea ice and was maximum in December. Future wave resource assessments in Alaska should include the effects of sea ice, especially for the northern parts of Alaska, where significant amounts and time periods of ice cover occur.
The SWAN model was validated using buoy data from one buoy near the mouth and one inside the sound. The model-data comparison improved at both sites when sea ice was added to SWAN, but the improvement was greater at the buoy near the mouth due to strong ice concentrations at that location and larger waves than inside the sound. Four sets of sea ice dampening coefficients were evaluated in SWAN and were found to have varying impacts on the waves. The larger coefficients dampened the waves more than the smaller coefficients and occasionally decreased the wave heights so much that the wave heights were underestimated relative to the buoy measurements. The coefficients are defined by ice type such as pancake or large floes, but no satellite information is available yet that describes ice type. The satellite information input to this implementation of SWAN was a 3.125 km grid of daily ice concentration values, but no information about the type of ice, its thickness, or its change over the 24 h period was available.
The current sea ice implementation in SWAN only considers one set of coefficients for the entire model domain. Future research on the implementation of spatially varying ice type should be considered to improve the prediction of sea ice/wave interactions for locations where the sea ice type changes within the model domain. Large coefficients may be needed close to the coast where the ice is thick and old, but small coefficients may be needed in open water that features freshly formed pancake or frazil ice.  Data Availability Statement: Wave buoy data are available from NDBC at https://www.ndbc.noaa. gov/station_page.php?station=46061 (accessed on 30 April 2020) and https://www.ndbc.noaa.gov/ station_page.php?station=46060 (accessed on 30 April 2020). Sea ice concentration data are available from ftp://ftp-projects.cen.uni-hamburg.de/seaice/AMSR2/3.125km/ and accessed on 30 April 2020. Model output is available by contacting zhaoqing.yang@pnnl.gov.

Acknowledgments:
The authors would like to thank Erick Rogers and Marcel Zijlema for help with the sea ice implementation in SWAN, Harry Stern for his help locating high-resolution sea ice concentration data, and Stefan Kern and Xiangshan Tian-Kunze for help obtaining the AMSR2 data.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Error Statistics
The error metrics used to evaluate model performance are the root mean square error (RMSE), mean percentage error (PE), scatter index (SI), bias, and the linear correlation coefficient (R): where N is the number of measurements, M i is the measured data, P i is the predicted results, and overbars represent time averages.