Numerical Analysis of Storm Surges on Canada’s Western Arctic Coastline

Abstract

A numerical study was conducted to characterize the probability and intensity of storm surge hazards in Canada’s western Arctic. The utility of the European Centre for Medium-Range Weather Forecasts Reanalysis 5th Generation (ERA5) dataset to force numerical simulations of storm surges was explored. Fifty historical storm surge events that were captured on a tide gauge near Tuktoyaktuk, Northwest Territories, were simulated using a two-dimensional (depth-averaged) hydrodynamic model accounting for the influence of sea ice on air-sea momentum transfer. The extent of sea ice and the duration of the ice season has been reducing in the Arctic region, which may contribute to increasing risk from storm surge-driven hazards. Comparisons between winter storm events under present-day ice concentrations and future open-water scenarios revealed that the decline in ice cover has potential to result in storm surges that are up to three times higher. The numerical model was also used to hindcast a significant surge event that was not recorded by the tide gauge, but for which driftwood lines along the coast provided insights to the high-water marks. Compared to measurements at proximate meteorological stations, the ERA5 reanalysis dataset provided reasonable estimates of atmospheric pressure but did not accurately capture peak wind speeds during storm surge events. By adjusting the wind drag coefficients to compensate, reasonably accurate predictions of storm surges were attained for most of the simulated events. The extreme value probability distributions (i.e., return periods and values) of the storm surges were significantly altered when events absent from the tide gauge record were included in the frequency analysis, demonstrating the value of non-conventional data sources, such as driftwood line surveys, in supporting coastal hazard assessments in remote regions.

1. Introduction

Extreme weather events along Canada’s Arctic seaboard can generate storm surges, defined as changes in water levels associated with atmospheric pressure and wind effects, which have led to flooding of coastal communities. Numerical modelling of storm surges in the Arctic Ocean has generally been scarce due to the lack of measurement data and the unique challenges that arise when modelling wave–ice interactions [1]. Additionally, the Arctic is greatly affected by climate change impacts [2]. The typical open-water season of June to October along the north coast of Canada is predicted to lengthen to a season spanning April to December by 2100 [3]. Ignoring potential climate change influences on storminess, which are highly uncertain [1], winter storm surges that would in the past have been attenuated by sea ice cover have the potential to become more severe in the future. Other projected climate change effects in the Beaufort Sea include relative sea-level rise, thawing permafrost, increasingly extreme wave climate, and accelerating coastal erosion [4], which will compound declining sea ice cover and hence exacerbate the existing hazards and vulnerabilities associated with extreme storm surges. In addition to the changes in ice conditions, projections indicate that relative sea level may rise by around 1 m by 2100 (95% confidence limit relative to 1986–2005) under the RCP8.5 scenario for the Tuktoyaktuk region [5]. The effects of relative sea-level rise on coastal hazards are intensified by the increasing permafrost temperatures [6], which can contribute to erosion and coastline retreat [7]. These climate change impacts are expected to increase the exposure of coastal communities and maritime infrastructure to storm surge hazards, exacerbating existing vulnerabilities, and causing a pressing need to study and understand the development of storm surges in the Arctic [6].

During the 1970s and 1980s, a series of studies were initiated by the Beaufort Sea Project focused on assessing storm surge hazards in the Canadian Arctic. These early studies characterized the weather patterns in the Beaufort Sea and the development of numerical storm surge models forced by surface winds calculated from estimated geostrophic winds [8,9,10,11], with spatial resolutions of roughly 20 km. During that time period, validation of storm surge models was not possible due to the lack of continuous, multi-year tide gauge records. Kowalik [12] simulated three storm surge events in the Beaufort and Chukchi Seas that occurred during ice-free periods using interpolated surface winds from 6-hourly calculations of geostrophic winds.

High-resolution climatic forcing data have since become available in the form of global reanalysis datasets, which combine models with observations. Of particular interest is the ERA5 reanalysis [13], recently released by the European Centre for Medium-Range Weather Forecasts (ECMWF), and superseding ERA-Interim [14]. ERA5 promises a higher spatial and temporal resolution and improved storm physics compared to ERA-Interim and other global reanalyses, with direct implications for application to storm surge modelling. While there have been some studies that assess the performance of ERA5 in approximating surface winds and pressures [15,16,17], these studies did not investigate, in detail, ERA5’s performance in the Arctic region. The spatial coverage of weather station observations to support the development of the reanalysis is limited in the Beaufort Sea region. Thus, the quality of ERA5 as the forcing input for Arctic storm surge models is still unknown.

