2.1 Forecast ocean forcing from TOPAZ4

The official European forecast for the Arctic is developed and run by the CMEMS Arc-MFC (Arctic Monitoring and Forecasting Centre)⁵. This uses the TOPAZ system (Simonsen et al., 2018; Sakov et al., 2012), which uses version 2.2.37 of the Hybrid Coordinate Ocean Model (HYCOM) (Bleck, 2002). In the current version (4) of TOPAZ, HYCOM is coupled to a sea ice model derived from version 4.1 of the Community Ice CodE (CICE: Hunke and Lipscomb, 2010); ice thermodynamics are described in Drange and Simonsen (1996), while the dynamics are based on the viscous–plastic (VP) sea ice rheology (implemented with the elastic–viscous–plastic (EVP) solver of Hunke and Dukowicz, 1997). The model's native grid covers the Arctic and North Atlantic oceans and has a horizontal resolution of between 11 and 16 km. TOPAZ4 uses the ensemble Kalman filter method (EnKF; Sakov and Oke, 2008) to assimilate remotely sensed sea level anomalies, sea surface temperature, sea ice concentration, sea ice thickness and Lagrangian sea ice velocities (the latter two in winter only), as well as temperature and salinity profiles from Argo floats and ice-tethered profilers. Data assimilation is performed weekly.

To force neXtSIM, we use the following daily-averaged variables from TOPAZ: the sea surface (0–3 m) ocean velocity, temperature and salinity (SST and SSS, respectively), and the mixed layer depth (MLD). We give more details of how they are used in Sect. 3.1 below.

2.2 Forecast atmospheric forcing from ECMWF

For our forecast demonstration, we use the latest version (Cycle 45r1) of the Integrated Forecast System from ECMWF (IFS: Owens and Hewson, 2018) to provide atmospheric forcing fields to neXtSIM. It consists of an atmospheric model coupled to the NEMO 3.4 ocean model (Nucleus for European Modelling of the Ocean), the LIM2 (Louvain-la-Neuve Sea Ice Model) sea ice model, the ECWAM (ECMWF Wave Model) wave model and a land surface model (HTESSEL). Its spatial resolution is 0.1^∘, and while its temporal resolution is 1 h we update our forcing fields less frequently (every 6 h).

The variables we use are the 10 m wind velocity, the 2 m air and dew point temperatures (the latter is used to determine the specific humidity of air for the latent heat flux calculation), the mean sea level pressure, the long- and short-wave downwelling radiation, and the total precipitation (this becomes snow if the 2 m air temperature is below 0 ^∘C).

2.3 Sea ice concentration products from OSI SAF

OSI SAF provides estimates of sea ice concentration derived from the Special Sensor Microwave Imager Sounder (SSMIS) radiometer (Tonboe et al., 2016; Tonboe and Lavelle, 2016; Lavelle et al., 2017) and from the Advanced Microwave Scanning Radiometer 2 (AMSR2: Lavelle et al., 2016a, b; Tonboe and Lavelle, 2015). The SSMIS algorithm uses the 19 GHz frequency (vertically polarised, footprint size about 56 km) and the 37 GHz frequency (both vertically and horizontally polarised, footprint size about 33 km). The AMSR2 algorithm uses three frequencies: 18.7, 36.5 and 89 GHz (also in vertical and horizontal polarisations with footprints from 22 to 5 km). The AMSR2 data are presented on a 10 km grid, and we chose this product over the higher-resolution (3.25 km) ASI product as we found it less noisy near the ice edge. These products are available daily within 12 h after acquisition and processing, so it is possible to assimilate these data in operational forecasts. However, the file for the day before the bulletin date (the day the model is run) does not arrive early enough to be assimilated in our daily run, which is launched at 03:00 (central European time). Therefore, we use the file from 2 d before the bulletin date.

As specified in the validation reports cited above the SSMIS has lower-resolution ice concentration but has the advantage of higher accuracy, while the AMSR2 algorithm has higher resolution but also higher uncertainties. In order to combine the advantages of these products we generated a blended product that was used for assimilation during the forecasts. Blending was performed with a weighted average of the two products (using the errors in the products to calculate the weights):

$$\begin{array}{}\text{(1)}& c_{\text{osisaf}}={\displaystyle \frac{c_{\text{ssmis}}{\mathit{\sigma }}_{\text{ssmis}}^{-\mathrm{2}}+c_{\text{ amsr2}}{\mathit{\sigma }}_{\text{amsr2}}^{-\mathrm{2}}}{{\mathit{\sigma }}_{\text{ssmis}}^{-\mathrm{2}}+{\mathit{\sigma }}_{\text{ amsr2}}^{-\mathrm{2}}}},\end{array}$$

where c denotes sea ice concentration and σ denotes the concentration uncertainty.

Sea ice extent, used as an evaluation metric of the model, was calculated from the concentration product as a sum of areas of all pixels within the model domain with concentration above 15 %. Sea ice extent uncertainty was calculated as a difference between the extents calculated from the sum of concentration and uncertainty and concentration alone.

We use both SSMIS and AMSR2 products for assimilation by the forecasts but only SSMIS for evaluation of the free run and our forecasts. This was because we found the AMSR2 product somewhat inconvenient due to missing sections of data, which made our evaluation statistics quite noisy. (OSI SAF SSMIS is therefore not an independent validation dataset for the forecasts.)

2.4 Sea ice drift from OSI SAF

