Impact of internal wave drag on Arctic sea ice

Abstract

A parameterization of the impact of internal waves on momentum transfer at the sea-ice–ocean interface based on previous work by McPhee has been implemented in a sea-ice model for the first time. The ice–ocean drag from internal waves is relevant for shallow mixed layer depth and the presence of a density jump at the pycnocline and is also a function of the strength of the stratification beneath the ocean mixed layer and geometry of the ice interface. We present results from a coupled sea-ice–ocean model where the parameterization of internal wave drag has been implemented. We conducted simulations spanning the years from 2000 to 2017. We find a deceleration of ice drift by 5–8% in both winter and summer, but with significant spatial and temporal variation reaching seasonal average values of ~10%. The spatial variation of ice transport leads to local impacts on deformed ice of magnitude ~0.05 m (2–5%), and reductions in ocean-to-ice heat fluxes of ~1 W m⁻², and a decrease in bottom melt of ~0.02–0.04 cm d⁻¹. There is an increase of up to 15% in thickness and ice concentration in the Canadian Arctic and a 10% overall impact on the total sea-ice volume.

1. Introduction

Internal waves (IWs) generated by an object moving over the stratified ocean have long been studied because of their effect on ship's movement: sailors have observed ships stalling in relatively calm conditions. This phenomenon, known as ‘dead water’, was observed and named by Nansen in 1893 (Nansen, 1902). During an oceanographic expedition on the polar schooner Fram in the Nordenskiöld Archipelago, Nansen observed that the speed of the ship would decrease to a fifth, even with the engines at full power. Nansen stated ‘The Fram appeared to be held back, as if by some mysterious force, and she did not always answer the helm. We made loops in our course, turned sometimes right around, tried all sorts of antics to get clear of it, but to very little purpose.’ The hull of the Fram was exciting IWs along the ocean pycnocline that were acting as a source of drag.

Several experiments have been carried out to investigate the dead water effect. The first was conducted by Ekman in 1904, who simulated the dead water phenomenon by filling a tank with salty water topped by a layer of fresh water and towing a model ship over it. He found that the drag generated by IWs followed a quadratic law depending on the fluid's density profile, the ship's hull size and the towing speed (Medjdoub and others, 2020). More recent experiments confirm Ekman's main conclusions (Fourdrinoy and others, 2020; Medjdoub and others, 2020), and that, while IWs can become nonlinear, linear analytical calculations well simulate the behaviour of the model ship in a stratified fluid.

A naturally occurring analogue of an object moving over a stratified fluid is sea ice moving across a stratified ocean: large sea-ice protrusions into the mixed layer such as keels can generate IWs that propagate momentum away from the ice and act to increase the drag coefficient at the ice–ocean interface (Dosser and Rainville, 2016; Cole and others, 2018; Kawaguchi and others, 2019).

The drag induced by IWs is additional to that from form drag, which is caused by redirection of fluid flow around an irregularly shaped ice–ocean interface. It has been demonstrated that the strength of the drag component associated with IWs depends upon the depth and spacing of under-ice topographic features, the magnitude of the relative velocity between the ice and ocean currents at the pycnocline depth, the density jump across the pycnocline, the depth of the mixed layer and the strength of the stratification of the underlying ocean (McPhee, 1987; McPhee and Kantha, 1989). For stratified flows, laboratory experiments have demonstrated that the ice–ocean drag induced by keels has a prominent impact with respect to the overall surface stress (Pite and others, 1995).

An increase of the magnitude of the ice–ocean drag by IWs has the potential to affect ocean and ice velocity and the spatial distribution (thickness, area fraction) of the sea-ice cover, with broad impacts on climate and weather. Here, we implement an IW drag parameterization based on McPhee's theory (McPhee and Kantha, 1989), into a local (CPOM) version of the sea-ice model CICE, accounting for its impact on the momentum balance, and couple this to the ocean model NEMO (Nucleus for European Modelling of the Ocean, www.nemo-ocean.eu/). Using a series of numerical simulations of our sea-ice–ocean model, we study the impact of IW drag on sea-ice drift, extent and thickness, on basin and local scales in the Arctic Ocean.

Our paper describes the first implementation of ice–ocean IW drag into a climate model set up, by modifying code in the CICE and NEMO models, to include a parameterization of a known physical process whose impact on the sea-ice state was unknown. Inclusion of a representation of this fundamental process into these climate model components and exploring its impacts on sea-ice state is the main goal of our study.

The paper is structured as follows. In section 2 we briefly present linear theory describing the impact of IWs on ice–ocean drag and describe how we have implemented this into the CICE sea-ice model. In section 3 we present numerical simulations evaluating the impact of IW drag on the motion of sea ice and its mass balance. In section 4 we discuss our results, and present concluding remarks in section 5.

2. Model development to bring internal wave drag into the NEMO-CICE sea-ice model

2.1 Ice–ocean drag generated by internal waves

The total ice–ocean drag is given by the sum of parasitic drag and the IW drag. Parasitic drag is defined as the sum of form drag, related to flow around larger obstacles such as keels or floe edges, and skin drag, due to the interaction between the ocean and the surface roughness (‘skin’) of the ice. IW drag is produced by the motion of large sea-ice roughness elements, i.e. keels, over the stratified ocean; the keels generate IWs in the ocean that radiate momentum away from the ice–ocean interface.

