Evaluation of Arctic sea ice drift and its dependency on near-surface wind and sea ice conditions in the coupled regional climate model HIRHAM–NAOSIM

We examine the simulated Arctic sea ice drift speed for the period 2003–2014 in the coupled Arctic regional climate model HIRHAM–NAOSIM 2.0. In particular, we evaluate the dependency of the drift speed on the near-surface wind speed and sea ice conditions. Considering the seasonal cycle of the Arctic basin averaged drift speed, the model reproduces the summer–autumn drift speed well but significantly overestimates the winter–spring drift speed, compared to satellite-derived observations. Also, the model does not capture the observed seasonal phase lag between drift and wind speed, but the simulated drift speed is more in phase with the near-surface wind. The model calculates a realistic negative correlation between drift speed and ice thickness and between drift speed and ice concentration during summer–autumn when the ice concentration is relatively low, but the correlation is weaker than observed. A daily grid-scale diagnostic indicates that the model reproduces the observed positive correlation between drift and wind speed. The strongest impact of wind changes on drift speed occurs for high and moderate wind speeds, with a low impact for rather calm conditions. The correlation under low-wind conditions is overestimated in the simulations compared to observation/reanalysis data. A sensitivity experiment demonstrates the significant effects of sea ice form drag from floe edges included by an improved parameterization of the transfer coefficients for momentum and heat over sea ice. However, this does not improve the agreement of the modeled drift speed / wind speed ratio with observations based on reanalysis data for wind and remote sensing data for sea ice drift. An improvement might be achieved by tuning parameters that are not well established by observations.

1 Introduction

Arctic sea ice has experienced a rapid decrease in recent decades (e.g., Serreze and Stroeve, 2015; Stroeve et al., 2012). On the one hand, the observed Arctic sea ice decline is caused by thermodynamic processes, such as increasing air and ocean temperatures, e.g., due to ocean heat transports into the Arctic (Serreze et al., 2007). On the other hand, dynamical processes, such as a changed sea ice drift in response to changing wind, ocean currents, and sea ice conditions (e.g., reduction of sea ice concentration and thickness), play an important role in the redistribution of sea ice (Rampal et al., 2011; Serreze et al., 2007; Spreen et al., 2011). Except for the fast-ice and shear zones, sea ice moves largely in response to local winds and ocean currents. More than 70 % of the variance of the central Arctic pack ice motion is explained by the geostrophic wind alone on the timescale of days to months (Thorndike and Colony, 1982). For compact sea ice in the inner Arctic or sea ice near the coasts, the internal friction can be as large as the forces due to the winds and ocean currents (Leppäranta, 2011). The seasonal Arctic basin-wide sea ice drift speed, however, is additionally correlated with the sea ice concentration and thickness: sea ice drift speed decreases with increasing ice concentration when ice concentration is low, and sea ice drift speed decreases with increasing ice thickness when ice concentration is high (Docquier et al., 2017; Olason and Notz, 2014). In order to understand and project Arctic sea ice changes, it is vital for climate models that they can realistically capture the observed sea ice drift and its dependency on the atmospheric and oceanic forcing and on the sea ice conditions at different timescales. Tandon et al. (2018) showed that only a few CMIP5 models capture the observed seasonal cycle of sea ice drift speed and that some models show a seasonal cycle of sea ice drift in phase with the near-surface wind speed.

The first aim of this paper is to evaluate the simulated sea ice drift speed and its relation with sea ice concentration, thickness, and near-surface wind speed in the coupled Arctic regional climate model HIRHAM–NAOSIM 2.0 (Dorn et al., 2019). We first evaluate the simulated Arctic basin-wide monthly mean drift, then we evaluate the relationship between sea ice drift speed and sea ice conditions/wind speed both on the Arctic basin-wide, multiyear monthly mean scale and on the daily grid scale.

The second aim of this paper is to explore the sensitivity of the simulated ice drift to the parameterization of the atmospheric near-surface transfer coefficients for heat and momentum over sea ice, in which the effect of sea ice form drag is included. The air–ice drag controls the sea ice drift under the influence of the wind forcing. It can be separated into the frictional skin drag due to microscale roughness elements on the sea ice surface and the form drag caused by large structures like pressure ice ridges and floe edges (Arya, 1973, 1975). Most climate models account for the sea ice skin drag and consider the form drag only by tuning the spatially and temporally invariant air–ice drag coefficient (Castellani et al., 2018; Tsamados et al., 2014). This approach is poorly constrained by the observations and fails to describe the variability of the air–ice drag processes correctly (Andreas et al., 2010; Lüpkes et al., 2012a; Lüpkes and Gryanik, 2015; Tsamados et al., 2014).

Several model studies that include the sea ice form drag were carried out (e.g., Castellani et al., 2018; Renfrew et al., 2019; Tsamados et al., 2014). Tsamados et al. (2014) implemented a complex sea ice form drag parameterization based on many sea ice cover properties (e.g., sea ice concentration, vertical extent and area of ridges, freeboard and floe draft, and the size of floes and melt ponds) into the stand-alone sea ice model CICE. Castellani et al. (2018) implemented a simpler sea ice form drag parameterization that only relies on sea ice deformation energy and concentration into the coupled ocean–sea ice model MITgcm. Both studies showed improvement in sea ice drift after the form drag had been included. Recently, Renfrew et al. (2019) implemented an observation-based parameterization of atmospheric form drag caused by floe edges based on suggestions by Lüpkes et al. (2012a), Lüpkes and Gryanik (2015), and Elvidge et al. (2016) into a stand-alone atmosphere model. The simulation results show an improved agreement of mean atmospheric variables and turbulent fluxes with measurements in cold-air outbreak situations over the Fram Strait when form drag is included.

The physical parameterizations suggested by Lüpkes and Gryanik (2015) were constrained by summertime observations over the sea ice pack and by aircraft observations over the marginal ice zone (MIZ) during winter. Later the parameterizations were once more validated using a larger and independent set of aircraft data obtained from campaigns during different seasons (Elvidge et al., 2016). This validation work concerned the momentum fluxes. The assumptions of Lüpkes and Gryanik (2015) about heat and moisture fluxes over the MIZ could not yet be evaluated by measurements. Thus, further research is necessary on this issue.

Despite the potential uncertainties in terms of the heat and moisture fluxes, the strategy of the present study is to investigate the performance of the parametrization by Lüpkes and Gryanik (2015) without further modification in a coupled atmosphere–ocean–sea ice model. The parametrization by Lüpkes and Gryanik (2015) contains constants that could be fine-tuned in the future to eventually improve the model results. We leave this optimization, however, to future work since it requires a large number of model runs and extensive analysis work that is beyond the scope of this paper focusing on the principal impact of sea ice form drag on the model results.

The organization of this paper is as follows: Sect. 2 describes the model simulations and the observational data as well as the analysis methods used. In Sect. 3, we present the simulated long-term average Arctic sea ice drift speed and its dependency on near-surface wind and sea ice conditions based on both the multiyear, Arctic basin-wide scale and the daily grid scale. In Sect. 4, we discuss the sea ice form drag impact on the atmospheric near-surface turbulent fluxes and on the sea ice drift speed. Finally, a summary and conclusions are provided in Sect. 5.

2 Data and analysis

2.1 Model and simulations

2.1.1 Coupled atmosphere–ice–ocean regional model