We use this product for evaluation of both our free run and our forecasts. It is not assimilated. To produce it, low-resolution ice drift datasets are computed on a daily basis from aggregated maps of passive microwave (e.g. SSMIS, AMSR-E) or scatterometer (e.g. ASCAT) signals (all channels are used) using the continuous maximum cross-correlation method (CMCC: Lavergne et al., 2010; Lavergne and Eastwood, 2010; Lavergne, 2010). Daily 48 h ice drift vectors can be obtained at a spatial resolution of 62.5 km. As part of our evaluation we apply a filter on the uncertainty given in the product to remove less accurate observations. We usually take the maximum allowed 2 d drift uncertainty to be 20 km, which allows a reasonable sample size of vectors to compare the model to, and right throughout the year. However, it is useful to sometimes focus on the more precise observations when considering winter drift. Therefore in those cases we apply a stricter filter, where the 2 d drift uncertainty is less than 2.5 km. This completely excludes the summer period of May to September (rms uncertainty is about 12 km), since surface melting and a denser atmosphere preclude the retrieval of precise information. From October to April, we can still retain about 75 % of the observation vectors after using this threshold, with the removed vectors generally being close to the ice edge, the coast or the north pole. The error is higher in these regions as the sub-images on which the CMCC method is applied must be reduced to limit them to being inside regions where there actually is ice (Lavergne and Eastwood, 2010). In the case of the north pole, there are fewer observations there, while the vectors in the marginal ice zone (MIZ) have especially high uncertainties (sometimes up to 12 km) due to the high velocities in those regions, combined with the relatively long time interval over which the drift is calculated.

In order to compare neXtSIM drift to this product more accurately, every day at 12:00 we seed synthetic Lagrangian drifters at the grid points of the OSI SAF drift product and advect them for 48 h according to the ice velocity in the model. That is, their drift from the original position is updated every model time step (Δt)

$${\mathit{d}}_{\text{mod}}(\mathit{x},t^n+\mathrm{1})={\mathit{d}}_{\text{mod}}(\mathit{x},t^n)+\mathit{u}(\mathit{x},t^n)\mathrm{\Delta }t,$$

where u is the ice velocity in the model. (Note that the drift d_mod is a global field.) At output time, and when the model mesh needs regridding due to it becoming too deformed, the drifter positions are updated with

$${\mathit{x}}_{\text{mod},i}\left(t^{n+\mathrm{1}}\right)={\mathit{x}}_{\text{mod},i}\left(t^n\right)+{\mathit{d}}_{\text{mod}}\left({\mathit{x}}_{\text{mod},i}\left(t^n\right),t^{n+\mathrm{1}}\right),$$

and d_mod is reset to zero. This step, which requires d_mod to be interpolated, is done as little as possible for the sake of model performance and to avoid the accumulation of interpolation error. The total drift after 48 h is then compared to the OSI SAF drift product.

2.5 Sea ice thickness from CS2–SMOS

We use the CS2–SMOS sea ice thickness product (version 2.3: Ricker et al., 2017) to initialise our free run and forecast (this is done once only and is not to be confused with our daily assimilation, which corrects the initial conditions of individual forecast bulletins) and for long-term evaluation of our free run. This product is a daily hybrid product which combines thickness estimates from the CryoSat-2 (CS2) altimeter (more reliable for thicker ice, ≳ 0.5 m) and from SMOS (Soil Moisture and Ocean Salinity: better for thin ice). The CS2 altimeter tracks are somewhat sparse, so optimal interpolation (OI) is used to fill the gaps between them and the areas of thin ice. For this reason each file covers a 7 d period. The OI method requires a background field to be created with full coverage, and it requires that the background field is independent of the target week so it is created from the CS2 values for the 2 weeks ahead and behind the target week, as well as from the SMOS values for the day before the target week.

The errors from this approach can be particularly high in coastal areas of thick ice (for example north of Greenland and Canada), which are too thick to be measured by SMOS, but may be only covered by altimeter tracks every 2–3 weeks. For these gaps in coverage, the product uses the background field.

3.1 The neXtSIM model

neXtSIM is a stand-alone sea-ice model which can use winds and currents from a variety of atmospheric and oceanic models (hindcasts or forecasts). This makes it quite flexible and light to run and therefore ideal for a forecasting context. Its dynamical core is the new brittle Bingham–Maxwell (BBM) rheology (Ólason et al., 2021). (The version of the BBM rheology corresponding to the results in this paper is also summarised in Appendix A.) Rampal et al. (2016) presented results using a previous version of neXtSIM including an elasto-brittle (EB) rheology as described in Bouillon and Rampal (2015), showing good agreement with observed drift and concentration in particular. More recently, Rampal et al. (2019) showed the ability to reproduce the characteristic multi-fractal scalings of deformation when using the previous version of neXtSIM, which included the Maxwell elasto-brittle (MEB) rheology (Dansereau et al., 2016). The key contribution of the MEB rheology was the addition of a viscous dissipation of stress in areas where the ice is damaged, allowing it to move more freely. However, in longer-term simulations the MEB showed unrealistic pile-up of ice particularly along the northeast coast of Greenland and the northwest coast of Svalbard. This led to the further addition of a frictional element to the rheology (Ólason et al., 2021), which provides some resistance to compression (up to a threshold). This framework, of a spring in series with a composite element made up a dashpot and fraction element in parallel, is called BBM.

