Articles | Volume 16, issue 11
Development and technical paper
14 Jun 2023
Development and technical paper |  | 14 Jun 2023

Adding sea ice effects to a global operational model (NEMO v3.6) for forecasting total water level: approach and impact

Pengcheng Wang and Natacha B. Bernier

In operational flood forecast systems, the effect of sea ice is typically neglected or parameterized solely in terms of ice concentration. In this study, an efficient way of adding ice effects to the global total water level prediction systems, via the ice–ocean stress, is described and evaluated. The approach features a novel, consistent representation of the tidal relative ice–ocean velocities, based on a transfer function derived from ice and ocean tidal ellipses given by an external ice–ocean model. The approach and its impact are demonstrated over four ice seasons in the Northern Hemisphere, using in situ observations and model predictions. We show that adding ice effects helps the model reproduce most of the observed seasonal modulations in tides (up to 40 % in amplitude and 50 in phase for M2) in the Arctic and Hudson Bay. The dominant driving mechanism for the seasonal modulations is shown to be the under-ice friction, acting in areas of shallow water (less than 100 m) and its accompanied large shifts in the amphidromes (up to 125 km). Important contributions from baroclinicity and tide–surge interaction due to ice–ocean stress are also found in the Arctic. Both mechanisms generally reinforce the seasonal modulations induced by the under-ice friction. In forecast systems that neglect or rely on simple ice concentration parameterizations, storm surges tend to be overestimated. With the inclusion of ice–ocean stress, surfaces stresses are significantly reduced (up to 100 % in landfast ice areas). Over the four ice seasons covered by this study, corrections up to 1.0 m to the overestimation of surges are achieved. Remaining limitations regarding the overestimated amphidrome shifts and insufficient ice break-up during large storms are discussed. Finally, the anticipated trend of increasing risk of coastal flooding in the Arctic, associated with decreasing ice and its profound impact on tides and storm surges, is briefly discussed.

1 Introduction

The ice conditions in the Arctic are changing rapidly. As the onset of the ice season is delayed and the return to the ice-free season is advanced (Johnson and Eicken2016; Parkinson2022), the period of exposure to coastal flooding is lengthened. The provision of an accurate and timely forecast of total water level (TWL) in ice-infested waters is thus becoming increasingly important. Under global warming, the increasing effects of the receding ice that protects the shorelines, combined with permafrost thawing that leads to coastal erosion, are resulting in increased exposure to coastal hazards. Many coastal communities in the Arctic and nearby bays and seas are already affected by larger storm surges and rising sea level (Pörtner et al.2022). For example, Shishmaref, a village on an island off the coast of northern Alaska, is facing the prospect of relocation. Tuktoyaktuk, the major port of the western Canadian Arctic, is experiencing severe coastal erosion (Whalen et al.2022), and its shoreline protection structures have been rapidly destroyed by storm surges and accompanying waves.

Sea ice affects both tides and storm surges, the dominant components of TWL, by adjusting the air–sea momentum flux and providing additional friction to the underlying ocean flow. In situ observations made by tide and bottom pressure gauges have shown remarkable seasonal variability in the M2 tidal amplitude in many parts of the Canadian Arctic, including the Beaufort Sea and the Amundsen Gulf (up to 50 %; Henry and Foreman1977; Godin and Barber1980), the Kitikmeot Sea (50 %–60 %; Rotermund et al.2021), and the Hudson Bay (HB) system (8 %–40 %; Prinsenberg1988; St-Laurent et al.2008). Large variability in the M2 amplitude was also reported in the Russian Arctic (Kulikov et al.2018, 2020), with up to 63 % in the Chukchi Sea (CS) and 9 % in the White Sea (Fig. 1). (We note, however, that the last two amplitude changes are calculated with respect to the annual mean and are thus larger compared to this and other studies that calculate changes using the maxima as the reference.) Significant delay or advance in the winter M2 phase was also observed, with up to 40 in the CS (Kulikov et al.2018) and eastern HB (Prinsenberg1988). On the Ross Ice Shelf of the Antarctic, analyses of global positioning system solutions together with tide gauge data (Ray et al.2021) reveal a counterintuitive M2 seasonal cycle, associated with a suppressed amplitude (10 %) and retarded phase during the ice-free season. Recently, altimeter-derived data at high latitudes were also used to study the M2 seasonality for the Arctic and connected regional seas (Bij de Vaate et al.2021). Although hampered by low temporal resolution and the presence of ice cover, Bij de Vaate et al. (2021) showed opposing responses to winter ice condition, with the M2 phase delayed in most of the Arctic but advanced in the HB.

Figure 1Tide gauges (circles) and in situ moorings (triangles) used in the present study. Red/green symbols indicate that data are available/unavailable during our study period from November 2018 to April 2022 (see Fig. 2 for data availability). The contour map shows bathymetry features in meters. Abbreviations are used for the White Sea (WS), East Siberian Sea (ESS), Chukchi Sea (CS), Beaufort Sea (BS), Canadian Arctic Archipelago (CAA), Hudson Bay (HB), and Gulf of St. Lawrence (GSL). The three subregions, their abbreviations, and stations numbers are as follows: (1) Arctic (Arctic, 1–27), (2) North Atlantic and Hudson Bay (NAHB, 28–43), and (3) northwestern Atlantic (NWA, 44–58).

To understand the underlying physics leading to the seasonal modulation of tides, it is desirable to isolate processes at play. Modeling studies can help separate ice effects from other relevant processes, such as the nonlinear tide–surge interaction (TSI; Bernier and Thompson2007) and baroclinicity (Müller et al.2014). Using a coupled ice–ocean model, St-Laurent et al. (2008), and later Müller et al. (2014), showed that the observed seasonal M2 modulation in the HB, derived from bottom pressure records, can be largely accounted for by the under-ice friction. As both studies focused on ice processes only, TSI and baroclinicity were not examined. Other studies are based on tide-only models, with under-ice friction expressed as additional bottom friction parameterized solely in terms of ice concentration (e.g., Dunphy et al.2005; Kleptsova and Pietrzak2018) or applied over landfast ice only (Bij de Vaate et al.2021; Rotermund et al.2021). These simple methods help produce the ice-induced modulation of tides over particular regions and periods but cannot account for its complex spatial and temporal variability.

For storm surges, ice-induced attenuation has been observed in the Baltic Sea (Lisitzin1974) and Beaufort Sea (Henry1975). Efforts have been made to include such effects in storm surge modeling. Kowalik (1984) and Danard et al. (1989) applied models that include ice–ocean interactions in the Beaufort and Chukchi seas, but neither one verified the ice effects on surges for winter storm events. Zhang and Leppäranta (1995) applied an ice–ocean model in the Baltic Sea and found that the sea surface slope in ice-covered cases may get down to one-third of the ice-free value. More recently, Joyce et al. (2019) incorporated the ice effects on surges through parameterizations of the wind drag coefficient and showed improvements on the coast of Alaska over particular periods. Kim et al. (2021) adopted the method of Joyce et al. (2019) and showed improvements for simulated peak winter surges at Tuktoyaktuk. However, one particular challenge with such a parameterization is that ice concentration alone cannot fully represent the internal ice stress or the ice strength, which is important for the ice drift response to winds (e.g., Fissel and Tang1991; Heil and Hibler2002) and the subsequent ice–ocean momentum transfer. The ice strength is usually a function of both ice concentration and ice thickness (e.g., Heil and Hibler2002); for example, higher ice concentration and thicker ice can enhance the ice strength and reduce the ice–ocean momentum transfer.

In Canada, sea ice effects on TWL forecasts are a major concern; sea ice is a prominent feature in the Canadian Arctic and Hudson Bay and to a much lesser extent on the east coast of Canada. This process is missing in the recently developed global high-resolution (1/12) TWL system (Wang et al.2021, 2022) running operationally at Environment and Climate Change Canada (ECCC). The system is currently under active development, by addressing important physical processes, whilst keeping the system easy to maintain and computationally efficient, so that ensemble forecasts can be performed and made available with sufficient lead time to allow maximum response time for the authorities and the public. Following this principle, we have developed effective and efficient methods to address TWL contributions from tides, storm surges, baroclinicity, and their interactions (Kodaira et al.2016a, b; Wang et al.2021, 2022). In the present study, we attempt to further address sea ice effects following the same principle.

