Background and Motivation

The changing climate is increasingly impacting coastlines around the world (for example, Vitousek and others, 2017; Storlazzi and others, 2018; Barnard and others, 2019; Oppenheimer and others, 2019; Reguero and others, 2019). One of the most affected regions is the Arctic, which is warming at a faster rate than any other region (Serreze and others, 2009; Cohen and others; 2014; Acosta Navarro and others, 2016). In the Arctic, the warming climate leads to an increase in the duration of the open-water season (the period during which the ocean is not covered by ice) and a delay in ice formation in the fall and winter (Serreze and others, 2009), which is typically the stormy season. Studies have found that the extended open-water season correlates with an increase in maximum and average wave heights (Erikson and others, 2020a; Nederhoff and others, 2022) at offshore locations, whereas in the nearshore the total wave energy over the year increases primarily due to an increase in the extent of the open-water season (Erikson and others, 2020a; Gibbs and others, 2021; Nederhoff and others, 2022).

Increases in wave energy will likely cause an increase in erosion, which in Arctic regions is driven by a combination of wave activity and thermal processes (Erikson and others, 2011; Stopa and others, 2016; Denali Commission, 2019; Gibbs and others, 2019a; Casas‐Prat and Wang, 2020; Erikson and others, 2020a; Gibbs and others, 2021). The American and Canadian coastlines exhibit the highest erosion rates in the Arctic (Lantuit and others, 2012). The United States parts of the Beaufort Sea and Chukchi Sea coastlines have been primarily erosional since the 1940s (Gibbs and Richmond, 2017; Gibbs and others, 2019b). At some locations, coastal erosion rates have more than doubled since the mid-20th century (Frederick and others, 2016; Jones and others, 2020; Gibbs and others, 2021).

To be able to assess past and anticipate future changes to the Alaskan coast, reliable nearshore wave data are needed. However, wave observations are sparse and difficult to obtain due to the harsh environmental conditions and remoteness of the region. Although reanalysis data are now available for the satellite era (1979 to present), the resolution is still too coarse (0.25 degrees for atmospheric and 0.5 degrees for oceanic processes) to capture nearshore processes such as refraction, diffraction, dissipation, and non-linear energy transfer processes, which are important to account for when considering erosion and inundation of coastal areas.

Here, we use a numerical wave model to downscale available reanalysis wave data for the Alaskan region from the United States-Canada border to the Bering Sea (fig. 1). The generated nearshore time series allows for a detailed analysis of changes in the wave climate along the entire coast. We caution, however, that this dataset is based on outdated bathymetry (see the “Bathymetry” section) and should be updated when modern bathymetric nearshore data are made available. The lack of modern nearshore bathymetric data and its importance in modeling coastal hazards along Alaska’s shorelines is well recognized (Denali Commission, 2019; Hamilton, 2021; Williams and Erikson, 2021). Extensive efforts are underway to update the well-outdated (greater than [>] 50 year) nearshore bathymetric data along the roughly 4,000-kilometer (km) stretch of coast that this study addresses, but the data are not anticipated to be available for another 5 to 10 years (Alaska Mapping Executive Committee, 2020). Simultaneously, coastal hazards are of immediate and increasing concern to vulnerable coastal ecosystems and villages, necessitating an understanding of both past and potential future coastal storms. To that end, the scientific community-at-large is at a temporary impasse of either postponing research or moving forward to advance the quantification of nearshore wave conditions, but with the latter action subject to uncertainties in accuracy. This study by the U.S. Geological Survey (USGS) was structured to address these acute needs but with the recognition that uncertainties in the results exist and that updates should be made when modern coastal bathymetry data are obtained.

Meteorological and Oceanic Reanalysis Forcing Data