IWs, a fundamental component of ocean dynamics, are gravity waves that propagate within stratified fluid layers beneath the ocean surface (Gill, 1982). These waves arise from the restoring force of buoyancy acting on perturbations in the density field, leading to oscillations in the vertical displacement of isopycnals. Mathematically, the governing equations for IWs can be described by the nonlinear, dispersive and rotating Navier–Stokes equations, coupled with the equation of state and continuity equation. The dispersion relation for IWs in a stratified fluid is given by the Taylor–Goldstein equation (LeBlond and Mysak, 2014), which relates the wave frequency, horizontal wavenumber and vertical structure of the wave. The energy flux associated with IWs can be quantified using the energy equation, considering the conversion between potential and kinetic energy as the waves propagate. Flow past keels cause a perturbation to the pycnocline even if the keels do not penetrate the pycnocline, resulting in an IW both at the pycnocline and in the stratified fluid below it. The radiation of momentum caused by these waves is the source of the IW drag.

In the polar regions, the presence of sea ice creates the conditions for the origination of IWs under the ice excited by the sea-ice bottom topography. The presence of strong pycnoclines observed during summer melt together with the presence of keels under the ice depicts the typical ‘dead water’ scenario, as described by Nansen in 1902. Shallow mixed layer depths lead to higher speeds of IWs between the pycnocline surface and the overlying ice cover (Guthrie and others, 2013). This increase in wave velocity would consequently cause greater dissipation or loss of wave energy at that boundary (Morison and others, 1985).

A description of a model of linear IWs generated by the sea-ice bottom topography, and their impact on ice–ocean drag, is presented by McPhee (1987), McPhee and Kantha (1989). Here we summarize the main points pertinent to our study. The under-ice surface is characterized by the depth of ice features H _k (keels) protruding into the ocean relative to the surrounding level ice, and the wavenumber k ₀ characteristic of the feature spacing (k ₀ = 2π/L where L is the spacing between the topographic features, see Fig. 1a). The assumption of linearity is that k ₀ H _k ≪ 1, i.e. the features are of small aspect ratio compared to their distance, and allows the development of analytical expressions for the IW drag (McPhee and Kantha, 1989). An illustrative sketch is shown in Figure 1a, where the variables involved in the calculations of the IW drag are labelled. McPhee showed that the IW drag coefficient C _IW is given by

(1)$$C_{{\rm IW}} = {\rm \Gamma }C_{{\rm DNW}}, \;$$

where C _DNW is the drag coefficient that considers the ocean stratified up to the ice–ocean interface (no mixed layer below the ice), while Γ is a dampening factor, related to the buoyancy jump across the pycnocline at the bottom of the mixed layer and the mixed layer depth. The drag coefficient C _DNW is given by

(2)$$\matrix{ {C_{{\rm DNW}} = \displaystyle{{\nu {( {k_{\rm c}H_{\rm k}} ) }^2} \over 2}\sqrt {1-\nu ^2} } \cr {0 < \nu = \displaystyle{{k_0} \over {k_{\rm c}}} < 1} \cr } , \;$$

where the critical wavenumber k _c is defined as k _c = N/u ₀, where N is the buoyancy frequency and u ₀ is the speed of the ocean current relative to the sea ice directly below the pycnocline. If k ₀ > k _c then the IWs generated are evanescent and do not contribute to drag. The dampening factor Γ has values between 0 and 1, but it can reach values over the unity for some combination of the input variables (McPhee and Kantha, 1989). The factor Γ has the effect of reducing the IW drag as the pycnocline density jump increases and is given by

(3)$$\Gamma = \left\{{1 + \left({\displaystyle{1 \over {v^2}} + R_{\rm b}^2 } \right)\sin h^2( k_0H) -R_{\rm b}\sin h( 2k_0H) } \right\}^{{-}1}, \;$$

were H is the depth of the mixed layer,

(4)$$R_{{\rm b} } = \displaystyle{{\Delta b} \over {k_0u_0^2 }}$$

is the Richardson number, and the reduced gravity is given by

(5)$$\Delta b = \displaystyle{{\Delta \rho } \over \rho }g.$$

Fig. 1. (a) Schematic of sea-ice moving relative to the ocean indicating variables used to calculate internal wave drag, (b) dependence of the internal wave drag coefficient (C _IW) on mixed layer depth and ν = k ₀/k _c keeping Δb = 0.002 m s⁻², (c) dependence of C _IW on mixed layer depth and Δb keeping ν = 0.2, and (d) dependence of C _IW on Δb and ν keeping the mixed layer depth at 10 m.

The value of Γ ranges between 0 and 1 and is larger for small values of Δb. McPhee and Kantha (1989) describes the interplay between the variables involved in the calculation of the IW drag coefficient.

