Edge displacement scores

As a consequence of a diminishing sea ice cover in the Arctic, activity is on the rise. The position of the sea ice edge, which is generally taken to define the extent of the ice cover, changes in response to dynamic and thermodynamic processes. Forecasts for sea ice expansion on synoptic timescales due to an advancing ice edge will provide information that can be of significance for open ocean operations in polar regions. However, the value of this information depends on the quality of the forecasts. Here, we present methods for examining the quality of forecasted sea ice expansion on sub-seasonal timescales and the geographic location where the largest expansions are expected from the forecast results. The algorithm is simple to implement, and an examination of 2 years of model results and accompanying observations demonstrates the usefulness of the analysis.

1 Introduction

Due to climate change the sea ice extent is in decline in the Arctic (Parkinson, 2014). This change has led to increased activity in the region, and commercial shipping in open waters via Arctic sea routes will become increasingly economically viable in the 21st century (Aksenov et al., 2017). Thus, data sets for monitoring and forecasting sea ice conditions are receiving growing attention.

The past years have seen a flurry of activity related to assessing the quality of sea ice data sets. Dukhovskoy et al. (2015) presented a review and comparison of various traditional metrics for the assessments of the skill of sea ice models. Goessling et al. (2016) introduced the integrated ice edge error (IIEE), a quantity for describing mismatching sea ice extents from two data sets to analyze the predictability of the sea ice edge. Melsom et al. (2019) took advantage of the IIEE in their examination of various metrics for the assessment of the quality of forecasts for the sea ice edge position. Methods for examining the quality of probabilistic results for sea ice conditions have been introduced by Goessling and Jung (2018) and Palerme et al. (2019). Recently, Cheng et al. (2020) examined the accuracy of a visually estimated ice concentration monitoring product.

The changing position of the sea ice edge is generally not only shifted by dynamic advection but can be significantly affected by the thermodynamics as well (Bitz et al., 2005). Thus, the temporal displacement of the sea ice edge will be affected by freezing along the perimeter of the sea ice extent in winter and melting in summer. Hence, pattern-recognition algorithms for displacements using maximum cross-correlation (MCC) methods such as those introduced by Leese et al. (1971) for wind vectors, and later for ocean surface currents (Tokmakian et al., 1990) and sea ice vectors (Lavergne et al., 2010), are not ideal for tracking displacements of the sea ice edge.

Ebert and McBride (2000) examined the position error of the contiguous rain area in weather forecasts. They determined the error vector from a total squared error minimization method when shifting the forecasted rain region to match the corresponding observations. Their preference of applying an error minimization algorithm rather than an MCC approach was motivated by the former having better representations of displacement of rainfall maxima. An object-based approach with a focus on assessing the quality of forecasts for highly localized and episodic phenomena was introduced by Davis et al. (2006). Like in Ebert and McBride (2000), their focus when evaluating displacement is on the objects' centroid. Displacement of the perimeter of the contiguous rain area was not addressed in either investigation.

We begin this study by presenting a new algorithm for assessing the quality of representations of the sea ice edge by comparing results for two different data sets. We examine displacements over time of the sea ice edge and compare model results for displacement distances vs. observed displacements. The method is described in Sect. 2 with an idealized case study. In Sect. 3 we apply the algorithm to analyze displacements of the sea ice edge in the Barents Sea. Finally, we provide our concluding remarks in Sect. 4. Technical details and extensions are provided in two appendices.

2 Methods

In order to illustrate the validation metrics that are introduced in this section, a set of idealized ice edges is introduced, as depicted in Fig. 1. The domain is divided into 1000×500 square grid cells, and we set the length of the side of a grid cell to 1. We denote the curve that separates regions with binary values 0 and 1 as an edge curve. Let L^O(t) and L^M(t) denote observed and modeled edges, respectively, at time t. Idealized examples with edges for L^O and L^M at two different times, t₀ and t₀+Δt, are displayed. In the context of forecasting, L^M(t₀) may be taken to represent the model initialization at t₀, and L^M(t₀+Δt) is then the forecast at a temporal range of Δt. The other binary fields can represent observations at the same times.

Figure 1 Binary fields with values of 1 (ice) and 0 (no ice/ocean) are displayed by white and blue color shading, respectively. Light shades of blue indicate regions with a non-overlapping ice cover, as indicated by the color legend. The derived modeled and observed ice edges L^M and L^O at $t=t_{\mathrm{0}}+\mathrm{\Delta }t$ are drawn as red and black curves, respectively. The corresponding ice edges that are taken to represent the situation at t₀ are drawn as light red and gray curves. The full black circle indicates the position on the observed ice edge at t₀+Δt which has the largest distance to the ice edge at t₀ (the full gray circle). The largest displacement of the model ice edge is marked by full diamonds. The full red circle is the position along the model ice edge at $t=t_{\mathrm{0}}+\mathrm{\Delta }t$ closest to the full black circle, while the full light red circle is the position of the observed ice edge at t₀ closest to the full red circle. All dashed lines represent displacements as defined by Eq. (3). Open circles indicate a random selection of displacement positions for the model results; see the text for details.

2.1 Validation metrics for ice edge displacement in one data set