The model used in this study is HIRHAM–NAOSIM 2.0 (Dorn et al., 2019), which consists of the regional atmospheric model HIRHAM5 and the regional ocean–sea ice model NAOSIM (North Atlantic/Arctic Ocean Sea Ice Model). HIRHAM5 is based on the numerical weather forecast model HIRLAM 7.0 (Undén et al., 2002) and applies the physical parameterizations of the atmospheric general circulation model ECHAM 5.4.00 (Roeckner et al., 2003). Köberle and Gerdes (2003) give a basic description of NAOSIM, while Fieg et al. (2010) describe NAOSIM's fine-resolution model version, which is used in HIRHAM–NAOSIM 2.0. NAOSIM's ocean component is based on the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model, MOM-2 (Pacanowski, 1996). The sea ice component is based on the dynamic–thermodynamic sea ice model described by Harder et al. (1998). The internal ice stress is described by a viscous-plastic rheology according to Hibler (1979). Thermodynamic snow-ice processes are handled using the zero-layer approach by Semtner (1976). Detailed information about HIRHAM–NAOSIM 2.0 is given by Dorn et al. (2019). The model is applied over a circum-Arctic domain at a horizontal resolution of 1∕4^∘ (∼27 km) in the atmosphere (HIRHAM5) and 1∕12^∘ (∼9 km) in the ocean (NAOSIM). The ocean–sea ice domain corresponds to the domain shown in Fig. 1.

Figure 1 Mean spatial pattern of sea ice drift speed (km d⁻¹) in the model (ensemble mean of HIRHAM–NAOSIM 2.0) for the 2003–2014 (a) winter (DJFM) and (c) summer (JJAS). Panels (b) and (d) are the model differences to the observation (model − KIMURA) for winter and summer respectively. The purple line in each panel indicates the study domain used for the basin-wide analysis.

2.1.2 Long-term simulations of the base configuration

A 10-member ensemble of multidecadal climate simulations for the period 1979 to 2016 was carried out by Dorn et al. (2019) with the base configuration of HIRHAM–NAOSIM 2.0. These ensemble simulations represent the basis for the present study and are referred to as BASE hereafter. The individual ensemble members used the same atmospheric initialization but applied different ice–ocean initial conditions, which were taken from 1 January of the last 10 years of a preceding 22-year-long coupled spin-up run for the period 1979–2000 (see Dorn et al., 2019, for more details). The regional model simulations were driven by ERA-Interim reanalysis data (Dee et al., 2011; referred to as ERA-I hereafter). ERA-I provided 6-hourly lateral atmospheric boundary conditions as well as daily surface boundary conditions where required (outside the coupling domain). The lateral ocean boundary conditions in the northern North Atlantic were taken from the Levitus climatology (Levitus and Boyer, 1994; Levitus et al., 1994).

2.1.3 Sensitivity experiments

To investigate the model's sensitivity to the momentum and heat transfer coefficients which explicitly include sea ice form drag contributions, both a control run (CTRL) and a sensitivity experiment (SENS) for 1 year (year 2007) were performed. CTRL applies the default atmospheric boundary layer parameterization of ECHAM 5.4 (Roeckner et al., 2003) with air–ice momentum and heat transfer coefficients depending only on atmospheric stability. The sea ice form drag effect is not included. SENS includes the improved parameterization of air–ice momentum and heat transfer coefficients adapted based on Lüpkes and Gryanik (2015). There, the transfer coefficients depend on sea ice concentration and include both the sea ice skin drag and form drag effects caused by the edges of ice floes in the marginal sea ice zone and of leads in the inner Arctic. The parameterization could also account for the effect by melt pond edges, but this additional effect would require the knowledge of the pond fraction on top of the ice, which is not available in HIRHAM–NAOSIM 2.0. The model applies a flux averaging method, which means that the total flux over a surface covered with sea ice of concentration A and with open water (1−A) is the concentration-weighted mean of both contributions. We describe here only the fluxes over sea ice, in which form drag is included in SENS, and refer to Roeckner et al. (2003) for the parameterization of fluxes over open water. In CTRL (as in BASE), however, the air–ice momentum transfer coefficient C_d,i and the heat transfer coefficient C_h,i include only the sea ice skin drag effect. They are calculated as

$$\begin{array}{}\text{(1)}& {\displaystyle }{C}_{\mathrm{d},\mathrm{i}}=C_{\mathrm{dn},\mathrm{i}}{f}_{\mathrm{m},\mathrm{i}},\text{(2)}& {\displaystyle }{C}_{\mathrm{h},\mathrm{i}}=C_{\mathrm{hn},\mathrm{i}}{f}_{\mathrm{h},\mathrm{i}},\end{array}$$

where C_dn,i (C_hn,i) is the drag (heat transfer) coefficient under neutral atmospheric stratification over ice and f_m,i (f_h,i) is the stability correction function over ice to adjust C_dn,i (C_hn,i) based on atmosphere stability. C_dn,i is calculated as

$$\begin{array}{}\text{(3)}& C_{\mathrm{dn},\mathrm{i}}={\displaystyle \frac{k^{\mathrm{2}}}{{\left[\ln \left(z_{\mathrm{L}}/z_{\mathrm{0},\mathrm{i}}+\mathrm{1}\right)\right]}^{\mathrm{2}}}},\end{array}$$

where k=0.4 represents the von Karman constant, z_L is the height of the lowest atmospheric model level, and z_0,i represents the skin drag roughness length over sea ice. C_hn,i is calculated as

$$\begin{array}{}\text{(4)}& C_{\mathrm{hn},\mathrm{i}}={\displaystyle \frac{k^{\mathrm{2}}}{\ln \left(z_{\mathrm{L}}/z_{\mathrm{0},\mathrm{i}}+\mathrm{1}\right)\ln \left(z_{\mathrm{L}}/z_{\mathrm{t},\mathrm{i}}+\mathrm{1}\right)}},\end{array}$$

where z_t,i represents the scalar roughness length over ice. Equations (1) to (4) are common descriptions of air–ice momentum and heat transfer coefficients except that “+1” was added to both $z_{\mathrm{L}}/z_{\mathrm{0},\mathrm{i}}$ and $z_{\mathrm{L}}/z_{\mathrm{t},\mathrm{i}}$ in Eqs. (3) and (4). This is done in the model to prevent the argument of the logarithm from going to zero, for which C_d,i (C_h,i) would go to infinity (see also Giorgetta et al., 2013).

In SENS, the new momentum transfer coefficient ${\widehat{C}}_{\mathrm{d},\mathrm{i}}$ and the new heat transfer coefficient ${\widehat{C}}_{\mathrm{h},\mathrm{i}}$ over ice, which both include skin drag and form drag effects, are calculated as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\widehat{C}}_{\mathrm{d},\mathrm{i}}=C_{\mathrm{dn},\mathrm{i}}{f}_{\mathrm{m},\mathrm{i}}+C_{\mathrm{dn},\mathrm{f}}\left[f_{\mathrm{m},\mathrm{i}}A+f_{\mathrm{m},\mathrm{w}}\left(\mathrm{1}-A\right)\right]/A,\text{(6)}& {\displaystyle }{\widehat{C}}_{\mathrm{h},\mathrm{i}}=C_{\mathrm{hn},\mathrm{i}}{f}_{\mathrm{m},\mathrm{i}}+C_{\mathrm{hn},\mathrm{f}}\left[f_{\mathrm{h},\mathrm{i}}A+f_{\mathrm{h},\mathrm{w}}\left(\mathrm{1}-A\right)\right]/A,\end{array}$$