The dynamical equations are solved with a finite element method on a Lagrangian (moving) triangular mesh. The code is currently a parallelised C++ code, used by Rampal et al. (2019), and presented by Samaké et al. (2017). Momentum input comes from the wind and ocean stresses (a turning angle of 25^∘ is applied to the ocean velocity from the ocean forcing), the Coriolis force, and sea surface slope, and there is also a basal stress applied at the bottom of the ice when it becomes grounded. For this basal stress we follow the scheme of Lemieux et al. (2016b), using the parameters k₁=10, k₂=15 Pa m⁻¹, α_b=20, and $u_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹.

There is also a thermodynamic component of the code, and beneath the ice is a slab ocean with three variables: temperature, salinity and thickness. The temperature and salinity are modified by the heat and salinity fluxes determined by the thermodynamical model as ice melts and freezes and as the model interacts with the atmosphere. They are relaxed towards the SST and SSS from TOPAZ over a timescale of about 1 month, while the thickness of the slab ocean is taken directly to be the MLD of TOPAZ. This varies spatially and evolves with time according to the forecast from TOPAZ. The thermodynamical model is a three-category model (detailed in Rampal et al., 2019, Appendix A): open water, newly formed ice (treated as one ice layer and one snow layer; Semtner, 1976) and older ice (treated as two ice layers and a snow layer; Winton, 2000).

The older ice is characterised in the model by concentration, c, and thickness averaged over the entire cell (effective thickness or, in other words, volume per unit area), h. The absolute thickness of ice can be computed as the ratio: $H=h/c$; there is also an effective snow thickness, s. The young ice has concentration c_y, effective thickness h_y, and snow thickness s_y. The absolute thickness of young ice $H_y=h_y/c_y$ is constrained so that $H_{\text{min}}\le H_y\le H_{\text{max}}$. If H_y<H_min, c_y is reduced to ${c^{\prime }}_y=\left(H_y{c}_y\right)/H_{\text{min}}$; if H_y>H_max, some ice is moved to the older ice category. During this process H_y is reduced to ${H^{\prime }}_y=H_{\text{max}}$, while c_y is also reduced to ${c^{\prime }}_y=c_y(H_{\text{max}}-H_{\text{min}})/(H_y-H_{\text{min}})$. This concentration reduction is intended to give some lateral decrease in young ice volume and not just vertical. The corrected values for h_y and s_y are ${h^{\prime }}_y={c^{\prime }}_y{H}_{\text{max}}$ and ${s^{\prime }}_y={c^{\prime }}_y(s_y/c_y$); the corresponding properties for the older ice h, c and s are then increased in a ice-and-snow-volume-conserving manner. The values that we used for the thresholds on the absolute thin ice thickness are H_min=0.05 m and H_max=0.275 m.

The domain we present simulations for is a pan-Arctic one, with a resolution of about 7.5 km. A 7 d forecast for this domain runs in about 30 min using 16 processors.

3.2 Initialisation of the model fields

Before the model can be run it has to be initialised. (Note this is only done once and is not to be confused with the daily assimilation, which modifies the initial conditions for each new forecast). First, the triangular mesh for the destination domain is generated with Unref (a component of the open-source mesh-generation library GMSH: Geuzaine and Remacle, 2009) and using the Global Self-consistent, Hierarchical, High-resolution Shoreline Database (GSHHS) (Wessel and Smith, 1996). Then model variables like the ice concentration, the ice and snow thicknesses, and the temperature and salinity of the slab ocean are initialised. We use the CS2SMOS product (Sect. 2.5) to set the ice concentration and thickness, and we use the simulated sea surface temperature and salinity from TOPAZ4 (Sect. 2.1). The ice concentration in the CS2SMOS product comes from the low-resolution OSI SAF product (SSMIS), which is known to be too low in areas with low thickness (Ivanova et al., 2015). Therefore total sea ice concentration (c_tot) is calculated by increasing the observed SIC according to the empirical formula:

$$\begin{array}{}\text{(2)}& c_{\text{tot}}={\displaystyle \frac{\mathrm{1}}{a_{\mathrm{0}}}}\text{csch}\left({\displaystyle \frac{h_{\text{cs2smos}}-a_{\mathrm{1}}{a}_{\mathrm{2}}}{\mathrm{2}{a}_{\mathrm{1}}}}\right)\times c_{\text{cs2smos}},\end{array}$$

where c_cs2smos and h_cs2smos are respectively the concentration and effective thickness in the CS2SMOS product, and a₀=0.9569, a₁=0.06787 and a₂=0.4255 are parameters fitted from the observations (Thomas Lavergne, personal communication, 2019).

As mentioned above the model has two ice categories – young ice and older ice with different rheological behaviour. At the initialisation step young ice c_y is set to comprise 20 % and the older ice c 80 % of the total SIC. If the observed absolute thickness $h_{\text{cs2smos}}/c_{\text{tot}}$ is below the young ice upper thickness limit H_max, then thickness is distributed between young and older ice in the same proportion; otherwise young ice thickness is calculated as $h_y=c_y{H}_{max}$ and for the older ice: $h=h_{\text{cs2smos}}-h_y$. It was identified that the model has little sensitivity to the fraction of young ice, and it may vary within 5 %–30 %.

The temperature and salinity of the slab ocean is taken to be equal to the TOPAZ4 surface forecast, while the ice velocity and damage are set to zero.

3.3 Assimilation of concentration