We aim at defining metrics that describe differences in maximum sea ice expansion from sea ice edge displacements between two data sets. In order to do so, we must first introduce a quantity that properly measures the maximum displacement in one data set. Here a definition is provided which is a gridded, signed, one-sided variation of the Hausdorff distance (Dukhovskoy et al., 2015).

For the remainder of this investigation we will take the binary fields to be representations of sea ice, with values assigned to 0 and 1 for conditions of no ice and ice, respectively. We will here associate the presence of ice (value 1) with sea ice concentration c exceeding c_edge=0.15. In a gridded representation the ice edge can then be taken to be constituted by the grid cells $e=[i,j]$ that meet the condition

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }c[i,j]\ge c_{\text{edge}}\wedge \\ \text{(1)}& {\displaystyle }& {\displaystyle }min\left(c\right[i-\mathrm{1},j],c[i+\mathrm{1},j],c[i,j-\mathrm{1}],c[i,j+\mathrm{1}\left]\right)<c_{\text{edge}},\end{array}$$

where ∧ is the logical AND operator. Denoting the N(t) grid cells that satisfy this condition for time t by $e_{\mathrm{1}}^{\left(t\right)},e_{\mathrm{2}}^{\left(t\right)},\mathrm{\dots },e_N^{\left(t\right)}$, the ice edge for time t is then the curve

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }L\left(t\right)=\left\{e_{\mathrm{1}}^{\left(t\right)},e_{\mathrm{2}}^{\left(t\right)},\mathrm{\dots },e_{N\left(t\right)}^{\left(t\right)}\right\}.\end{array}$$

This follows the algorithm presented in Melsom et al. (2019). Let L(t₀) and L(t₀+Δt) denote the sea ice edges at times t₀ and t₀+Δt, respectively. Furthermore, let $d_n^{\mathrm{\Delta }t}$ be the displacement distance between a grid cell $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ in L(t₀+Δt) and curve L(t₀), i.e.,

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{d}_n^{\mathrm{\Delta }t}=s_{n}min\left|\left|e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}-L\left(t_{\mathrm{0}}\right)\right|\right|.\end{array}$$

Here, $\left|\right|z\left|\right|$ is the Euclidean distance of z, and the minimum is considered for the distances between grid cell $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ and each of the grid cells belonging to L(t₀). Furthermore, s_n is +1 or −1 when $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ is on the no ice or ice side of L(t₀), respectively, i.e.,

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{s}_n=\left\{\begin{array}{rl}-\mathrm{1}& \text{if }c\left[e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}\right]\left(t_{\mathrm{0}}\right)\ge c_{\text{edge}}\\ & \\ +\mathrm{1}& \text{if }c\left[e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}\right]\left(t_{\mathrm{0}}\right)<c_{\text{edge}}\end{array}\right.,\end{array}$$

where $c\left[e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}\right]\left(t_{\mathrm{0}}\right)$ denotes the sea ice concentration at the time t₀.

Figure 1 shows an idealized example in which a modeled sea ice edge and an observed sea ice edge are displaced. Presently, we consider metrics for one product, and as an illustrative example, we focus on the ice edges derived from the model product (light red and red lines). The length of dashed lines in Fig. 1 then corresponds to model displacements $min\left|\right|e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}-L\left(t_{\mathrm{0}}\right)\left|\right|$ for selected cells $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$. We introduce the maximum expansion displacement as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{d}_{\text{max}}^{\mathrm{\Delta }t}=max\left(d_n^{\mathrm{\Delta }t}\right),n=\mathrm{1},\mathrm{2},\mathrm{\dots },N(t_{\mathrm{0}}+\mathrm{\Delta }t).\end{array}$$

Note that the definition of the sign s in Eq. (3) has been chosen so that Eq. (5) will return the largest positive value among $d_n^{\mathrm{\Delta }t}$. If all values of $d_n^{\mathrm{\Delta }t}$ are negative, the result is the negative value distance with the lowest magnitude. The definition of s was designed so that $d_{\text{max}}^{\mathrm{\Delta }t}$ will represent the displacement of the largest sea ice advance from L(t₀) to L(t₀+Δt).

If we briefly introduce $d_m^{-\mathrm{\Delta }t}$ as the shortest distance from a grid cell e_m of L(t₀) to the curve L(t₀+Δt), we note that the Hausdorff distance d_H between curves L(t₀+Δt) and L(t₀) is

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{d}_{\mathrm{H}}=max\left(\left|d_n^{\mathrm{\Delta }t}\right|,\left|d_m^{-\mathrm{\Delta }t}\right|\right).\end{array}$$

Here, d_H is symmetric with respect to the distance when swapping the two compared data sets (usually an observed and a modeled feature). The corresponding distance measure introduced in Eq. (5), on the other hand, is asymmetric by design: distances are computed from the grid cells of the updated ice edge at t₀+Δt to the ice edge curve from the initialization at t₀.

For the various quantities we reference model results and observations by superscripts M and O. For the set of binary fields depicted in Fig. 1, we find that $d_{\text{max}}^{\mathrm{M};\mathrm{\Delta }t}=\mathrm{113.2}$ (the distance between the red and light red diamonds in the figure), while $d_{\text{max}}^{\mathrm{O};\mathrm{\Delta }t}=\mathrm{97.9}$ (the distance between the black and gray full circles).

