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

. In operational ﬂood forecast systems, the effect of sea ice is typically neglected or parameterized solely in terms of ice concentration. In this study, an efﬁcient 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 M 2 ) in the Arctic and Hud-son 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


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 Eicken, 2016;Parkinson, 2022), 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 M 2 tidal amplitude in many parts of the Canadian Arctic, including the Beaufort Sea and the Amundsen Gulf (up to 50 %; Henry and Foreman, 1977;Godin and Barber, 1980), the Kitikmeot Sea (50 %-60 %; Rotermund et al., 2021), and the Hudson Bay (HB) system (8 %-40 %; Prinsenberg, 1988;St-Laurent et al., 2008).Large variability in the M 2 amplitude was also reported in the Russian Arctic (Kulikov et al., 2018(Kulikov et al., , 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 M 2 phase was also observed, with up to 40 • in the CS (Kulikov et al., 2018) and eastern HB (Prinsenberg, 1988).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 M 2 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 M 2 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 M 2 phase delayed in most of the Arctic but advanced in the HB.
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 Thompson, 2007) 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 M 2 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 Pietrzak, 2018) or applied over landfast ice only (Bij de Vaate et al., 2021;Rotermund et al., 2021).These simple methods help produce the iceinduced 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 (Lisitzin, 1974) and Beaufort Sea (Henry, 1975).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 Tang, 1991;Heil and Hibler, 2002) 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 Hibler, 2002); 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(Wang et al., , 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., 2021Wang et al., , 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(Smith et al., , 2018;;Lemieux et al., 2015Lemieux et al., , 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.

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.

The ocean model
The NEMO modeling framework (Madec, 2008) 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 u h represents the horizontal velocity vector (u, v), u denotes the complete velocity vector in three dimensions (u, v, w), f is the Coriolis parameter, p a 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 selfattraction and loading (Stepanov and Hughes, 2004).The lateral eddy viscosity and diffusivity coefficients are set to constant (A h = 100 m 2 s −1 and K h = 10 m 2 s −1 ), and the vertical eddy viscosity and diffusivity coefficients (A z and K z ) 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 u h towards the observed current u obs calculated using the tidal amplitude and phase of eight major tidal constituents (M 2 , S 2 , N 2 , K 2 , O 1 , K 1 , P 1 , and Q 1 ) provided by TPXO8 (Egbert and Erofeeva, 2002).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  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 (T f and S f ) provided by a coarser-resolution (1/4 • ) data-assimilative system -the Global Ice Ocean Prediction System (GIOPS; Smith et al., 2016Smith et al., , 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 highfrequency variability to evolve freely (Wang et al., 2022).
At the surface, the boundary condition for Eq. ( 1) is given by where ρ a denotes the air density, and u 10 is the wind velocity at 10 m height.In the absence of sea ice (i.e., ice concen-tration α equals zero), the surface stress τ s equals the airocean stress τ ao .The air-ocean drag coefficient C ao equals 1.2 × 10 −3 for |u 10 | < 8 m s −1 and then increases linearly with |u 10 |, with a slope of 0.065 × 10 −3 for every 1 m s −1 increase in |u 10 | (Bernier and Thompson, 2007).Hourly fields of u 10 and p a 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 where τ b refers to the bottom stress, u b denotes the current velocity at the bottom, and C db 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.

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: τ io can be parameterized by a quadratic drag law in terms of the relative velocity between ice and surface currents (u ice − u surf ) as follows: where C io 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 elasticviscous-plastic solver for ice dynamics (Hunke and Dukowicz, 1997), 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 u ice − u surf 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.

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 iceocean 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: As u T ice is mainly forced by u T surf , we introduce a transfer function, such that where and the spatially varying scale factor a T and the rotation angle ϕ can be inferred from u T * ice and u T * surf provided by external ice-ocean models (the asterisk * denotes quantities from external models).Specifically, a T 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 where I is the identity matrix.Unlike periodic tides, storm surges are sporadic and occur on a local scale driven by the atmospheric forcing.Applying a transfer function to u S surf , similar to Eq. ( 10), is not feasible, since u S surf is forced by u S ice .Instead, we expect that the https://doi.org/10.5194/gmd-16-3335-2023Geosci.Model Dev., 16, 3335-3354, 2023  2020).Middle and bottom panels show the derived monthly a T and ϕ for the M 2 tide for the same period.Note that for areas with very weak u T * ice (major axis velocity magnitude less than 5 × 10 −3 m s −1 in this study), ϕ is irrelevant, its estimation is also not reliable, and so it is set to 0. 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, a S , can be used to adjust the residual relative velocities to the surface levels of the TWL model as follows: For example, if the surface level in the external model is shallower than the TWL model, then u S * surf will be stronger than u S surf , and a S should be larger than unity.In practice, a S can be tuned to best reproduce the observed residual water level in the presence of ice.
Finally, the total relative velocity is In practice, u T surf = u surf 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).

Ice-ocean stress
Gridded fields of hourly ice concentration (α), ice velocity (u T * ice ; u S * ice ), and surface current (u T * surf ; u S * surf ) were obtained from the assimilation component of GIOPS (Smith et al., 2016(Smith et al., , 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 NEMObased 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 u T * ice and u T * surf , 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 (M 2 , S 2 , K 1 , and O 1 ) in GIOPS, as they can be adequately resolved with monthly harmonic analyses.
Figure 3 (middle panels) illustrates the monthly estimates of a T for M 2 from December 2020 to March 2021.Note that a T captures the mobility of sea ice, with landfast for a T = 0, non-free drift for 0 < a T < 1, and free drift for a T = 1.In theory, only landfast ice and non-free drift ice exert friction on the underlying ocean flow.Regions of a T → 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 M 2 (not shown), although over some regions (e.g., ESS and CS) K 1 and O 1 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 M 2 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.

Ice-ocean drag coefficient
The parameterization of the ice-ocean drag coefficient C io is a complex issue, as C io 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 C io 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 C io 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 C io 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 C io to 1.00 × 10 −2 , which produces reasonable agreement between observed and predicted seasonal M 2 modulations.
The empirical scale factor for residual relative velocity, a S 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 (a S ) 2 C io for the residual stress based on u S * ice − u S * surf 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.

Experimental design and analysis
Two basic runs, Run AO and Run AIO , were conducted to examine the impact of adding ice effects on predicted water levels.In Run AO , ice effects are not considered, and τ s equals τ ao (see Eq. 5).In Run AIO , τ 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 Run AIO , including the baroclinic effects (by using constant T and S; Run1), TSI due to τ io (Run2), TSI due to τ b (Run3), underice friction (by setting u T ice − u T surf to zero; Run4).Specifically, for Run2, TSI due to τ io was removed by setting , where u T io and u S io are, respectively, the tidal and residual relative ice-ocean velocities given by Eqs. ( 12 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. https://doi.org/10.5194/gmd-16-3335-2023 Geosci.Model Dev., 16, 3335-3354, 2023 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 Amin, 1986).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 K 2 (P 1 ) was inferred from S 2 (K 1 ), 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 (A Sept , φ Sept ), when sea ice has the minimum cover in the Arctic.Thus, where the subscript i denotes a particular month.The iceinduced 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; Ray, 2022).In the main text, we focus on the seasonality of the dominant M 2 constituent, which has negligible astronomical contribution.Other minor constituents, including S 2 , K 1 , and O 1 , are briefly discussed in the Supplement.

Total water level
We first compare the model skill for Run AO and Run AIO 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.3cm).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 (Ray, 2016).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.

Seasonal variability
Figure 5 shows the M 2 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 Run AO and Run AIO (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 M 2 modulations in Run AO , 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 (S 2 , K 1 , and O 1 ; see the Supplement for additional details).For M 2 , one exception is station 18, where observations show an anomalous positive ÃMar (Fig. 5e), while Run AIO 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 Run AIO is roughly parallel to that coast.We note, however, that Run AIO 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.
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 M 2 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 (Run AO ), 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 (Run AIO ), forecasts of Ãi and φ i are greatly improved at most stations across the season.
We next examine the interannual variability in the seasonal M 2 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 Run AO , while they are reasonably captured by Run AIO .One exception is Churchill (station 39), where observations show almost no modulations in amplitude, while Run AIO generates 5 %-10 % modulations.We note, as before, that observed tides at Churchill may have quality or drift issues (Ray, 2016).

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 M 2 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 M 2 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 (Taylor, 1922) 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 Teubner, 1980;Prinsenberg, 1988;Roos and Schuttelaars, 2011).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 M 2 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.
With ice conditions changing rapidly in the Arctic, it is important to know where the presence of ice or the underice 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 Run AIO .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 https://doi.org/10.5194/gmd-16-3335-2023 Geosci.Model Dev., 16, 3335-3354, 2023   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).
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.

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 M 2 (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 M 2 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.

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 highfrequency 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 highfrequency contributions, we applied a low-pass filter with a https://doi.org/10.5194/gmd-16-3335-2023 Geosci.Model Dev., 16, 3335-3354, 2023  cutoff period of 10 h to tidal residuals to obtain η S .We further decompose η S into a wind and pressure-driven component so that 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 barometeronly prediction (i.e., a run driven with surface air pressures only).We then removed the predicted η P from both observed 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.
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 Run AIO (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 Run AIO from mid-March to mid-April (Fig. 13d).
In contrast, the successful attenuation of η W predicted by Run AIO 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 Run AO , up to 1.3 m at Tuktoyaktuk, and 0.4 m at Alert.In Run AIO , 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 Run AO and Run AIO 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). https://doi.org/10.5194/gmd-16-3335-2023 Geosci.Model Dev., 16, 3335-3354, 2023

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) iceinduced 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 highresolution (1/12 • ) global operational TWL forecast system.The external model is a coarser-resolution (1/4 • ), dataassimilative 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 M 2 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., S 2 , K 1 , and O 1 ; see the Supplement).
The dominant mechanism for seasonal modulations is the under-ice friction due to the presence of landfast or nonfree 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 M 2 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 https://doi.org/10.5281/zenodo.7662916(Wang, 2023).The original code was modified to include the new parameterization of the ice-ocean stress.Caldwell et al., 2015), and EMODnet (https://emodnet.ec.europa.eu/geoviewer/; EMODnet, 2022).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 "CO-DAmoorings_masterfile.mat" at http://hdl.handle.net/1773/47139(Thomson et al., 2021).The observed weekly fast ice extent can be accessed at https://doi.org/10.7265/46cc-3952(U.S. National Ice Center, 2020).
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.
Disclaimer.Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

Figure 3 .
Figure 3. Top 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 a T and ϕ for the M 2 tide for the same period.Note that for areas with very weak u T * ice (major axis velocity magnitude less than 5 × 10 −3 m s −1 in this study), ϕ is irrelevant, its estimation is also not reliable, and so it is set to 0.
) and (13).For Run3, TSI due to τ b was removed by setting τ b = ρ 0 C db (|u T b |u T b + |u S b |u S b ), where u T b and u S b are isolated using the online tidal filter of Wang et al. (2021), with u T b = u b ; u S b = u b − u b .The six runs are summarized in Table

Figure 4 .
Figure 4.The rms of observations (black) and RMSE for Run AO (blue) and Run AIO (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 5 .
Figure 5. Modulation of the M 2 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 Run AO (a, b) and Run AIO (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.

Figure 6 .
Figure 6.Normalized monthly M 2 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 Run AO (blue lines) and Run AIO (red lines).The title of each subplot gives the station code, station name, and averaged M 2 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 7 .
Figure 7. Same as Fig. 6 but for the monthly M 2 phase anomaly ( φ i , • ) relative to September.

Figure 8 .
Figure 8. Monthly normalized M 2 amplitude anomaly ( Ãi ; left panels) and phase anomaly ( φ i ; right panels) relative to September 2019 at five permanent gauges observed (black), Run AO (blue), and Run AIO (red) for the period from November 2018 to April 2022.

Figure 9 .
Figure 9. (a) M 2 amphidromes in March predicted by Run AO , with amplitude (m) in the contour map and phase (every 30 • ) shown with white lines.(b) Ice-induced shifts of M 2 amphidromes in March, taken as the difference in the predicted amphidromic points between Run AIO and Run AO .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).

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

Figure 13 .
Figure 13.Time 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.

Figure 14 .
Figure 14.Snapshots of surface stress (τ s ) used in Run AO (left) and Run AIO (middle) and differences, η W (Run AO -Run AIO ), 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).

Table 1 .
Summary of water level observations collected from various institutes and publications (see Fig.1for station code and Fig.2for 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).