Figure 1 illustrates how the IW drag coefficient C _IW depends upon the mixed layer depth, the k ₀/k _c ratio (ν) and reduced gravity Δb. To do so, we keep in turn one parameter fixed and study the behaviour of C _IW with respect to the others. We set the velocity shear at 0.05 m s⁻¹ and the size of the bottom topography feature at 7 m. Figure 1a lists graphically the variables used in the IW drag parameterization. Figure 1b shows the relationship between C _IW and the mixed layer depth and k ₀/k _c for Δb = 0.002 m s⁻², corresponding to a salinity jump of 0.25 PSU (McPhee and Kantha, 1989); Figure 1c shows the relationship between C _IW and k ₀/k _c and Δb for a mixed layer depth of 10 m; and Figure 1d shows the relationship between C _IW and mixed layer depth and Δb for k ₀/k _c = 0.2. The IW drag increases by a factor of ten for mixed layers shallower than 10 m, peaks at ν ~ 0.7, and is greater for smaller density jumps at the pycnocline interface. Larger density jumps have been observed in the Arctic more recently (Toole and others, 2010), but the complex interactions between all the inputs to the IW drag do not lead to the conclusion that larger density jump will undermine the importance of the inclusion of this parameterization into the CICE model. Conditions under which IW drag may become notable frequently arise in the central Arctic during summer (McPhee and Kantha, 1989), where mixed layers between 5 and 7 m have been observed.

Keels moving in shallow mixed layer create a larger displacement of the stratified ocean beneath the mixed layer, therefore creating more powerful IWs and consequently a larger IW drag. This is due to the presence of a layer of cold and relatively fresh water resulting from sea-ice melt.

Furthermore, the existence of a density jump implies, with suitable forcing (i.e. relative flow past keels), an IW whether there is or is not stratification beneath the mixed layer. A stronger pycnocline implies a weaker coupling of the mixed layer from the ocean below the mixed layer that leads to weaker IWs. Smaller density jumps at the pycnocline cause stronger IW drag.

Stronger/weaker stratification below the mixed layer implies smaller/larger amplitude of the IWs; this has a nonlinear and not very easily interpretable impact on the drag, calculated from Eqns (1–5). As can be seen from Eqn (2), weaker stratification below the ML (N approaches zero) results in a decrease of the IW drag.

2.2 Adjustments of drag formulation in CICE to include internal wave drag

The IW drag described above has been brought into the CPOM implementation of the sea-ice model CICE version 5.1 (Hunke and others, 2015), which is coupled to the ocean model NEMO version 3.6 (Nucleus for European Modelling of the Ocean) (Madec, 2008). The NEMO 3.6 model is modified to pass variables to CICE needed to calculate the IW drag. Other than the changes related to IW drag, we use the configuration of CICE 5.1.2 described in Stroeve and others (2018) and the global ocean GO6 configuration of NEMO 3.6 described in Storkey and others (2018). The CICE configuration mirrors that outlined in Schröder and others (2014), a predictive melt pond model (Flocco and others, 2012), and an elastic anisotropic-plastic rheology (Tsamados and others, 2014). CICE 5.1.2 includes a form drag parameterization (Tsamados and others, 2014).

The sea-ice–ocean simulations are performed on a 1-degree tripolar grid (~40 km grid resolution for the Arctic Ocean) over the period 2000–2017 using NCEP-DOE-2 Reanalyses data (Kanamitsu and others, 2002, updated, 2017) as atmospheric forcing. The IW drag is calculated in CICE using variables calculated in the form drag routine (Tsamados and others, 2014) and several variables calculated in NEMO and passed to CICE. Here we describe details of our implementation.

In the absence of IW drag, the ice–ocean drag coefficient C _DW is given (Tsamados and others, 2014) by

(6)$$C_{{\rm DW}} = C_{{\rm dwf}} + C_{{\rm dwk}} + C_{{\rm dws}}, $$

where C _dwf is the form drag contribution from floe edge draft, C _dwk is the form drag contribution from keels, and C _dws is the bottom surface skin drag contribution. The relative importance of the drag contributions varies spatially and temporally with form drag typically dominating. In their numerical simulations, Tsamados and others (2014) found that during winter the main contribution to the total ice–ocean drag coefficient came from form drag from keels in the heavily ridged regions in the Lincoln Sea and north of the Arctic Canadian Archipelago, while in the summer, as the sea-ice concentration drops, the contribution from floe edges becomes meaningful.

We have calculated the IW drag contribution from floe edges to be, for typical values of floe draft (~1 m), an order of magnitude smaller than that from keels (~7 m depth).

To account for IW drag, an extra term is added to the total ice–ocean drag C _DW, so that Eqn (6) becomes

(7)$$C_{{\rm DW}} = C_{{\rm dwf}} + C_{{\rm dwk}} + C_{{\rm dws}} + C_{{\rm IW}}, \;$$

where C _IW is given by Eqn (1).

The ice state parameters entering the equation for C _IW are the keel depth H _k and wavenumber, determined by keel spacing D _k. The keel depth and spacing are diagnosed in the form drag parameterization (Tsamados and others, 2014) from the CICE tracers for the deformed sea-ice volume v _ardg and area a _ardg, with the assumption of triangular keels and sails with constant angles of repose (α _r = α _k = 22°) and a constant ratio between the keels and sails depths, D _s, and spacing, H _s.

The slope angle of the keel has been set to 22° as in Tsamados and others (2014), an average value taken from several observational studies (Wadhams and Doble, 2008; Wadhams and Toberg, 2012; Ekeberg and others, 2015; Salganik and others, 2023).

Following Tsamados and others (2014), we set H _k/H _s = R _h = 4 and D _k/D _s = R _d = 1 to obtain