The maximum distance in Eq. (5) provides a single measure to examine the ice edge displacement. However, it can be more informative to analyze the whole distribution of the displacements $d_n^{\mathrm{\Delta }t}$ defined by Eq. (3) rather than their maximum only. This can be done by inspecting a histogram of the displacements $d_n^{\mathrm{\Delta }t}$ (Fig. 2). Another option is to present the cumulative probability distribution of $d_n^{\mathrm{\Delta }t}$ (Fig. 3).

To avoid inflating the sample size beyond its degrees of freedom, displacement results can be subsampled at the spatial decorrelation length along the ice edge. For the distribution of $d_n^{\mathrm{\Delta }t}$ a proper decorrelation length can be computed if the edge cells $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ are in sequence along L(t₀+Δt). A detailed description for this procedure is given in Appendix A. In Sect. 3.3 we will present results for the time series $d_{\text{max}}^{\mathrm{\Delta }t}\left(t\right)$ derived from subsampling based on decorrelation.

Figure 2 Histograms for the distribution of displacement distances computed from Eq. (5) for the ice edges displayed in Fig. 1. The mean displacement distances for model results and observations are 82 grids and 61 grids, respectively. The corresponding median values are 88 grids and 70 grids, respectively.

Figure 3 Cumulative distributions of the separation (ice edge displacement) distances from t₀ to t₀+Δt for model results (red curve) and observations (black curve). Shown here are results for the idealized example displayed in Fig. 1, with distances subsampled at intervals of the decorrelation lengths, which are 42 and 38 grid cells along the ice edge for the model results and observations, respectively. The mean separation distance difference for the present subsample of ice edge grid cells is the integral of the area between the curves, here displayed by gray shading. In this case, the mean difference is 20.9 in grid cell units with larger displacement values in model results than from observations.

2.2 Comparison of the displacements of modeled and observed ice edges

In the previous section we focused on metrics which describe the displacement of a single (modeled or observed) sea ice edge. In this section we extend these to assess the differences in the displacements of the modeled vs. observed ice edge. For this purpose, binary fields that are taken to represent observations, as well as model results, are analyzed, as displayed in Fig. 1.

We aim to compare model displacement distances with the corresponding displacements from observational data. Here, we discuss results for the idealized case depicted in Fig. 1.

In Fig. 2 the histograms for all displacements $d_n^{\mathrm{\Delta }t}$ are shown, with the histograms to the left and right representing model results and observations, respectively. The results have been binned into categories that each span distances of 20 grid cell units. Displacement distances are not evenly distributed as they have about half of the distance values in one of six categories. This specific category is shifted by one category from observations to model results, indicating an overprediction of the sea ice displacement.

Next, the cumulative distributions of displacement distances are displayed in Fig. 3. Here, the results were subsampled by the decorrelation length along the edge curves prior to the analysis. We find that the model displacements are shifted approximately 20 grid units higher for the entire distribution as the two curves are nearly parallel.

Hence, Figs. 2 and 3 both reflect the distributions' near uniform shift of 20 grid units and show that this shift is not qualitatively impacted from subsampling at the decorrelation length.

From the perspective of an observer, a useful attribute is the quality of the forecasted maximum displacement of the binary field over the forecast period. A simple metric which provides this type of information is the difference in the maximum displacement as given by Eq. (5), i.e.,

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{d}_{\text{max}}^{\mathrm{\Delta }t}=d_{\text{max}}^{\mathrm{M};\mathrm{\Delta }t}-d_{\text{max}}^{\mathrm{O};\mathrm{\Delta }t},\end{array}$$

where $d_{\text{max}}^{\mathrm{O};\mathrm{\Delta }t}$ is computed from observed ice edges at t₀ and t₀+Δt (black and gray curves in Fig. 1, respectively), and $d_{\text{max}}^{\mathrm{M};\mathrm{\Delta }t}$ is computed from the corresponding model results. For the results in the idealized example that was introduced in Sect. 2.1, the model is overestimating the maximum displacement by $\mathrm{\Delta }{d}_{\text{max}}^{\mathrm{\Delta }t}=\mathrm{15.3}$ grid cell units.