where C_dn,f (C_hn,f) represents the form drag (heat transfer) coefficient related to neutral conditions over ice, f_m,w (f_h,w) is the stability correction function over water to adjust C_dn,f and C_hn,f to atmospheric stability. Equation (5) is obtained by combining the Eqs. (6), (52), and (70) by Lüpkes and Gryanik (2015). Equations (6) is obtained by combining the Eqs. (9), (64), and (74) by Lüpkes and Gryanik (2015). After adding “+1” to both 10∕z₀ and z_L∕z₀ and replacing z₀ with z_0,f in Eq. (65) by Lüpkes and Gryanik (2015), C_dn,f is calculated as

$$\begin{array}{}\text{(7)}& C_{\mathrm{dn},\mathrm{f}}=C_{\mathrm{e}\mathrm{10}}{\left[{\displaystyle \frac{\ln \left(\mathrm{10}/z_{\mathrm{0},\mathrm{f}}+\mathrm{1}\right)}{\ln \left(z_{\mathrm{L}}/z_{\mathrm{0},\mathrm{f}}+\mathrm{1}\right)}}\right]}^{\mathrm{2}}A{\left(\mathrm{1}-A\right)}^{\mathit{\beta }},\end{array}$$

where C_e10 represents the effective resistance coefficient (both the aerodynamic resistance coefficient of individual floes and the shape factor), z_0,f represents the (form) roughness length, and β is a constant exponent describing the dependence of the crosswind dimension of a floe on A. The values of C_e10, z_0,f, and β are $\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{3}}$, $\mathrm{0.57}\times {\mathrm{10}}^{-\mathrm{3}}$ m, and 1.1 respectively. The value of C_e10 is the average given in Eqs. (48) and (49) by Lüpkes and Gryanik (2015). The value of z_0,f is an average resulting from measured roughness lengths by various campaigns considered by Andreas et al. (2010), Lüpkes et al. (2012a), and Castellani et al. (2014). Note that this value is not critical for the parametrization. The value of β comes from Eq. (59) by Lüpkes et al. (2012a). C_hn,f is calculated as

$$\begin{array}{}\text{(8)}& C_{\mathrm{hn},\mathrm{f}}={\displaystyle \frac{C_{\mathrm{dn},\mathrm{f}}}{\mathrm{1}+C_{\mathrm{a},\mathrm{f}}\sqrt{C_{\mathrm{dn},\mathrm{f}}}}},\end{array}$$

where

$$\begin{array}{}\text{(9)}& C_{\mathrm{a},\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{k}}\ln \left({\displaystyle \frac{z_{\mathrm{0},\mathrm{i}}}{z_{\mathrm{t},\mathrm{i}}}}\right).\end{array}$$

Equation (8) represents a simple algebraic transformation of Eq. (60) by Lüpkes and Gryanik (2015), making use of their Eqs. (59) and (61) with α_f=α. The skin drag roughness length over sea ice is set to

$$\begin{array}{}\text{(10)}& {\displaystyle }{z}_{\mathrm{0},\mathrm{i}}=\mathrm{0.69}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m},\end{array}$$

and the scalar roughness length over sea ice is parameterized as mostly done in the literature as

$$\begin{array}{}\text{(11)}& z_{\mathrm{t},\mathrm{i}}=\mathit{\alpha }{z}_{\mathrm{0},\mathrm{i}},\end{array}$$

with

$$\begin{array}{}\text{(12)}& \mathit{\alpha }=\exp \left[\mathrm{3.0}-\mathrm{29.53}{z}_{\mathrm{0},\mathrm{i}}^{\mathrm{0.25}}\right].\end{array}$$

The CTRL and SENS simulations each comprise an ensemble of 10 members, which only differ in their ice–ocean initial state. The ice–ocean initial conditions for CTRL and SENS were produced in exactly the same way as for BASE (see Sect. 2.1.2). ERA-I provided the boundary forcing also as in the BASE simulations.

2.2 Datasets for evaluation

2.2.1 Sea ice drift

For the evaluation of the sea ice drift speed, satellite-based daily sea ice drift observations from Kimura et al. (2013; referred to as KIMURA hereafter) are used. There, the improved maximum cross-correlation method (Kimura and Wakatsuchi, 2000, 2004) was applied to detect ice motions based on AMSR-E (Advanced Microwave Scanning Radiometer for EOS) brightness temperature. The horizontal resolution of the KIMURA dataset is 60 km×60 km. The uncertainty of the KIMURA data over the Arctic in summer is from 1.12 to 1.47 km d⁻¹ and depends on the drift speed (Sumata et al., 2015a). In winter, the uncertainty is at least 50 % smaller than in summer and depends on the drift speed too (Sumata et al., 2015b). Although the accuracy of KIMURA drift speed is lower than that of buoy data, it has a much wider spatial and temporal coverage and is therefore appropriate for regional model evaluation (Sumata et al., 2015a). Another advantage of the KIMURA product is that it provides ice drift data both in winter and summer. More details are given by Kimura et al. (2013) and Sumata et al. (2015a).

In addition, daily sea ice drift speed from the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS; Zhang and Rothrock, 2003) is used. This enables a consistent and simultaneous evaluation of sea ice drift speed, concentration and thickness, and near-surface wind speed. The PIOMAS data are downloaded from ftp://pscftp.apl.washington.edu/zhang/PIOMAS/data/v2.1/ (last access: 1 June 2019). The mean horizontal resolution in the Arctic is approximately 22 km according to the PIOMAS description given in http://psc.apl.uw.edu/research/projects/projections-of-an-ice-diminished-arctic-ocean/model/ (last access: 30 April 2020). Detailed information about the PIOMAS dataset is given by Schweiger et al. (2011).

2.2.2 Sea ice concentration