ERA5 (Hersbach and others, 2020) is a detailed reanalysis of the global atmosphere, land surface, and ocean waves from 1950 onwards produced by the ECMWF. This dataset provides, among other variables, estimates of atmospheric parameters such as air temperature, pressure, and wind on a 0.25-degree grid and sea ice concentrations and information on waves over the global oceans at a 0.5-degree resolution. The reanalysis combines model data with observations from across the world into a globally complete and consistent dataset. ERA5 has been shown to perform well in capturing observed weather and climate variability in Alaska and the Arctic (Graham and others, 2019). Here, offshore significant wave height (H_s), mean period (T_m), and direction (D_m) are used as input to the SWAN model. Wind conditions (u₁₀, v₁₀) from this reanalysis dataset are additionally applied across all model domains. ERA5 does not provide wave conditions if the sea ice cover is greater than 30 percent, therefore the reconstructed time series assumes no waves inshore if there are no waves in ERA5. Additionally, the 0.25-degree sea ice resolution does not resolve landfast ice, which could potentially cause an overestimate of the wave climate (Hošeková and others, 2021).

Nearshore Wave Model

Model hindcasts were run with the third generation, physics-based, phase-averaging wind-wave model SWAN (version 3.07.00.63652). SWAN is based on the wave action balance and includes physical processes such as wind-wave generation, shoaling, breaking, and nonlinear wave-wave interactions (see Booij and others, 1999, for details). Simulations were carried out in a stationary mode for 2,500 sea states (see “Boundary Conditions and Downscaled Wave Database” section of this report) over eight two-dimensional (2D) curvilinear and one 2D rectilinear domains for a total of nine domains (fig. 1). The rectilinear domain was chosen for the Kotzebue Sound area to ensure numerical stability. The offshore extent for the model domains was roughly defined by the 20-m isobath, and the domain length and width varied, depending on the offshore extent of each grid and local shoreline curvature. Selection of the approximately 20-m isobath as the offshore extent was based on (1) the need to reduce computation cost by optimizing the spatial extent in the cross-shore direction while accounting for the shallow and wide continental shelf that abuts the north and west Alaska coastline, and (2) comparisons between ERA5 reanalysis waves and five temporary buoy deployments along the Beaufort Sea part of the Alaska Coast which showed an acceptable accuracy (average root mean square error [RMSE] and bias =0.35 m and 0.16 m respectively; Kasper and others, 2023) in 22- to 30-m water depth, supporting the use of these data as open boundary conditions for the nearshore model. Nearshore grid size resolution was on average less than or equal to [≤] 200×200 m. In the offshore regions, cross-shore grid resolution varied between about 300–1,000 m and 200–300 m in the alongshore direction. Since grid resolution was coarse offshore, due to the large simulation domain, SWAN’s “obstacle” option was used to ensure that barrier islands were properly resolved in areas with larger cell sizes. Obstacles, implemented either as option “dam” or “sheet” with a polyline, can be located between grid cells and interrupt the wave propagation (Deltares, 2020). In this study, dams were used to represent barrier islands, with transmission coefficients depending on (1) the wave conditions at the obstacle and (2) obstacle configuration. Empirical coefficients $\alpha$ and $\beta$ that describe obstacle shape, were kept at default values of 2.6 and 0.15, respectively. Wave reflection off the obstacle was set to zero, and the dam height was set to 2 m, on the basis of 2010–2011 light detection and ranging (lidar) elevation measurements along Alaska’s north coast (Gibbs and Richmond, 2017; Hamilton and others, 2021).

SWAN was run with the white capping formulation of van der Westhuysen and others (2007) using 72 directional bins (5-degree resolution) and 46 frequency bands from 0.03 to 2.5 hertz (Hz). The friction formulation of Collins (1972) with a coefficient of 0.02 was used after calibration by Nederhoff and others (2022) for the Alaskan coast. Wave conditions at the boundaries and wind forcing, applied homogenously over each domain, were determined by the ERA5 sea states (see the “Boundary Conditions and Downscaled Wave Database” section in this report). The presence of sea ice was not accounted for in these SWAN model runs. Whereas landfast sea ice conditions are not explicitly accounted for in these runs, the database could be expanded to include nearshore sea ice by running the same conditions with different levels of sea ice cover.

Bathymetry