Field measurements or observations of storm surges are needed to enable development, calibration, and validation of predictive numerical models. For the Arctic region, there are few reliable, long-term, continuous tide gauge records, reflecting the sparsity of human settlement and the challenges for deployment, maintenance, and survival of instrumentation in this harsh environment. In the absence of conventional water-level records, evidence inferred from seasonal field survey campaigns provides an alternative means to gain insight into historical storm surges. In 1977, Reimnitz and Maurer [18] conducted field investigations along the coast of northern Alaska, USA, and their observations provided insight into a significant storm surge event that occurred on 13 September 1970. Reimnitz and Maurer [18] deduced sea-level elevations associated with the 1970 event from a field survey of driftwood debris lines deposited by the storm. Similarly, Harper et al. [19] surveyed driftwood debris lines at Tuktoyaktuk, on Canada’s Beaufort Sea coast, and inferred estimates of the maximum storm surge elevations for the area. In addition, the survey produced detailed storm surge elevation estimates for an event that occurred on 23 August 1986. Records from a tide gauge installed at Tuktoyaktuk in 1961, though affected by several gaps and periods when the gauge was not in operation, suggest that storm surge elevations have not since exceeded the values inferred by Harper et al. [19].

Accurate simulation of storm surges in Arctic regions can be challenging owing to the presence and dynamics of sea ice, which can amplify or attenuate storm surges [20]. Kowalik’s [12] model incorporated coupled interactions between ice and water. However, the model’s performance in capturing the influence of ice on storm surges was not tested, as only events coinciding with ice-free conditions were simulated. Danard et al. [21] simulated eight storm surge events in the Beaufort Sea region that occurred during the months of August, September, and October in the 1970s and 1980s. These events coincided with low ice concentrations and an absence of fast ice. The study showed some agreement between the measured and modelled storm surges, but the timing of the modelled peak storm surge water level was delayed by up to 14 h. The authors noted that the large grid size (18 km) and a lack of sufficient spin-up time were limitations to their study. Additionally, the work by Danard et al. [21] was not verified for winter events with high ice concentrations. Joyce et al. [22] demonstrated the performance of a new set of formulas that improve on the work of Chapman et al. [23] to account for the influence of ice cover on storm surges on the western Alaskan coast. Four simulated storm surge events yielded root-mean-square errors in total water levels of less than 20 cm.

In this study, the utility of the ERA5 reanalysis dataset as atmospheric forcing input for storm surge modelling in Canada’s western Arctic region was examined. A two-dimensional (depth-averaged) hydrodynamic model of the region was developed using the TELEMAC-2D (v7p3r1) shallow water (Saint-Venant) equations solver [24]. Customized subroutines were employed to incorporate the impacts of sea ice on storm surges, using the formulation of Joyce et al. [22]. Model performance was evaluated by simulating fifty historical storm surge events that spanned a range of sea ice conditions, ranging from full ice cover to open-water conditions. Subsequently, a hindcast of the significant storm surge event that occurred on 23 August 1986 was conducted, and the results were compared to peak storm surges inferred from driftwood lines surveyed by Harper et al. [19] more than 30 years ago. The comparison provided insight into the veracity and potential utility of storm surge records derived from non-conventional observation methods. To the knowledge of the authors, this is the first attempt to conduct a large-scale numerical investigation of storm surge hazards in the Beaufort Sea since the work of Danard et al. [21]. Significant advances in atmospheric datasets and computational hydrodynamic models have occurred in the three decades since, which are leveraged in this study to provide an improved understanding of the probability distribution of extreme storm surges in the Beaufort Sea.

2. Methodology

An overview of the study methodology is shown in Figure 1. Each step is described in further detail in the proceeding sections—the study area and computational domain is defined in Section 3; the analysis of tide gauge records is presented in Section 4; the model setup is described in Section 5; and the model calibration, validation, and simulation results are shown in Section 6.

3. Study Area

The hydrodynamic model’s computational domain (Figure 2) was selected to encompass the approximate historical maximum extent of open water (i.e., coinciding with the minimum sea ice cover). As the storm surge simulations performed in this study span a range of historical ice cover extents, this methodology ensures sufficient fetch to capture the development of storm surges due to wind shear stresses and barometric pressure fluctuations associated with the passage of storm systems through the Beaufort Sea. Monthly ice records for the Beaufort Sea during the period of 1979 to 2019 from the U.S. National Snow and Ice Data Centre (NSIDC) [25] were used to determine the minimum ice extent that occurred during that time period. The computational domain spans roughly 1000 km in the east–west direction and 800 km in the north–south direction, and touches multiple coastal communities adjacent to the Beaufort Sea and Amundsen Gulf. The domain extends to offshore water depths of more than 3700 m.