(8)$$D_{\rm k} = 2R_{\rm d}\displaystyle{{v_{{\rm rdg}}} \over {a_{{\rm rdg}}}}\displaystyle{{\alpha \tan ( \alpha _{\rm k}) R_{\rm d} + \beta \tan ( \alpha _{\rm r}) R_{\rm h}} \over {\phi _{\rm s}\tan ( \alpha _{\rm k}) R_{\rm d} + \phi _{\rm k}\tan ( \alpha _{\rm r}) R_{\rm h}^2 }}$$

and

(9)$$H_{\rm k} = 2R_{\rm h}H_{\rm s}\displaystyle{{a_{\rm i}} \over {a_{{\rm rdg}}}}\left({\displaystyle{\alpha \over {\tan ( \alpha_{\rm r}) }} + \displaystyle{\beta \over {\tan ( \alpha_{\rm k}) }}\displaystyle{{R_{\rm h}} \over {R_{\rm d}}}} \right), \;$$

where ϕ _s and ϕ _k are the sail and keels porosity and have a constant value of 0.8, while α and β are weight functions that define the amount of ridged area assumed to be covered by sails and keels. We use the default values of α  = 0 and β  = 0.75 following Tsamados and others (2014). Setting α  = 0 implies that the ridged area is determined by the area of keels; the choice of β  = 0.75, along with some other parameter choices, resulted in realistic sail/keel spacing and height/depth when compared to field observations. Further details are in Tsamados and others (2014). It follows that in the absence of ridged ice (therefore of keels), the IW drag is zero.

In particular, the parameters present in the internal drag formula are calculated in the ocean model and passed to CICE. In the NEMO ocean model, the reduced gravity Δb is calculated using the ocean density jump across the pycnocline (one grid spacing above minus one grid spacing below). Mixed layer depth is not a prognostic variable in NEMO (or most ocean models). We have used the common approach of defining the mixed layer depth by the depth of the pycnocline, which is diagnosed from a change in relative density gradient of 0.01 m^–1. Shallower than 10 m, limited spatial and temporal resolution can result in this density gradient depth being quite variable leading to inconsistent pycnocline depths. As a consequence, in our IW model we only use pycnocline/mixed layer depths of 10 m and deeper. The buoyancy frequency is calculated at the pycnocline.

There are regions where the keel depth H _k exceeds the mixed layer depth diagnosed from NEMO. Observations show that the mixed layer deepens around drifting keels (Mercier and others, 2011, Morison and Goldberg, 2012; Grue and others, 2016), and so in these cases we set the depth of the mixed layer to be equal to 110% of the keel depth (the modified value of the mixed layer depth was only used for the IW drag calculation and did not affect the rest of the simulation results). In our simulations, the keels only exceed the ML depth in winter and in limited areas (not shown). Scouring can occur when keels are deep enough to reach the ocean floor (Barnes and Reimnitz, 1997): in these cases, we set the keel depth equal to 99% of the ocean depth.

Finally, the calculation of the heat fluxes at the ice–ocean interface is described in the following parameterization:

(10)$$Q_{{\rm ocean \ndash ice}} = - cp_{{\rm ocn}}\rho _{\rm w}CH_{{\rm ocn}}{\rm \Delta }Tu^\ast , \;$$

where cp _ocn (3980 J kg⁻¹) is the specific heat for sea water near freezing, ρ _w is the sea water density, CH _ocn is the ice–ocean heat transfer coefficient, ΔT is the ice–ocean temperature difference and u* is defined by

(11)$$u^\ast{ = } \sqrt {\displaystyle{{str_{{\rm ocn}}} \over {\rho _{\rm w}}}} , \;$$

where str _ocn is the ice–ocean stress (Maykut and McPhee, 1995). We use CH _OCN equal to the default ice–ocean drag coefficient (not altered by the IW drag). The IW drag will indirectly affect the heat transfer rates by altering the ice–ocean speed difference/shear.

3. Results

3.1 Sea-ice–ocean (CICE-NEMO) model set up and reference simulation

We are interested in the seasonal and regional impact of the new IW drag parameterization on sea-ice concentration, thickness and motion. We perform atmosphere-forced simulations of the sea-ice and ocean system using the CPOM implementation of CICE 5.1 coupled to the ocean model NEMO 3.6. Here we give essential details of the model configuration, pertinent to the present study; the full description of the ocean modelling component can be found, e.g. in Kelly and others (2018) and Storkey and others (2018). For the purpose of this study, it is key to resolve surface and subsurface stratification and pycnocline. The model has 75 vertical levels, with resolution from 1 m at the surface to under 2 m at 10 m depth, and to 204 m at 6000 m (there are 31 model levels in the upper 200 m) (Kelly and others, 2018). This suffices for resolving upper pycnocline through the year. Partial steps in the model bottom bathymetry allow for improving model approximation of the steep seabed relief near the continental shelves (e.g. Bacon and others, 2015; Hordoir and others, 2022).

The CPOM CICE model activates form drag (Tsamados and others, 2014), and the IW drag parameterization described in the previous section.

We present results from a reference run without IW drag (REF run hereafter), and a run including IW drag (hereafter referred to as the IW run).