As in other operational centers, ECCC has recently developed advanced operational ice–ocean systems with data assimilation and more realistic representations of ice physics and its interaction with the ocean (Smith et al.2016, 2018; Lemieux et al.2015, 2016; Roy et al.2015). These systems are generally not suitable for accurate and timely water level forecast, as their horizontal resolution is too coarse (1/4), their computational cost is not sufficiently low, and/or tides or other processes critical to TWL forecasts are not considered. However, they can offer the information necessary to account for ice effects in higher-resolution, computationally efficient systems optimized to forecast TWLs. In this study, we aim to address the following questions: (1) can we design a new parameterization to include ice effects in a global ocean model for forecasting TWL and improve forecast skill in polar regions? (2) Can we isolate and explain the contribution of dominant physical processes (e.g., under-ice friction, baroclinicity, and nonlinear tide–surge interaction) to the seasonal modulation of tides?

The structure of the paper is as follows. The observations of coastal TWL are described in Sect. 2. The ocean model is introduced in Sect. 3. The new parameterization of ice effects, via the ice–ocean stress, is described in Sect. 4. The experimental design and analysis are presented in Sect. 5. The impact of adding ice effects on the forecast skill and the underlying physics are examined in Sect. 6. The results are summarized and discussed in the final section.

St-Laurent et al. (2008)Kulikov et al. (2018)

Table 1Summary of water level observations collected from various institutes and publications (see Fig. 1 for station code and Fig. 2 for data availability). Abbreviations are used for Marine Environmental Data Service (MEDS), National Oceanic and Atmospheric Administration (NOAA), University of Hawaii Sea Level Center (UHSLC), and European Marine Observation and Data Network (EMODnet).

Download Print Version | Download XLSX

2 Observations

The present study uses 58 stations grouped into three subregions (Fig. 1). Permanent tide gauges (red circles) in the Arctic are very sparse and primarily located around the Beaufort Sea, Chukchi Sea, and northern Norway. In an effort to maximize observations available for verification, we collected data, including tide gauge records, bottom pressure records, and monthly tidal constants from various institutes and publications (Table 1), dating back as early as 1957. Our initial criteria were that stations have records of at least 12 continuous months. More details about data availability, for each station, are given in Fig. 2.

Data quality control was conducted by removing isolated and clustered spikes in TWL records and tidal residuals following careful visual inspection. In addition, the tide gauge record at station 2 before 1969 was not used, as it shows significantly different statistical properties (e.g., variance, seasonal cycle, and datum) than the remaining data. The bottom pressure record at station 40 from August 2003 to August 2004 was discarded, as it has a much coarser temporal resolution than the remainder of the record.

Figure 2Availability of observed water levels as a function of station code (see Fig. 1) from January 1957 to April 2022.


3 The ocean model

The NEMO modeling framework (Madec2008) is used to solve the governing equations (i.e., the momentum equation, continuity equation, and equations for heat and salt transport). They are as follows:


where uh represents the horizontal velocity vector (u,v), u denotes the complete velocity vector in three dimensions (u,v,w), f is the Coriolis parameter, pa denotes atmospheric pressure at the sea level, ρ0 denotes the reference density (1025 kg m−3), η denotes the sea surface height, and ηA represents the gravitational tidal potential. The depth-dependent coefficient αs is used to parameterize the impact of self-attraction and loading (Stepanov and Hughes2004). The lateral eddy viscosity and diffusivity coefficients are set to constant (Ah= 100 m2 s−1 and Kh= 10 m2 s−1), and the vertical eddy viscosity and diffusivity coefficients (Az and Kz) are determined using the turbulent kinetic energy scheme introduced in Gaspar et al. (1990).

In Eq. (1), the last term on the right-hand side represents the tidal nudging technique introduced in Wang et al. (2021). It nudges the model's depth-averaged current uh towards the observed current uobs calculated using the tidal amplitude and phase of eight major tidal constituents (M2, S2, N2, K2, O1, K1, P1, and Q1) provided by TPXO8 (Egbert and Erofeeva2002). The angle brackets indicate that the nudging is filtered temporally to isolate the variability in tidal frequency bands. The strength of the nudging is determined by a spatially varying coefficient λ(x). South of 66 N its global distribution is given by Wang et al. (2021). North of 66 N, we set λ(x) to zero because the nudging would dampen the ice-induced seasonal modulation of tides. Recall that TPXO8 does not take ice effects into account, and so nudging is not desired when ice effects are considered.

In Eqs. (3) and (4), the last terms on the right-hand side correspond to the nudging of the model's temperature and salinity (T and S), towards operational forecasts (Tf and Sf) provided by a coarser-resolution (1/4) data-assimilative system – the Global Ice Ocean Prediction System (GIOPS; Smith et al.2016, 2018). The strength of the nudging is controlled by a spatially uniform coefficient r. We set r= 0.2 d−1, and this adds the high-quality, low-frequency variability (with periods exceeding about 15 d) provided by the 1/4 model to our TWL model, while allowing high-frequency variability to evolve freely (Wang et al.2022).

At the surface, the boundary condition for Eq. (1) is given by

(5) A z u h z z = 0 = τ s ρ 0 , τ s = τ ao = ρ a C ao | u 10 | u 10 ( α = 0 ) ,

where ρa denotes the air density, and u10 is the wind velocity at 10 m height. In the absence of sea ice (i.e., ice concentration α equals zero), the surface stress τs equals the air–ocean stress τao. The air–ocean drag coefficient Cao equals 1.2 × 10−3 for |u10|< 8 m s−1 and then increases linearly with |u10|, with a slope of 0.065 × 10−3 for every 1 m s−1 increase in |u10| (Bernier and Thompson2007). Hourly fields of u10 and pa were obtained from the assimilation component of ECCC's operational Global Deterministic Prediction System (GDPS; Buehner et al.2015), with a grid spacing of roughly 15 km.

At the bottom, the boundary condition for Eq. (1) is

(6) A z u h z z = - H = τ b ρ 0 , τ b = ρ 0 C db | u b | u b ,

where τb refers to the bottom stress, ub denotes the current velocity at the bottom, and Cdb is the bottom drag coefficient, which is set equal to 2.5 × 10−3.

The extended version of a tripolar ORCA grid (eORCA12) is used as the model grid. It covers the Antarctic ice shelves and has a horizontal grid spacing of 1/12. The bathymetry is obtained from GEBCO_2014 (Weatherall et al.2015), with local adjustments made in the HB and on the Labrador and Newfoundland shelves, based on bathymetric data provided by Florent H. Lyard (personal communication). Following Wang et al. (2022), the vertical grid is nine z levels, which is able to capture baroclinic variability via T and S nudging, whilst maintaining low computational cost. Partial steps are employed for the bottom layer to achieve a more accurate representation of the bathymetry. Mode splitting was used, with time steps of 240 and 6 s for the internal and external modes, respectively.

4 Surface stress in the presence of sea ice

In this section, we describe the surface stress τs in the presence of sea ice. We review parameterizations or models previously developed and introduce a new, cost-efficient, method for its parameterization in TWL systems.

In the presence of sea ice, τs is generally approximated by a combination of the air–ocean stress τao (see Eq. 5) and the ice–ocean stress τio weighted by the ice concentration α as follows:

(7) τ s = ( 1 - α ) τ ao + α τ io .

τio can be parameterized by a quadratic drag law in terms of the relative velocity between ice and surface currents (uiceusurf) as follows:

(8) τ io = ρ 0 C io | u ice - u surf | ( u ice - u surf ) ,

where Cio is the ice–ocean drag coefficient.

To address τio, there are several options with different levels of complexity. The most complex option is to couple the ocean model with a sophisticated ice model that simulates ice thermodynamics, dynamics, transport, and ridging (Hunke et al.2010). Unfortunately, such an option is not suitable for relatively high-resolution systems built with computational efficiency in mind. A simpler option is to solve the ice momentum equation with prescribed ice concentration and ice thickness. However, the major time-consuming part in most modern ice modeling, the sub-cycling of the standard elastic–viscous–plastic solver for ice dynamics (Hunke and Dukowicz1997), is still required. In addition, the simplified coupling neglects mass transport, and so it has deleterious effects known as “artificial inertial resonance” in the presence of both tidal and wind forcing (Hibler et al.2006). Another option is to take uiceusurf directly from an external ice–ocean model, which requires only a negligible additional computational cost. The main issue with this option is the potentially large inconsistency in predicted tides between the TWL and external ice–ocean models. These differences can have various sources such as model resolution, bottom topography, open boundary conditions, and parameterization of dissipation. In the following section, we propose and evaluate a new approach to address these inconsistencies.

4.1 Mapping ice effects on currents