A similar quantity that provides local information is the local difference in displacement of the model ice edge in proximity of the maximum displacement found in the observations. Let $e_{\mathrm{0}}^{\mathrm{O};(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ be the position in L^O(t₀+Δt) to which the maximum edge displacement is found in the observations. Then, determine ${\mathit{ϵ}}_{\mathrm{0}}^{\mathrm{M};(t_{\mathrm{0}}+\mathrm{\Delta }t)}$, the model edge grid cell positioned closest to $e_{\mathrm{0}}^{\mathrm{O};(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ at the same time. In Fig. 1, the positions $e_{\mathrm{0}}^{\mathrm{O};(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ and ${\mathit{ϵ}}_{\mathrm{0}}^{\mathrm{M};(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ are indicated by the full black and full red circles, respectively. Following Eq. (3) the corresponding local edge displacement in the model results is

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}=smin\left|\left|{\mathit{ϵ}}_{\mathrm{0}}^{\mathrm{M};(t_{\mathrm{0}}+\mathrm{\Delta }t)}-L^{\mathrm{M}}\left(t_{\mathrm{0}}\right)\right|\right|.\end{array}$$

For the idealized example, we find that ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}=\mathrm{83.9}$. The local difference in displacement between model and observations, with reference to the position $e_{\mathrm{0}}^{\mathrm{O}(t_{\mathrm{0}}+\mathrm{\Delta }t)}$, becomes

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{\mathit{\delta }}_{\text{max}}^{\mathrm{\Delta }t}={\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}-d_{\text{max}}^{\mathrm{O};\mathrm{\Delta }t}.\end{array}$$

We recall from Sect. 2.1 that $d_{\text{max}}^{\mathrm{O};\mathrm{\Delta }t}=\mathrm{97.9}$, so, for the idealized example, we have $\mathrm{\Delta }{\mathit{\delta }}_{\text{max}}^{\mathrm{\Delta }t}=-\mathrm{14}$, i.e., a local underestimation of the displacement in the model results.

One aspect which is not disclosed by the metrics introduced thus far is to what degree forecasts manage to reproduce the geographical location of the observed maximum displacements. In order to examine such a relation, we first compute the decorrelation length of displacements given by Eq. (3). If we denote this grid distance by Δn, we restrict the analysis of grid cells and corresponding displacements to

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\left\{\mathrm{\dots },{\mathit{ϵ}}_{\mathrm{0}-\mathrm{2}\mathrm{\Delta }n}^{\mathrm{M}},{\mathit{ϵ}}_{\mathrm{0}-\mathrm{\Delta }n}^{\mathrm{M}},{\mathit{ϵ}}_{\mathrm{0}+\mathrm{\Delta }n}^{\mathrm{M}},{\mathit{ϵ}}_{\mathrm{0}+\mathrm{2}\mathrm{\Delta }n}^{\mathrm{M}},\mathrm{\dots }\right\}(t_{\mathrm{0}}+\mathrm{\Delta }t),\text{(11)}& {\displaystyle }& {\displaystyle }\left\{\mathrm{\dots },{\mathit{\delta }}_{\mathrm{0}-\mathrm{2}\mathrm{\Delta }n}^{\mathrm{M};\mathrm{\Delta }t},{\mathit{\delta }}_{\mathrm{0}-\mathrm{\Delta }n}^{\mathrm{M};\mathrm{\Delta }t},{\mathit{\delta }}_{\mathrm{0}+\mathrm{\Delta }n}^{\mathrm{M};\mathrm{\Delta }t},{\mathit{\delta }}_{\mathrm{0}+\mathrm{2}\mathrm{\Delta }n}^{\mathrm{M};\mathrm{\Delta }t},\mathrm{\dots }\right\},\end{array}$$

respectively, limited by the first and last cells along the curve L^M(t₀+Δt). Next, we construct bins analogous to the method used for producing rank histograms (Talagrand diagrams) for ensemble forecasts (Hamill, 2001; Talagrand et al., 1997; Anderson, 1996): first, distances listed in Eq. (11) are sorted by increasing values, and then bins are introduced for values smaller than the minimum distance, for the intervals between the sorted distances, and for values larger than the maximum distance. The bin placement of ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}$ then gives the rank of this displacement. With a perfect model, the maximum of δ^M;Δt will occur at the geographical position corresponding to the maximum displacement in the observations ($d_{\text{max}}^{\mathrm{O};\mathrm{\Delta }t}$). This will place ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}$ in the highest (rightmost) bin. A flat histogram indicates a model with no skill, since each bin is equiprobable for random draws.

In the present idealized example we find that Δn=42, and the rank of ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}$ in the 24 resulting bins is 9. When multiple forecasts are examined, the decorrelation length will generally change, as will the length of the edges. Thus, we randomly subsample a fixed number of intervals from Eq. (11) so that the number of bins is equal across different cases and results can be aggregated.

For the idealized example, a set of nine randomly subsampled edge positions from those given by Eq. (10) for the model results at $t=t_{\mathrm{0}}+\mathrm{\Delta }t$ is displayed by open circles in Fig. 1. For this particular case, in the range from 1 to 10 the rank of the displacement ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}$ is 3. So, in this idealized, synthetic example, the model exhibits a poor positioning of the maximum observed displacement.

Figure 4 Map of the full SVIM simulation domain. The Barents Sea analysis region in the present study is shown as a highlighted region where a sample sea ice concentration distribution is depicted. The 40^∘ E meridian which will subsequently be used for dividing the domain into two parts is displayed by the red curve. The shading of ice concentration values is given in the label bar, where c is the sea ice concentration fraction. This sample shows the model results for 15 April 2000 with the horizontal resolution from the SVIM experiment. See the text for details about the SVIM archive.

3 Application of the new validation method

3.1 Description of sea ice data sets

To illustrate the methodology introduced in Sect. 2, we examine model results from a coupled ocean–sea ice model and compare them with relevant observational data. The model results are taken from the SVIM hindcast archive (SVIM, 2015). For the present illustrative purpose we limit the analysis to the 2-year period from 1 January 2000 to 31 December 2001. Results are available as daily means on the model configuration's native 4 km stereographic grid projection (Lien et al., 2013).

The ocean module of the coupled model used for the regional simulation is the Regional Ocean Modeling System (ROMS), described in Haidvogel et al. (2008) and references therein. The sea ice module was developed by Budgell (2005). The ice model dynamics are based on the elastic–viscous–plastic (EVP) rheology following Hunke and Dukowicz (1997) and Hunke (2001), and the ice thermodynamics are based on Mellor and Kantha (1989) and Hakkinen and Mellor (1992).

The model results for sea ice concentration are somewhat noisy on the grid cell scale owing to the dispersiveness of the numerical scheme. In some regions, the grid cells that constitute the ice edge as defined by Eq. (1) can then appear as a mesh-like collection of cells. In order to reduce the impact of this issue, we applied the second order checkerboard suppression algorithm (Li et al., 2001) to the model results before conducting the present analysis. The sea ice concentrations from observations do not suffer from this type of noise; thus such an algorithm was not applied to the observational data set.

We compare model results with observations from the Arctic Ocean Sea Ice Concentration Charts Svalbard which is a multi-sensor data set that uses data from synthetic aperture radar (SAR) instruments as its primary source of information (WMO, 2017). This observational data set will be referred to as the ice chart data hereafter.

The ice chart data cover the northern Nordic Seas, the Barents Sea and adjacent ocean regions. The ice chart data deviate from a passive microwave product in this region, particularly in the final months of the melting season (e.g., Sect. 6 in Melsom et al., 2019). Hence the sea ice edge observations are estimated with an uncertainty that is not known.

The ice chart data are available on a stereographic grid projection with a resolution of 1 km. Data availability is restricted to working days. During a regular week, we then have 4 d with 24 h displacement results. The data set is also slightly reduced due to holidays, and a total of 354 d with 24 h ice edge displacement results were available from the present 2-year period.

The present study will be restricted to results and data for the Barents Sea. The SVIM simulation domain is displayed in Fig. 4, in which the Barents Sea analysis region is highlighted. Ice chart results are integrated into the SVIM domain using a mass-conserving Riemann integral approach. All grid cells inside the Barents Sea region which lack proper values (usually due to the presence of land) in at least one of the data sets are masked prior to the analysis. The analysis region is then composed of 80 399 wet grid cells, which represent an area of 1.29×10⁶ km².

Figure 5 Sample scene displaying the changes in model sea ice extent from 23 October 2001 (day 1) to 24 October 2001 (day 2). The black line indicates the maximum displacement distance ($d_{\text{max}}^{\mathrm{\Delta }t}$, given by Eq. 5) with the original algorithm, while the red line shows the result when grid cells along the open boundaries and coastlines are included ($\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{L}}\left(t_{\mathrm{0}}\right)$ from Eq. B7). The color coding is given by the label bar, and note that only the northern part of the Barents Sea analysis region is displayed.

3.2 Open boundaries and coasts

For some cases, the algorithm described in Sect. 2.1 does not describe properly the true displacement. This particularly affects the maximum value as defined by Eq. (5). We illustrate this issue with a case study displayed in Fig. 5, showing the 24 h change in the ice edge position from 23 to 24 October 2001. In this case the ice edge was displaced into the validation domain across an open boundary to the north. The general algorithm in Sect. 2.1 mismatches the ice edge grid cell (since ice beyond the open boundary is unseen) and leads to an unrealistic maximum ice edge displacement of 285 km, as given by the thick black line close to the subdomain's north border in the figure.

In order to address this issue, a modification of the algorithm was implemented in which ocean open boundaries are considered as continuations of the ice edge. The modified algorithm is described in full detail in Appendix B. For the case illustrated in Fig. 5, the maximum ice edge displacement calculated with the modified algorithm becomes 79 km. This displacement is indicated by the red line in the eastern part of the validation domain.

It must be noted that if the ice is advected into the domain, the distances associated with such a displacement will be underestimated since the real position of the ice edge outside of the analysis domain at t₀ is unknown. Another situation in which unrealistic representations for displacements may arise is when ice freezes along the coast, for example, due to colder air in the vicinity of continents, possibly in combination with less salty water masses close to the coastline. This issue is treated analogously to advection across an open boundary by considering coastlines as continuations of the ice edge (see Appendix B for details).

Figure 6 Histograms for the distribution of daily maximum displacement distances $\mathrm{\Delta }{d}_{\text{max}}^{\mathrm{\Delta }t}$, defined by Eq. (5). Horizontal bars pointing left and right correspond to results from SVIM model simulation and ice chart data, respectively. Light colored bars display results from the original algorithm in Sect. 2.1, while bars with regular colors (labeled “OB&C restricted”) result after the extension in Appendix B for open boundaries and coastlines is applied. Results from 354 d of 24 h edge displacements have been analyzed, see the text for further details.

Figure 7 Rank histogram for model results for the local ice edge displacement corresponding to the position of the maximum observed displacement. (a) Sets of nine alternative model displacements were derived for each of 235 d with 24 h displacements results. The nine displacement values were ordered from lowest to highest, and the local displacement was given a rank from the slot in which this value belonged; see the text for details. The blue curve shows the average local model displacement distances for results belonging to each of the ranks, with negative numbers corresponding to local sea ice retreat in the model results. The average maximum observed displacement is 39 km. (b) Results obtained for a subdivision as indicated by the red line in Fig. 4. Here, sets of seven alternative displacements results were derived for each subdomain. Then, results were available for 179 and 224 d for the western and eastern subdomains, respectively. The average maximum observed displacements are 36 km in both subdomains. The use of gray shading and line colors for the subdomains are indicated by the inset label. The y axes in (a) and (b) have been specified so that the black horizontal 0-axis lines for the right-side axes correspond to the frequency level of a flat distribution of frequencies (left-side y axes).

3.3 Validation results

We first examine the distribution of daily maximum ice edge displacements. From Fig. 4 we note that this examination will be performed for a domain in which advection of sea ice across the open boundaries is relevant, as well as freezing along coastlines of the continent and archipelagos. To address this issue, the analysis based on the algorithms in Sect. 2 will be extended following the outline in Sect. 3.2 and detailed in Appendix B.

In Fig. 6 results from the 354 d with 24 h maximum displacements from both data sets are displayed. We note that about two thirds of the maximum displacements in model results are in the range of 10–30 km. The corresponding distribution of results from the ice chart data has two maxima: one for the range of 20–40 km which accounts for nearly half of the cases and a secondary maximum for short (0–10 km) maximum displacements. The medians of the daily maximum displacement distances are 23 and 32 km for the SVIM results and the ice chart data, respectively. We also note that the distribution frequencies for the 2-year period change only moderately when the adjustments for open ocean boundaries and coastlines that were described in Sect. 3.2 are included in the analysis.

A conclusion that can be drawn from these results is that the largest expansions of sea ice extent in the model (SVIM) are underestimated when compared with observations (ice chart data). This is generally the case as the SVIM median is in the range of 20–30 km, while the median of the ice chart data is in the range of 30–40 km. The underestimation is also seen for extreme cases as the frequency of maximum expansion exceeding 60 km is about 5 times as high for the ice chart data.

In order to examine the degree to which SVIM results reproduce the geographical location of the observed maximum displacement, we apply the ranking method described in Sect. 2.2. We consider a fixed number of 10 bins for the present investigation. Hence, for each set of 24 h results for displacement distances, nine values are randomly selected from the displacements in Eq. (11). The requirement of at least nine additional ice edge positions separated by the decorrelation length scale restricts the cases that can be considered in this analysis: from the full set of 354 cases with 24 h displacement results, only 235 cases could then be kept in the analysis of ranked displacements as these cases met the requirement of at least 10 statistically independent displacement distances. The size of the set of independent values is restricted by the temporally varying degrees of freedom, as given by the ice edge length and the decorrelation length scale.

The resulting frequency distribution for each of the 10 ranks is displayed as gray vertical bars in Fig. 7a with rank values from 0 to 9. The highest rank (9) results when the model displacement close to the site with maximum displacement in the observations (the reference displacement, ${\mathit{\delta }}_{\mathrm{0}}^{\mathrm{M};\mathrm{\Delta }t}$) is larger than all displacements from the nine subsampled ice edge positions. The next rank (8) corresponds to cases when one and only one of the subsampled positions has a larger displacement than the reference displacement, and so on. In other words, high ranks indicate situations in which the position of the maximum displacement is described with a relatively high quality.

The histogram in Fig. 7a has nearly twice as many entries in ranks 5–9 than ranks 0–4. This mode for higher ranks indicate some skill for SVIM in detecting the location of the maximum displacements in the observations. The average rank in the present analysis is 5.48. For a random distribution of 235 integer numbers in the range of 0–9 percentiles 0.5 and 99.5 of the average rank are 4.02 and 4.98, respectively. Thus, the analysis reveals that while the model results are far from perfect, the average rank of 5.48 is significantly higher than results from random spatial distributions of ice edge displacements. We have also applied the more traditional χ² test for rank flatness (Wilks, 2019) and find that for the present histogram we have χ²=52. This value is nearly twice the magnitude of ${\mathit{\chi }}_{\text{crit}}^{\mathrm{2}}$ when this is set to reject the null hypothesis of a flat distribution at an α level of 0.001.

We supplement this analysis by dividing the region into two subdomains separated by the 40^∘ E meridian, as indicated by the red line in Fig. 4. In order to retain the majority of the days in the analysis, the number of randomly chosen displacements was reduced to seven values. This was due to the reduced degrees of freedom when the same decorrelation length scale was applied in smaller domains with shorter ice edges.

The results for the eastern and western Barents Sea subdomains are displayed in Fig. 7b. Contrasts between the frequency distributions between panels (a) and (b) arise for several reasons. First, the domain split leads to a set of two time series of maximum displacements in which the maxima from the full domain will be distributed between the two subdomain time series, and new maxima are introduced for the alternative subdomain. Next, the separation line between the two subdomains is manifested in the analysis as a new, shared open boundary. Hence, the shapes of the subdomain distributions may to some degree deviate from the full domain distribution. Note that the distribution peak at ranks 4–5 for the full domain disappears when considering the rank frequencies for the two subdomains.

The average ranks of the model displacements corresponding to the largest observed displacements are 3.95 and 4.13 in the western and eastern subdomains, respectively, whereas the ranges spanned between percentiles 0.5 and 99.5 for random distributions are [3.06,3.94] and [3.11,3.89] for the western (179 d) and eastern (224 d) subdomains, respectively.

4 Concluding remarks

In this study we present a new algorithm for the examination of the displacement of the edge (or the front) of a binary field. The algorithm computes different attributes of the ice edge displacement, such as the maximum and the distribution of the distances. Then, different methodologies for comparing these attributes are introduced. The method introduced enables the assessment of the forecast quality for the displacement of the ice edge and expands on existing validation metrics such as, for example, the integrated ice edge error (Goessling et al., 2016) and the various ice edge metrics considered by Melsom et al. (2019): the methods presented here provide summary statistics for the quality of model results for ice edge displacements in the presence of an expanding sea ice cover, as exemplified by Figs. 6 and 7, that are not provided with existing metrics. Such quality assessments are of high relevance for planned or ongoing site-specific activities in regions which can potentially become ice infested.

The present study has been framed in the context of results for displacements of the sea ice edge. Thus, the investigation in Sect. 3 was based on data for the sea ice edge from satellite observations and simulation results from a coupled ocean–sea ice model. However, the algorithm that was introduced in Sect. 2 can be applied to the displacement of the perimeter of any physical variable or feature that can be represented by a spatially continuous binary field. Stratiform precipitation is an example of another physical variable for which the method presented here could be suitable.

Note that we have used the term displacement rather than advection. The reason for this is that displacements need not be purely of an advective nature. In the case of sea ice, the displacement of the initial edge will generally be affected by freezing or melting along the perimeter of the sea ice extent. Analogously, displacement of the area affected by stratiform precipitation can be affected by new condensation or partial depletion of the cloud.

As demonstrated in the example depicted in Fig. 5, the original algorithm described in Sect. 2.1 and 2.2 can sometimes mismatch ice edge grid cells and hence diagnose unrealistic displacements, which may yield misleading results. Here, we have amended situations in which the sea ice enters a limited area domain across an open boundary and situations when freezing takes place next to a physical boundary (the coast). Modifications of the algorithm which include distances from ice edges to coasts and open boundaries are described in Sect. 3.2 and detailed in Appendix B. This approach eliminates unphysical edge displacement distance values, as shown by the case illustrated in Fig. 5. However, there may be other issues that can distort results that are produced by the analysis presented here. One example is cases when features are seen to arise seemingly spontaneously from one time of analysis to another: the algorithm in Sect. 2 can, for example, if applied to precipitation data, give rise to unrealistic results for displacements when convective precipitation cells develop.

Results from the algorithms that are introduced in the present study give valuable information regarding the changing extent of sea ice and how well the displacements of the observed and modeled sea ice edges agree. These algorithms have proven to provide simple yet robust and informative assessments for the quality of ice edge forecasts both with respect to the largest displacements from one time to another and with respect to the reproduction of the geographical position where the largest displacement occurs.

Appendix A: Decorrelation length of displacements

Assume that we have a set of N edge grid cells e_n (i.e., satisfying Eq. 1) that form a curve,

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }L=\mathit{\{}{e}_{\mathrm{1}},e_{\mathrm{2}},\mathrm{\dots },e_N\mathit{\}},\end{array}$$