The scope of our paper is the assessment of the impact that the new implementation has on the sea-ice state variables, not to achieve the perfect sea-ice state forecast, but we still need to ensure that our coupled model simulation produce realistic results. Figure 2a shows an ice thickness climatology calculated using our REF configuration to be compared with mean ice thickness calculated from the Pan-Arctic Ice-Ocean Modelling and Assimilation System (PIOMAS), and Figure 2b shows a simulated September sea-ice extent to be compared with that retrieved using SSMI data using the Bootstrap algorithm (Comiso, 2017). PIOMAS is a reanalysis dataset that produces daily and monthly estimates of Arctic sea-ice thickness using a coupled ice and ocean model from January 1979 to almost real-time, developed by the Polar Science Centre (PSC) of the University of Washington's Applied Physics Laboratory (APL), responsible for the model's development and upkeep.

Fig. 2. The top panel shows the thickness climatology in our NEMO-CICE reference run REF in comparison to PIOMAS between 2000 and 2017, while the lower panel shows the mean September sea-ice extent comparison derived from SSMI sea-ice concentration data and our REF run from 2000 to 2017.

Our CICE simulations underestimate the September sea-ice extent in some years (e.g. ~2 × 10⁶ km² in September 2003), and slightly overestimate the annual cycle of sea-ice volume. However, our REF run is realistic enough to explore the sensitivity and to assess the impact of our new parameterization of IW drag.

3.2 Impact of internal wave drag on the total ice–ocean drag coefficient

A climatology of the total sea-ice–ocean drag is shown in Figure 3 together with the IW drag coefficient. The IW run shows higher values throughout the year, with the IW drag coefficient presenting values in the order of 10⁻⁴ and maxima during summer. The maximum calculated value of the IW drag is actually observed in September, due to the shallower mixed-layer and smaller density jump in summer, but this is less representative than the impact in July and August, as the ice is at its minimum in September. The difference between the sea-ice–ocean drag between the two runs can be considered always present during the year.

Fig. 3. Climatology of the ice–ocean drag C _DW showing the impact of the IW drag (C _IW) throughout the year from 2000 to 2017. In yellow the IW drag climatological value.

Climatological maps averaged over the 18 years of simulation of the total ice–ocean drag and IW drag coefficient are shown in Figures 4a and b, together with the relative impact in percentage of IW drag over the total ice–ocean drag (Fig. 4c). The IW drag ranges between 0 and 1.5 × 10⁻³, with an overall increase of the total ice–ocean drag of about 10%, with lower values in the Russian Arctic and peaks of up to 20% in the Canadian Arctic.

Fig. 4. The top three panels (a, b, c) show the 2000–2017 climatology of the total ice–ocean drag coefficient, the internal wave drag coefficient, and the IW impact on the total drag coefficient as a percentage (IW drag/total ice–ocean drag). The bottom panels (d, f) show the climatology of mixed layer depth, density jump at the mixed layer depth and the keel depth.

The strongest IW drag occurs in the Fram Strait, north of Queen Elizabeth Island and in the western the Beaufort Sea (orange in Figs 4b and c). These areas are characterized by a shallow mixed layer (green in Fig. 4d), and by the presence of a density jump (blue/green in Fig. 4e) and by large keel depths, in the order of ~15 m (orange in Fig. 4f). The maps show the average of these variables over the 18 years of simulation (e.g. the average mixed layer depth is calculated by averaging over the 18 years of simulation the mixed layer variable diagnosed by the NEMO-CICE, 5 d output). The pattern of the IW drag coefficient observed suggests a link with the Transpolar Drift. IW drag has a stronger impact for slowly moving ice, since fast moving ice would exceed the IW's wake, dispersing energy locally. The Transpolar Drift has accelerated, leading to faster ice drift by up to 6 km d⁻¹ compared to the mean of the previous decade (Dethloff and others, 2022). Furthermore, sea ice is thinner due to the fast export related to the Transpolar Drift, compared to the Canadian Arctic, where thicker multi-year ice is still present; therefore, we expect the observed sharp decrease of the IW drag component in the Central Arctic compared to the Canadian Arctic.

While the principal interplay of the factors controlling IW drag can be seen from the annual climatology, more insight is gained by considering seasonal climatologies. We present results for January and August climatologies in Figures 5 and 6, representing winter and summer conditions for the main sea-ice state and oceanic variables impacted by our parameterization. In winter, the ice–ocean drag presents lower values than in summer. In winter (January climatological average), the IW drag shows values of up to 30% of the total drag along Greenland and in the Kara Sea, while in the Central Arctic we mostly see values ranging between 15 and 20% of the total drag (Figs 5a–c). Figures 5d and e show the mixed layer depth and the density jump in the IW run, and the keel depth: we can see that increases in the mixed layer depth correspond to low values of IW drag as expected. In Figures 5g–i, we plot k ₀, k _c and k ₀/k _c to show the order of magnitude that these parameters present in our simulations. As expected, the regions where the strongest impact of IW drag is observed correspond to those with values of 0.6 < k ₀/k _c < 0.8.

Fig. 5. The top three panels (a, b, c) show the 2000–2017 January climatology of the total ice–ocean drag coefficient, the internal wave drag coefficient, and the IW impact on the total drag coefficient as a percentage. The middle panels (d, e, f) show the January climatology of mixed layer depth, density jump at the mixed layer depth and the keel depth. The bottom panels (g, h, i) show the January climatological values for k ₀, k _c and k ₀/k _c.