For sea ice concentration (SIC), the NSIDC bootstrap daily SIC over the Northern Hemisphere, version 3 (https://nsidc.org/data/nsidc-0079, last access: 25 May 2019), is used. This SIC dataset is based on Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data and has been generated using a bootstrap algorithm with daily varying tie points. It is gridded on the 25 km×25 km polar stereographic grid and provides an overall retrieval accuracy of approximately 5 %–10 %, but the uncertainty could be larger due to factors such as the proximity to land or ice edge, melting during summer, the presence of a large fraction of thin ice or melt ponds within a pixel, and the presence of stormy weather conditions (Comiso, 2017).

2.2.3 Sea ice thickness

Since Arctic basin-wide long-term sea ice thickness (SIT) observations are not available, SIT data from PIOMAS are used in this study as a substitute for observational SIT as done in previous studies (Docquier et al., 2017; Johnson et al., 2012; Shu et al., 2015; Stroeve et al., 2014). However, we have to recall that PIOMAS ice thicknesses represent simulation results of a coupled ice–ocean model, even though this model is constrained through the assimilation of observed sea ice concentrations and sea surface temperatures. Schweiger et al. (2011) showed that the PIOMAS ice thickness agrees with ICESat ice thickness retrievals (in the order of 0.1 m mean difference) and that the spatial thickness patterns agree with each other (pattern correlation coefficients >0.8). However, PIOMAS appears to overestimate thin ice thickness and to underestimate thick ice (Schweiger et al., 2011).

2.2.4 Near-surface wind

For the near-surface wind speed (WS), daily 10 m wind speed from ERA-I is used. The ERA-I data were downloaded from the MARS archive at ECMWF and interpolated to the same $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid as used in the model's atmosphere component HIRHAM5. More information about this dataset is given by Berrisford et al. (2011). The NCEP/NCAR 10 m wind speed with $\mathrm{1.875}{}^{\circ }\times \mathrm{1.9}{}^{\circ }$ horizontal resolution is used to accompany the PIOMAS sea ice data as this data were used as the wind forcing for PIOMAS. More information about the NCEP/NCAR dataset is given by Kalnay et al. (1996).

2.3 Analysis methods

As the different evaluation datasets for sea ice have different spatial resolution, a bilinear interpolation method is used to remap them onto the NAOSIM grid. The common analysis period used in this study is the 12-year-long period 2003–2014 (only December is included for 2012). This limitation is related to the KIMURA data, which are only available since October 2002, and January to November data are missing for 2012 due to the transition from AMSR-E to AMSR-2. We focus on extended summer (JJAS) and winter (DJFM) in this study.

The domain for the basin-wide analysis covers the Arctic Ocean (referred to as the study domain hereafter; enclosed by the purple line in Fig. 1). The study domain is defined following Tandon et al. (2018) and excludes the grid points within a distance of 150 km from each coastline.

The simulations are evaluated by means of the commonly used climatological approach. For this, we present multiyear seasonal mean spatial maps and results spatially averaged over the study domain. The model data represent an average of the 10 ensemble members, which were first spatially averaged over the study domain and then time averaged over the period 2003–2014. Following Olason and Notz (2014) and Docquier et al. (2017), we use scatter plots showing Arctic basin-wide and multiyear averaged monthly mean sea ice drift speed against sea ice conditions (sea ice concentration and thickness) to evaluate the relationships between sea ice drift speed and sea ice conditions. The linear fit lines are added in the scatter plots to assist the comparison of the relationship in the model and in the observation/reanalysis. Furthermore, we present evaluation results based on daily data on the grid scale, i.e., for all grid points within the study domain. With this we aim to statistically evaluate the high spatial and temporal variability in the domain, and we represent this by means of box–whisker plots. Therein, the horizontal bar represents the median, the notch represents the 95 % confidence interval of the median, the dot represents the mean, the top and bottom of the box represent the 75th and 25th percentiles, and the upper/lower whiskers represent the maximum/minimum values within 1.5 times the interquartile range (IQR) to 75th/25th percentiles.

Normalized ensemble mean differences between SENS and CTRL are used to investigate the influence of sea ice form drag on atmosphere–ice momentum and heat fluxes, sea ice states, and motion and were calculated by dividing ensemble mean differences of SENS minus CTRL with ensemble standard deviations of these differences. Assuming that the differences between two random simulations are normally distributed around zero, the normalized differences enable a rough estimate of the statistical significance of the differences. Normalized differences greater than 2 (3) or lower than −2 (−3) indicate that the difference is significant on the 95 % (99.7  %) level.

3 Evaluation of simulated sea ice drift speed

First, the skill of the BASE simulated mean sea ice drift speed (SID) is quantified (Sect. 3.1). SID is mainly governed by near-surface wind, sea ice conditions, and ocean currents. The dependency of SID on near-surface wind speed (WS; Sect. 3.2) and on sea ice conditions (Sect. 3.3 and 3.4) is evaluated in each case in terms of both the climatological and the daily grid-scale views.

3.1 Multiyear mean sea ice drift speed (SID)

The simulated mean SID shows a distinct spatial pattern with highest drift speed near the ice edges in the Barents Sea, Greenland Sea, and Labrador Sea in winter and over the Alaskan coast in summer (Fig. 1a). Compared to the KIMURA dataset, the study-domain-mean bias (RMSE) is 1.72 km d⁻¹ (2.12 km d⁻¹) in winter and −0.03 km d⁻¹ (0.91 km d⁻¹) in summer. The model generally overestimates SID in the ice edge zone and north of the Canadian Arctic Archipelago, the region of thickest ice, with a maximum bias of ca. 6 km d⁻¹ in winter and a smaller bias in summer (Fig. 1b and d). Compared with the uncertainty in the KIMURA sea ice drift speed provided by Sumata et al. (2015a, b), the model bias in summer is close to or slightly smaller than the uncertainty of the KIMURA data. This indicates that the sign of the model bias in summer SID is uncertain. In winter, however, the model clearly overestimates the SID over the central Arctic and north of the Canadian Arctic Archipelago and Greenland, even if considering the uncertainty of the KIMURA data. This overestimation of SID in the thick-ice region may be linked to the underestimation of SIT and SIC and overestimation of WS in that region (Supplement Figs. S1–S3). The SIT differences between HIRHAM–NAOSIM 2.0 and PIOMAS can partly be attributed to the feedback between atmosphere and ocean–sea ice in the fully coupled atmosphere–ocean–sea ice model, which is not present in an ocean–sea ice model like PIOMAS. Analysis of the SIT differences between HIRHAM–NAOSIM 2.0 and CryoSat-2 during winter 2010–2014 (Fig. S2; see Hendricks and Ricker, 2019, for details about the CryoSat-2 SIT data used here) confirms that HIRHAM–NAOSIM 2.0 underestimates the SIT over the central Arctic and north of the Canadian Arctic Archipelago and Greenland, at least in winter. The underestimation of SIT and sea ice volume compared to PIOMAS has been discussed by Dorn et al. (2019) for the longer time period of 1979–2016.

The model reproduces the basic mean seasonal cycle of SID of the KIMURA data (Fig. 2) with respect to timing and absolute values between May and December. The average bias between June and November is almost zero (−0.0003 km d⁻¹), and the corresponding RMSE is 0.28 km d⁻¹. The maximum ice drift speed (10 km d⁻¹) occurs in October both in the model and in the KIMURA data. However, the model substantially overestimates the SID speed in winter and early spring. The averaged bias (RMSE) between December and May reaches 1.52 km d⁻¹ (1.64 km d⁻¹). In addition, the modeled minimum SID occurs in May, 2 months delayed compared to the KIMURA data (minimum SID in March). This reflects a longer-lasting and much weaker SID reduction from autumn to spring in the model compared to the KIMURA data. Furthermore, the interannual variation of SID in the model is lower than in the KIMURA data (Fig. 2). For example, the interannual amplitude is 0.70 km d⁻¹ (1.26 km d⁻¹) in the model (KIMURA data) in March and 0.36 km d⁻¹ (0.74 km d⁻¹) in September, respectively.

Figure 2 Mean annual cycle of sea ice drift speed (km d⁻¹, solid lines) and 10 m wind speed (m s⁻¹, dashed lines), based on the model (ensemble mean of HIRHAM–NAOSIM 2.0; blue lines) and observation/reanalysis (KIMURA ice drift, ERA-I wind; black lines) for 2003–2014 over the study domain (indicated in Fig. 1). The across-ensemble scatter (standard deviation) of the simulations is included as shaded area. The interannual variation is shown by error bars.

The overestimation of winter SID could be related to the underestimation of winter SIT (Fig. S2). Besides, the prescribed values of the ice–ocean drag coefficient and the ice strength parameter could also play a role. The ice–ocean drag coefficient, C_dw, in the base configuration of HIRHAM–NAOSIM 2.0 ($C_{\mathrm{dw}}=\mathrm{5.5}\times {\mathrm{10}}^{-\mathrm{3}}$) is comparable to other CMIP5 models, even though a few models use 3–4 times higher coefficients (see Tandon et al., 2018). Higher ice–ocean drag might damp the SID and its strong dependency on the wind speed. Docquier et al. (2017) showed that the higher the ice strength parameter, the lower the winter SID and the lower SIT. The ice strength parameter, P^∗, in the base configuration of HIRHAM–NAOSIM 2.0 ($P^{\ast }=\mathrm{30}$ 000 N m⁻²) is already slightly higher than in all CMIP5 models (see Tandon et al., 2018). This value has been established for stand-alone ice–ocean simulations with daily wind forcing but might still be too low considering the hourly wind forcing from the interactively coupled atmosphere.

3.2 Dependency of sea ice drift speed (SID) on near-surface wind speed (WS)

3.2.1 Climatological view

As shown in previous studies (Docquier et al., 2017; Kushner et al., 2018; Olason and Notz, 2014; Tandon et al., 2018), the observed distinct mean seasonal cycle of SID (maximum in autumn, minimum in spring) is obviously not solely controlled by the wind speed, which is strongest in winter and weakest in summer (Fig. 2). The phase lag between the seasonal cycle of SID and WS is about 3–4 months in observations/reanalysis (KIMURA ice drift/ERA-I wind). The modeled seasonal cycle and magnitude of the WS agrees well with ERA-I. According to the delayed SID minimum (Sect. 3.1), the phase lag between the simulated seasonal cycle of SID and WS is reduced to about 1 month, like in many CMIP3 and CMIP5 models (Rampal et al., 2011; Tandon et al., 2018), leading to a higher correlation between SID and WS. This indicates that the modeled SID is much more strongly controlled by the wind speed than the observed SID.

Another metric to quantify the mean relationship between SID and WS is the wind factor, which is the ratio of SID to WS. Averaged over the study domain, the simulated wind factor is 1.77 % in winter and 1.87 % in summer, which agrees with the observations/reanalysis (KIMURA ice drift/ERA-I wind) in the sense that the averaged wind factor is smaller in winter (1.42 %) than in summer (1.96 %). However, compared to the observations, the simulated wind factor is too large in winter. As mentioned in Sect.3.1, too-high sensitivity of the SID to the wind in winter may be related to the underestimated SIT and model parameters governing the sea ice dynamics.

Figure 3 shows the spatial patterns of the multiyear mean wind factor over the Arctic. In winter, the simulated wind factor has its maximum associated with strong wind and low ice thickness over the Baffin Bay and Greenland Sea. Within the study domain, the maximum occurs along the Transpolar Drift Stream (∼2 %) and the minimum over the north of Greenland and the Canadian Arctic Archipelago where sea ice is thickest (<0.5 %). The simulated pattern roughly agrees with the corresponding pattern of KIMURA ice drift/ERA-I wind (Figs. 3 and S4) and previous results using SSM/I ice drift/NCEP wind (Spreen et al., 2011). In winter, however, the simulated wind factor is overestimated compared to the KIMURA/ERA-I data almost everywhere over the study domain, with the maximum bias reaching 1 % over the thick ice north of the Canadian Arctic Archipelago (Fig. 3). In summer, the modeled wind factor peaks (∼3 %) along the marginal ice zone, such as in the coastal Beaufort Sea. This could be the result of a dynamical coupling between sea ice and the coastal ocean as suggested by Nakayama et al. (2012): in a coastal ocean covered with sea ice, wind-forced sea ice drift excites coastal trapped waves and generates fluctuating ocean currents. These ocean currents can enhance the sea ice drift when the current direction is the same as the wind-driven drift direction. In contrast to winter, the modeled wind factor in summer is underestimated over the study domain. These seasonally different bias patterns in the wind factor are associated with those of the SID (Fig. 1b and d).

Figure 3 Mean spatial pattern of wind factor (%) in the model (ensemble mean of HIRHAM–NAOSIM 2.0) for the 2003–2014 (a) winter (DJFM) and (c) summer (JJAS). Panels (b) and (d) are the model differences to the observation/reanalysis (model − KIMURA/ERA-I) for winter and summer respectively. The purple line indicates the study domain used for the basin-wide analysis.

3.2.2 Daily and grid-scale view

Previous studies reported on the ice drift's nearly linear increase with increasing wind for high and moderate WS (Thorndike and Colony, 1982) but no clear relationship for low WS (Leppäranta, 2011; Rossby and Montgomery, 1935). To investigate this in detail, Figs. 4 and 5 present the analysis of the dependency of SID on WS based on daily data and on the grid scale. For the observation/reanalysis reference we use KIMURA ice drift and ERA-I wind.

Figure 4 Box–whisker plots of the relationship between sea ice drift speed and sea ice concentration for different near-surface wind speed classes (different colors) for the 2003–2014 (a) winter (DJFM) and (c) summer (JJAS) in the model (HIRHAM–NAOSIM 2.0). Panels (b) and (d) are the same as (a) and (c), respectively, but based on the observation/reanalysis data. For the model, all 10 ensemble members are included. The plot is based on daily data and on all grid points within the study domain indicated in Fig. 1. The horizontal bar represents the median, the notch represents the 95 % confidence interval of the median, the dot represents the mean, the top and bottom of the box represent the 75th and 25th percentiles, and the upper/lower whiskers represent the maximum/minimum value within 1.5 times the interquartile range (IQR) to 75th/25th percentiles. The numbers above the boxplots represent the slopes of near-surface wind and sea ice drift speed fit lines (unit: km d⁻¹ per 1 m s⁻¹ wind speed change; font colors as for the wind speed classes). The numbers right of the boxplots represent the slopes of sea ice concentration and sea ice drift speed fit lines (unit: km d⁻¹ per 10 % sea ice concentration change). A bold and asterisked number indicates that the slope of the fit line is significant at the 95 % level. In the labels of different sea ice concentration and 10 m wind speed classes, “(” means exclusive and “]” means inclusive.

Figure 5 Relationship between sea ice drift speed and sea ice concentration for different near-surface wind speed classes (different colors) for the 2003–2014 (a) winter (DJFM) and (d) summer (JJAS) in the model (HIRHAM–NAOSIM 2.0). Panels (b) and (e) are the same as (a) and (d), respectively, but based on the observation/reanalysis data. Panels (c) and (f) are the same as (a) and (d), respectively, but based on PIOMAS data. The points in the plot are the median value of all the daily data and on all grid points for a certain wind speed and sea ice concentration, within the study domain indicated in Fig. 1. In the labels of different sea ice concentration and 10 m wind speed classes, “(” means exclusive and “]” means inclusive.

During winter, the simulated SID significantly increases with increasing WS. According to the median values of SID and the associated best fit line slope of SID vs. WS, the increase in SID varies from 1.29 to 1.49 km d⁻¹ per 1 m s⁻¹ WS increase under different SIC classes (Figs. 4a and 5a). The strongest SID increase per WS increase occurs for SIC of 50 %–70 % (Fig. 4a). In the observations/reanalysis, the SID only consistently increases with increasing WS when SIC is higher than 90 % (Figs. 4b and 5b). The corresponding slope indicates 0.93 km d⁻¹ per 1 m s⁻¹ WS increase, which is an approximately 40 % smaller increase than the modeled value. Both the model (for all SIC classes) and the observation/reanalysis (for SIC >90 %) show a linear increase in SID with increasing WS for strong winds (WS >4 m s⁻¹) but a weak dependency of SID on WS for lower winds.

During summer, both the model and the observation/reanalysis show a consistent general increase in SID with increasing WS for all SIC classes when WS >4 m s⁻¹ (Figs. 4c, 4d, 5d and 5e). The simulated magnitude of SID increase per 1 m s⁻¹ WS increase is similar to that in winter. Again, the simulated SID increases faster with increasing WS than the observed SID, by a factor of about 2–2.5. Another striking difference between the simulations and the observation/reanalysis is that the simulated SID–WS relation is rather linear (only slightly attenuated when WS ≤2 m s⁻¹). In contrast, this relation is highly nonlinear in the observation/reanalysis (Figs. 4d, 5e). The observed increase in SID with increasing WS is much stronger when WS passes a threshold (of ca. 4 m s⁻¹), compared to lower WS. For weak winds (WS lower than 4 m s⁻¹), SID is not much affected by WS changes.

Generally, the observed SID (for a certain WS and SIC class) shows an about 2 times higher spatiotemporal variability (indicated by the larger IQR; Fig. 4) than the modeled one. This discrepancy might be linked to a data source difference or inconsistency between the used KIMURA SID observations and ERA-I reanalysis WS data, but it might also hint to a model weakness.

The PIOMAS-based SID–WS relation (PIOMAS ice drift, NCEP/NCAR wind) is very similar to that in the HIRHAM–NAOSIM 2.0 simulations: SID increases with increasing WS for all SIC classes, both in winter and summer (Fig. 5c and f). However, the SID–WS relation in PIOMAS is stronger than in HIRHAM–NAOSIM 2.0 during winter and summer (Fig. 5). This may reflect a higher air–ice drag in PIOMAS than in HIRHAM–NAOSIM 2.0.

3.3 Dependency of sea ice drift speed (SID) on sea ice concentration (SIC)

3.3.1 Climatological view

Figure 6a shows the relationship between SID and SIC in terms of the mean seasonal cycle. An inverse correlation between SID and SIC is expected because a reduction in SIC reduces the internal friction, brings the SID closer to free drift conditions, and increases SID. We find such an inverse correlation between SID and SIC both in the model and in the observations (KIMURA ice drift, NSIDC bootstrap ice concentration), as well as in PIOMAS from June to September when SIC is low. The corresponding fit line slope in the model is 4.04 km d⁻¹ SID increase per 10 % SIC decrease, which agrees with the slopes in the observations and in PIOMAS. During the winter months (January to May) when SIC >95 %, such a relationship is not found both in the observations and in the models. The SID of the 12 months can be categorized into two groups based on SIC values: group one for SIC smaller than 95 % and group two for SIC larger than 95 %. The observed mean SID value in group one (9.23 km d⁻¹) is obviously larger than in group two (7.28 km d⁻¹). This is also found in PIOMAS. The model reproduces the SID in group one but largely overestimates the SID in group two. Therefore, the simulated mean SID values in group one (9.24 km d⁻¹) and two (8.95 km d⁻¹) are quite close to each other. This model deficit may result from the too-strong coupling between the SID and the WS in the model (Fig. 2; Sect. 3.2.1). Following Docquier et al. (2017), we also calculate the slope of the SID–SIC best fit line based on data from 12 months, even though there is no clear SID–SIC linear relation over the full year data. The slope in the model (0.66 km d⁻¹ SID increase per 10 % SIC decrease) is much weaker than in the observations (6.93 km d⁻¹ SID increase per 10 % SIC decrease) and in PIOMAS (3.16 km d⁻¹ SID increase per 10 % SIC decrease). The simulated relationship is also much weaker than in the ocean model NEMO–LIM3.6 as reported by Docquier et al. (2017). The underestimated relationship in the model can be explained by the overestimated SID in winter (Fig. 2).

Figure 6 Scatter plots of monthly mean sea ice drift speed against (a) sea ice concentration and (b) sea ice thickness, averaged over the period of 2003–2014 and the study domain (indicated in Fig. 1). Numbers denote the months. Results are shown for the model (HIRHAM–NAOSIM 2.0) ensemble simulations (red), KIMURA for ice drift speed plus NSIDC bootstrap for ice concentration (black), and for PIOMAS data (blue).

3.3.2 Daily and grid-scale view

On daily grid scale, the model shows an inverse correlation of SID with SIC during winter (Figs. 4a, 5a), differently to the climatological view (Sect. 3.3.1). SID increases with decreasing SIC when SIC is larger than ca. 30 %. The slopes of the simulated SID–SIC relation across all SIC and WS classes vary from 3.41 to 5.64 km d⁻¹ SID increase per 10 % SIC decrease. In the observation, there is no clear SID–SIC relation because the median SID uncertainty is very high when SIC is lower than 90 % (Figs. 4b, 5b). The small sample size is one of the reasons for the high uncertainty (Table 1). However, considering the uncertainty of the median value (represented by the notch), the statistically lower SID for high SIC (>90 %), compared to the SID for lower SIC, stands out both in the observations and in the simulations. The observations indicate an abrupt drop of SID at 90 % SIC. This is in accordance with observations by Shirokov (1977), who showed that SID starts to decrease linearly with increasing SIC but drops abruptly at a SIC of 90 % when internal friction forcing starts its significant influence on SID. This is not reproduced by the simulations and indicates too-low internal friction in the model.

Table 1 The sample sizes of sea ice drift speed data under different 10 m wind speed and sea ice concentration classes in Fig. 4. In the labels of different sea ice concentration and 10 m wind speed classes, “(” means exclusive and “]” means inclusive.

The model also shows an inverse correlation of SID with SIC for WS >2 m s⁻¹ during summer (Figs. 4c and 5d). Under different WS classes, the modeled SID–SIC relation varies from 2.12 to 3.16 km d⁻¹ SID increase per 10 % SIC decrease. The higher WS is, the stronger the SID–SIC relation is. However, these findings from the simulations are not confirmed by the observations. The observations (KIMURA drift/NSIDC bootstrap ice concentration) do not show an inverse correlation between SID and SIC for all SIC and WS classes in summer (Figs. 4d and 5e).

The SID–SIC relation in PIOMAS is just the opposite compared to HIRHAM–NAOSIM 2.0, both in winter and in summer: SID increases with increasing SIC when SIC is lower than 90 % (Fig. 5c and f). PIOMAS shows a SID–SIC relation inconsistent with the observed relation. Thus, the PIOMAS relation might violate physical consistency due to the optimal interpolation method used to obtain a realistic sea ice field (concentration). This assimilation method contains addition/subtraction of sea ice into the system at every assimilation time step, when the modeled sea ice concentration differs from the observed one. Due to the addition/subtraction of sea ice (called increment or innovation in the terminology of data assimilation), PIOMAS does not necessarily preserve the physical relations described in the underlying sea ice–ocean model. Such an inconsistency is one of the drawbacks of the optimal interpolation method. Therefore, relations to assimilated physical variables should be examined with caution.

Generally, the obvious difference between the SID–SIC relation on the multiyear Arctic mean scale and on the daily grid scale emphasizes its strong dependency on the temporal and spatial scale.

3.4 Dependency of sea ice drift speed (SID) on sea ice thickness (SIT)

3.4.1 Climatological view

All datasets (HIRHAM–NAOSIM 2.0 simulations, KIMURA drift/PIOMAS thickness, and PIOMAS drift/PIOMAS thickness) show that SID decreases with increasing SIT (Fig. 6b). The SID–SIT fit line slope calculated based on HIRHAM–NAOSIM 2.0 is 1.10 km d⁻¹ SID decrease per 1 m SIT increase. The observed SID–SIT relation is stronger. Data of KIMURA drift/PIOMAS thickness and PIOMAS drift/PIOMAS thickness indicate 2.17 and 1.83 km d⁻¹ SID decrease per 1 m SIT increase, respectively. The simulated weaker dependency of SID on SIT is mainly due to the overestimated SID during winter and spring (December to May), when the ice is thick. As pointed out by Olason and Notz (2014), the inverse correlation between drift speed and thickness in winter, when the ice concentration is high, is physically plausible, but the inverse correlation in summer, when the ice concentration is lower, is probably only of statistical nature.

3.4.2 Daily and grid-scale view

As in the climatological view, the simulations show generally that SID decreases with increasing SIT (Figs. 7a and c and S5a, c), both in winter and in summer. In winter, this relation occurs only for SIT smaller than approximately 4 m. The simulated SID–SIT relation that calculated based on all thickness categories varies from 0.38 to 0.95 km d⁻¹ SID decrease per 1 m SIT increase across different wind categories (Supplement Fig. S5a). PIOMAS simulates a stronger SID–SIT relation (Supplement Fig. S5b); the slope of the SID–SIT fit line varies from 0.65 to 2.64 km d⁻¹ SID decrease per 1 m SIT increase, which is larger by a factor of up to 3.5 compared to the simulations. In summer, the simulated inverse SID–SIT correlation is stronger than in winter, particularly for high WS (Figs. 7a and c and S5a, c). The SID–SIT slope varies from 0.71 to 1.69 km d⁻¹ per 1m SIT increase across different wind categories. In PIOMAS, the SID–SIT slopes are not significant when WS >4 m s⁻¹. This is mainly influenced by the increase in SID when SIC increased from (0.1,0.3] to (0.3, 0.5]. Furthermore, Fig. 7 confirms the nonlinearity of the SID–WS relation, with a strong impact of WS changes on SID for high and moderate wind speeds, but with a lower impact for low wind speeds (WS ≤4 m s⁻¹).

Figure 7 Relationship between sea ice drift speed and sea ice thickness for different near-surface wind speed classes (different colors) for the 2003–2014 (a) winter (DJFM) and (c) summer (JJAS) in the model (HIRHAM–NAOSIM 2.0). Panels (b) and (d) are the same as (a) and (c), respectively, but based on PIOMAS data. The points in the plot are the median value of all the daily data and on all grid points within the study domain based on Fig. 1. In the labels of different sea ice thickness and 10 m wind speed classes, “(” means exclusive and “]” means inclusive.

4 Sensitivity to the parameterization including sea ice form drag

4.1 Model results

In the previous section, we investigated the dependency of SID on WS and sea ice conditions in the BASE simulations, which do not account for the sea ice form drag. In this section, we analyze how an additional sea ice form drag influences the near-surface atmospheric fluxes, SID, and its dependency on WS and sea ice condition. For this we compare the 1-year ensemble simulations CTRL (without sea ice form drag) and SENS (with sea ice form drag; see Sect. 2.1.3 for details). Tsamados et al. (2014) showed that an additional form drag mainly improves the sea ice drift in summer. Therefore, our analysis is focused on summer (JJAS).

Figure 8a indicates that the increased ice roughness due to the additional form drag leads to an increased surface momentum flux over most regions of the Arctic Ocean, even though this increase is not significant from a statistical point of view. In a few regions the momentum flux does not increase, which indicates that either the neutral drag coefficient over ice is very low or the wind speed and/or the atmospheric stratification has changed. SID mainly increases over the MIZ along the Russian coasts (Fig. 8b). The spatial pattern of SID changes is similar to that of the momentum flux changes, with a pattern correlation of 0.66 (0.76) between their normalized (not normalized) differences (normalized differences were calculated by dividing ensemble mean differences with ensemble standard deviations of these differences; see Sect. 2.3). This reveals that the SID increase is mainly associated with an increased momentum flux from the atmosphere to the ice. Related with this, WS decreases over most parts of the Arctic Ocean when the sea ice becomes rougher (Fig. 8c). Changes in turbulent heat fluxes, SIC, and 2 m air temperature are small and mixed with regional increases and decreases (Fig. 8d–f). Overall, the ensemble mean differences between SENS and CTRL are statistically insignificant, indicating a large across-ensemble scatter due to high internally generated model variability.

Figure 8 Normalized ensemble mean differences (SENS minus CTRL) of (a) momentum flux, (b) sea ice drift speed, (c) 10 m wind speed, (d) turbulent heat flux, (e) sea ice concentration, and (f) 2 m air temperature for the 2007 summer (JJAS). The ensemble mean difference is normalized by the cross-ensemble standard deviation of the differences. The purple contours represent the ensemble mean sea ice concentration in the CTRL experiment.

Nonetheless, the increase in SID for low SIC in the MIZ and the smaller effect on the more consolidated ice in the central Arctic should imply a stronger Arctic-wide SID–SIC relation. Figures 9 and 10a confirm the stronger increase in SID with decreasing SIC in SENS. The increased SIC dependency of SID can be explained by the fact that the form drag contribution to the air–ice drag is effectively proportional to (1−A)^β (see Eqs. 5 and 7). This means that the lower the SIC, the higher the form drag and the higher the air–ice drag. When SIC is around 90 %, the sum of form drag and skin drag in SENS is of the same order of magnitude as the skin drag in CTRL.

Figure 9 Box–whisker plots of the relationship between sea ice drift speed and sea ice concentration for different near-surface wind speed classes (different colors) for the 2007 summer (JJAS) in the (a) CTRL and (b) SENS experiment. The plot is based on daily data and on all grid points within the study domain indicated in Fig. 1. The horizontal bar represents the median, the notch represents the 95 % confidence interval of the median, the dot represents the mean, the top and bottom of the box represent the 75th and 25th percentiles, and the upper/lower whiskers represent the maximum/minimum value within 1.5 times the interquartile range (IQR) to 75th/25th percentiles. The numbers above the boxplots represent the slopes of near-surface wind and sea ice drift speed fit lines (unit: km d⁻¹ per 1 m s⁻¹ wind speed change; font colors as for the wind speed classes). The numbers right of the boxplots represent the slopes of sea ice concentration and sea ice drift speed fit lines (unit: km d⁻¹ per 10 % sea ice concentration change). A bold and asterisked number indicates that the slope of the fit line is significant at the 95 % level. In the labels of different sea ice concentration and 10 m wind speed classes, “(” means exclusive and “]” means inclusive.

Figure 10 The (a) simulated relationship between sea ice drift speed and sea ice concentration for different near-surface wind speed classes (different colors) for the 2007 summer (JJAS) in the CTRL (circle marker and solid line) and SENS (cross marker and dashed line) experiment. The relationship based on KIMURA sea ice drift speed, NSIDC bootstrap sea ice concentration, and ERA-Interim 10 m wind speed is shown in (b). The points in the plot are the median value of all the daily data and on all grid points within the study domain indicated in Fig. 1 under certain wind speed and sea ice concentration classes. In the labels of different sea ice concentration and 10 m wind speed classes, “(” means exclusive and “]” means inclusive.

Furthermore, there is an increased WS dependency; i.e., the SID change is larger for high WS. For example, the SID–SIC relation for low-wind conditions (WS of 2–4 m s⁻¹) in SENS is 5.26 km d⁻¹ SID increase per 10 % SIC decrease, which is about 2.5 times higher than in CTRL. The differences between SENS and CTRL increase even more for high wind speeds. For WS within 8–10 m s⁻¹, the SID–SIC relation in SENS is 11.79 km d⁻¹ SID increase per 10 % SIC decrease, which is about 4.6 times higher than in CTRL. The dependence of SID on wind is complex. Two external forces act on the flow, the ocean drag and atmospheric drag. When the ocean drag is much larger than the atmospheric drag (low wind speed), the drift is mainly governed by the ocean drag even with a doubled atmospheric drag coefficient. However, for high wind speeds, a doubling of the atmospheric drag coefficient has a large effect on the drift speed, also due to the proportionality of the air–ice stress to the square of WS.

4.2 Model versus observation

The increased dependency of SID on WS and SIC in SENS compared to CTRL does not reduce the deviation to the observation/reanalysis. In contrast, Fig. 10 shows that the overestimation of the dependency of SID on WS and SIC in SENS is larger than in CTRL.

Figure 11 shows the spatial distribution of the summer 2007 wind factor, WS, and near-surface air temperature from observation/reanalysis data, as well as the deviations of these three variables in CTRL and SENS from observation/reanalysis. It is obvious that both the bias patterns and the bias magnitudes in CTRL and in SENS are quite similar. Considering the ensemble mean bias and taking the internal model variability into account, it is hard to detect significant changes in SENS compared to CTRL, as discussed in Sect. 4.1.

Figure 11 The 2007 summer (JJAS) (a) wind factor, (d) 10 m wind speed, and (g) 2 m air temperature from ERA-Interim (KIMURA sea ice drift is used for wind factor calculation). (b) Wind factor, (e) 10 m wind speed, and (h) 2 m air temperature differences between the CTRL experiment and ERA-Interim (KIMURA and ERA-Interim for the wind factor). Panels (c), (f), and (i) are the same as (b), (e), and (h) but show the differences between the SENS experiment and ERA-Interim (KIMURA and ERA-Interim for the wind factor).

Although the ensemble mean of Arctic basin-wide mean SID from July to September in SENS is larger than in CTRL, the differences are not statistically significant due to a large ensemble spread (Fig. 12). Actually, there are no statistically significant differences in the Arctic basin-wide mean SID between CTRL and SENS in all months. From January to May, the simulated Arctic basin-wide mean SID (both in CTRL and in SENS) is higher than that in KIMURA. With respect to the summer months (June to September), the August SID in CTRL is lower than in KIMURA, while the July and September SID in SENS is higher than in KIMURA. For the Arctic basin-wide mean WS, there is no significant difference between CTRL and SENS as well as between model and reanalysis, except for May, when both model simulations significantly overestimate the WS.

Figure 12 Mean annual cycle of sea ice drift speed (km d⁻¹, solid lines) and 10 m wind speed (m s⁻¹, dashed lines), based on the CTRL experiment (ensemble mean; blue lines), SENS experiment (ensemble mean; red lines) and the observation/reanalysis (KIMURA ice drift, ERA-I wind; black lines) for 2007 over the study domain (indicated in Fig. 1). The across-ensemble scatters (standard deviation) of the simulations are included as shaded area (light blue for CTRL, orange for SENS).

4.3 Follow-up study

Although the new parameterization does not improve the simulated dependency of SID on WS and sea ice conditions compared to observations/reanalysis, the sensitivity study clearly shows that the new parameterization does increase the SID due to the added form drag. In a follow-up study, we are going to focus on several aspects to improve the simulations. First, tunable parameters of the new parameterization, such as z_0,i, z_t,i, z_0,f, C_e10, and β, represent an opportunity to better adapt the form drag parameterization itself to the observations. A first step could be the use of values found by Elvidge et al. (2016). A large effect can be expected by a modification of the skin drag coefficient, since a large region would be affected, and large variations in the drag due to pressure ridges allow a wide range of values. Second, model parameters outside the new parameterization, which have direct impact on SID, like the ice strength and ice–ocean drag coefficient, need to be harmonized with the new parameterization, since their values were chosen empirically in terms of adequately balanced performance of the model. A key is probably the oceanic form drag. Its effect is accounted for in the present study only indirectly via the constant oceanic drag coefficient. Such a parametrization is probably too simple, especially when atmospheric form drag is included (see also Tsamados et al., 2014). Birnbaum and Lüpkes (2002) as well as Lüpkes et al. (2012b) found in a mesoscale modeling study that oceanic form drag can have a strong decelerating effect on SID, especially when the sea ice concentration is low so that the discussed drawbacks for small sea ice fraction would be reduced or even removed. This effect of form drag on SID was discussed also by Steele et al. (1989). The parametrizations are evidently not balanced anymore after improving one key process of the SID-related atmosphere–ocean–ice interaction. A previous study on the surface–albedo feedback by Dorn et al. (2009) showed that an improved simulation can only be achieved by a harmonized combination of more sophisticated parameterizations of the related subprocesses. It can be assumed that this holds true for the SID-related subprocesses.

5 Summary and conclusions

We evaluated the sea ice drift speed (SID) and its dependency on the near-surface wind and the sea ice conditions (ice concentration and thickness) on the multiyear, Arctic-wide mean scale during 2003–2014. Compared with observations, the model does not fully capture the observed SID seasonal cycle but overestimates SID in winter–spring. Regardless, the model realistically describes the main drivers of the seasonal and long-term variations of Arctic SID: when the sea ice concentration (SIC) is lower than 95 % in summer–autumn, SID increases with decreasing SIC. However, when SIC is higher than 95 % in winter, the sea ice thickness (SIT) is the main factor for SID changes (higher SID for lower SIT). As the simulated SID is overestimated during the cold seasons, the modeled strength of the SID–SIC and SID–SIT relations is underestimated compared with observational data. The SID overestimation during winter in the model cannot simply be attributed to the underestimation of SIT or too-strong coupling between SID and WS. Further in-depth analyses of the wintertime sea ice dynamics and thermodynamics and the atmospheric and oceanic forcing are needed in the future to identify the possible cause of the SID overestimation.

The analysis on the daily and grid scale revealed that the SID relations with SIC and SIT are complex due to the large spatiotemporal variation of the sea ice. Based on observations, it is difficult to find a clear relation between SID and SIC. In the model, when SIC is larger than 30 %, a higher SIC is accompanied by a significantly lower median SID compared to the SID for lower SIC. In agreement with the multiyear and Arctic-wide findings, the SID decreases with increasing SIT, both in winter and in summer. However, in winter, this relation occurs only for high wind speeds. The simulated SID–SIT relation is stronger (weaker) in summer (winter) compared to PIOMAS.

We also evaluated the dependency of SID on near-surface wind speed (WS) on the daily and grid scale. The simulated increase in SID with increasing WS is consistent with the observation/reanalysis. Our analysis supports the earlier discussed strong nonlinearity of the SID–WS relation, with a strong impact of WS changes on SID for high and moderate wind speeds but with only a low impact for lower wind speeds. The weak dependency of SID on WS, when WS is low, was also shown by Lund et al. (2018), based on shipboard marine radar sea ice drift measurement over the MIZ. This changing relation in the low-wind regime may be caused by the increasing importance of upper-ocean currents for the drift compared to the wind (Lund et al., 2018). Another explanation was given by Leppäranta (2011) and Rossby and Montgomery (1935), who argued that the large-scale wind dominates over the local wind effect.

Finally, we investigated the impact of the changed parameterization of the transfer coefficients for heat and momentum over sea ice on the SID and its dependency on WS and sea ice conditions. The consideration of sea ice form drag effects increases the air-to-ice momentum flux and accordingly the SID over most of the Arctic. The largest effects appear for low SIC in the MIZ under high-wind conditions. The reason is that in the new parameterization the lower SIC is associated with a potentially larger atmosphere–ice momentum flux, since the form drag contribution is completely transferred to the ice. As a consequence, the increases in SID with decreasing SIC is stronger when the form drag is included.

The inclusion of the melt pond effect on the atmospheric form drag in the model might be beneficial. In the current version, form drag was only considered at the edges of ice floes, mainly in the marginal sea ice zone but not on top of the ice, where melt ponds cause form drag also during summer (Andreas et al., 2010; Lüpkes et al., 2012a). Additional form drag at the ice–ocean interface may further improve the simulated SID–WS relation, because the oceanic form drag has normally the opposite effect on the ice motion as the atmospheric form drag (Steele et al., 1989; Lüpkes et al., 2012b). Systematic biases in the reanalysis used for the calculation of the observed wind factor cannot be excluded, since form drag is not taken into account in the underlying atmospheric model. Therefore, the increased deviation of the simulated SID–WS relation from the observations/reanalysis does not necessarily mean that the implemented new parameterization worsens the SID–WS relation.