The assimilation is performed before each forecast run using the data insertion method – an updated (analysis) concentration ($c_{\text{tot}}^{\text{a}}$) is calculated as a function of the forecast variable ($c_{\text{tot}}^{\text{f}}$) and observations (c_osisaf), the blended SSMIS/AMSR2 concentration from Eq. (1). Other variables (particularly ice and snow thicknesses and the SST of the slab ocean) are adjusted for consistency, the simulation is then restarted using the updated variables, and the model is run for 8 d to provide 1 d of hindcast and a 7 d forecast. The assimilation is performed on the model mesh – the satellite observations, originally provided on a regular grid, are linearly interpolated onto the centres of the mesh elements so that they can be directly compared with the neXtSIM prognostic variables.

The concentration update is done as follows. A target concentration is calculated using a weighted average approach:

$$\begin{array}{}\text{(3)}& c_{\text{target}}={\displaystyle \frac{c_{\text{tot}}^{\text{f}}{\mathit{\sigma }}_{\text{model}}^{-\mathrm{2}}+c_{\text{osisaf}}{\mathit{\sigma }}_{\text{osisaf}}^{-\mathrm{2}}}{{\mathit{\sigma }}_{\text{model}}^{-\mathrm{2}}+{\mathit{\sigma }}_{\text{osisaf}}^{-\mathrm{2}}}}.\end{array}$$

The uncertainties of the observed concentration (σ_osisaf) are obtained from the input products (the root mean square of the SSMIS and AMSR2 errors), while the value for the uncertainty in the forecast concentration is set to 0.3 and can also be thought of as introducing a timescale (in days) of $\mathrm{1}+{\mathit{\sigma }}_{\text{osisaf}}^{\mathrm{2}}/{\mathit{\sigma }}_{\text{model}}^{\mathrm{2}}$ for relaxation towards the observations that depends on the relative uncertainties. Assimilation of concentration is performed in all elements that have valid observations.

In the pack, the total (young and old ice) model concentration c_tot is generally close to 100 % (except inside leads and cracks), while the OSI SAF concentration can be around 90 %–95 %. We found that using c_target as the actual update $c_{\text{tot}}^{\text{a}}$ produced too much of a drop in the pack concentration, lowering the internal stress to near zero and allowing too much drift. This quickly led to rapid build-up of very thick ice in unusual places. Therefore we had to take a quite conservative approach to our correction and only change the model to make sure the ice mask (the part of the domain where the concentration is higher than 15 %) was correct (this is a kind of assimilation of extent):

$$\begin{array}{}\text{(4)}& {\displaystyle }{c}_{\text{tot}}^{\text{a}}=\left\{\begin{array}{l}\mathrm{0}{c}_{\text{target}}<\mathrm{0.15},\\ c_{\text{target}}{c}_{\text{tot}}^{\text{f}}<\mathrm{0.15}\text{and}{c}_{\text{target}}\ge \mathrm{0.15},\\ c_{\text{tot}}^{\text{f}}\text{otherwise.}\end{array}\right.\end{array}$$

Table 1Accuracy of the free run. Concentration and extent are evaluated against OSI SAF SSMIS; thickness is evaluated against CS2-SMOS; drift is evaluated against OSI SAF drift, where observations with reported error up to 10 km d⁻¹ are considered. Results are 2-monthly-averaged.

Download Print Version | Download XLSX

After calculating the updated variables as specified in Eq. (4) the fractions and thicknesses of young and older ice are calculated in the same way as in the initialisation procedure (see Sect. 3.2).

Once all the assimilation steps have been performed the model fields are checked for consistency with each other and corrected if necessary. First, the concentration of young ice is reduced in the elements where the total concentration exceeds 100 %. Second, the volume of ridged ice, damage, and the ice and snow thicknesses are set to zero in the added ice. Lastly, we need to correct the SST of the slab ocean. Where new ice is added during assimilation it is set to the freezing point, but if ice is removed then the situation is more complicated. Then we proceed as in the following section.

3.4 Compensating heat flux

One of the side effects of assimilation is that the heat balance in the model is disturbed: reducing the concentration opens more water, and depending on the temperatures, atmospheric humidity and ocean salinity provided by the forcing, the heat flux out of the ocean can increase dramatically, causing the ice to freeze up again very fast. This effect is strongest when the atmosphere is very cold. Therefore a compensating heat flux is added to the total ocean heat flux in order to keep the heat balance and prevent refreezing of ice where it was removed by assimilation, thus prolonging the effect of assimilation. If $c_{\text{tot}}^{\text{f}}$ and $c_{\text{tot}}^{\text{a}}$ are respectively the total concentrations before and after assimilation, a compensating heat flux (Q_comp) is calculated if ice was removed – i.e. if $c_{\text{tot}}^{\text{a}}<c_{\text{tot}}^{\text{f}}$ – according to the following formula:

$$\begin{array}{}\text{(5)}& Q_{\text{comp}}=Q_{\text{ocean}}\left({\left(c_{\text{tot}}^{\text{a}}/c_{\text{tot}}^{\text{f}}\right)}^n-\mathrm{1}\right),\end{array}$$

where Q_ocean is the total heat flux from the ocean (sum of flux from the ocean to sea ice and to the atmosphere), and n is a parameter controlling the strength of correction. The function was chosen so that Q_comp is zero when the concentration update is zero, and $Q_{\text{comp}}=-Q_{\text{ocean}}$ when $c_{\text{tot}}^{\text{a}}=\mathrm{0}$; i.e. all the ice was removed by assimilation. With n=1 the reduction of Q_comp from 0 to −Q_ocean is linear, and with n>1 it becomes steeper for values of $c_{\text{tot}}^{\text{a}}$ closer to $c_{\text{tot}}^f$. We use n=4 for the runs presented here, as this gives a suitably strong heat flux for a modest reduction in concentration.

4.1 Evaluation of free model run

In this section we demonstrate that the model is generally able to reproduce the overall drift, concentration and thickness patterns in observations. For all comparisons we average the model fields in time over an appropriate time window (in practice 1, 2 or 7 d), apply some spatial smoothing (being guided by the size of the satellite footprint) and interpolate onto the observation grid. For scalar variables, we define bias as 〈V_mod−V_obs〉 (with 〈⋅〉 defining the spatial mean over pixels where both model and observation are defined, and where either model or observations have ice). We also define RMSE as $\langle (V_{\text{mod}}-V_{\text{obs}}{)}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}$. For the ice extent, we define bias in extent as A₁₀−A₀₁, where A₁₀ is the total area of pixels where neXtSIM predicts the presence of ice (total concentration greater than 15 %) but the observation has no ice, while A₀₁ is the total area of pixels where neXtSIM predicts no ice but the observation does have ice. Instead of an RMSE, for the ice extent we define IIEE (integrated ice edge error, Goessling et al., 2016) as A₁₀+A₀₁. Thus the IIEE is always positive, and the bias is positive if neXtSIM is overestimating the total extent and negative if it is underestimating it. For the ice velocity we define the bias and RMSE as the mean and rms values of the difference in speed respectively – i.e. $\langle \left|{\mathit{u}}_{\text{mod}}\right|-\left|{\mathit{u}}_{\text{obs}}\right|\rangle$ and $\langle \left(\right|{\mathit{u}}_{\text{mod}}|-|{\mathit{u}}_{\text{obs}}|{)}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}$. We also define the vector RMSE (VRMSE) as the rms of the vector difference, $\langle |{\mathit{u}}_{\text{mod}}-{\mathit{u}}_{\text{obs}}{|}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}$. Note that these drift errors can only be calculated where both model and observation have ice.