where L is continuous in the sense that grid cells e_n and e_n+1 are neighbors. Furthermore, associate displacement distances d_n to each e_n as defined in Sect. 2.1. Then, the spatial autocorrelation of displacements can be estimated using a sample Pearson correlation coefficient approach:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }r\left(\mathit{\eta }\right)=\sum _{n=\mathrm{1}}^{N-\mathit{\eta }}\left[\left(d_n-\stackrel{\mathrm{‾}}{d_n}\right)\left(d_{n+\mathit{\eta }}-\stackrel{\mathrm{‾}}{d_{n+\mathit{\eta }}}\right)\right]/\\ \text{(A2)}& {\displaystyle }& {\displaystyle }{\left[\sum _{n=\mathrm{1}}^{N-\mathit{\eta }}{\left(d_n-\stackrel{\mathrm{‾}}{d_n}\right)}^{\mathrm{2}}\sum _{n=\mathrm{1}}^{N-\mathit{\eta }}{\left(d_{n+\mathit{\eta }}-\stackrel{\mathrm{‾}}{d_{n+\mathit{\eta }}}\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

We have r(0)=1 and we define the decorrelation length of the displacements, Δn, as

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }n=\underset{\forall \mathit{\eta }}{min}\left(\mathit{\eta }\right|r\left(\mathit{\eta }\right)<\mathrm{1}/\text{e}),\end{array}$$

where e is Euler's number. If, for a given time, the ice edge is discontinuous, each continuous curve segment is treated separately, and the weighted mean value of the results for Δn from each segment is used. In that case, weights are applied according to the number of edge grid cells in each curve segment.

Appendix B: Open boundaries and coasts

As discussed by Melsom et al. (2019), open boundaries and coastlines can potentially have significant impacts on the results for the metrics for the position of the sea ice edge. Here, we introduce a method which will give more physically meaningful results if ice either freezes near a coastline or enters into a domain across an open boundary. Moreover, the method affects the results modestly or not at all for an edge that is displaced inside of the domain.

First, set the open boundary grid lines to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{L}^{\text{OB}}\left(t_{\mathrm{0}}\right)=\mathit{\{}{e}_{{\mathrm{1}}_{\text{OB}}},e_{{\mathrm{2}}_{\text{OB}}},\mathrm{\dots },e_{N_{\text{OB}}}\mathit{\}},\\ \text{(B1)}& {\displaystyle }& {\displaystyle }c\left[e_{n_{\text{OB}}}\right]\left(t_{\mathrm{0}}\right)<c_{\text{edge}},\end{array}$$

where $e_{n_{\text{OB}}}$ is any ocean grid cell along the boundary of the domain which was on the open ocean side of the ice edge at t=t₀. Then L(t₀) in Eq. (3) can be replaced by

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{L}\left(t_{\mathrm{0}}\right)=L\left(t_{\mathrm{0}}\right)\cup L^{\text{OB}}\left(t_{\mathrm{0}}\right),\end{array}$$

and for the corresponding distances we introduce the notation $\stackrel{\mathrm{̃}}{d}$, so Eq. (3) becomes

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{̃}}{d_n}}^{\mathrm{\Delta }t}=min\left|\left|e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}-\stackrel{\mathrm{̃}}{L}\left(t_{\mathrm{0}}\right)\right|\right|,\end{array}$$