Fig. 6. The top three panels (a, b, c) show the 2000–2017 August climatology of the total ice–ocean drag coefficient, the internal wave drag coefficient and the IW impact on the total drag coefficient as a percentage. The middle panels (d, e, f) show the climatology of mixed layer depth, density jump at the mixed layer depth and the keel depth.

In summer (August climatological average), the total ice–ocean drag coefficient is larger than in winter therefore even though the IW drag is larger in absolute value, its contribution to the total sea-ice–ocean drag is smaller in percentage than in winter (Figs 6a and b), with an impact of up to 15% in the Canadian Arctic (Fig. 6c). In Figures 6d and e, we can see a different situation compared to winter: the mixed layer is shallower (acting to strengthen IW drag) and the small density jump at the pycnocline overcomes this with an overall increase of the impact of the IW drag in the Canadian Arctic, where also the depth of the keels is larger.

Note that our results set a lower bar on the impact of IW on drag and should be considered conservative, considering that our coupled NEMO-CICE model does not allow mixed layers shallower than 10 m.

3.3 Impact of IW drag on sea-ice motion

Figures 7a, c and e show respectively the annual average, the winter (January) and summer (August) climatologies of ice drift between 2000 and 2017 in the IW run. Figures 7b, d and f show the corresponding differences with the REF run (no IW drag). Figure 7b shows that IW drag leads to an overall slowdown of the ice drift, in the range of 0.05 m s⁻¹, the equivalent of ~5% slowdown in speed with a 10% peak in the Russian Arctic. We see a larger impact of IW drag on ice drift in winter (Fig. 7d) relative to the annual average. The magnitude of the decrease in ice drift due to IW drag is less pronounced and is more varied in summer (Fig. 7f), with the strongest impact in the Canadian Arctic and along Greenland. Where the IW drag is negligible (Fig. 6b), we also observe a negligible difference in ice drift anomaly (Fig. 7f).

Fig. 7. Maps of 2000–2017 climatology of sea-ice drift. The left panels (a, d, g) show values from our IW run for the annual average, January and August. The middle panels (b, e, h) show the respective differences between the IW run and the reference REF run. The right panels (c, f, i) show the climatological differences between the REF run and the Pathfinder's observations.

To estimate the improvement in the sea-ice state representation, we compared our sea-ice drift with available Pathfinder observations (Tschudi and others, 2019), from 2000 to 2017. The coupled NEMO-CICE REF run we used overestimates the ice drift overall by values in the order of 0.1 m s⁻¹ (Figs 7c, f and i). By implementing the IW drag, we observe a slowdown in the order of 0.01 m s⁻¹ (see middle panels of Fig. 7), which correspond to an improvement in the ice drift bias of around 10%. It is interesting to notice that in regions where the REF run shows stronger overestimates of the ice drift, e.g. in the Fram Strait and the Beaufort Sea, the IW run shows larger differences (in particular a larger slowdown), compared to the REF run, therefore leading to a larger compensation of the discrepancy between model results and observations.

3.4 Impact of IW drag on ice concentration and thickness

The area averaged ice concentration in the IW run compared to the REF run is shown in Figure 8. When looking at the areal changes we only observe a small increase of annual mean ice concentration (Figs 8a and b), but we can see meaningful local increases in the summer months ranging between 5 and 20%, particularly in the Russian Arctic and along the coast in the Beaufort Sea: this can be seen in the August climatology difference (Fig. 8d). It is interesting to highlight the localized increase in sea-ice concentration in the Kara Sea in August, related to areas of higher values of IW drag observed in Figures 6b and c.

Fig. 8. Maps of sea-ice concentration climatology: the left panels show absolute values for the annual average and summer (a and c); the right panels show the respective differences between the IW run and the reference run (b–d).

Figure 9a shows the annually averaged ice thickness climatology in our IW run. Figure 9b reveals that inclusion of IW drag leads to spatially variable changes in ice thickness with an overall small increase in thickness. Figures 9d and f show that inclusion of IW drag leads to a pronounced increase by up to 10 cm thickness in the Canadian Arctic and north of Greenland during winter and summer.

Fig. 9. Maps of sea-ice thickness climatology: the left panels show absolute values for the annual average, January and August (a, c and e); the right panels show the respective differences between the IW run and the reference run (b, d and f).

4. Discussion

The implementation of IW drag in a coupled sea-ice–ocean model (NEMO-CICE) has seasonal and regional effects on the sea-ice thickness, concentration and motion. As it would be expected, the increase in the ice–ocean drag mainly impacts the dynamics of sea ice without noticeably altering the climatological sea-ice concentration, nonetheless affecting the sea-ice thickness. The impact of IW drag we simulate is likely an underestimate since in our parameterization only minimum mixed layer depths of 10 m are considered, whereas shallower pycnoclines have been observed in the Central Arctic during summer (McPhee and Kantha, 1989; Randelhoff and others, 2014; Rynders, 2017), and these are expected to be associated with increased IW drag as we have shown in our results. Toole and others (2010) presented data showing summer mixed layers of 16 m (overestimated, as they mention in their paper), and 24 m in winter, values in agreement with our simulation results.

Observations on bottom sea-ice topography is not easy to obtain, even though statistically available thanks to the combination of numerous datasets (Martin, 2007). Further observations on bottom topography would be desirable to improve its characterization and therefore its depiction in sea-ice models.