4. Tide Gauge Analysis

Water-level measurements were available for several tide gauges in the region operated by the Canadian Hydrographic Service (CHS) of Fisheries and Oceans Canada, with some records starting in the 1960s. Many of the tide gauge records are associated with temporary monitoring stations, and there are few continuous, long-term (i.e., multi-decadal) water-level records at locations in the Beaufort Sea. Of the 42 tide gauges in the numerical model domain, the station at Tuktoyaktuk (CHS gauge 6485) is the only permanent tide gauge and provides the longest-standing and most complete record of water levels. The Tuktoyaktuk gauge was therefore used as the primary source of water-level data for model calibration and validation. The water-level records are referenced vertically to local chart datum (CD). Based on the tide gauge records between August 2008 and February 2018, CD is approximately 0.454 m below the mean sea level (MSL) at the Tuktoyaktuk gauge for the mean epoch of December 2010. All total water levels and elevations reported in this paper are referenced vertically to local CD at Tuktoyaktuk.

Following the procedures outlined in Foreman et al. [27], astronomical tidal constituents for the Tuktoyaktuk station provided by CHS were used to develop a continuous time series of tide predictions for the period 1961–2019. The predicted tidal elevations were subtracted from the total water-level measurements for the Tuktoyaktuk gauge, to obtain a time series of water-level residuals. The residuals were assumed to be representative of storm surges arising from wind and atmospheric pressure set-up or set-down. This assumption may be a source of uncertainty, in case the water level records are contaminated by wave effects, which are projected to become more frequent and intense in the Arctic region [28,29]. The water-level residual record was linearly de-trended to remove the long-term trend associated with relative sea-level rise.

Figure 3a,b show the measured water levels and the predicted tidal elevations, while Figure 3c shows the residuals calculated from the de-trended sea-level elevations.

Identification of Storm Surge Events and Extreme Value Analysis

A peaks-over threshold (POT) approach [30] was applied to identify extreme storm surge events. The analysis was conducted using the Wave Analysis and Fatigue and Oceanography (WAFO) toolbox [31]. Extreme values occurring within 2 days of each other were counted as a single event to ensure event independence, a prerequisite for the extreme value analysis. An initial threshold based on the average residual plus two standard deviations was used for preliminary screening and identification of extreme events. A final threshold value of 0.96 m was selected using the ‘mean residual life’ plotting method [32]. A Generalized Pareto Distribution (GPD) was fit to the identified peak water level residuals to estimate return periods associated with storm surge events. Figure 4 shows the GPD fit to the data with 95% confidence intervals.

Although the Tuktoyaktuk gauge is the most complete tide gauge in the studied area, there are long periods where the gauge was not operational and was not able to capture large surge events (Figure 3). In particular, based on interviews with the local residents of Tuktoyaktuk, two separate storm surge events on 1 September 1944 and 14 September 1970 were not captured. These events reportedly led to water levels of approximately 3 m (MSL) [33]. These estimates are likely subject to some uncertainty, as they were based solely on the witness reports and recollections of residents. Harper et al. [19] inferred peak water levels closer to 2.5 m MSL (2.95 m CD) for both events based on field surveys of driftwood deposits. If these estimates are accurate, they exceed all high water-level events ever captured in the Tuktoyaktuk tide gauge record (

Table 1. The ten highest water levels recorded at the Tuktoyaktuk tide gauge station and the two highest water-level events inferred by Harper et al. [19] from the driftwood line field surveys (shown in bold text).

Rank	Date	Water Level (m CD)
1	1 September 1944	2.95
1	14 September 1970	2.95
3	4 October 1963	2.23
4	4 September 1962	2.15
5	4 November 2017	2.02
6	1 September 2013	1.96
7	17 August 2018	1.95
8	21 July 2019	1.94
9	1 September 2018	1.89
10	30 July 1963	1.85
11	27 August 2015	1.82
12	1 September 1962	1.78

). Harper et al. [19] noted that their driftwood line elevation observations are bounded by a ±0.3 m measurement error.