where $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ is a grid cell on L(t₀+Δt), as before. The set of grid cells $e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}$ is not affected, so the number of displacement distances considered in Eq. (5), N(t₀+Δt), is unchanged. Note that here, the additional curve L^OB is only added to L(t₀). Otherwise, if it was also added to L(t₀+Δt), the trivial score from perfect matching segments of edge curves would significantly impact the results as they are, for example, displayed in Figs. 2 and 3.

A sample grid cell to which the displacement distance is significantly affected by this modification is displayed as e_n3 in Fig. B1. It must be noted that if the ice is imported into the domain, the distances $\stackrel{\mathrm{̃}}{d}$ associated with such a displacement will be underestimated since the real position of the ice edge outside of the analysis domain at t₀ is unknown. Moreover, for regular displacement of ice inside the domain, results will be affected slightly when occurring in the vicinity of the open boundary (e.g., e_n2 in Fig. B1).

Similarly, there can be cases where freezing of ice occurs along the coastline, for example, as an effect of colder air in the vicinity of continents or less salty water masses close to the coastline. This is another case where the algorithm above can yield grossly exaggerated displacement distances. Again, the problem can be alleviated by including additional grid lines.

Figure B1 Binary fields with values of 1 (ice) and 0 (no ice/ocean) are displayed by white and blue color shading, respectively. Land is indicated as a black region. Light shades of blue indicate regions with a non-overlapping ice cover, as indicated by the color legend. Open boundary grid cells are depicted as gray lines. Ice edges for t₀ and t₀+Δt are drawn as light red and red lines, respectively. Dashed black lines show the edge displacements as defined in Sect. 2.1 for a selection of labeled ice edge grid cells from t₀+Δt, marked by white circles. Full black lines display the displacements $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{d}}$ that result from the modifications described in Appendix B; see Eq. (B8). For the set of grid cells that are highlighted here, only e_n1 is unaffected by the modified definitions.