Our IW drag parameterization has the potential to reduce known biases in sea-ice state variables in climate models. CMIP6 models are known to overestimate ice drift by up to a factor of 2 in the Arctic compared to observations (Wang and others, 2023), and the inclusion of the IW drag parameterization decreases the ice drift in our simulations.

When comparing our new simulation results with Pathfinder observations, we can clearly see an improvement in sea-ice drift difference compared to observations of about 10%, with the REF run overestimating the observations by values in the order of 0.1 m s⁻¹, compensated by a reduction of about 1 cm s⁻¹ when implementing the IW drag parameterization (see details in Fig. 7).

It is also important to underline that the form drag parameterization produces a larger drag coefficient where the ice is more heavily ridged, and IW drag is also largest in these regions. Our study contributes to bring wider awareness of IW drag to the climate modelling community. The IW drag coefficient is also critically affected by mixed layer depth, unlike form drag, and the seasonal variation of this cannot be accounted for by tuning of the form drag scheme.

The location of the increased sea-ice concentration and thickness in our simulations acts to reduce the observed bias in sea-ice models producing too thin ice in the Beaufort Sea (Ridley and others, 2018). Though in practise the level of uncertainty in observations, e.g. 10% in sea-ice concentration rising to 30% in the summer seasonal ice cover (e.g. Comiso, 2017) and 0.5 m rising to 33–40% error in summer monthly mean sea-ice thickness (Landy and others, 2022), is too great to reliably distinguish the impact of IW drag in large-scale model-observation comparisons.

Our simulation is run globally, therefore we could assess the impact of the IW drag in the Antarctic, which shows a similar impact during summer, while is found not to be as strong as during the Antarctic winter. We suggest that the reason for the observed winter difference is due to the deeper mixed layer observed in the Antarctic.

4.1 Impact of IW drag on thermodynamic and dynamic thickness tendencies

The total ice–ocean drag is dependent on the IW drag, directly via the ice–ocean drag coefficient C _DW, and indirectly through the slowdown of the ice drift relative to the ocean. There is an impact on ice thickness and to gain insight into the causes of this we consider the dynamic or thermodynamic origin of the changes. To do so, the difference between IW and REF in thermodynamic (melt/freeze) and dynamic (convergence and advection) tendency variables are analysed here. This assessment is presented in Figure 10, focusing on winter (January, panels 10a–e) and summer (August, panels 10f–j) climatologies. In panels 10a and 10f we show the changes in volume of ridged ice, since the keel's distribution is retrieved from the volume of ridged ice in the CICE model. The ridged ice volume changes are in the order of 5 cm and present different patterns for summer and winter, with a decrease in Beaufort Sea and along the Canadian coast in winter, and increased values in the Eurasian Basin with maxima along the path of the Transpolar Drift. In summer, the pattern shows a decrease of ridged ice volume in the Beaufort Sea and an increase in the north of Greenland.

Fig. 10. Maps of differences of IW run minus REF run climatologies of volume of ridged ice, dynamic and thermodynamic tendencies bottom melt and ocean-to-ice heat fluxes for winter (top panels, a–e), and summer (bottom panels, f–j).

The dynamic tendency difference presents a spatial variability with both positive and negative anomaly values, varying between −0.2 and 0.2 cm d⁻¹ both in winter and summer (respectively Fig. 10b and 10g), indicating the measure of the impact of the IW drag implementation on the sea-ice dynamics. The dynamic tendency reacts in a variable way to the new implementation throughout the year, while the thermodynamic tendency shows a more homogeneous behaviour with mostly positive anomalies during summer (~0.04 cm d⁻¹, Fig. 10h) and small negative anomalies during winter (~−0.02 cm d⁻¹, Fig. 10c). The bottom melt decreases by 0.05 cm d⁻¹ on average, with a decrease in ocean-to-ice heat flux in the order of 1 W m⁻² (Fig. 10i), which is the equivalent of up to 10% of the observed value of ocean to ice heat fluxes (Lin and Zhao, 2019), and corresponding to 0.05 cm d⁻¹ of sea-ice melt/freeze (Fig. 10j).

Since the IW drag coefficient does not alter the ice–ocean heat transfer coefficient we can infer that any change in heat flux from ice to ocean is related to the consequences of the increased drag, namely to the slower ice drift. Increased IW drag acts to equalize ice drift and surface ocean current, reducing the ice to ocean heat flux and reducing sea-ice melt at the bottom of the ice.

4.2 Case study: April 2002

The annual, January and August climatologies we have presented above, necessarily average over interannual variability. To illustrate the typical spatial patterns of impact of IW drag on sea ice we present results for April 2002. We have chosen April 2002 because this year/month does not present any particular climatic behaviour that could overshadow or exalt the influence of the IW drag. In Figures 11a and b we present the total ocean drag and the IW drag, while in Figure 11c we can see that the impact of the IW drag in most of the Russian Arctic reaches values of up to 30% of the total ice–ocean drag. The pattern of high IW drag along the north of Greenland is clearly visible in the total drag coefficient in Figure 11c.

Fig. 11. Case study April 2002: maps of ice–ocean drag coefficient (a), internal wave drag coefficient (b) and percentage of impact of IW drag over total drag (c).