To address the inconsistency issue in tides associated with the use of an external ice–ocean model to define the ice–ocean stress while maintaining computational efficiency, we decompose the relative velocity into a tidal component (denoted with the superscript T) and a residual/surge component (superscript S) as follows:

(9) u ice - u surf = u ice T - u surf T + u ice S - u surf S .

As uiceT is mainly forced by usurfT, we introduce a transfer function, such that

(10) u ice T a T ( x ) R ( φ ( x ) ) u surf T ,


(11) R ( φ ( x ) ) = cos φ ( x ) - sin φ ( x ) sin φ ( x ) cos φ ( x ) ,

and the spatially varying scale factor aT and the rotation angle φ can be inferred from uiceT and usurfT provided by external ice–ocean models (the asterisk  denotes quantities from external models). Specifically, aT and φ are derived by scaling and rotating the ice and ocean tidal ellipses so that their semi-major axes are equal. The tidal relative velocity is thus written as

(12) u ice T - u surf T = a T ( x ) R ( φ ( x ) ) - I u surf T ,

where I is the identity matrix.

Figure 3Top panels show the observed frequency of the landfast ice occurrence for December 2020 to March 2021, calculated based on the weekly fast ice extent provided by the U.S. National Ice Center (2020). Middle and bottom panels show the derived monthly aT and φ for the M2 tide for the same period. Note that for areas with very weak uiceT (major axis velocity magnitude less than 5 × 10−3m s−1 in this study), φ is irrelevant, its estimation is also not reliable, and so it is set to 0.

Unlike periodic tides, storm surges are sporadic and occur on a local scale driven by the atmospheric forcing. Applying a transfer function to usurfS, similar to Eq. (10), is not feasible, since usurfS is forced by uiceS. Instead, we expect that the inconsistency in surges between our model and the external ice–ocean model is acceptable, given that their atmospheric forcings are similar. The residual relative velocity can thus be taken directly from the external model. Since the tidal and residual relative velocities are calculated based on surface currents coming from the TWL and external models, respectively, their valid surface levels could be different. To be consistent, an empirical scale factor, aS, can be used to adjust the residual relative velocities to the surface levels of the TWL model as follows:

(13) u ice S - u surf S = a S u ice S - u surf S .

For example, if the surface level in the external model is shallower than the TWL model, then usurfS will be stronger than usurfS, and aS should be larger than unity. In practice, aS can be tuned to best reproduce the observed residual water level in the presence of ice.

Finally, the total relative velocity is

(14) u ice - u surf = a T ( x ) R ( φ ( x ) ) - I u surf T + a S u ice S - u surf S .

In practice, usurfT=usurf can be obtained using an efficient online tidal filter (Wang et al.2021) denoted by the angle brackets. Note that the same filter is also used for the tidal nudging shown in the last term on the right-hand side of Eq. (1).

4.2 Ice–ocean stress

Gridded fields of hourly ice concentration (α), ice velocity (uiceT; uiceS), and surface current (usurfT; usurfS) were obtained from the assimilation component of GIOPS (Smith et al.2016, 2018) developed and run operationally at ECCC. In GIOPS, the CICE-based (Hunke et al.2010) ice component has 10 categories of ice thickness, and the NEMO-based ocean component has 50 vertical levels and a horizontal resolution of 1/4. The initialization of GIOPS involves using analyses created by Mercator Océan's System d'Assimilation Mercator, version 2 (SAM2; Tranchant et al.2008). Further information regarding the initialization procedures can be found in Smith et al. (2018). We note that although GIOPS is a global system, its model grid has a southern limit at about 77 S, which excludes ice cavities in the Ross Sea and Weddell Sea. In the present study, we thus focus on ice-infested waters of the Northern Hemisphere.

As in many global ice–ocean systems, the operational GIOPS does not include tides. To obtain uiceT and usurfT, we reran GIOPS by activating the astronomical tidal potential forcing. In the present study, the transfer function is updated monthly, which is sufficient to capture its seasonality. We focus on four major tidal constituents (M2, S2, K1, and O1) in GIOPS, as they can be adequately resolved with monthly harmonic analyses.

Figure 3 (middle panels) illustrates the monthly estimates of aT for M2 from December 2020 to March 2021. Note that aT captures the mobility of sea ice, with landfast for aT=0, non-free drift for 0<aT<1, and free drift for aT=1. In theory, only landfast ice and non-free drift ice exert friction on the underlying ocean flow. Regions of aT→0 are seen along the Arctic coast, in parts of the Canadian Arctic Archipelago (CAA), and the East Siberian Sea (ESS). Identified regions are in reasonable agreement with observed landfast ice occurrences (top panels). Non-free drift ice is found further away from the coast in the Arctic and HB, and in particular, it covers broad areas of the ESS and CS. We note that not all the identified landfast or non-free drift ice are relevant for the seasonality of tides. We return to this point later (see the end of Sect. 6.2.2). The rotation angle φ is only relevant for drift ice. Its main impact is found in the ESS and CS, where the tidal ice and ocean velocity vectors can be nearly 180 out of phase (bottom panels of Fig. 3), effectively enhancing the under-ice friction. Elsewhere, the absolute value of φ is relatively small (within 20), which has minimal effects on the calculation of the under-ice friction.

Results for the other three constituents are roughly similar to M2 (not shown), although over some regions (e.g., ESS and CS) K1 and O1 are too weak to derive reliable estimates. This similarity between constituents indicates that the transfer function, or the response of sea ice relative to tidal current, is largely determined by ice characteristics. Thus each monthly transfer function for M2 was applied to other constituents in the present study. We note that the main advantage of this new approach is that it is not sensitive to differences in predicted tides between the ice–ocean and TWL models, as the mapping of ice effects is achieved via a transfer function. Therefore, in theory, the approach can be used with any ice–ocean systems, regardless of the model skill in tides – as long as they are realistic.

4.3 Ice–ocean drag coefficient

The parameterization of the ice–ocean drag coefficient Cio is a complex issue, as Cio depends on various ice characteristics, such as surface roughness, floe size, ridge height, and keel depth (Lu et al.2011; Tsamados et al.2014). Instead, constant values are commonly used and determined by matching model predictions with observations (e.g., St-Laurent et al.2008; Rotermund et al.2021). Typical measured Cio values range from 1.05 × 10−3 to 4.70 × 10−2 in field investigations (see Table 1 in Lu et al.2011). GIOPS has 50 vertical levels, and its Cio was set to 2.32 × 10−2, based on a log layer assumption, using its first-layer currents at 0.5 m and an undersurface sea ice roughness length scale of 0.030 m (Roy et al.2015). In the TWL system, only nine vertical levels are used (Wang et al.2022), and so our first-layer currents are valid at 7.5 m. The Cio value for our model is thus expected to be smaller than that used by GIOPS, but the log layer assumption used in GIOPS is not suitable for our first layer. Based on sensitivity tests, we set Cio to 1.00 × 10−2, which produces reasonable agreement between observed and predicted seasonal M2 modulations.

The empirical scale factor for residual relative velocity, aS in Eq. (14), was tuned to 1.64, which produces reasonable agreement between observed and predicted residuals. To offer an alternative interpretation of this value, we note that the resulting drag coefficient (aS)2Cio for the residual stress based on uiceS-usurfS is 2.69 × 10−2, indicating a slightly higher roughness length scale of about 0.036 m, compared to 0.030 m used in GIOPS (Roy et al.2015), which is well within the range given in the literature. As an example, McPhee (2008) provides a mean value of 0.049 m, with a standard deviation ranging between 0.016 and 0.146 m, based on estimates for a typical multiyear sea ice floe.

5 Experimental design and analysis

Two basic runs, RunAO and RunAIO, were conducted to examine the impact of adding ice effects on predicted water levels. In RunAO, ice effects are not considered, and τs equals τao (see Eq. 5). In RunAIO, τs is computed as the combination of τao and τio (see Eq. 7). In order to quantify the contribution of individual physical processes on the seasonality of tide, four process-oriented runs were also conducted by gradually removing the relevant processes from RunAIO, including the baroclinic effects (by using constant T and S; Run1), TSI due to τio (Run2), TSI due to τb (Run3), under-ice friction (by setting uiceT-usurfT to zero; Run4). Specifically, for Run2, TSI due to τio was removed by setting τio=ρ0Cio(|uioT|uioT+|uioS|uioS), where uioT and uioS are, respectively, the tidal and residual relative ice–ocean velocities given by Eqs. (12) and (13). For Run3, TSI due to τb was removed by setting τb=ρ0Cdb(|ubT|ubT+|ubS|ubS), where ubT and ubS are isolated using the online tidal filter of Wang et al. (2021), with ubT=ub; ubS=ub-ub.

