the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# 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.

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.

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*_{0} and *t*_{0}+Δ*t*, are displayed. In the context of forecasting, *L*^{M}(*t*_{0}) may be taken to represent the model initialization at *t*_{0}, and *L*^{M}(*t*_{0}+Δ*t*) is then the forecast at a temporal range of Δ*t*. The other binary fields can represent observations at the same times.

## 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

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

This follows the algorithm presented in Melsom et al. (2019). Let *L*(*t*_{0}) and *L*(*t*_{0}+Δ*t*) denote the sea ice edges at times *t*_{0} and *t*_{0}+Δ*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*_{0}+Δ*t*) and curve *L*(*t*_{0}), i.e.,

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*_{0}). 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*_{0}), respectively, i.e.,

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*_{0}.

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

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*_{0}) to *L*(*t*_{0}+Δ*t*).

If we briefly introduce ${d}_{m}^{-\mathrm{\Delta}t}$ as the shortest distance from a grid cell *e*_{m} of *L*(*t*_{0}) to the curve *L*(*t*_{0}+Δ*t*), we note that the Hausdorff distance *d*_{H} between curves *L*(*t*_{0}+Δ*t*) and *L*(*t*_{0}) is

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*_{0}+Δ*t* to the ice edge curve from the initialization at *t*_{0}.

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*_{0}+Δ*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.

## 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.,