While the impact of IW drag on sea-ice concentration (not shown) is not substantial (<10%), it leads to sea-ice thickness increases by up to 25 cm in the Central Arctic (Figs 12a and b). IW drag tends to reduce ice drift, by up to −0.01 m s⁻¹, which corresponds to ~7–8% of the total ice drift (Figs 12c and d). Slower, less mobile ice is more likely to dwell in the Arctic without moving southward. By dwelling in the Arctic cold waters, it is also less likely to melt quickly, inducing a positive feedback mechanism.

Fig. 12. Case study April 2002: maps of sea-ice thickness (12a), ice drift (12c) for the IW run with the corresponding differences between the IW and the REF run (b, d). Panels e and f show respectively the difference in dynamical tendency and thermodynamic tendency between the IW run and the REF run.

As expected, the dynamic and thermodynamic tendencies show larger values compared to the ones shown in the climatologies, while following the same behaviour with very variable values in the dynamical tendency patterns and an overall decrease in the thermodynamic tendency with values reaching 0.1 cm d⁻¹ (Figs 12e and f).

4.3 Monthly climatology of IW impact

We here show the annual behaviour of the sea-ice volume during the period of our simulation for each month, in order to give a quantitative assessment of the impact of the new parameterization implemented in the CICE model.

From Figure 13 we observe a noticeable seasonality of the IW drag parameterization throughout the simulation time.

Fig. 13. Monthly climatology of the Pan-Arctic integrated sea-ice volume difference between the IW run and the REF run.

The Pan-Arctic climatology of volume anomaly (the volume difference between the REF and the IW) presents a positive trend in the first three years of our simulations in every month that can be attributed to model spin-up. It also shows a marked seasonality with the highest values in summer months (i.e. July, August, September), of 100–550 km³ representing a 10% impact on the total sea-ice volume calculated in our simulations. This is consistent with our previous results, which show larger climatological values of the IW drag coefficient during summer (see Fig. 3), and a greater spatial impact in the summer climatology (Figs 6b and c).

The positive summer anomaly is consistent with the seasonality of the IW drag, which is known to be stronger in low sea-ice concentration and shallower mixed layer depths. Further work that falls out of the scope of this study would be needed to understand the possible relationships between the observed anomalies in sea-ice volume and the Arctic sea-ice minima, in particular focussing on the spatial distribution of the observed IW impact.

5. Concluding remarks

IWs form because of the interaction between flow past the ice bottom topography and the stratified ocean underneath. These interactions create an additional drag (IW drag C _IW), which increases in the presence of shallow pycnoclines with a density jumps across them. IW drag is responsible for the ‘dead water’ effect first noted by Nansen (1902).

Our work describes the first implementation of ice–ocean IW drag into a climate model set up, by modifying code in the CICE and NEMO models, to include a parameterization of a known physical process whose impact on the sea-ice state was unknown. The theory of IW drag suggests a priori the potential for a material impact on the sea-ice state, but this depends upon a complex, nonlinear interplay of various factors including keel depth and spacing, mixed layer depth, pycnocline strength, upper ocean stability and shear under the ice. To quantify the impact of IW drag, we included a parameterization of IWs drag at the ice–ocean interface into a sea-ice–ocean model (NEMO – CICE), and here present results of 18 years of simulation with and without IW drag. Our results show that including IW drag results in an overall increase of the ice–ocean drag coefficient, presenting regional stronger effects in the Canadian Arctic, where the ice–ocean drag coefficient is impacted by the IW drag in the range of 10–20% in the climatological average, and by up to 30% in the month of April 2002, used in this work as a case study. We also find an overall decrease in sea-ice drift between 5 and 8% with stronger impact in winter and a sea-ice concentration increase by up to 10% in the Central Arctic and in the Russian Arctic during the summer months. This effect improves the model's results by about 10% when comparing our results with Pathfinder observations of sea-ice drift. We also observe an increase in ice thickness associated with inclusion of IW drag, by up to 20 cm during winter in the Canadian Arctic: this thickening is due to the redistribution of sea ice due to emerging new dynamics linked to a slower ice drift which allows more ridging (anomalies in sea-ice ridged area in the order of 5 cm). The sea-ice thickness increase is also a consequence of the observed decrease in sea-ice basal melt, which is due to the smaller ice-to-ocean heat flux (~1 W m⁻²). The new IW drag parameterization leads to an average increase in the Pan-Arctic July sea-ice volume of 350 km³, with a maximum of 550 km³.

There are considerable uncertainties about the future evolution of the Arctic Ocean and its sea-ice cover, but climate projections such as those in the 6th Coupled Model Intercomparison Project (CMIP6) consistently project further loss of Arctic sea ice in the coming decades (Notz and SIMIP Community, 2020). The impact of altered thermohaline and mechanical forcing on the upper ocean is currently unknown but will almost certainly alter the spatial and temporal evolution of the mixed layer depth, strength of pycnocline and upper ocean stratification. The IW drag parameterization presented here is capable of dealing with all such future conditions; we might expect IW drag to become more (or less) important depending on, for example, a decrease in local mixed layer depth or strengthening of the pycnocline. Current climate models use constant neutral drag coefficients for ice–atmosphere and ice–ocean drag, but future model generations are expected to include form drag parameterizations (e.g. Lüpkes and others, 2012, Tsamados and others, 2014; Lüpkes and Gryanik 2015; Whalen and others, 2020; Sterlin and others, 2023). Our results show that IW drag has a sensible impact on ice state and should be included in such ice–ocean drag parameterizations.