To examine the effects of the 1944 and 1970 event omissions from the tide gauge record on the estimated storm surge probability distributions, the extreme value analysis was repeated with the two events added to the record, based on the peak water levels inferred by Harper et al. [19]. As Harper et al.’s [19] study inferred peak water levels from the maximum elevation of driftwood deposits, the time and astronomical tidal stage associated with the peak water-level events in 1944 and 1970 could not be determined from the surveys. The uncertainty surrounding the tide level during the 1944 and 1970 peak water-level events compounds the measurement error in estimating the storm surge contribution to the total water level. Separate extreme value analyses were therefore conducted for different possible storm surge magnitudes associated with the 1944 and 1970 events, taking into account measurement error and the range of predicted astronomical tide stages on each date. Nine scenarios were tested to examine the impact of including the two peak water levels from driftwood lines in the extreme value analysis. The scenarios consisted of combinations of the survey error (i.e., +0.3 m, 0 m, and −0.3 m) and the predicted tide level for the day (i.e., minimum, average, and maximum tide level) for the two events. For example, the maximum possible storm surge for the 1944 and 1970 event would correspond to a +0.3 m measurement error and the minimum predicted tide level on that day. The upper and lower bounds of the GPD fit to storm surge events, after considering all possible combinations and estimates of the 1944 and 1970 events, are shown compared to the fit relying only on tide gauge records in Figure 5. The results show that the inclusion of these two events in the analysis significantly alters the inferred probability distribution of extreme storm surges at Tuktoyaktuk.

The baseline (best fit) return value of the storm surge associated with the 1 in 100 year return period was 2.10 m, accounting only for the surge events captured by the gauge (dotted line in Figure 5). The inferred estimates of the peak storm surges associated with the 1944 and 1970 storm surge events increased the 1 in 100-year return value to between 2.51 m (lower bound) and 3.50 m (upper bound). The range reflects uncertainty associated with the measurement errors and the timing of the peak storm surge relative to astronomical tides. The results of this analysis demonstrate the potential significance of considering non-conventional observation methods (such as driftwood surveys) to support estimates of the probability distribution of extreme storm surges in remote regions with limited continuous, long-term water-level records, with implications for flood hazard and risk assessment, and coastal engineering design.

5. Model Setup

A hydrodynamic model developed in TELEMAC-2D [24] was used to simulate storm surges in the Beaufort Sea. TELEMAC-2D solves the two-dimensional depth averaged Navier–Stokes equations (Saint-Venant shallow water equations). Bathymetry data were obtained from the General Bathymetric Chart of the Oceans (GEBCO) [34], which is derived from multiple data sources, interpolated to a grid with a resolution of approximately 500 m. The 2019 version of the GEBCO gridded dataset was linearly interpolated to the unstructured computational mesh using Blue Kenue [35] and is shown in Figure 6.

GEBCO strives to deliver a comprehensive bathymetry tied to the mean sea level (MSL) datum; however, it is an amalgamation of multiple bathymetry datasets, some of which may not be tied to MSL [36] and this is noted as a potential source of error. The shoreline of the model was based on digital shorelines obtained from Open Street Maps [37]. Closed boundaries were used for the entire domain since the effect of artificial boundaries was minimized by using a large domain and by having the open-water boundary being located in deep water. The model was calibrated by adjusting the wind drag coefficient and comparing the simulated storm surges (during open-water/ice-free periods) to the water-level residuals derived from the Tuktoyaktuk tide gauge during storm surge events. The model was forced by spatially and temporally varying wind and atmospheric pressure fields. Other model parameter selections included a constant bed friction (n = 0.025 s/m^1/3), the k-Epsilon turbulence model, and inclusion of Coriolis force effects. Each event had a simulated duration of 6 days, including 3 days prior to the peak storm surge to allow sufficient time for the model to spin up from a zero-elevation initial water-level condition.

5.1. Reanalysis Dataset

The ERA5 reanalysis dataset [13] was downloaded from the Copernicus Climate Data Store, implemented by the European Centre for Medium-Range Weather Forecasts (ECMWF). Downloaded data included hourly wind and surface pressure fields on a 0.25° (roughly 30 km) grid encompassing the study area for the period of January 1979 to the present day. As the total duration of the time period spanned by the reanalysis dataset is shorter than the tide gauge record at Tuktoyaktuk, only surge events that occurred after 1979 (i.e., within the ERA5 data time span) were simulated using the numerical model, to enable direct comparisons between the modelled storm surges and water level residuals derived from the tide gauge measurements. For context, the 50th largest storm simulated in this study is equivalent to the 89th largest storm captured on the gauge.