Table 2Design of the model experiments.

Download Print Version | Download XLSX

The six runs are summarized in Table 2. Note that we chose to remove processes gradually instead of the more traditional removal of a process at a time because it is not possible to completely isolate the under-ice friction which is a prerequisite for the TSI due to τio. Sensitivity studies (not shown) confirm that the impact of the combined removal approach on other processes that can be isolated (i.e., the TSI due to τio and τb) is negligible. Each model run starts on 21 September 2018 and finishes on 30 April 2022. The first 40 d of a run are discarded to allow for model spin-up, which is mainly determined by the spin-up of the tidal nudging (Wang et al.2021).

We use the root mean square error (RMSE) to evaluate the model performance for individual stations. To facilitate the comparison with different scales, we also compute the root mean square (rms) of the observations. Both metrics were calculated for TWL, tides, and tidal residuals at each station. Tides were reconstructed based on the eight major tidal constituents used in this study (see Sect. 3), using the T_TIDE package of Pawlowicz et al. (2002).

Monthly harmonic analyses of observed and predicted TWL were conducted to examine the ice-induced seasonality of tides. It is noted that for diurnal tides, there is substantial variability in the standard monthly analysis, possibly due in part to the contamination from non-tidal energy (Cartwright and Amin1986). To minimize such an effect and focus on the seasonal variability, we conducted another set of monthly analyses, using a sliding window of 90 d to obtain the estimates for diurnal tides only. The unresolvable constituent K2 (P1) was inferred from S2 (K1), and the inference parameters including amplitude ratios and phase differences were taken from the yearly analysis. Nodal corrections were performed. Estimates for stations with a signal-to-noise ratio (SNR; see Pawlowicz et al., 2002, for details) lower than 2 are not used.

For each year, we calculate the normalized amplitude anomaly (ΔÃi) and phase anomaly (Δϕi) relative to their corresponding values in September (ASept, ϕSept), when sea ice has the minimum cover in the Arctic. Thus,


where the subscript i denotes a particular month. The ice-induced maximum modulation occurs in March (ΔÃMar; ϕMar) when sea ice reaches its maximum. We note that the seasonality of tides can be caused by a variety of mechanisms, including astronomical motions, frictional/advective interactions, and climate processes (e.g., baroclinicity, sea ice, and river discharge; Ray2022). In the main text, we focus on the seasonality of the dominant M2 constituent, which has negligible astronomical contribution. Other minor constituents, including S2, K1, and O1, are briefly discussed in the Supplement.

Figure 4The rms of observations (black) and RMSE for RunAO (blue) and RunAIO (red) for TWL (a), tides (b), and storm surges (c). All rms and RMSE values are in centimeters. Note that the missing stations are those for which only historical observations are available, i.e., data do not overlap with the study period.


Figure 5Modulation of the M2 amplitude (ΔÃMar; panels a, c, and e) and phase (ΔϕMar; panels b, d, and f) in March relative to September. The contour map shows results predicted by RunAO (a, b) and RunAIO (c, d) averaged over 2019–2021. Filled circles show results taken from available observations during 1957–2021. (e, f) Comparison between observation and prediction as a function of station code (Fig. 1). Shaded areas indicate the 10th–90th percentile range. Only stations with a SNR greater than 2 are plotted.

6 Results

6.1 Total water level

We first compare the model skill for RunAO and RunAIO in predicting TWL in terms of RMSE at 34 permanent tide gauges (Fig. 4a). Improvements from the addition of ice effects are seen in the Arctic, most particularly in the Canadian Arctic (stations 2, 9, and 12). The reductions in RMSE are relatively small (i.e., 1.0–3.3 cm). This is expected, considering that sea ice matters mostly in winter and, in particular, during winter storms which are relatively rare in parts of the Arctic. For example, over the four ice seasons of the study period, there are only two storms at station 2 and no storms at station 9, both of which are located in the CAA. We note, however, that the impact on peak water levels can be large (up to 1.0 m; see Sect. 6.3). The impact at stations in the northern North Atlantic (28–31) and Gulf of St. Lawrence (44–58) is negligible, indicating that the predicted ice is mostly in free drift over these regions. A slight increase in RMSE value of about 4.0 cm is noted at Churchill (station 39) in the HB. This is largely due to the existing bias in the predicted tides under ice-free conditions, thus adding the ice-induced modulation increases the bias in the ice season (Fig. 4b). However, we note that observations at Churchill have possible quality or drift issues, as observed tides have undergone large changes since 1998 (Ray2016). Finally, we note that this evaluation is limited to permanent gauges which are very sparse in the Arctic and HB. Next, we focus on the ice effects on the seasonality of tides at all available stations and storm surges during large storm events.

6.2 Tides

6.2.1 Seasonal variability

Figure 5 shows the M2 modulation in March relative to September (ΔÃMar, ΔϕMar). The ice effects on tides can be understood as a combination of two processes, namely (1) direct frictional effects that result in amplitude reductions (i.e., negative ΔÃMar) and phase delays (i.e., positive ΔϕMar) and (2) indirect effects through the amphidrome shift that leads to both positive and negative changes in amplitude and phase. Comparison of predictions given by RunAO and RunAIO (Fig. 5a–d) reveals large ice-induced modulations in the Arctic and HB, with the largest modulations occurring around amphidromic points. Adding ice effects generally reduces the amplitudes at the coast (Fig. 5c), while the opposite also occurs due to the amphidrome shift. Adding ice effects also leads to phase delays in most of the Arctic and the CAA and phase advances in parts of the CAA and most of the HB (Fig. 5d) due to contributions from both the direct and indirect effects of friction. We return to the point regarding the amphidrome shift later (see Sect. 6.2.2). Figure 5a and b also reveal non-negligible M2 modulations in RunAO, indicating that nonlinear TSI and/or baroclinicity also contribute to seasonal variability (see Sect. 6.2.3 for details).

Observations (filled circles) are also plotted on top of predictions in Fig. 5a–d. Their comparison is further plotted as function of station code in the bottom panel. Adding ice effects significantly improves the model skill in predicting ΔÃMar (up to 40 %) and ΔϕMar (up to 50) at most stations in the Arctic and HB. In the Arctic, improvements are also found for other smaller constituents (S2, K1, and O1; see the Supplement for additional details). For M2, one exception is station 18, where observations show an anomalous positive ΔÃMar (Fig. 5e), while RunAIO generates a negative ΔÃMar. The only possible scenario for positive ΔÃMar is the indirect effect of friction via a shift of the local amphidrome away from the coast where station 18 is located. Unfortunately, the shift predicted by RunAIO is roughly parallel to that coast. We note, however, that RunAIO reproduces the observed positive ΔϕMar at station 18 well because the direct frictional effects on ΔϕMar over this area are much stronger than the indirect effect associated with the amphidrome shift.

Figure 6Normalized monthly M2 amplitude anomaly (ΔÃi; %) relative to September at 14 selected stations. Observed mean (black line) and 10th–90th percentile range (gray shading) are presented, based on available data from 1957–2022. Model predictions from 2019–2022 are provided by RunAO (blue lines) and RunAIO (red lines). The title of each subplot gives the station code, station name, and averaged M2 amplitude (in cm) in September (number in square brackets). The locations of selected stations and their numbering are shown in the bottom right panel.

Figure 7Same as Fig. 6 but for the monthly M2 phase anomaly (Δϕi) relative to September.

Figures 6 and 7 show monthly time series of ΔÃi and Δϕi at 14 stations. Selected stations cover various geographical areas, where noticeable amplitude or phase modulations are observed, and they all have a M2 amplitude of at least 9 cm and a phase modulation of at least 5. Observations show large-amplitude reductions (up to 40 %–50 %; Fig. 6) in the Canadian Arctic (stations 6, 7, and 12), CS (stations 14 and 17), HB (stations 42 and 43), and Northumberland Strait (station 53) and large phase modulations (up to 40–50; Fig. 7) at Tuktoyaktuk (station 12), in the CS (stations 14 and 17) and the eastern HB (station 42). Large modulations can last up to 8 months of the year (e.g., station 12). Model results show that in the absence of ice-induced stress (RunAO), the predicted ΔÃi and Δϕi are pretty flat, except in the CS (stations 14 and 17), where other processes (e.g., TSI and baroclinicity; see Sect. 6.2.3 for details) contribute up to 20 % to the amplitude modulation. When ice-induced stress is included (RunAIO), forecasts of ΔÃi and Δϕi are greatly improved at most stations across the season.