Table 1 shows 2-monthly summary statistics for observations that are available year-round (although they are much less reliable between May and September). Errors which are particularly high are the RMSE in ice concentration and the IIEE for extent between May and August – we will discuss reasons behind this below. The bias in these quantities is quite low, being highest between November and December, reflecting an ice advance which is too fast. The drift on the other hand compares quite well throughout the period, with the VRMSE being generally less than 4.5 km d⁻¹ and the RMSE in speed less than 3.3 km d⁻¹ (with the exception of the first 2 months). The bias is quite low, staying between ±0.8 km d⁻¹. The first 2 months have a bias of zero, despite their high VRMSE and RMSE, with a Beaufort Gyre that is too slow cancelling a drift that is too fast in the triangle between the north pole, Svalbard and Franz Josef Land.

Table 2Accuracy of the free run for the 2018–2019 and 2019–2020 winters. Thickness is evaluated against CS2-SMOS and drift against OSI SAF drift, where only observations with reported error less than 1.25 km d⁻¹ are considered. Results are 2-monthly-averaged.

Download Print Version | Download XLSX

Table 2 shows the effect of restricting the drift evaluation to regions where the observational error is less than 1.25 km d⁻¹. As discussed in Sect. 2.4 this removes the entire period from May to September, as well as the MIZ, the area around the north pole and areas close to the coast. Without the less accurate observations, the drop in VRMSE is between 0.3 and 1 km d⁻¹, with their values now between 3 and 4.1 km d⁻¹. However, the less accurate regions are also regions where the model has difficulties with other variables like thickness and extent, so this could also contribute to the reduction in error. Notwithstanding this, the November–December 2018 period is now much less of an outlier, now being quite close to November–December 2019. The RMSE in speed is now less than 2.6 km d⁻¹, while the bias is between ±0.3 km d⁻¹. There is a clear deterioration going from the first to the second winter.

Table 2 also summarises the comparison of model thickness with the CS2-SMOS product. RMSE is lower (0.26 m) in November–December 2018 since thickness was initialised using CS2-SMOS and grows to 0.51 m over the first winter. The RMSE grows from 0.4 to 0.56 m over the second winter. The bias is quite low in the first winter but about 0.2 m from January–April 2020.

Figure 2Temporal comparison of model and OSI SAF SSMIS concentrations. The shaded area in the lower plot shows the rms uncertainty of the concentration for reference when considering the errors. The mean concentration is the mean over the ocean points in the domain, while the bias is defined as the mean error (model minus observation) over the region where either the model or observations have ice.

Download

Figures 2 and 3 compare the mean concentration and extent from neXtSIM and the OSI SAF concentration product (see Sect. 2.3). Since we do not have enough information to create an error model for the concentration observations (for example it cannot be Gaussian due to its restriction to being between 0 % and 100 %) we could not generate confidence intervals for the mean concentration but instead plot a shaded area corresponding to $\pm \langle {\mathit{\sigma }}_{\text{OSISAF}}{\rangle }^{\mathrm{1}/\mathrm{2}}$ for reference when considering the bias and RMSE in concentration. The seasonal cycle for mean concentration and extent is captured reasonably well, but the bias is clearly positive in winter and clearly negative in summer. As noted from Table 1 the RMSE in concentration and the IIEE are very high between mid-April and August. The RMSE in concentration is generally comparable to the rms error level outside these periods.

Figure 3Temporal comparison of extents (region where the concentration exceeds 15 %) from model and OSI SAF SSMIS concentrations. The bias is defined as modelled extent–observed extent.

Download