Time series of surface wind and pressure from the ERA5 dataset, extracted at grid points approximately co-located to weather stations near Tuktoyaktuk, were compared to the weather station records, to assess the quality of the reanalysis in the study area. Five weather stations (i.e., Stations 1699, 1700, 26987, 53581, and 53582), all in close proximity at Tuktoyaktuk, have been operated by Environment and Climate Change Canada (ECCC) since 1958. The available records at these stations have varying sampling frequencies: records were available for Station 1699 at six-hourly intevals from 1958 to 1993 at an elevation of 18.3 m (MSL), and for Station 26987 at hourly intervals from 1994 onwards at an elevation of 4.6 m (MSL). Station 1700 was operational between 1970 and 2015 and the archives for this station contain hourly measurements between 6 a.m. and 10 p.m. (Mountain Standard Time) at an elevation of 4.3 m (MSL). The archives for stations 53581 and 53582 contain hourly wind speeds and surface pressures measured at an elevation of 4.3 m (MSL), starting from 2015. Due to the overlaps in measurement periods between the weather stations, the comparisons between the ERA5 reanalysis and the measured weather data were performed using data from the weather station that had the highest measurement frequency during the selected storm surge event.

A comparison between the atmospheric pressure recorded at the weather stations and the 10 m surface pressure time series at the nearest ERA5 reanalysis grid point is shown in Figure 7 (a time series plot for a 7-day period in 2019 coinciding with a storm surge event, and a correlation plot showing measured and reanalysis pressure minima associated with 50 historical storm surge events).

ERA5 surface pressure minima showed good agreement with the atmospheric pressure minima recorded at the weather station during the fifty identified historical storm surge events, with a correlation coefficient of 0.996, a root-mean-squared error (RMSE) of 77 Pa, and less than a 1% mean difference in pressure minima across all fifty events. The coefficient of variation (defined as the ratio of the standard deviation and the mean) was less than 1 for the aforementioned values. These comparison metrics demonstrate that the ERA5 surface pressure data exhibit reasonable accuracy for application to storm surge modelling in the region, considering that a 100 Pa difference in atmospheric pressure roughly corresponds to a change in sea level of 1 cm [38].

Figure 8 shows the comparison between the wind speed recorded at the weather station and the surface wind speed estimates provided in ERA5.

Although the wind speed data from the ERA5 reanalysis grid point closest to the Tuktoyaktuk meteorological stations generally capture the magnitude, timing, trend, and variability of measured wind speeds, the storm peak wind speeds are consistently under-estimated by ERA5. For the 50 storm events investigated, ERA5 underpredicted the peak wind speeds by an average of 27% (Figure 8b). This under prediction of peak values by ERA5 held true for comparisons to both the hourly and six-hourly wind speed measurements.

5.2. Effect of Ice Presence on Storm Surges

A unique aspect of modelling storm surges in the Arctic is the need to consider the presence of sea ice, which is observed year-round under varying levels of concentration. Sea ice has been known to dampen the magnitude of storm surges by acting as a barrier between the water surface and winds [20,39]. Conversely, marginal ice concentrations (defined as the area fraction of which is covered with sea ice) can amplify the wind-driven components of storm surges by increasing the form drag, sea surface roughness and associated shear stresses [40]. Joyce et al. [22] performed a numerical investigation of storm surges in western Alaska and proposed the following modified formulation of the drag coefficient that characterizes the transfer of momentum from air and ice to the sea surface:

$$C_{d\_mod}=\left(1-AF\right)C_d+\left(AF\right)C_{d-is}+C_{d-if}$$

(1)

where AF is the area fraction of ice coverage, which can range from 0 (open water conditions) to 1 (complete ice cover); C_{d_mod} is the air–sea drag coefficient after the modification for the presence of sea ice; C_d is the air–sea drag coefficient in open water; C_d-is is the contribution of ice skin drag, set to a constant value of 0.0015 based on Lüpkes et al. [41]; and C_d-if is the form drag contribution from the ice, which depends on AF [41]:

$$C_{d-if}=4C_{d-if,max}\left(AF\right)\left(1-AF\right)$$

(2)

where C_d-if,max is the maximum value of C_d-if that occurs at 50% ice coverage. C_d-if,max was set to 0.0025 [22]. This modified drag formulation was added to TELEMAC-2D by the authors by modifying the model source code. AF was determined for each of the 50 events by importing weekly ice charts issued by CIS.

5.3. Sensitivity Study

Mesh sensitivity tests were conducted to assess the spatial resolution required to achieve a reasonable balance between computation times and model accuracy. A 1.4 million element mesh with characteristic element edge lengths in the range 500 m (nearshore) to 50 km (offshore) was finally selected. A variable time step within the range 0–30 s was applied to satisfy a Courant number less than 0.9 for numerical stability during the simulation. Computational time was roughly 1 h per event (6 day simulation period) using six cores of an Intel Core i7-8850H processor. Given these relatively short computational times, the numerical model has potential for future use as an operational storm surge forecasting model.