We note that there remains room for improvement. For example, ΔÃi in RunAIO (Fig. 6) shows slight overestimations ( 10 %) at Resolute and Churchill (stations 4 and 39), moderate underestimations (20 %) at Tuktoyaktuk and Inukjuak (stations 12 and 42), and large underestimations (40 %) at Shediac Bay (station 53). We note that Shediac Bay is located near the narrow Northumberland Strait (width of about 13 km) that the 1/4 external model cannot resolve. For Δϕi (Fig. 7), the results show moderate (10–20) underestimation at Tuktoyaktuk, Wrangle, and La Grande Rivière (stations 12, 17, and 43) and overestimation at Red Dock dock (station 14). Some of these discrepancies can be explained by the overpredicted amphidrome shifts. We return to this point later (see Sect. 6.2.2).

We next examine the interannual variability in the seasonal M2 modulation at five permanent gauges during the four ice seasons of the study period (Fig. 8). Observations show interannual variability in both the duration and magnitude of the maximum modulation; changes in duration are up to 2 months (e.g., station 14), while changes in magnitude are up to 10 % in amplitude and 10 in phase. These features are apparently missed by RunAO, while they are reasonably captured by RunAIO. One exception is Churchill (station 39), where observations show almost no modulations in amplitude, while RunAIO generates 5 %–10 % modulations. We note, as before, that observed tides at Churchill may have quality or drift issues (Ray2016).

Figure 8Monthly normalized M2 amplitude anomaly (ΔÃi; left panels) and phase anomaly (Δϕi; right panels) relative to September 2019 at five permanent gauges observed (black), RunAO (blue), and RunAIO (red) for the period from November 2018 to April 2022.


6.2.2 Ice-induced amphidrome shift

Ice-induced tidal modulations are associated with shifts of amphidromes. In this section, we examine these shifts. We focus on March, at the peak of the ice cover, and examine the M2 amphidrome shifts averaged over the study period (Fig. 9). We note that interannual variabilities in shifts are relatively small, except for several small amphidromes in the CS. In general, facing the direction of the shift, amplitudes decrease in the front, while they increase in the back. Still facing the direction of the shift, phase delays and phase advances occur on the left and right sides, respectively. The largest shift is found in the ESS where the amphidrome (marked “A”) moves towards the coast by 90–125 km (note that the arrows in Fig. 9b do not scale with the background distances) due to the ice-induced strong tidal dissipation on the onshore side of the system (see Fig. 5c). The second largest occurs close to the center of the Arctic, where the system (marked “B”) moves towards the Canadian Arctic by 70–90 km. The two shifts are clearly responsible for the dominant large-scale features of M2 modulations in the Arctic (see Fig. 5). There are also small to moderate shifts (10–50 km) of numerous systems in the Russian Arctic and across the Bering Strait, which affect regional, small-scale modulations.

In the CAA and HB, the frictional effects in Taylor's problem of reflection of Kelvin waves in semi-enclosed basins (Taylor1922) explain most of the shifts. The primary effect is the exponential decay of Kelvin wave amplitude in the direction of wave propagation, which causes the amphidromes to shift towards the coast where the reflected Kelvin waves travel (Rienecker and Teubner1980; Prinsenberg1988; Roos and Schuttelaars2011). This behavior applies to both real amphidromes (i.e., amphidrome over the ocean, as marked by “D”–“G”) and virtual amphidromes (i.e., amphidrome over land, as marked by “C” and “H”). We note that real amphidromes may become virtual. This is the case for “D” and “E” in the CAA, which shift over land due to the frictional effects.

These shifts explain the observed and predicted M2 modulations in the CAA and HB shown in Figs. 6 and 7. For example, the shift of “D” leads to the phase advance at station 7 located to its right side (recall that the relative direction is referred to the direction facing the amphidrome shift). The shift of “C” causes the phase delay at station 4 located to its left side. The shift of “E” is found to be responsible for the overestimated amplitude modulation at station 4, suggesting this shift is overestimated. In the HB, the shifts of “F” and “G” led to the phase advance in most of the HB. In the southern extension of HB, the shift of “H” led to a phase delay on its left side, where station 43 is located, which counters the phase advance induced by the shifts of “F” and “G”. The overall effect is an underestimated phase advance at station 43, indicating that the shift of “H” is overestimated.

Figure 9(a) M2 amphidromes in March predicted by RunAO, with amplitude (m) in the contour map and phase (every 30) shown with white lines. (b) Ice-induced shifts of M2 amphidromes in March, taken as the difference in the predicted amphidromic points between RunAIO and RunAO. The unfilled circles denote real amphidromes. Their shifts, averaged for 2019–2022, are denoted by colored arrows. The filled circles denote virtual amphidromes or amphidromes that become virtual after the shift. Their positions and shifts (denoted by black arrows) are for illustration purpose only (not to scale).

With ice conditions changing rapidly in the Arctic, it is important to know where the presence of ice or the under-ice friction is most relevant for the amphidrome shifts shown above. Similar to bottom friction, we expect that the impact of the under-ice friction is most relevant over shallow waters where tidal dissipation is significant. We thus conducted four additional sensitivity runs by applying the under-ice friction over regions with water depths of less than 50, 100, 150, and 200 m, respectively. We found that applying the friction over a water depth less than 100 m can reproduce almost all modulations and the associated amphidrome shifts of RunAIO. These important regions (see Fig. 1 for bathymetry), combined with significant presence of landfast or non-free drift ice (see the middle panels of Fig. 3), cover the bulk of the ESS and CS and the shallow waters (less than 100 m) of the Canadian Arctic and the HB system. This also indicates that although there are large volumes of landfast ice in the deeper waters (greater than 100 m; middle panels of Fig. 3) of the CAA, due primarily to arch formation between islands, their impact on tides is insignificant.

Figure 10Differences in ΔÃMar (top panels) and ΔϕMar (bottom panels) for M2 tides between process-oriented runs, including Run3–Run4 (left), Run2–Run3 (middle left) and Run1–Run2 (middle right), and RunAIO–Run1 (right), corresponding to effects of the under-ice friction, TSI due to τb, TSI due to τio, and baroclinicity.

Figure 11Time series of ηW observed (black) and predicted by RunAO (blue) and RunAIO (red) at five permanent gauges during February 2019–April 2020. All values are in meters.


6.2.3 Impact of tide–surge interaction and baroclinicity

We next examine the individual contributions from TSI, baroclinicity, and under-ice friction to ΔÃMar and ΔϕMar for M2 (Fig. 10), based on process-oriented runs (Runs 1–4; Table 2). As expected, the under-ice friction (left panels) has the largest influence, while the effect of TSI due to bottom stress is negligible (middle left panels), due to weak bottom currents. The effects of TSI due to τio (middle right panels) and the effects of baroclinicity (right panels) are non-negligible. The two mechanisms act in different ways, leading to different spatial and temporal signatures. Strong TSI due to τio occurs predominantly over the ESS and CS in March, due to the combination of large ice cover and strong wind-driven surface currents induced by frequent winter storms. Its main effect is to reduce the local M2 amplitude (> 10 %), resulting in shifts in many small to moderate local amphidromes. It also drives a phase delay of 10–40 over most of the ESS, CS, and Beaufort Sea. Overall, over most of the affected areas, this mechanism reinforces the modulations induced by the under-ice friction (compare panels to the left and middle right).

In contrast, we found that modulations induced by baroclinicity occur mainly in September. This is consistent with Müller et al. (2014), who found that the annual maximum tide occurs in summer in the western Yellow Sea and North Sea. This can be explained by baroclinic effects on the vertical profile of eddy viscosity (Müller et al.2014); the presence of the pycnocline leads to a stabilized water column and thus reduced tidal dissipation through turbulent processes. Baroclinic effects appear to have a larger scale, mainly affecting several relatively large amphidromes in the Arctic (right panels of Fig. 10). In September, relative to March, this mechanism leads to increased amplitudes (up to 10 %) in the ESS, CS, and north of Norway. The corresponding amphidrome shifts are responsible for most of the phase modulations (about 10). Compared to the left panels of Fig. 10, baroclinicity also reinforces the modulations induced by the under-ice friction.

Figure 12Same as Fig. 11 but for 2021–2022.


6.3 Storm surges

Storm surges are primarily driven by surface winds and air pressure, and they are usually represented as the tidal residuals, the differences between observed water levels, and predicted tides. However, tidal residuals also contain high-frequency contributions, such as instrument errors and seiches (see Fig. 5 of Wang et al.2022), which could interfere with the comparison of surges. To attenuate such high-frequency contributions, we applied a low-pass filter with a cutoff period of 10 h to tidal residuals to obtain ηS. We further decompose ηS into a wind and pressure-driven component so that