where ${d}_{\text{max}}^{\mathrm{O};\mathrm{\Delta}t}$ is computed from observed ice edges at *t*_{0} and *t*_{0}+Δ*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*_{0}+Δ*t*) to which the maximum edge displacement is found in the observations. Then, determine ${\mathit{\u03f5}}_{\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{\u03f5}}_{\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

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

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

respectively, limited by the first and last cells along the curve *L*^{M}(*t*_{0}+Δ*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.

## 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^{6} km^{2}.

## 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*_{0} 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).

## 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 *χ*^{2} test for rank flatness (Wilks, 2019) and find that for the present histogram we have *χ*^{2}=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.

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.

Assume that we have a set of *N* edge grid cells *e*_{n} (i.e., satisfying Eq. 1) that form a curve,

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:

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

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.

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

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*_{0}. Then *L*(*t*_{0}) in Eq. (3) can be replaced by

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

where ${e}_{n}^{({t}_{\mathrm{0}}+\mathrm{\Delta}t)}$ is a grid cell on *L*(*t*_{0}+Δ*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*_{0}+Δ*t*), is unchanged. Note that here, the additional curve *L*^{OB} is only added to *L*(*t*_{0}). Otherwise, if it was also added to *L*(*t*_{0}+Δ*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{\u0303}}{d}$ associated with such a displacement will be underestimated since the real position of the ice edge outside of the analysis domain at *t*_{0} 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.

Set the coastal grid lines as

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*_{0}. Then *L*(*t*_{0}) can be replaced by

Here, Eq. (3) will be replaced by

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

and Eq. (3) can be written

The idealized distributions of concentrations that are depicted in Fig. 1 and source code for computing results displayed in Figs. 2 and 3 are available from https://doi.org/10.5281/zenodo.4545686 (Melsom, 2021). Model results that were analyzed in Sect. 3 are available from https://archive.norstore.no/pages/public/datasetDetail.jsf?id=10.11582/2015.00014 (SVIM, 2015), while the observations are available from https://thredds.met.no/thredds/catalog/arcticdata/met.no/iceChartSat/catalog.html (last access: 28 October 2020, WMO, 2017).

The author declares that there is no conflict of interest.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This research has been supported by the Norges Forskningsråd (grant no. 276730). Additional funding was provided by the Copernicus Programme via Mercator Océan (grant no. 2015/S 009-011301).

This paper was edited by Lars Kaleschke and reviewed by two anonymous referees.

Aksenov, Y., Popova, E. E., Yool, A., Nurser, A. G., Williams, T. D., and Bertino, L.: On the future navigability of Arctic sea routes: High-resolution projections of the Arctic Ocean and sea ice, Mar. Policy, 75, 300–317, https://doi.org/10.1016/j.marpol.2015.12.027, 2017. a

Anderson, J. L., A method for producing and evaluating probabilistic forecasts from ensemble model integrations, J. Climate, 9, 1518–1530, https://doi.org/10.1175/1520-0442(1996)009<1518:AMFPAE>2.0.CO;2, 1996. a

Bitz, C. M., Holland, M. M., and Hunke, E. C.: Maintenance of the Sea-Ice Edge, J. Climate, 18, 2903–2921, https://doi.org/10.1175/JCLI3428.1, 2005. a

Budgell, P. W.: Numerical simulation of ice-ocean variability in the Barents Sea region: Towards dynamical downscaling, Ocean Dynam., 55, 370–387, https://doi.org/10.1007/s10236-005-0008-3, 2005. a

Cheng, A., Casati, B., Tivy, A., Zagon, T., Lemieux, J.-F., and Tremblay, L. B.: Accuracy and inter-analyst agreement of visually estimated sea ice concentrations in Canadian Ice Service ice charts using single-polarization RADARSAT-2, The Cryosphere, 14, 1289–1310, https://doi.org/10.5194/tc-14-1289-2020, 2020. a

Davis, C. A., Brown, B. G., and Bullock, R. G.: Object-based verification of precipitation forecasts, Part I: Methodology and application to mesoscale rain areas, Mon. Weather Rev., 134, 1772–1784, https://doi.org/10.1175/MWR3145.1, 2006. a

Dukhovskoy, D. S., Ubnoske, J., Blanchard-Wrigglesworth, E., Hiester, H. R., and Proshutinsky, A.: Skill metrics for evaluation and comparison of sea ice models, J. Geophys. Res.-Oceans, 120, 5910–5931, https://doi.org/10.1002/2015JC010989, 2015. a, b

Ebert, E. E. and McBride, J. L.: Verification of precipitation in weather systems: determination of systematic errors, J. Hydrol., 239, 179–202, https://doi.org/10.1016/S0022-1694(00)00343-7, 2000. a, b

Goessling, H. F. and Jung, T.: A probabilistic verification score for contours: Methodology and application to Arctic ice-edge forecast, Q. J. Roy. Meteor. Soc., 144, 735–743, https://doi.org/10.1002/qj.3242, 2018. a

Goessling, H. F., Tietsche, S., Day, J. J., Hawkins, E., and Jung, T.: Predictability of the Arctic sea ice edge, Geophys. Res. Lett., 43, 1642–1650, https://doi.org/10.1002/2015GL067232, 2016. a, b

Haidvogel, D. B., Arango, H. G., Budgell, P. W., and 17 coauthors: Ocean forecasting in terrain-following coordinates: Formulation and skill assessment of the Regional Ocean Modeling System, J. Comput. Phys., 227, 3595–3624, https://doi.org/10.1016/j.jcp.2007.06.016, 2008. a

Häkkinen, S., and Mellor, G. L.: Modelling the seasonal variability of a coupled arctic ice-ocean system, J. Geophys. Res., 97, 20.285–20.304, https://doi.org/10.1029/92JC02037, 1992. a

Hamill, T. M.: Interpretation of Rank Histograms for Verifying Ensemble Forecasts, Mon. Weather Rev., 129, 550–560, https://doi.org/10.1175/1520-0493(2001)129<0550:IORHFV>2.0.CO;2, 2001. a

Hunke, E. C.: Viscous-plastic sea ice dynamics with the EVP model: linearization issues, J. Comput. Phys., 170, 18–38, https://doi.org/10.1006/jcph.2001.6710, 2001. a

Hunke, E. C. and Dukowicz, J. K.: An elastic-viscous-plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, https://doi.org/10.1175/1520-0485(1997)027<1849:AEVPMF>2.0.CO;2, 1997. a

Lavergne, T., Eastwood, S., Teffah, Z., and Schyberg, H.: Sea ice motion from low resolution satellite sensors: An alternative method and its validation in the Arctic, J. Geophys. Res., 115, C10032, https://doi.org/10.1029/2009JC005958, 2010. a

Leese, J. A., Novak, C. S., and Clark, B. B.: An Automated Technique for Obtaining Cloud Motion from Geosynchronous Satellite Data Using Cross Correlation, J. Appl. Meteorol., 10, 118–132, https://doi.org/10.1175/1520-0450(1971)010<0118:AATFOC>2.0.CO;2, 1971. a

Li, Q., Steven, G. P., and Xie, Y. M.: A simple checkerboard suppression algorithm for evolutionary structural optimization, Struct. Optimization, 22, 230–239, https://doi.org/10.1007/s001580100140, 2001. a

Lien, V. S., Gusdal, Y., Albretsen, J., Melsom, A., and Vikebø, F. B.: Evaluation of a Nordic Seas 4 km numerical ocean model hindcast archive (SVIM), Fisken og Havet, 2013–2017, Bergen, Norway, 80 pp., 2011. a

Mellor, G. L. and Kantha, L.: An ice ocean coupled model, J. Geophys. Res., 94, 10.937–10.954, https://doi.org/10.1029/JC094iC08p10937, 1989. a

Melsom, A.: Edge metrics software v.1.0 [Software], Zenodo [code], https://doi.org/10.5281/zenodo.4545686, 2021. a

Melsom, A., Palerme, C., and Müller, M.: Validation metrics for ice edge position forecasts, Ocean Sci., 15, 615–630, https://doi.org/10.5194/os-15-615-2019, 2019. a, b, c, d, e

The NCAR Command Language (Version 6.4.0) [Software], UCAR/NCAR/CISL/TDD [code], Boulder, Colorado, https://doi.org/10.5065/D6WD3XH5, 2017. a

Norwegian Meteorological Institute, Institute of Marine Research: SVIM ocean hindcast archive, National Infrastructure for Research Data [data set], https://doi.org/10.11582/2015.00014, 2015. a, b

Palerme, C., Müller, M., and Melsom, A.: An intercomparison of skill scores for evaluating the sea ice edge position in seasonal forecasts, Geophys. Res. Lett., 46, 4757–4763, https://doi.org/10.1029/2019GL082482, 2019. a

Parkinson, C. L.: Spatially mapped reductions in the length of the Arctic sea ice season, Geophys. Res. Lett., 41, 4316–4322, https://doi.org/10.1002/2014GL060434, 2014. a

Talagrand, O., Vautard R., and Strauss B.: Evaluation of probabilistic prediction systems, Proc. Workshop on Predictability, ECMWF, Reading, UK, 20–22 October 1997, 25 pp., 1997. a

Tokmakian, R., Strub, P. T., and McClean-Padman, J.: Evaluation of the Maximum Cross-Correlation Method of Estimating Sea Surface Velocities from Sequential Satellite Images, J. Atmos. Ocean. Tech., 7, 852–865, https://doi.org/10.1175/1520-0426(1990)007<0852:EOTMCC>2.0.CO;2, 1990. a

Wilks, D.: Indices of Rank Histogram Flatness and Their Sampling Properties, Mon. Weather Rev., 147, 763–769, https://doi.org/10.1175/MWR-D-18-0369.1, 2019. a

WMO: Sea-Ice Information Services in the World, World Meteo. No. 574, World Meteorological Organization, Geneva, Switzerland, 103 pp., 2017. a, b