6. Results

6.1. Calibration and Validation

The storm surge model was calibrated by adjusting the wind drag coefficient (C_d). Three large storm surge events that occurred during the summer months (i.e., July 2019, August 2018, and September 2013) were investigated as part of the calibration exercise. The selected events corresponded to ice-free conditions in the vicinity of Tuktoyaktuk, based on weekly ice charts from the Canadian Ice Service (CIS) [42]. This was an important consideration to ensure confidence in the model’s ability to simulate storm surges under ice-free conditions, prior to modifying the source code to include the effects of ice cover on air-sea momentum transfer. The TELEMAC-2D default relationship between C_d and wind speed, based on Flather [43], was adjusted through trial and error until modelled peak storm surges were within 5% of the measured values for the three events. As ERA5 systematically underpredicts the peak wind speeds, the threshold wind speed for transition to the maximum drag coefficient was reduced. The final drag coefficient dependence on wind speed, as applied in the calibrated storm surge model, is shown in Figure 9 compared to the default Flather [43] formulation, and other popular parametrizations from Sheppard [44], Wu [45], Large and Pond [46], and Yelland et al. [47].

The time series of the measured water-level residuals and simulated (modelled) storm surges using the calibrated wind drag coefficient formula is shown in Figure 10, for three events.

For all three events, the simulated peak storm surge was within 3% of the measured water-level residual and the timing and duration of the flooding was well captured. Following model calibration, the remaining 47 historical storm surge events were then simulated (with calibrated model input parameters applied) under open-water conditions to validate the performance of the model. Figure 11a shows the performance of all fifty storm surge events that were simulated with the calibrated model.

While the majority of the fifty events were adequately simulated, events occurring during the months of December and January, when fast ice is observed on the shorelines along the north coast of Canada, showed modelled surges much greater than the measured surge when assuming open-water conditions.

Figure 12 shows a comparison between measured water-level residuals and modelled storm surges during a January 2005 winter event where high ice concentrations and fast ice were found near Tuktoyaktuk, for simulated scenarios with and without the effects of ice on air-sea momentum transfer incorporated. The latter scenario represents the model-predicted storm surge for this event under a hypothetical, ice-free condition.

With the standard wind–sea drag formulation, the hydrodynamic model over-predicted the peak storm surge by 171%. Following modification of the drag coefficient to incorporate sea ice effects, the peak predicted storm surge was within 12% of the measured peak water level residual. These findings suggest that, under certain conditions, sea ice plays a significant role in reducing peak storm surges in the Beaufort Sea during winter months. If, as climate model projections suggest, sea ice extent and ice season duration continue to decline in the Arctic [4,48] throughout the 21st century, future winter storms could generate more severe storm surges than comparable storms today. Figure 11b shows a comparison of the measured and modelled storm surges after accounting for the presence of ice—the extreme outliers present in the correlation plot for simulations without consideration for ice are eliminated. In addition, Figure 11c shows that the adjustment of the model to incorporate ice effects yielded higher correlation coefficients, lower RMSE, and smaller absolute differences for peak storm surge magnitudes.

6.2. Hindcast of the 23 August 1986 Storm Surge Event

Harper et al. [19] estimated the maximum water levels based on field observations of the driftwood lines for a significant high water-level event that occurred on 23 August 1986. This event was not captured on the Tuktoyaktuk (6485) tide gauge and was simulated to numerically validate the estimated water levels from driftwood line surveys. The simulation included the influence of sea ice as described in Section 5.2. Near Tuktoyaktuk, Harper et al. [19] estimated a peak local water level of approximately 2.1 m CD (1.6 m MSL). This provides an estimated storm surge level ranging from 1.53 m to 1.90 m, depending on the tide level during the storm, and considering the ±0.3 m measurement error noted by Harper et al. [19]. The 23 August event was hindcast in the numerical model and Figure 13 shows the simulated storm surge elevations for areas covered by Harper et al.’s [19] survey.

The simulated local maximum surge near Tuktoyaktuk was approximately 1.40 m, which falls within the bounds of the estimates by Harper et al. [19], with allowances for measurement error (i.e., 1.53 m to 1.90 m ±0.3 m). However, the model does not simulate overland flooding or wave effects, which would likely have influenced driftwood elevations, and thus the high water-level estimates inferred from driftwood line measurements cannot be precisely verified. However, the results confirm that storm surges of magnitude similar to those inferred from the survey measurements are possible, and could reasonably have occurred during the 23 August 1986 storm. This numerical exercise suggests that driftwood lines can be used to provide reasonable estimates of elevated water levels due to coastal storm events in places where tide gauge records are short, intermittent, or non-existent. Driftwood lines may also provide useful information on high water-levels in areas where gauges are operational, as they can malfunction during extreme events. The numerical analysis lends confidence to Harper et al.’s [19] estimates of peak water levels reaching up to 2.95 m, which are thought to exceed all storm surge events in the region within the past 100 years. Thus, driftwood lines may provide additional information when calculating the probability distribution of the return periods and magnitudes of storm surge events, with Figure 5 demonstrating the potentially significant implications if large events are omitted.