(17) η S η W + η P ,

where ηW is the isostatic wind adjustment part of ηS, and ηP is the inverse barometer effect. A significant contribution to the variability in surges comes from ηP, and this part of the surge signal is almost unaffected by τio. This causes difficulties in visualizing and analyzing the ice effects. For this reason, in this section we remove ηP from the surge level (ηW=ηS-ηP). To obtain ηP, we produced an inverse barometer-only prediction (i.e., a run driven with surface air pressures only). We then removed the predicted ηP from both observed and predicted ηS. Figures 11 and 12 show time series of observed and predicted ηW at five permanent gauges where the ice effects are sufficiently large to affect ηW. Adding ice effects in RunAIO significantly improves the model skill during storm events. For example, it attenuates the peaks in ηW by up to 1.0 m at Tuktoyaktuk (station 12) and 0.25 m at Alert (station 2). Low-frequency variability is also improved (e.g., at Alert) due to the persistent presence of sea ice. However, some over-attenuated peaks in ηW are also found, particularly at Prudhoe Bay (station 13) with, for instance, the largest negative ηW in 2019 and largest positive ηW in 2020 (third row of Fig. 11). Further investigation shows that increasing the ice–ocean drag coefficient even to unrealistically large values does not help. This leads us to speculate that ice velocities are underestimated by GIOPS, possibly as a result of insufficient ice break-up during strong storms.

Figure 13Time series of (a) wind stress provided by the GDPS and (b, c) observed and predicted zonal ice and surface current velocities at an offshore mooring near Jones Islands, located about 50 km west of Prudhoe Bay (station 13), from January to April 2020. (d) Time series of observed and predicted ηW at Prudhoe Bay.


To verify this speculation, we use in situ measurements of ice and surface current velocities collected at an offshore mooring (station S2 offshore in Hošeková et al.2021) located only about 50 km west of Prudhoe Bay. The wind stress, ice, and current velocities are primarily zonal and parallel to the coast of Prudhoe Bay (Fig. 13). Prior to 15 March, the observed ice was mainly landfast (i.e., ice velocity close to zero), except during a storm event in mid-January. The currents predicted by GIOPS agree well with observations and are consistent with reasonable attenuations of ηW in RunAIO (Fig. 13d). Around 15 March, during the passage of another storm, observations indicate an ice break-up event associated with strong ice and current velocities (up to 0.8 m s−1). These observed values are greatly underestimated (about 60 %) by GIOPS. After 15 March, the observed ice appears to keep drifting most of the time. This feature was poorly modeled by GIOPS, leading to systematically underestimated ice and current velocities and results in an over-attenuation of ηW in RunAIO from mid-March to mid-April (Fig. 13d).

Figure 14Snapshots of surface stress (τs) used in RunAO (left) and RunAIO (middle) and differences, ΔηW (RunAO-RunAIO), during two storm events on 15 March 2020 (a–c) and 2 April 2020 (d–f). Arrows in the left two columns show the wind stress vectors. Circles show the location of three tide gauges (clockwise from the bottom are station 2 at Alert, station 12 at Tuktoyaktuk, and station 13 at Prudhoe Bay).

In contrast, the successful attenuation of ηW predicted by RunAIO at Tuktoyaktuk and Alert throughout February–April suggests that sea ice over the two regions has much stronger resistance to large storms. To further illustrate the impact of ice, we examine the ice-induced changes in τs during two large storm events on 15 March 2020 at Tuktoyaktuk and 2 April 2020 at Alert (Fig. 14). In both cases, winds blow parallel to the coast (Fig. 14a and d), generating the Ekman setup associated with large ηW predicted by RunAO, up to 1.3 m at Tuktoyaktuk, and 0.4 m at Alert. In RunAIO, sea ice associated with strong internal stress greatly reduces τs. Reductions in τs in the Tuktoyaktuk region are concentrated closer to the coast, and in particular, τs is completely shut down at the coast by landfast ice. Reductions in τs in the Alert region reach 70 %–80 % for the entire storm. This leads to significant ηW attenuation, spanning about 1000 km along the coast for each storm (Fig. 14c and f), including the western Canadian Arctic (up to 1.0 m attenuation), northeastern Canadian Arctic, and northern Greenland (up to 0.5 m attenuation).

Finally, we note that the comparison of results from RunAO and RunAIO also shows large attenuations of ηW up to 1.0 m in the Russian Arctic. Although observations are not available, the estimation is expected to be reasonable considering the large volume of landfast ice in the ESS (see Fig. 3).

7 Summary and conclusions

The present study outlines, and evaluates, a novel approach for adding sea ice effects to a global TWL forecast model. The following two overriding questions are addressed: (1) can we design an efficient parameterization to include ice effects in a global ocean model for forecasting TWL and improve forecast skill in polar regions? (2) Can we isolate and explain the contribution of dominant physical processes (e.g., under-ice friction, baroclinicity, and nonlinear tide–surge interaction) to the seasonal modulation of tides?

The approach incorporates the total (tide plus surge) ice-induced ocean stress by taking advantage of already available external forecast fields (i.e., ice concentration, ice velocity, and surface ocean currents) produced by operational ice–ocean systems. The new method's novel feature is a consistent representation of the tidal relative ice–ocean velocity based on a transfer function derived from ice and ocean tidal ellipses given by external ice–ocean models. This effectively helps circumvent inconsistencies in tides among different models. The approach was applied to ECCC's high-resolution (1/12) global operational TWL forecast system. The external model is a coarser-resolution (1/4), data-assimilative global ice–ocean prediction system also running operationally at ECCC. Model predictions of TWL were generated for the period from November 2018 to April 2022, covering four ice seasons.

The impact of adding ice effects was quantified using observed hourly sea level at 58 tide gauges and moorings in ice-infested waters of the Northern Hemisphere. Adding ice effects is shown to help reproduce most of the observed seasonal modulations in the dominant M2 tide (up to 40 % in amplitude and 50 in phase) in the CS, Canadian Arctic, and HB. The observed interannual variability in the modulation (up to 10 % in amplitude and 10 in phase) during the four ice seasons is also captured by the model, with the addition of ice effects. Improvements are also found, mostly in the Arctic, for other smaller constituents (i.e., S2, K1, and O1; see the Supplement).

The dominant mechanism for seasonal modulations is the under-ice friction due to the presence of landfast or non-free drift ice. It is mostly relevant in areas of shallow waters (less than 100 m) with strong tidal dissipation. These areas cover most of the ESS and CS and parts of the Canadian Arctic and HB system. The under-ice friction leads to amplitude reductions and phase delays. In turn, they can drive large shifts of amphidromes (up to 125 km), resulting in opposite responses (i.e., amplitude enhancement and phase advance). Remote effects of overpredicted shifts help explain some of the discrepancies between observations and model predictions. In addition to the under-ice friction, important contributions from baroclinicity and TSI due to ice–ocean stress were found. The impact of TSI due to ice–ocean stress is found predominantly in March, in shallow areas (i.e., ESS and CS), where large ice cover and strong winter storms occur. In contrast, baroclinic effects are prominent in September, owing to the presence of a pycnocline. The baroclinic mechanism also affect amphidromes and, in particular, the relatively large amphidromes located in the Arctic Ocean. Both mechanisms, TSI and baroclinic effects, generally reinforce the seasonal modulations induced by under-ice friction.

Adding ice effects also greatly improves the model skill in predicting storm surges. This is achieved via ice-induced surge attenuation (up to 1.0 m) over regions (e.g., western and northeastern Canadian Arctic) associated with weak ice mobility and thus strong ice strength. The attenuation is due to considerable ice-induced reductions (70 %–100 %) of the surface stress. Large attenuations, up to 1.0 m, are also predicted along the coast of ESS associated with strong ice strength. Forecast challenges remain in regions with intermediate ice strength (i.e., the coast of northern Alaska), where the inclusion of ice effects leads to over-attenuated surges. Insufficient ice break-up during strong storms in the external model was shown to lead to underestimated ice and current velocities (about 60 %) compared to in situ measurements, and they carry across in the form of over-attenuated surges in the TWL system.