The bathymetry for the domains was primarily based on products provided by the International Bathymetric Chart of the Arctic Ocean (IBCAO; Jakobsson and others, 2020), which have a resolution of 200×200 m. IBCAO is based on all available bathymetric datasets north of 64° N. latitude. In the nearshore region, these are mostly National Oceanic and Atmospheric Administration (NOAA), National Ocean Service (NOS) hydrographic surveys, which were conducted between the 1940s and 1970s. It is assumed that the vertical IBCAO datum is the chart datum of the NOS samples (Mean Lower Low Water, MLLW), which were converted to mean sea level (MSL) by using datum information from NOAA’s Prudhoe Bay AK station (https://tidesandcurrents.noaa.gov/datums.html?id=9497645). Changes to the bathymetry were made in the vicinity of Cross Island, where a NOS dataset (H07760 from 1950) appeared to have been horizontally misplaced and which was manually rectified. Sea level rise was considered as an increase in water levels of 2 millimeters per year (mm/yr) in the region (Sultan and others, 2011) between the year of the surveys and today (roughly 0.15 m), which was added to the bathymetry. The land boundary is based on data extracted from USGS topographic maps from the 1950s to 1990s at scales between 1:63,360 and 1:250,000.

Inconsistencies between the land boundary and IBCAO products were adjusted by removing bathymetric points along the boundary in a roughly 200- to 400-m-wide section, setting the immediate coastline to zero and using a smoothing process (internal diffusion, see Deltares, 2021) to interpolate between the coastline and nearest bathymetric depths. Some uncertainty in bathymetry and land boundaries exists. In particular, many of the NOS hydrographic surveys in the nearshore region were conducted between the 1940s and 1970s and likely differ in places from current and recent past bathymetry over the simulation time period from 1979 to 2019. Erosion of the coastline and barrier island migration have led to substantial change in the nearshore bathymetry since data were last collected (for example, Bristol and others, 2021; Gibbs and others, 2021; Hamilton and others, 2021). Erosion is especially severe along the United States part of the Beaufort Sea coast, where, for example, shoreline erosion rates measured in the region around Drew Point are as high as 17.2 m/yr in the years between 2007 and 2016 (Jones and others, 2018). In addition, nearshore surveys in the southern Chukchi region are generally sparse (https://www.ncei.noaa.gov/maps/bathymetry/). Overall, uncertainties related to bathymetry and shoreline position are substantial and cannot be specified owing to a lack of knowledge and local differences in erosion and accretion rate. Therefore, the data have been manually inspected and improved where needed and accepted as is for the remainder of this study.

Boundary Conditions and Downscaled Wave Database

Carrying out model computations over the period 1979–2019 along roughly 4,000 km of coastline would be computationally intensive and require considerable time to complete; in other words, it would be very computationally expensive. To improve this, the model was forced with a reduced set of likely combinations of wind and wave parameters, hereafter termed “sea states,” which are used as boundary conditions to the SWAN model domains. These parameters were derived from ERA5 reanalysis time series during the open season (when sea ice cover is less than 30 percent) covering the period 1979 to 2019. For each grid, appropriate ERA5 output locations were chosen along the length of the grid (which ranged from approximately [~] 150 to 280 km in the alongshore direction), and all ERA5 time series were extracted. The number of ERA5 locations varied with grid length and coastline alignment. Those time series were combined for each grid, and representative sea states were established with a multivariant maximum-dissimilarity algorithm (MDA), which determined representative combinations of significant wave height (H_s), mean wave period (T_m), mean wave direction (D_m), wind speed (W_spd), and wind directions (W_dir), resulting in 2,500 different combinations. The MDA method allows for a full representation of the marine climate because the determined sea states are uniformly distributed across all the data (fig. 3), including extreme events (Camus and others, 2011).

SWAN was forced with all 2,500 wind and sea states to create a DWDB at each grid point. The wind was assumed to be homogenous for each sea state over the model domain. To construct nearshore hindcast time series, the 40-year record of ERA5 hindcasts along the boundary was first matched with sea states. For each timestep, the closest combination, as described below, of significant wave height, mean wave period, mean wave direction, wind speeds, and wind directions to ERA5 was found. The algorithm started with a combination of small differences (H_s less than or equal to [≤] 0.05 m, T_m≤0.5 second [s], D_m≤5 degrees, wind speed <2 meters per second [m/s] and wind direction <5 degrees), and if no sea state could be found, a scheme was deployed that allowed for a gradual widening of the bins up to a difference in H_s of 1 m, in T_m of 2 s, in D_m of 20 degrees, in wind speed of 3 m/s, and in wind direction of 20 degrees. Values for which the difference in H_s is greater than 0.15 m are flagged to inform the user of this difference. Flags are set for differences of 0.20 m, 0.25 m, 0.50 m, 0.75 m, and 1 m. While the MDA method allows for the inclusion of extreme events, values for some timesteps could not be found (less than 0.6 percent of the open season values and mostly coming from about 100 degrees), owing to the limit of 2,500 wind and sea states that were chosen to save computational time. Most missing values were caused by an underrepresentation of the wave direction or wind speed (particularly by wind speeds smaller than 4 m/s) in the sea states. Out of the 52 ERA5 locations that were used for the reconstruction, six locations are missing one to three of the highest 10 wave heights. For all values that could not be found, as well as for times when no ERA5 wave data were available (during the ice-covered season) values were set to −9999. Time steps for which the differences between ERA5 wind or wave directions and the found sea state was equal or larger than 15 degrees were flagged (Flag D) in the time series files (available from Engelstad and others, 2024). Wind and wave directions can considerably affect the wave evolution, especially when varying in the alongshore direction between onshore and offshore, so the user is advised to use caution under these circumstances. Although not all conditions are a perfect match at the boundary (fig. 4) owing to the constraints of using a limited set of sea states instead of running the model brute force, extremes for all parameters are captured well and the difference between the original ERA5 time series and the reconstructed values at the boundary are minimal.

For each of the sea states found at the boundary, the corresponding values were extracted from the DWDB at nearshore locations. Those nearshore locations were determined from shore-normal transects derived by the USGS for shoreline change studies (Gibbs and Richmond, 2017) using the Digital Shoreline Analysis System (DSAS) extension for ArcGIS (Thieler and others, 2017). These transects cover the open-ocean and barrier coast as well as the lagoon coast. Approximately every 8th transect of the roughly 50-m alongshore spaced transects was chosen for an alongshore resolution of roughly 400 m. Data were extracted from the DWDB along the 5-m and 10-m isobath closest to each selected transect (fig. 5). For the reconstruction of the time series, the time series of the closest ERA5 locations to each nearshore location was used. Since ERA5 locations are roughly 55 km apart, wave and wind time series between ERA5 locations can differ substantially. This, in turn, can cause inconsistencies between nearshore time series for adjacent output locations that utilize time series from two different ERA5 locations (see fig. 5). To make this apparent to users, adjacent nearshore locations using different ERA5 locations are marked by a flag, Flag S, in the time series files. In addition, Flag S is also used when adjacent nearshore locations were located in separate numerical domains, which can cause inconsistencies even with a wide grid overlap (dual transect locations were removed roughly midpoint of the grid overlap).

Although the method of reconstructing the time series was described above for the 40-year hindcast, data can be extracted from the DWDB in a short amount of time at every location on the grid for any period within the DWDB time frame.

Foggy Island Bay

Observations in the summers of 2019 and 2020 (table 1) in the vicinity of Foggy Island Bay were collected by the University of Alaska Fairbanks, with Sofar Spotter devices (Raghukumar and others, 2019). These buoys measure 3D surface displacements at 2.5 Hz. Bulk wave statistics, such as significant wave heights, mean and peak wave periods, and mean and peak wave directions, are computed from the surface elevation variance density spectrum every 30 minutes in the frequency range from 0.0293 to 0.6543 Hz. Note that in 2020 the Foggy Island Bay Spotter buoy (SPOT 0519) was displaced during the measurement campaign, probably by an ice float (fig. 6, table 1). Wave measurements at Wainwright were obtained by the USGS with a Nortek 1-MHz acoustic wave and current profiler (AWAC), sampling during 8.53-minute bursts at a rate of 2 Hz, (fig. 6, table 1) from the end of August to the beginning of October 2009. Wave measurements at Arey Island were collected in August of 2011 with a Nortek 1-MHz AWAC, sampling at 1 Hz in 17-minute-long hourly bursts.

The comparison between SWAN results and observations in Foggy Island Bay in 2019 shows good agreement (fig. 7) for H_s, T_m, and D_m, suggesting that the model can reproduce the wave field in the area. Skill scores for wave heights are 0.12 m and 0.09 m for the RMSE and MAE, respectively (table 2), and the SCI is 24 percent. Wave periods have a RMSE of 0.8 seconds (s), a MAE of 0.5 s, and a SCI of 23 percent. Further, SWAN runs with the DWDB method have similar reproductive skill as the more elaborate and time-consuming brute-force method (fig. 7). Skill score differences for H_s between the two methods are 0.01 m for the RMSE, 0.01 m for MAE and a bias of 0.03 m, and the SCI is the same. For T_m, the difference between the two methods is 0.1 s for the RMSE, 3 percent for the SCI and a −0.1 s bias, and the MAE is the same. The largest discrepancies in H_s between model(s) and observations can be seen between August 31 and September 2, 2019. Although it is possible that SWAN was not able to reproduce the onset of the sudden large wave events, for example, owing to the stationary mode with hourly wind forcing, the strong increase in measured wave heights over short periods of time (for example, from H_s=0.85 m on September 2, 20:30 hour to H_s=2.34 m on September 2, 21:00 hour) suggests that the difference could also have been caused by measurement errors. Also note that the ERA5 wind speed used to force the model was higher by roughly 4 to 5 m/s during this period (fig. 7D). Smaller mismatches between modeled and observed H_s, such as on August 26, September 4, and September 8, are mostly caused by differences between ERA5 winds (which, as mentioned earlier, were forced homogenously over the SWAN domain) and local winds, as supported by a model run using the local wind forcing (not shown).

The wind also seems to affect the comparison of modeled and observed H_s at the SPOT 0518 and SPOT 0519 locations in 2020 (figs. 8, 9). Discrepancies between modeled and observed H_s can be seen for times of disagreements between ERA5 and local winds (see, for example, July 25, July 29–August 1, and August 7 in figs. 8, 9). On the other hand, the largest event with wave heights up to 1 m on August 2, for which ERA5 and local winds were similar, was captured well by the DWDB method, but it was underpredicted by the brute-force model run for SPOT 0518. Although the discrepancy in wind speeds was largest (~5 m/s) around July 27, the wind direction was from the southwest and did not affect wind-wave growth as much because the fetch (the area where the wind blows to generate waves) was limited. In general, modeled and observed T_m agreed well (figs. 8B, 9B), aside from some single higher than modeled periods, as between July 30–31 when T_m for SPOT 0518 was as high as 12 s, whereas modeled T_m was around 2 s. This was not observed for SPOT 0519. Wave directions were mostly aligned with the wind directions, and differences between modeled and observed D_m (figs. 8C, 9C) are apparent when local and ERA5 wind fields differ.

Wainwright

Observations are also available (Erikson and others, 2022) for a location (F1) in the Chukchi Sea fronting the community of Wainwright (table 1, fig. 6). Although for this location, small and medium H_s (~1 m) are generally captured well by the model (fig. 10), larger events in the period from September 2–8, 2009, were underestimated, with the largest difference on September 4 of about 0.9 m. Those differences are likely generated by variations in observed water levels (water levels were kept constant in the SWAN runs) owing to the wind direction (fig. 11), which are not captured by the model because wind speed and directions implemented in the model agree well with local wind observations (National Centers for Environmental Information, n.d.) at Wainwright airport at this time (fig. 11C, D). Generally, wind from the east tends to decrease water levels, whereas winds from the west tend to increase water levels along the North Slope (Reimnitz and Maurer, 1979; Sultan and others, 2011; Erikson and others, 2020a). This is caused by the Coriolis force, which deflects wind- and wave-forced currents to the right of the flow direction. Additionally, strong winds can push water masses downwind, which at the Wainwright location could also potentially increase water levels at times when the wind is from the west and decrease water levels when the wind is from the east. On the other hand, wave heights were overestimated on September 22 by about 0.5 m, which could have been caused by a lower than modeled wind speed and (or) by the lower water level owing to the easterly winds (fig. 11B, C). Water level fluctuations coinciding with a mismatch between model and observations are not noticed for the other locations in the vicinity of Foggy Island Bay, which are sheltered by barrier islands, whereas the Wainwright location is fully exposed to the open ocean but close to shore.

Arey Island

For all previously presented model-data comparisons, the agreement between observations was good or could be explained by differences between ERA5 offshore winds and locally observed winds. In contrast, modeled wave heights offshore of Arey Island at BI01 are consistently biased high by about 0.20 m (fig. 12A) even though DWDB, ERA5, and local winds compare well (fig. 12D, E). The low wave heights (0.50 m) observed during periods of high westerly and easterly sustained winds ≥12 m/s suggest that the measured H_s may be biased low, rather than the model estimating too high. In fact, once a uniform bias correction of 0.19 m is applied, the RMSE reduces from 0.26 to 0.17 m for H_s (table 2). Satellite data indicate that sea ice was more than 300 km offshore, yielding a fetch area larger than 90,000 square kilometers (km²) that is only marginally limiting to wave growth. Simple calculations of wave heights in approximately 5-m water depth using linear wave theory under conditions of 12-m/s winds and an assumed ice-free region indicate that wave heights should be on the order of 1 m, supporting the SWAN model results. Because the SWAN model and the DWDB method provided otherwise good results, this suggests that the difference could be caused by local factors such as changes in bathymetry that are not captured by the model for the particular simulation time.

The fact that the agreement between observations and model is better for some locations than others is also reflected in the model skills, where the largest spread away from the line of perfect agreement is found for the Arey Island comparison (fig. 13, table 2). Although the model-observation comparison is not perfect, the combined RMSE (table 2) for H_s is relatively small (0.17m) and mostly unbiased (0.09 m). The model-observation comparison shows a larger spread in the data for the combined T_m (SCI of 44 compared to 38 for H_s) for which the bias is small (−0.2 s), although the combined RMSE is 1.7 s. (For statistics regarding the model skill at single locations, see table 2).

Model-Observation Comparison for Historical Wave Observations

The historical observations from 1985 (Envirosphere, 1991) are separated here from the more recent observations introduced above and will not be subjected to a skill performance analysis. The reason is twofold. First, wave height measurements were taken once a day with a wave staff (comparable to a yard stick), whereas the wave direction was determined by measuring the travel direction with a hand-held compass and therefore measurements lack precision. Second, owing to erosional and (or) depositional changes in the area and the construction of Endicott Island in 1987, locations relative to barrier islands and the coastline have changed. However, owing to the sparsity of measurements on the North Slope, these observations are valuable, and we use them to confirm that model results are in the expected value range.

We used 4 out of a total of the 27 measurement locations (fig. 14) for which observed and modeled water depths were similar. Numbers in the location name refer to all historical observation locations and are led by the letter “F,” an abbreviation for the original name “fyke” (meaning a bag net for catching fish). Historical wave heights were measured in 0.1-m increments, whereas wave directions were recorded in 10-degree increments. Wave directions were converted to true north, assuming a magnetic declination of 30 degrees. Although observations are available from July to mid-September, the length of the comparison record is restricted by offshore ERA5 wave conditions, which were only available after September 3 due to offshore ice cover. Each of the daily observed values was compared to the average daily modeled value as well as the minimum and maximum daily values (fig. 15).

Overall, the model-observation comparison shows fair results (fig. 15) considering the measurement methods, uncertainties in bathymetry, and local changes in barrier island locations. While the comparison certainly shows some discrepancies between observations and model (such as for wave heights at F24 and wave directions at F27, see fig. 15), differences are still within an order of magnitude. F18, located furthest away from the barrier islands, shows the best agreement between observations and model. Additionally, the highest observed wave heights (0.6 m), measured at F22, are represented by the model.