7. Discussion

7.1. Quality of the ERA5 Dataset

The advancement of atmospheric modelling has led to the development of high-resolution reanalysis datasets such as ERA5. Among the many atmospheric variables that these datasets offer, surface pressure and wind fields are of direct relevance for driving numerical storm surge models. Previous studies have demonstrated the improvements offered by ERA5 compared to predecessor datasets, when used to force storm surge models focused on regions in low latitudes [15,49]. This study compared surface pressure and wind data from ERA5 to measurements at stations in the Beaufort Sea region, to assess its utility for driving numerical simulations of storm surge events. Compared to local measurements, the ERA5 dataset demonstrated a high level of skill in reproducing observed surface pressures, with mean errors within 1% of the pressure minima during fifty historical storm surge events. The ERA5 reanalysis generally agreed with the timing, variability, and magnitude of mean wind speeds at weather stations but underestimated peak storm wind speeds by 27%, on average. The wind drag coefficient resulting from the calibration process described resulted in drag coefficient values exceeding those typically found in the literature [50]. Although some variability in drag coefficients may be expected due to a variety of factors, including fetch, sea surface roughness, and the methodology employed to measure wind drag coefficients, the majority of the discrepancies in the calibrated drag coefficients from typical values in the literature can be ascribed to ERA5’s underestimation of the peak wind speeds. Since C_d is proportional to the square of wind speed, the 27% underestimation of peak wind speeds by ERA5 (on average for the 50 events investigated) would require almost a doubling of the drag coefficient (or 88% increase) to compensate, which is consistent with values shown in Figure 9. This was not a significant limitation for numerical simulation of storm surges, as scaling the wind speeds (e.g., by adjusting drag coefficients) gave accurate results when the model was calibrated to observations. Simulated peak storm surges associated with fifty-one historical surge events showed good agreement with the water-level residuals derived from the Tuktoyaktuk tide gauge record and inferred from driftwood surveys.

7.2. Quality of Historic Tide Gauge Records and Impacts on Uncertainty

Extreme value statistics are often employed to describe the annual exceedance probability or recurrence interval of a natural phenomenon for risk assessment or engineering design applications. The quality of the extreme value analyses is intrinsically linked to the length and quality of the data records on which they rely on.

In this study, return period estimates for peak storm surge magnitude were determined based on water level records available from the Tuktoyaktuk tide gauge. As discussed in Section 4, the tide gauge at Tuktoyaktuk was operational starting from 1961. However, the gauge record is interrupted by periods of inactivity where the gauge was not operational. Most notably, the tide gauge record does not include two high-water-level events that occurred in 1944 and 1970, which Harper et al. [19] estimated, based on driftwood line surveys, to exceed any that were captured in the tide gauge record. It is uncertain whether these high-water events coincided with a low, mid, or high astronomical tidal stage. Therefore, it is not possible to elucidate a single storm surge magnitude for the two events and rather, a range of possibilities exist. However, even low-end estimates of storm surge magnitudes associated with the 1944 and 1970 events (approximately 2.10 m), inferred from the driftwood surveys, exceed any other water level residual events in the tide gauge record (up to 1.90 m). This suggests that the 1944 and 1970 events are the highest known storm surge events to have occurred in Canada’s western Arctic region in recent history. As such, the inclusion, or omission, of these events in extreme value analyses was demonstrated to have a significant impact on the GPD fit.

7.3. Effects of Sea Ice Conditions on Storm Surges

The use of the formulas presented by Joyce et al. [22] to account for the effect of ice on storm surges in western Alaska significantly improved the accuracy of the simulated peak storm surges. The modified drag formulation correctly accounted for the damping effects of sea ice on storm surges for winter events where fast ice was present in the vicinity of Tuktoyaktuk (Figure 11). These results demonstrate the role played by sea ice in attenuating winter storm surges in the Beaufort Sea (Figure 12). These findings are significant in the context of the lengthening open-water season durations in the Arctic, decreasing ice concentrations and declining ice-cover extents [48]. These changing ice conditions have potential to result in higher storm surges and associated impacts, and future studies that forecast storm surges under various climate change scenarios and ice conditions may be useful for community planners along Arctic coasts.