Over the next 100 years, climate change is expected to accelerate, causing a general reduction in ice cover in the Arctic (Pörtner et al.2022). Our results imply that the reduction in ice concentration and strength over areas of shallow waters, in particular, will dramatically increase tidal amplitudes over most coastal areas. As effects of ice and baroclinicity are expected to decrease and increase, respectively, in winter, tidal amphidromic systems will be pushed towards their ice-free states. The reduced ice cover is also expected to enhance the intensification of winter storms (Crawford et al.2022), contributing to higher storm surges. Future changes in both tides and storm surges thus pose an increasing risk of coastal flooding and erosion, particularly for coastal areas in the Canadian Arctic currently fully or partially protected by the ice cover.

In terms of future work, we plan to extend the present study to the Antarctic once ice cavities in the Ross Sea and Weddell Sea are included in the external ice–ocean model. It will also be interesting to investigate the dynamical mechanisms behind the observed counterintuitive M2 modulation in the Ross Sea, as reported by Ray et al. (2021).

Code availability

Source code of NEMO v3.6 and its configuration for this study can be accessed at (Wang2023). The original code was modified to include the new parameterization of the ice–ocean stress.

Data availability

Hourly tide gauge records collected from four institutes are publicly available at MEDS (; DFO2022), NOAA (; NOAA2022), UHSLC (; Caldwell et al.2015), and EMODnet (; EMODnet2022). Data obtained from two publications, including the monthly tidal constants in the Russian Arctic and bottom pressure records in the Hudson Bay, will be available upon request to the corresponding author. In situ measurements of ice and surface current velocities collected at an offshore mooring can be accessed via a single MATLAB file “CODAmoorings_masterfile.mat” at (Thomson et al.2021). The observed weekly fast ice extent can be accessed at (U.S. National Ice Center2020).


The supplement related to this article is available online at:

Author contributions

Both authors contributed to the conceptualization, methodology, and the writing of the paper. PW coded the parameterization of ice–ocean stress and performed the simulation and evaluation. NBB provided supervision and managed the project.

Competing interests

The contact author has declared that neither of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We would like to offer our gratitude and pay our deepest respects to our recently deceased long-time friend, colleague, and mentor Keith R. Thompson. It was a privilege to brainstorm and share the love of our science with him. He made valuable comments and suggestions in the development of the new parameterization. We thank him for that and much more. The authors also thank Mikhail E. Kulikov for providing the monthly tidal constants in the Russian Arctic, and Pierre St-Laurent for providing the bottom pressure records in the Hudson Bay. We also thank Mathieu Plante for his detailed comments on an early draft of this paper. Last, the authors would like to thank William Pringle and one anonymous reviewer for their constructive comments.

Review statement

This paper was edited by Taesam Lee and reviewed by William Pringle and one anonymous referee.


Bernier, N. B. and Thompson, K. R.: Tide-surge interaction off the east coast of Canada and northeastern United States, J. Geophys. Res.-Oceans, 112, C06008,, 2007. a, b

Bij de Vaate, I., Vasulkar, A., Slobbe, D., and Verlaan, M.: The Influence of Arctic Landfast Ice on Seasonal Modulation of the M2 Tide, J. Geophys. Res.-Oceans, 126, e2020JC016630,, 2021. a, b, c

Buehner, M., McTaggart-Cowan, R., Beaulne, A., Charette, C., Garand, L., Heilliette, S., Lapalme, E., Laroche, S., Macpherson, S. R., Morneau, J., and Zadra, A.: Implementation of deterministic weather forecasting systems based on ensemble–variational data assimilation at Environment Canada. Part I: The global system, Mon. Weather Rev., 143, 2532–2559, 2015. a

Caldwell, P. C., Merrifield, M. A., and Thompson, P. R.: Sea level measured by tide gauges from global oceans – the Joint Archive for Sea Level holdings (NCEI Accession 0019568), Version 5.5, NOAA National Centers for Environmental Information [data set], (last access: 3 October 2022), 2015. a

Cartwright, D. E. and Amin, M.: The variances of tidal harmonics, Deutsche Hydrografische Zeitschrift, 39, 235–253, 1986. a

Crawford, A. D., Lukovich, J. V., McCrystall, M. R., Stroeve, J. C., and Barber, D. G.: Reduced Sea Ice Enhances Intensification of Winter Storms over the Arctic Ocean, J. Climate, 35, 3353–3370, 2022. a

Danard, M., Rasmussen, M., Murty, T., Henry, R., Kowalik, Z., and Venkatesh, S.: Inclusion of ice cover in a storm surge model for the Beaufort Sea, Nat. Hazards, 2, 153–171, 1989. a

DFO: Marine Environmental Data Section Archive, Ecosystem and Oceans Science, Department of Fisheries and Oceans Canada [data set], (last access: 3 October 2022), 2022. a

Dunphy, M., Dupont, F., Hannah, C. G., and Greenberg, D.: Validation of a modelling system for tides in the Canadian Arctic Archipelago, Can. Tech. Rep. Hydrogr. Ocean Sci. 243: vi + 70 pp., Fisheries and Oceans Canada, Dartmouth, N.S., Canada, 2005. a

Egbert, G. D. and Erofeeva, S. Y.: Efficient inverse modeling of barotropic ocean tides, J. Atmos. Ocean. Tech., 19, 183–204, 2002. a

EMODnet: In situ near real time sea level data, EMODnet [data set], (last access: 3 October 2022), 2022. a

Fissel, D. and Tang, C.: Response of sea ice drift to wind forcing on the northeastern Newfoundland shelf, J. Geophys. Res.-Oceans, 96, 18397–18409, 1991. a

Gaspar, P., Grégoris, Y., and Lefevre, J.-M.: A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station Papa and Long-Term Upper Ocean Study site, J. Geophys. Res.-Oceans, 95, 16179–16193, 1990. a

Godin, G. and Barber, F.: Variability of the tide at some sites in the Canadian Arctic, Arctic, 33, 30–37, 1980. a

Heil, P. and Hibler, W. D.: Modeling the high-frequency component of Arctic sea ice drift and deformation, J. Phys. Oceanogr., 32, 3039–3057, 2002. a, b

Henry, R.: Storm surges, Beaufort Sea Project Technical Report No. 19, Department of the Environment, Victoria, Canada, 1975. a

Henry, R. and Foreman, M.: Numerical Model Studies of Semi-diurnal Tides in the Southern Beaufort Sea, Pacific Marine Science Report 77-11, Institute of Ocean Sciences, Victoria, Canada, 1977. a

Hibler, W., Roberts, A., Heil, P., Proshutinsky, A. Y., Simmons, H., and Lovick, J.: Modeling M2 tidal variability in Arctic sea-ice drift and deformation, Ann. Glaciol., 44, 418–428, 2006. a

Hošeková, L., Eidam, E., Panteleev, G., Rainville, L., Rogers, W. E., and Thomson, J.: Landfast ice and coastal wave exposure in northern Alaska, Geophys. Res. Lett., 48, e2021GL095103,, 2021. a

Hunke, E. C. and Dukowicz, J. K.: An elastic–viscous–plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, 1997. a

Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE: the Los Alamos Sea Ice Model Documentation and Software User's Manual Version 4.1 LA-CC-06-012, T-3 Fluid Dynamics Group, Los Alamos National Laboratory, Los Alamos, NM, USA, 675, 500, 2010. a, b

Johnson, M. and Eicken, H.: Estimating Arctic sea-ice freeze-up and break-up from the satellite record: A comparison of different approaches in the Chukchi and Beaufort Seas, Elementa: Science of the Anthropocene, 4, 000124,, 2016. a

Joyce, B. R., Pringle, W. J., Wirasaet, D., Westerink, J. J., Van der Westhuysen, A. J., Grumbine, R., and Feyen, J.: High resolution modeling of western Alaskan tides and storm surge under varying sea ice conditions, Ocean Model., 141, 101421, 2019. a, b

Kim, J., Murphy, E., Nistor, I., Ferguson, S., and Provan, M.: Numerical analysis of storm surges on Canada's western Arctic coastline, Journal of Marine Science and Engineering, 9, 326,, 2021. a

Kleptsova, O. and Pietrzak, J.: High resolution tidal model of Canadian Arctic Archipelago, Baffin and Hudson Bay, Ocean Model., 128, 15–47, 2018. a

Kodaira, T., Thompson, K. R., and Bernier, N. B.: The effect of density stratification on the prediction of global storm surges, Ocean Dynam., 66, 1733–1743,, 2016a. a

Kodaira, T., Thompson, K. R., and Bernier, N. B.: Prediction of M2 tidal surface currents by a global baroclinic ocean model and evaluation using observed drifter trajectories, J. Geophys. Res.-Oceans, 121, 6159–6183,, 2016b. a