In January–February there is a general overestimation in the pack as the model is approximately 100 % with some leads, while the SSMIS field has large regions of lower concentration that contribute to the errors (85 %–90 %). There is also a general overestimation in ice extent, which is also the case in November–December, primarily located in the Greenland Sea and the Barents seas. The areas north of Novaya Zemlya (around the Santa Anna Trough), around Franz Josef Land, and northwest of Svalbard stood out as places where the ice edge was not being located correctly. This continues into March–April and can be seen in Fig. 4, which shows selected maps of the concentration bias from March–October 2019. The region of underestimation to the northwest of Svalbard is related to the development of an arch between there and the northeast tip of Greenland. (This is also related to relatively large thicknesses in these areas, which increases the capacity of undamaged ice to stay undamaged.) This arching also reduces the ice export through the Fram Strait, the effect of which is seen in the maps for May–June and July–August. In May–June we also see that the ice that is in the Greenland Sea and is displaced too far to the east, producing a double penalty effect in the RMSE score. This is also a result of the arching in this area – the ice at the corner of Greenland is too thick to be exported, and so the ice that is exported has detached from the arch away from the coast and has then travelled roughly parallel to the coast without being pushed back towards it. In July–August the thicker ice has mostly melted and there is less of a dipole situation and more of a clear-cut underestimation.

Another problem area is around the Novosibirsk Islands (also known as New Siberian Islands) and the Laptev Sea. April sees the model opening slightly too far to the north; in May and June there is a strong underestimation to the north of the islands and too much fast ice. The underestimation is partly from melting that is too fast and from the fact that the fast ice has not detached and flowed into this area.

There is a similar problem to the west of Wrangel Island in the Chukchi Sea – an opening that is too large compared to the observations and with too much land-fast ice close to the coast which should be flowing into this region. There is also a large region of overestimation to the east of Wrangel Island, which is actually an artefact of the boundary conditions that we are using at the open Bering Strait boundary. If there is inflow at any open boundaries the value of any tracers is taken to be its value in the nearest mesh element, and in this case we have too much ice being imported through the Bering Strait, which leads to a build-up at this location.

The last things to mention about the July–August map are the strong overestimation around Franz Josef Land and the dipole situation in the Beaufort sea, which is the result of a Beaufort Gyre that is too slow.

However, by September–October, the situation has improved substantially. The main disagreement is in the ice edge from Severnaya Zemlya round to the just past Svalbard. This is due to an ice advance that is slightly too fast.

Figure 4Selected maps of the 2-monthly-averaged concentration biases between the free model run and from OSI SAF SSMIS. The bias is defined as model minus observation.

Figure 5Temporal comparison of model and CS2-SMOS thickness for two winters, 2018–2019 (a, c) and 2019–2020 (b, d). The shaded regions show the rms uncertainty in the CS2-SMOS product for reference. The mean thickness is the mean over the ocean points in the domain, while the bias is defined as the mean error (model minus observation) over the region where either the model or observations have ice.

Download

Figure 6Comparisons of time-averaged model and CS2-SMOS thickness for early (a–c) and late (d–f) in the 2019–2020 winter. (a, d) Bias maps; (b, e) neXtSIM thickness; (c, f) CS2-SMOS thickness. The bias is defined as model minus observation.

Figure 5 shows time series of mean thicknesses and of thickness errors when compared to CS2-SMOS, which have already been discussed to some extent in the context of Table 2. The mean observed and modelled thicknesses show a steady increase after November, but the modelled value is increasing faster than the observed one. The timing of when the modelled increase starts is close to when the observed increase starts. The lower plots have the rms uncertainties plotted for reference, and the model bias is generally at a similar level to this uncertainty. The RMSE is about twice this uncertainty in the second winter. We note here that the error levels in the CS2-SMOS product are only the interpolation error and are thus a lower bound as they do not include uncertainties in the individual CS2 and SMOS products. CS2 in particular is sensitive to the ice and snow densities used or the snow thickness which affect the conversion from freeboard to thickness (Zygmuntowska et al., 2014).

Figure 6 also shows the spatial distribution of the errors. Throughout the winter there is overestimation off the north coast of Greenland round to Axel Heiberg Island. Further to the west there is a thinning (but not an opening) from Axel Heiberg Island to Ellef Ringnes Island, which also seems to be related to the westward drift along this coast being too high. This persists throughout the winter, but the affected area reduces with time. In October–December there is a dipole pattern where the ice is too thin in the Beaufort Sea and too thick in the triangle between the north pole, Svalbard and Franz Josef Land. However, in January–April, this dipole reverses.

The modelled thickness is also quite high around the north of Svalbard, in contrast with the CS2-SMOS product, which has quite thin ice. This build-up contributes to the arching discussed in relation to the concentration errors, reducing the ice export through the Fram Strait. This stems from damaged ice not having enough resistance to compression. In the BBM rheology there is a balance between resisting compression enough to stop the build-up of ice at the coasts and resisting it so much that drift becomes too slow.

Figure 7Temporal comparison of model and OSI SAF drift. The shaded region shows the rms uncertainty in the OSI SAF drift product. The bias is defined as the mean modelled speed minus observed speed. The mean drift and errors are averages over the regions where both model and observations have valid drift vectors. Valid observations here are those that have reported errors less than 10 km d⁻¹.

Download