7.4. Limitations and Research Needs

The formulas presented by Joyce et al. [22] consider only the effect of the presence of ice on storm surges. However, the type of ice can be a significant factor controlling its effect on the propagation of long waves. Land-fast ice attenuates waves more than floating ice, and the thickness of ice can also influence the degree of damping [38,51]. Ice floes can scatter wave energy, with the magnitude of scatter a function of wave frequency, ice floe distribution, and the ratio of wave length to ice floe length [39]. Differentiating the effects of fast ice and floating ice on the wind drag coefficient in the formulas by Joyce et al. [22] may improve the performance of the numerical model. Available ice charts from CIS would support advanced formulations as it already classifies the types of ice, thickness, and stage of development, in addition to ice concentrations. However, ice data at more frequent intervals (e.g., daily or sub-daily) would enable the influence of shifting ice conditions to be more accurately captured in storm surge models.

While this study relied on the recently released, high-resolution ERA5 global reanalysis as the source of atmospheric forcing data, there are potential alternatives. The quality of other reanalysis datasets, such as the Japanese 55-year reanalysis [52] and the MSC Beaufort Wind and Wave Reanalysis [53], were not investigated. The latter has a higher spatial resolution than ERA5 (2 km grid resolution compared to 30 km). However, the spatial coverage of the MSC data is limited compared to the ERA5 reanalysis, and surface pressure data is not available. A comparative analysis of surface wind fields from various reanalysis datasets and weather station records, would help to identify the strengths and limitations of atmospheric forcing data for storm surge modelling in the Beaufort Sea region.

The numerical storm surge model was calibrated using data from a single tide gauge at Tuktoyaktuk, which has the longest standing, and most complete records in the Beaufort Sea region. Future studies should examine other tide gauges in the region and simulate events where storm surges were captured on multiple gauges. This would lend insight to model bias outside of the Tuktoyaktuk region. Additionally, during the study, ERA5’s reanalysis period was from 1979 to present but, just recently, it was extended from 1950 to the present [54]. This provides the opportunity to simulate the September 1970 storm surge event, believed to be one of the highest storm surges that has occurred in the Beaufort Sea [19]. This exercise could further validate the use of driftwood lines as a non-conventional source for high water-level marks and provide information on the potential magnitude of extreme storm surges in the Beaufort Sea. Furthermore, a full, long-term water level or storm surge hindcast for the entire ERA5 duration could help to fill gaps in tide gauge records and identify historical storm surges that were not recorded. The regional model could then be used to feed higher-resolution inundation models for individual communities to enable risk assessment. The model also offers potential for application to investigating climate change impacts on storm surge frequency and severity. Forcing the storm surge model with Regional Climate Models (RCM) and/or simulating storm surge events under future hypothetical sea ice conditions [55] could help to quantify changes in storm surge risk, and identify implications to guide planning for coastal resiliency.

8. Conclusions

A numerical investigation of storm surges in the Beaufort Sea was conducted. A two-dimensional hydrodynamic model was developed and forced by surface pressure and wind fields to simulate storm surges in Canada’s western Arctic region. The model was calibrated and validated using tide gauge records and high water-level marks inferred from driftwood field surveys. The surface drag formulation was modified to incorporate the effects of sea ice on air–sea momentum transfer, significantly improving the accuracy of simulated peak storm surges during periods of ice cover. Surface pressure and wind data from the recently released ERA5 global reanalysis was compared to weather station records to assess its utility for driving numerical models of storm surges in the Beaufort Sea region. The following conclusions were drawn:

• Inferences based on driftwood surveys and astronomical tide predictions identified two historical storm surge events that exceeded all events captured by the regional tide gauges. One of the events was successfully verified by a hindcast using the numerical model.
• Inclusion of the two historical events in an extreme value analysis significantly altered the best fit probability distribution of storm surges in the region, with implications for risk assessment and coastal engineering design.
• The ERA5 reanalysis was shown to be an adequate source of surface pressures and wind speeds, except for a systematic underprediction of peak wind speeds at the Tuktoyaktuk weather station. Acceptable numerical modelling results were achieved through calibration of the wind drag coefficient; large wind drag coefficients were prescribed to compensate for the generally low peak wind speeds reported in ERA5.
• Sea ice has a significant damping impact on storm surges along the Beaufort Sea coast and should be considered in storm surge analyses for the region. Under open-water scenarios representative of future Arctic conditions, storm surges were found to be as much as three times higher than under present-day ice conditions.