Kowalik, Z.: Storm surges in the Beaufort and Chukchi seas, J. Geophys. Res.-Oceans, 89, 10570–10578, 1984. a

Kulikov, M., Medvedev, I., and Kondrin, A.: Seasonal variability of tides in the Arctic Seas, Russian Journal of Earth Sciences, 18, 1–14, 2018. a, b, c

Kulikov, M., Medvedev, I., and Kondrin, A.: Features of Seasonal Variability of Tidal Sea-level Oscillations in the Russian Arctic Seas, Russ. Meteorol. Hydrol., 45, 411–421, 2020. a

Lemieux, J.-F., Tremblay, L. B., Dupont, F., Plante, M., Smith, G. C., and Dumont, D.: A basal stress parameterization for modeling landfast ice, J. Geophys. Res.-Oceans, 120, 3157–3173, 2015. a

Lemieux, J.-F., Dupont, F., Blain, P., Roy, F., Smith, G. C., and Flato, G. M.: Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges, J. Geophys. Res.-Oceans, 121, 7354–7368, 2016. a

Lisitzin, E.: Sea-level changes, Elsevier, Amsterdam, ISBN 0-444-41157-7, 1974. a

Lu, P., Li, Z., Cheng, B., and Leppäranta, M.: A parameterization of the ice–ocean drag coefficient, J. Geophys. Res.-Oceans, 116, C07019,, 2011. a, b

Madec, G.: NEMO ocean engine, Note du Pôle de modélisation, de l'Institut Pierre-Simon Laplace (IPSL) No. 27, France, ISSN 1288-1619, 2008. a

McPhee, M.: Air–ice–ocean interaction: Turbulent ocean boundary layer exchange processes, Springer Science & Business Media, New York,, 2008. a

Müller, M., Cherniawsky, J. Y., Foreman, M. G., and von Storch, J.-S.: Seasonal variation of the M2 tide, Ocean Dynam., 64, 159–177, 2014. a, b, c, d

NOAA: Tides & water levels, NOAA [data set], (last access: 24 May 2022), 2022. a

Parkinson, C. L.: Arctic Sea Ice Coverage from 43 Years of Satellite Passive-Microwave Observations, Front. Remote Sens., 3, 1021781,, 2022. a

Pawlowicz, R., Beardsley, B., and Lentz, S.: Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE, Comput. Geosci., 28, 929–937,, 2002. a

Pörtner, H.-O., Roberts, D. C., Adams, H., Adler, C., Aldunce, P., Ali, E., Begum, R. A., Betts, R., Kerr, R. B., Biesbroek, R., and Birkmann, J.: Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK and New York, NY, USA, 2022. a, b

Prinsenberg, S.: Damping and phase advance of the tide in western Hudson Bay by the annual ice cover, J. Phys. Oceanogr., 18, 1744–1751, 1988. a, b, c

Ray, R. D.: On measurements of the tide at Churchill, Hudson Bay, Atmos. Ocean, 54, 108–116, 2016. a, b

Ray, R. D.: Technical note: On seasonal variability of the M2 tide, Ocean Sci., 18, 1073–1079,, 2022. a

Ray, R. D., Larson, K. M., and Haines, B. J.: New determinations of tides on the north-western Ross Ice Shelf, Antarct. Sci., 33, 89–102, 2021. a, b

Rienecker, M. and Teubner, M.: A note on frictional effects in Taylor's problem, J. Mar. Res., 38, 183–191, 1980. a

Roos, P. C. and Schuttelaars, H. M.: Influence of topography on tide propagation and amplification in semi-enclosed basins, Ocean Dynam., 61, 21–38, 2011. a

Rotermund, L. M., Williams, W. J., Klymak, J. M., Wu, Y., Scharien, R. K., and Haas, C.: The effect of sea ice on tidal propagation in the Kitikmeot Sea, Canadian Arctic Archipelago, J. Geophys. Res.-Oceans, 126, e2020JC016786,, 2021. a, b, c

Roy, F., Chevallier, M., Smith, G. C., Dupont, F., Garric, G., Lemieux, J.-F., Lu, Y., and Davidson, F.: Arctic sea ice and freshwater sensitivity to the treatment of the atmosphere–ice–ocean surface layer, J. Geophys. Res.-Oceans, 120, 4392–4417, 2015. a, b, c

Smith, G. C., Roy, F., Reszka, M., Surcel Colan, D., He, Z., Deacu, D., Belanger, J.-M., Skachko, S., Liu, Y., Dupont, F., and Lemieux, J.-F.: Sea ice forecast verification in the Canadian global ice ocean prediction system, Q. J. Roy. Meteor. Soc., 142, 659–671, 2016. a, b, c

Smith, G. C., Bélanger, J.-M., Roy, F., Pellerin, P., Ritchie, H., Onu, K., Roch, M., Zadra, A., Colan, D. S., Winter, B., and Fontecilla, J. S.: Impact of coupling with an ice–ocean model on global medium-range NWP forecast skill, Mon. Weather Rev., 146, 1157–1180, 2018. a, b, c, d

St-Laurent, P., Saucier, F., and Dumais, J.-F.: On the modification of tides in a seasonally ice-covered sea, J. Geophys. Res.-Oceans, 113, C11014,, 2008. a, b, c, d

Stepanov, V. N. and Hughes, C. W.: Parameterization of ocean self-attraction and loading in numerical models of the ocean circulation, J. Geophys. Res.-Oceans, 109, C03037,, 2004. a

Taylor, G. I.: Tidal oscillations in gulfs and rectangular basins, P. Lond. Math. Soc., 2, 148–181, 1922. a

Thomson, J., Eidam, E., and Hosekova, L.: Mooring results from Coastal Ocean Dynamics in the Arctic (CODA), ResearchWorks Archive [data set], (last access: 21 February 2023), 2021. a

Tranchant, B., Testut, C.-E., Renault, L., Ferry, N., Birol, F., and Brasseur, P.: Expected impact of the future SMOS and Aquarius Ocean surface salinity missions in the Mercator Ocean operational systems: New perspectives to monitor ocean circulation, Remote Sens. Environ., 112, 1476–1487, 2008. a

Tsamados, M., Feltham, D. L., Schroeder, D., Flocco, D., Farrell, S. L., Kurtz, N., Laxon, S. W., and Bacon, S.: Impact of variable atmospheric and oceanic form drag on simulations of Arctic sea ice, J. Phys. Oceanogr., 44, 1329–1353, 2014.  a

U.S. National Ice Center: U.S. National Ice Center Arctic and Antarctic Sea Ice Concentration and Climatologies in Gridded Format, Version 1, Tech. rep., NSIDC: National Snow and Ice Data Center, Boulder, Colorado USA [data set], (last access: 3 October 2022), 2020. a, b

Wang, P.: NEMO v3.6 with the addition of the parameterized ice-ocean stress, Zenodo [code],, 2023. a

Wang, P., Bernier, N. B., Thompson, K. R., and Kodaira, T.: Evaluation of a global total water level model in the presence of radiational S2 tide, Ocean Model., 168, 101893,, 2021. a, b, c, d, e, f, g

Wang, P., Bernier, N. B., and Thompson, K. R.: Adding baroclinicity to a global operational model for forecasting total water level: Approach and impact, Ocean Model., 174, 102031,, 2022. a, b, c, d, e, f

Weatherall, P., Marks, K. M., Jakobsson, M., Schmitt, T., Tani, S., Arndt, J. E., Rovere, M., Chayes, D., Ferrini, V., and Wigley, R.: A new digital bathymetric model of the world's oceans, Earth Space Sci., 2, 331–345,, 2015. a

Whalen, D., Forbes, D. L., Kostylev, V., Lim, M., Fraser, P., Nedimović, M. R., and Stuckey, S.: Mechanisms, volumetric assessment, and prognosis for rapid coastal erosion of Tuktoyaktuk Island, an important natural barrier for the harbour and community, Can. J. Earth Sci., 59, 945–960, 2022. a

Zhang, Z.-H. and Leppäranta, M.: Modeling the influence of ice on sea level variations in the Baltic Sea, Geophysica, 31, 31–46, 1995. a

Short summary
Effects of sea ice are typically neglected in operational flood forecast systems. In this work, we capture these effects via the addition of a parameterized ice–ocean stress. The parameterization takes advantage of forecast fields from an advanced ice–ocean model and features a novel, consistent representation of the tidal relative ice–ocean velocity. The new parameterization leads to improved forecasts of tides and storm surges in polar regions. Associated physical processes are discussed.