Figure 7 shows the drift bias and RMSE of neXtSIM when compared to the OSI SAF drift product. The shaded area corresponds to $\langle \mathrm{2}{\mathit{\sigma }}_{\text{OSI SAF}}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}$, where the factor of 2 comes from the fact that drift has two components. We also note that if one cell has components (x_i,y_i) with (not-necessarily Gaussian) noise (ε_i,δ_i) added to it (each with zero mean and variance ${\mathit{\sigma }}_i^{\mathrm{2}}$), the mean-squared drift in the presence of noise,

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\langle {\stackrel{\mathrm{̃}}{d}}^{\mathrm{2}}\rangle & ={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\langle (x_i+{\mathit{\epsilon }}_i{)}^{\mathrm{2}}+(y_i+{\mathit{\delta }}_i{)}^{\mathrm{2}}\rangle \\ & ={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\langle x_i^{\mathrm{2}}+y_i^{\mathrm{2}}+{\mathit{\epsilon }}_i^{\mathrm{2}}+{\mathit{\delta }}_i^{\mathrm{2}}\rangle \\ & ={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left(x_i^{\mathrm{2}}+y_i^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}\right)\\ & =\langle d^{\mathrm{2}}\rangle +\mathrm{2}\langle {\mathit{\sigma }}_{\text{OSI SAF}}^{\mathrm{2}}\rangle ,\end{array}\end{array}$$

is higher than the noise-free (“true”) value 〈d²〉 by $\mathrm{2}\langle {\mathit{\sigma }}_{\text{OSI SAF}}^{\mathrm{2}}\rangle$. It is tempting to then just subtract this amount from the speed when comparing with the model, but we decided that the uncertainty in the quoted uncertainty could be too influential and also that the model drift has some unknown uncertainty associated with it, so we persisted with directly comparing the drift in the OSI SAF product with the model drift in order to judge when the model is too fast or too slow and plotted the error level $\langle \mathrm{2}{\mathit{\sigma }}_{\text{OSI SAF}}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}$ for reference. The difference in speed (model speed minus observed speed) fluctuates somewhat but stays largely within the error limits. The rms error for this product is approximately 1.8 km d⁻¹ from October to May but increases to about 4 km d⁻¹ from May to September. However, the model is starting to show signs of being too slow in April and May 2020. The RMSE in speed ranges between about 3–4.5 km d⁻¹, which is outside the estimated error for the more accurate colder months. The VRMSE (vector RMSE) is higher than the RMSE by about 1.5 km d⁻¹.

Figure 8Comparisons of 2-monthly-averaged model and OSI SAF drift from the 2019–2020 winter. The drift bias colour maps (left column) show the bias in speed, while the directions show the difference between the model and observation directions (arrows pointing up indicate the directions are the same). The central column shows ice velocity vectors from neXtSIM, while the right-hand column shows OSI SAF vectors. The concentration colour maps show the average neXtSIM concentration over the analysis period. The bias is defined as modelled speed minus observed speed.

Figure 9Two-monthly-averaged forecast errors grouped by lead time (lead times are indicated in the figure legends). The free run errors for the corresponding periods are plotted as dotted lines for reference on the bias and RMSE/IIEE plots. (a–c) Two-monthly-averaged bias, RMSE and forecast skill for concentration grouped by lead time. (d–f) Two-monthly-averaged bias, IIEE and forecast skill for extent grouped by lead time.

Download

Figure 8 shows the general spatial pattern in the drift errors. The maps for October–December 2019 show that the general circulation of the model and the observations are agreeing well, with the main feature being the gyre in the central Arctic Ocean. However, this gyre is slightly too fast and there is also too much westward drift along the Canadian Arctic Archipelago. The observations also show a small gyre in the East Siberian Sea which is not apparent in the model drift. The strongest positive bias is in the Kara Sea although the errors in the observations are high there, being close to the coast and the MIZ. The model is also showing strong positive biases in the Laptev and East Siberian seas where the ice is thinner. On the other hand, the drift in the Beaufort and Chukchi seas is too slow.

In January–March 2020, the general pattern in the observations is a strong flow from the Laptev over to the Greenland Sea, with a weaker transpolar drift from the Chukchi Sea joining this outflow. The model also produces the flow from the Laptev Sea, but the flow from the Chukchi Sea is more directed towards the Beaufort Sea than across to the Greenland Sea. This may be contributing to the biases seen in the thicknesses – too thick in the Beaufort Sea while too thin in the region between the Laptev and Greenland seas. The drift is also starting to become too slow north of Greenland and Svalbard, and this is most likely due to the ice build-up in those areas.

In April–June 2020 there is a stronger Beaufort Gyre and a similarly strong cyclonic gyre in the Laptev Sea. These connect to form quite a wide stream from the central Arctic to the Greenland Sea. While the model captures the gyre in the Laptev Sea, the Beaufort Gyre is less well defined in this period, possibly due to the ice being too thick in the Beaufort Sea at this time. As in January–March, drift is also too slow to the north of Greenland.

Figure 10Selected maps of daily-averaged biases in concentration compared to OSI SAF SSMIS for three examples of forecasts, starting on 3 January 2019 (a–c, positive skill) and the other on 7 March 2019 (d–f, negative skill) and 6 September 2019 (g–i, positive). Panels (a), (d) and (g) show the evaluation for the first day of simulation (lead time −1: analysis or hindcast), while panels (b), (e) and (h) show the evaluation for the eighth day of simulation (lead time 7). Panels (c), (f) and (i) show the bias in the persistence forecast (defined as the results for the first day of simulation) for lead time 7. The bias is defined as model minus observation.

4.2 Evaluation of forecasts with assimilation

The performance of the forecasting system with assimilation of concentration was evaluated over the same period as the free run was in Sect. 4.1 (that is the 20 months from November 2018 to the end of June 2020). In order to reduce computation time, the 7 d forecasts were launched approximately only every 7 d, with 1 d analyses being launched in between so that the assimilation was still performed daily.

The model performance was evaluated against satellite observations of the sea ice SSMIS concentration and drift from OSI SAF (see Sect. 2.3 and 2.4). The metrics we used in our assessment were bias and RMSE in concentration; bias and IIEE in extent; and bias, RMSE and VRMSE in drift. A persistence forecast of concentration (the initial concentration, defined for technical convenience as the average of the first hour) was used as a benchmark for concentration and extent, and forecast skill was defined as

$$\begin{array}{}\text{(7a)}& {\displaystyle }& {\displaystyle }\text{Concentration skill}=\mathrm{1}-{\left({\displaystyle \frac{\text{RMSE}}{{\text{RMSE}}_{\text{persistence}}}}\right)}^{\mathrm{2}},\text{(7b)}& {\displaystyle }& {\displaystyle }\text{Extent skill}=\mathrm{1}-{\displaystyle \frac{\text{IIEE}}{{\text{IIEE}}_{\text{persistence}}}}.\end{array}$$

Thus, a model that is agreeing perfectly with the observations has skill equal to 1, while it is negative (with no lower bound) if the persistence error is lower than the forecast error. We decided that given the strong dependence of drift on wind, its rapid variability in time would render a persistence forecast for drift quite easy to beat, so for this variable we just present the forecast errors. However, we can make some rough comparisons to the drift errors of other products, even though they use different observations for their evaluation so we cannot make a direct comparison. We note that the drift from the TOPAZ forecast generally has a bias in speed of about 2 km d⁻¹ and a VRMSE of about 5–8 km d⁻¹ Melsom et al. (2018), while Metzger et al. (2017) report an rms drift speed error of about 5–8 km d⁻¹ in the Arctic for the GOFS 3.1 system.

Figure 11Maps of daily-averaged biases (a, d) in concentration compared to OSI SAF SSMIS, and modelled (b, e) and observed (c, f) concentrations, for an example forecast starting on 6 June 2019. Panels (a)–(c) show the comparison for the first day of simulation, while panels (d)–(f) show the comparison for the eighth day of simulation. The bias is defined as model minus observation.

Figure 9 shows the 2-monthly averages for the different metrics applied to the concentration and extent. We only plot results for the first year and 2 months, since the error statistics from the 2020 results were nearly identical to the previous years. As with the free run, the concentration from the forecast is generally higher than OSI SAF in the winter months. As mentioned earlier, we do not try to correct this as reducing the concentration in the pack causes serious problem with the drift and thickness. The extent, which is adjusted during assimilation, is also lower than the observed one for the whole period, by differing amounts. Figure 10 shows some example forecasts from winter and autumn – one from January, March and September 2019 respectively. The first and third examples are forecasts with positive skill in extent (model does better than persistence), while the second has negative skill. However, they all have relatively low IIEE and RMSE in concentration. The January example shows that the rapid advance in the Kara Sea is captured relatively well by the model, and while the extent in the Greenland and Barents seas is too high, the ice edge is quite variable, and the model is doing better than the persistence at this time. The March example shows an example where there is quite a persistent negative bias to the northwest of Svalbard and in the Barents Sea. These errors probably originate in the ocean forcing, which we are unable to overcome even with our flux compensation. The persistence forecast does well in this example since the ice conditions are not changing very much. In the September example, the bias is similar to the free run, but it has been reduced by the assimilation. There is still an ice advance that is too fast in the European Arctic, but it is closer to the observations than the persistence forecast.

Figure 12Two-monthly-averaged bias, RMSE and VRMSE for drift grouped by lead time (lead times are indicated in the figure legends). The lead times refer to the start of each 2 d drift evaluation period – e.g. “0” covers the period from day 0 (12:00) to day 2 (12:00). The free run errors for the corresponding periods are plotted as dotted lines for reference. The bias is defined as modelled speed minus observed speed.

Download

In the summer, the forecasts are systematically lower in concentration, partly due to differences in the pack, which we do not correct for, and in the extent. We are not too concerned by differences in the pack, as the differences are still comparable to observation error. Differences in the extent however are more concerning. In the summer, the forecasts score particularly badly with this metric, largely due to the dynamical problems also seen in the free run: arching above the Fram Strait, as well as fast ice in the Laptev and East Siberian seas that is too slow to detach. Figure 11 shows a typical example from this time period. Despite the correction in extent during assimilation, the ice that is added in the Fram Strait is quite thin and it melts quite quickly.

Figure 12 shows the bias, RMSE and VRMSE for the forecast drift. For smaller lead times it has similar quality to the free run, which is to be expected since the assimilation does not modify fields in the pack, with minor differences that are probably just due to differences in extent determining which observations are included or not. There is a small negative bias in speed in the winter months and a slight positive bias in the summer. The bias is not showing a large variability with forecast lead time, although there are some exceptions. With the exception of the first 2 months, which is slightly higher than the others, the RMSE (RMSE in speed) is quite consistent throughout the whole evaluation period, with a steady increase with lead time. The RMSE for the final day is still mostly less than 5 km d⁻¹. The VRMSE, which evaluates the drift direction as well as its speed, shows a larger dependence on lead time, starting between 4.5–5 km d⁻¹ but getting up to between 7–8 km d⁻¹ by the final day (occasionally being under 7 km d⁻¹ or over 8 km d⁻¹). This dependence on lead time is almost certainly due to increasing inaccuracy in the forecast winds.