Set the coastal grid lines as

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\displaystyle }{L}^{\mathrm{C}}\left(t_{\mathrm{0}}\right)=\mathit{\{}{e}_{{\mathrm{1}}_{\mathrm{C}}},e_{{\mathrm{2}}_{\mathrm{C}}},\mathrm{\dots },e_{N_{\mathrm{C}}}\mathit{\}},c\left[e_{n_{\mathrm{C}}}\right]\left(t_{\mathrm{0}}\right)<c_{\text{edge}},\end{array}$$

where $e_{n_{\mathrm{C}}}$ is any ocean grid cell along the coastline which was ice free or had a sea ice concentration below c_edge at t=t₀. Then L(t₀) can be replaced by

$$\begin{array}{}\text{(B5)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{‾}}{L}\left(t_{\mathrm{0}}\right)=L\left(t_{\mathrm{0}}\right)\cup L^{\mathrm{C}}\left(t_{\mathrm{0}}\right).\end{array}$$

Here, Eq. (3) will be replaced by

$$\begin{array}{}\text{(B6)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{d_n}}^{\mathrm{\Delta }t}=min\left|\left|e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}-\stackrel{\mathrm{‾}}{L}\left(t_{\mathrm{0}}\right)\right|\right|,\end{array}$$

and the computation of the displacement distance to grid cell e_n4 in Fig. B1 is severely affected.

For a regional model, the typical situation is that there are both open boundaries and coastlines. In that case, we may combine the above modifications of the algorithm by adopting

$$\begin{array}{}\text{(B7)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{L}}\left(t_{\mathrm{0}}\right)=L\left(t_{\mathrm{0}}\right)\cup L^{\text{OB}}\left(t_{\mathrm{0}}\right)\cup L^{\mathrm{C}}\left(t_{\mathrm{0}}\right),\end{array}$$

and Eq. (3) can be written

$$\begin{array}{}\text{(B8)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{d_n}}}^{\mathrm{\Delta }t}=min\left|\left|e_n^{(t_{\mathrm{0}}+\mathrm{\Delta }t)}-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{L}}\left(t_{\mathrm{0}}\right)\right|\right|.\end{array}$$