## Abstract

Sea ice models have become essential components of weather, climate, and ocean models. A realistic representation of sea ice affects the reliability of process representation, environmental forecast, and climate projections. Realistic simulations of sea ice kinematics require the consideration of both large-scale and finescale geomorphological structures such as linear kinematic features (LKF). We propose a multiscale directional analysis (MDA) that diagnoses the spatial characteristics of LKFs. The MDA is different from previous analyses in that it (i) does not detect LKFs as objects, (ii) takes into account the width of LKFs, and (iii) estimates scale-dependent orientation and intersection angles. The MDA is applied to pairs of deformation fields derived from satellite remote sensing data and from a numerical model simulation with a horizontal grid spacing of ~4.5 km. The orientation and intersection angles of LKFs agree with the observations and confirm the visual impression that the intersection angles tend to be smaller in the satellite data compared to the model data. The MDA distributions can be used to compare satellite data and numerical model fields using conventional metrics such as a Euclidean distance, the Bhattacharyya coefficient, or the Earth mover’s distance. The latter is found to be the most meaningful metric to compare distributions of LKF orientations and intersection angles. The MDA proposed here provides a tool to diagnose if modified sea ice rheologies lead to more realistic simulations of LKFs.

## 1. Introduction

The fast retreat of the Arctic summer sea ice and the continually thinning multiyear ice facilitate an increase of marine transport and offshore operations, the exploitation of natural resources, and more marine industries in general in the Arctic (e.g., Melia et al. 2017; Smith and Stephenson 2013). The growing number of activities require accurate predictions not only of weather but also of the sea ice state for marine hazard identification, risk assessment, and minimizing possible but irreparable environmental impacts.

Recent work on sea ice forecast verification has focused mostly on the location of the sea ice edge (e.g., Dukhovskoy et al. 2015; Goessling et al. 2016; Goessling and Jung 2018; Melsom et al. 2019; Palerme et al. 2019). Sea ice deformation is another critically important variable as it is directly related to sea ice strength, ice pressure, and openings of leads [see Weiss (2013) for a definition of sea ice deformation]. Thus, to meet requirements of marine forecasts users, a verification method needs to evaluate quantitatively the spatial characteristics of complex structures of geomorphological features (e.g., leads) in high-resolution sea ice deformation fields. These features are important both for climate and weather research because they affect crucial physical properties such as surface turbulent heat transfer (e.g., Marcq and Weiss 2012; Alam and Curry 1995) or momentum and kinetic energy balances (Parmerter and Coon 1972).

Furthermore, properly representing the orientation and patterns of linear features (leads and pressure ridges) in models is thought to be directly related to the sea ice rheology of the model (e.g., Lindsay and Rothrock 1995; Hutchings et al. 2005). It is therefore desirable to evaluate linear kinematic features (LKF; Kwok 2001) emerging in high-resolution sea ice models against observations (Coon et al. 2007; Hutter and Losch 2020). In the development process of new rheology models, the generation of leads and their orientation are important parameters against which the success of a parameterization is gauged (Schreyer et al. 2006; Schulson 2004). In addition, the width and orientation of features are important for determining momentum transfer across sea ice (e.g., Vihma 1995).

In this paper, we present a new evaluation method that can be used to assesses the skills of sea ice models in generating linear features. We introduce a multiscale directional analysis (MDA) as a spatial verification method. It estimates and compares the main (i.e., predominant) local orientation, intersection angles of the sea ice structures on one or more spatial scales. Envisaged applications of the MDA are the assessment of simulated sea ice deformation fields in general and the estimation of the actual predictive skills of sea ice forecasts on short-term time scales in particular.

The methods introduced here build in part on recent work that investigated potential predictability of sea ice linear features in twin experiments (Mohammadi-Aragh et al. 2018). An important part of our method is the detection of geomorphological features. Finding a unique definition for a linear kinematic feature that is applicable to both modeled and observed fields is a crucial part of the MDA. For various reasons, deformation fields from numerical models are generally smoother than observed fields, so that the existing algorithms for detecting LKF need to be tuned separately for model simulations and for observations. As a result, the verification depends largely on specific parameter choices especially when a frequency-based metric is used. We define an LKF as the locus of points that have specific geomorphological properties (section 3a). Following Levy et al. (2008), an LKF is a one-dimensional curve in a two-dimensional sea ice deformation field. We use an algorithm that is able to separate the linear features from other types of detected geometrical features.

Sea ice models at coarse resolution tend to produce smooth concentration and thickness fields with very little detail, but at higher spatial resolution it is possible to model sea ice failure zones (Hutchings et al. 2005) using a viscous–plastic sea ice model and an isotropic rheology. With increasing resolution, the number of resolved discontinuities increases (e.g., Spreen et al. 2017; Mohammadi-Aragh et al. 2018; Hutter et al. 2018), but linear features in the simulated ice concentration cannot realistically represent leads because the width of sea ice fractures (leads and refrozen leads) is below 3 km (e.g., in Baffin Bay and the North Water area; Barry et al. 1989; Steffen 1987) and mostly even below 1 km (see also Feltham 2008). Even with the highest (published) grid spacing of less than 1 km (Hutter et al. 2018), the models are only beginning to represent these scales.

Instead of focusing on sea ice concentration and thickness, we follow Levy et al. (2008) and evaluate the geometrical features observed in sea ice deformation fields. The associated shear and divergence processes form ridges and open water areas in the sea ice cover with elongated, curvilinear, or piecewise linear features. More importantly, the spatial characteristics of these LKFs are highly correlated with the spatial pattern of leads. For example, the directional orientations of newly formed leads in the winter ice pack of the Beaufort Sea in data from the *European Earth Resources Satellite-1* (*ERS-1*) Synthetic Aperture Radar (SAR) were shown to be consistently related to the principal direction of sea ice shear (Cunningham et al. 1994). In addition, sea ice shear stress amplifies the tendency of divergence in the production of leads (Miles and Barry 1998; Stern et al. 1995). Sea ice shear and divergence are generally associated with leads (Barry et al. 1989). Thus, it is practical to detect the LKFs observed in the deformation field as indicative of leads and ridges.

Previously published verification methods of linear features in the sea ice cover are limited either to frequency-based indices of agreement (e.g., Levy et al. 2008) or to comparing time series of a global mean variable such as lead fraction, sea ice divergence, and shear (Wang et al. 2016). By contrast, we use and further develop a new multiscale directional orientation estimation method as the fundamental part of our MDA following a development in the context of neural imaging (Singh et al. 2016). Our orientation approximation method differs from previous approaches that measure the orientations of linear (kinematic) features in sea ice in at least four aspects: (i) Following Singh et al. (2016), we take into account the width of the individual LKF instead of only centerlines (as, e.g., Cunningham et al. 1994; Banfield 1992; Linow and Dierking 2017; Hutter et al. 2019). This reduces at least the uncertainties in measuring the orientation of wide and small-scale linear features. (ii) By applying the MDA we consider different spatial scales. This helps to evaluate how simulated deformation fields are affected, for example, by numerical schemes and parameterizations on different spatial scales. (iii) LKFs, for example, observed in satellite images are categorized as one object class (Linow and Dierking 2017; Hutter et al. 2019). The MDA works directly on gridded data without the need to isolate LKFs as objects, similar to neighborhood verification methods (Ebert 2008; Casati et al. 2008). (iv) The multiscale directional orientation estimation, which is one component of the MDA, also slightly differs from the original method of Singh et al. (2016) in three aspects: First, we find the main local orientation with higher accuracy by using a parabolic fit; second, we do not split the LKFs into individual segments; and third, our analysis is not limited to orientation but is extended to intersection angles.

The paper is structured as follows. In section 3 we review the components of the methodology and explain the algorithms underlying the MDA: feature detection, estimation of main orientation, feature classification, estimation of intersection angles, and distance and similarity measures. In section 4, we test and illustrate the MDA with a synthetic field of linear features (hence feature detection is not necessary). Section 4 is devoted to assessing the performance of the MDA in verifying pairs of actual modeled and observed deformation fields. In the last section we present the summary and conclusions.

## 2. Data–sea ice model and observation data source

We designed an Arctic-wide sea ice–ocean model experiment based on the configuration of Nguyen et al. (2011) of the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. 1997; Adcroft et al. 2019). Our configuration includes sea ice dynamics with a viscous–plastic rheology and 0-layer thermodynamics (Losch et al. 2010). The simulation spans the years 2001–15 and is forced with the ERA-Interim reanalysis data (Dee et al. 2011). The horizontal grid spacing of the experiment is ~4.5 km; the ocean model has 50 vertical ocean layers. We extract monthly mean open boundary conditions in the Pacific and the Atlantic from a global simulation (Menemenlis et al. 2008). For more details see Mohammadi-Aragh et al. (2018) or Spreen et al. (2017).

For model evaluation, we use the gridded dataset of the RADARSAT Geophysical Processor System (RGPS) sea ice divergence and shear products (available for the years between 1996 and 2008 at https://asf.alaska.edu/data-sets/derived-data-sets/seaice-measures/sea-ice-measures-data-products/). The gridded product is an interpolated version of the Lagrangian dataset, which originated from feature tracking (e.g., Bouillon and Rampal 2015). The data values cover 3-day periods with a grid spacing of 12.5 km. To allow the direct comparison with simulated data, we compute corresponding time averages of the modeled sea ice deformation fields and interpolate them spatially to the RGPS grid by nearest neighbor interpolation.

## 3. Methodology

In this section we introduce the components of the MDA. In each subsection, we first give a brief overview of the recent history of the algorithms and alternative approaches and then explain and discuss the logic behind our selection from (and modification of) the available methods for designing the components of the MDA.

### a. Detecting linear kinematic features: Edge detecting and image separation (EDIS)

Our main goal is to design an approach that detects comparable geomorphological features in both modeled and observed deformation fields. One approach for detecting geomorphological features in sea ice is to use a constant threshold of the considered quantity (in our case deformation) independent of the characteristics of the background field (e.g., Dierking and Dall 2007; Lindsay and Rothrock 1995; Levy et al. 2008; Cunningham et al. 1994; Banfield 1992; Barry et al. 1989). This means that variations in the magnitude of the background deformation field caused primarily by the structure of the large-scale wind (and ocean) forcing are not considered. When the threshold is too high, we detect only large linear features and little noise, but fail to detect many linear features whose signal is just above the noise level. In contrast, when the deformation threshold is too low, spurious features that are not LKF are identified as signal. Mohammadi-Aragh et al. (2018) suggested an approach that defines a threshold relative to the average total deformation of the local background. This method still requires different tuning parameters for the threshold and neighborhood radius to derive meaningful patterns of features from the observed and simulated deformation fields. Alternatively, Linow and Dierking (2017) applied an edge detection and feature enhancement method to an observed deformation field. Other sophisticated methods to detect geometrical features based on shearlets (e.g., Yi et al. 2009; Reisenhofer et al. 2016; Easley et al. 2008; Kayasandik et al. 2018) have to our knowledge not yet been applied to sea ice fields. Shearlets are a multiscale framework that allows the construction of efficient algorithms for denoising and inverse problems, such as image enhancement and edge detection.

Because we aim for an algorithm that does not require separate tuning for the different types of data, we apply an edge detection and feature enhancement method similar to Linow and Dierking (2017). Our approach includes the detection of zero-crossing regions in the field of second spatial derivatives of total deformation. In a second step, LKFs are recognized as the locus of the points in proximity of the edges with higher sea ice deformation than the surrounding detected edge points.

As will be shown in the results section 4, the detected features include some point-like features (which may consist of several pixels that form a circular cluster) that cannot be regarded as LKFs. We identify and remove point-like features with an image separation method (Kutyniok and Lim 2010) based on wavelets and compactly supported shearlets.

This algorithm is part of the ShearLab toolbox available online (http://www.shearlab.org). We developed a python version of the algorithm for our analysis.

In summary, the EDIS algorithm used here consists of the following three steps:

Edge detection (Fig. 1c): We apply a two-dimensional Marr–Hildreth operator (MHO) on sea ice deformation fields (Figs. 1a,b). The MHO is an edge-detecting operator and generates closed and connected contours (Mohammadi-Aragh et al. 2006). It does so by searching zero crossings in the field of second derivatives and suppresses noise before enhancement with a Gaussian filter [see Spontón and Cardelino (2015) for details and a review on classical edge detectors].

Filling algorithm (Fig. 1d): LKFs are recognized as the locus of the points in short distance of the edges that have higher sea ice deformation than edge points.

Image separation (Figs. 1e,f): Curve-like features are separated from point-like features. This step also helps to remove the known gridscale noise associated with the deformation rate derived from RGPS (Bouillon and Rampal 2015).

### b. Estimation of main orientation

The multiscale directional orientation estimation is the central component of the MDA and overcomes some deficiencies of previous methods. So far, measuring the orientation of linear features in the sea ice cover has been mostly restricted to estimating the orientation of centerlines (e.g., Banfield 1992; Linow and Dierking 2017; Hutter et al. 2019). A major advantage of considering centerlines instead of the whole extent of a linear structure is data reduction. It typically yields only single orientation estimates for entire features. To enhance efficiency, some of the existing methods eliminate spatial noise of centerlines using a smoothing function with tuning parameters for the degree and spatial scale of smoothing. Moreover, robustly estimating the orientation typically requires linear features to have low curvature (e.g., Farmer and Li 1995). This assumption is a source of error, as structures with high curvature are frequently found in observations (see, e.g., the LKFs in Fig. 3 of Wang et al. 2003), and LKFs usually reveal varying orientations along their extent.

More useful for our purpose (of evaluating model simulations against observations) is to estimate orientations locally along LKFs. One approach for this is to compute tangent vectors to the centerlines (e.g., Cunningham et al. 1994). Due to both intrinsic characteristics of the discretized domain and unrealistic irregularities in centerlines, however, linear features exhibit local nonphysical irregularities. Using tangent vectors thus adds another source of error and spatial-scale dependence (e.g., Yi et al. 2009). Alternative approaches toward local orientation estimation are based on discrete shearlets (e.g., section III.A in Yi et al. 2009) or curvelets (e.g., Ma and Plonka 2010) but they do not take into account the width of the LKFs either.

Here we overcome the issues by using multiscale directional local orientation estimation. The major advantage of this approach is that it provides local estimates of orientation so that segmentation and connection algorithms are not required. Our approach builds on a conceptual framework with a set of multiscale directional filters (Singh et al. 2016). The filters are constructed by dilating and rotating a rectangle. The convolution of the spatial filter at a given point along the centerline of a linear feature produces a histogram (Fig. 2b) that describes the directional response as a function of orientation (see section 4a for more details). In theory, the maximum directional response occurs when the filter is aligned with the linear feature. The histogram shape also depends on the dimensions of selected spatial filters (and hence the implied spatial scales). For example, a feature with wiggles has various local orientations while at the same time the large-scale orientation of the feature as a whole may be uniform. Accordingly, small filters capture the small-scale variations of orientations along LKFs, whereas large filters detect the uniform main orientation. The MDA can take these different scales into account by applying rectangular filters with different dimensions (i.e., length scales).

Single-orientation methods require isolating individual linear features to avoid interference from neighboring structures (Singh et al. 2016). Our method differs from the multiple orientation estimation of Singh et al. (2016) and from the methods that derive just one orientation per feature in that we do not mask adjacent features. In our analysis we identify and use the main orientation but consider the orientation signals of the surrounding features in an MDA component related to the determination of intersection angles, which is described below. In summary, our orientation estimation algorithm consists of the following three steps:

Thinning (centerline tracing): This algorithm generates a graph of points (i.e., centerlines) through the midlines of LKFs. We use the skeletonize subroutine (Zhang and Suen 1984) of scikit-image (available online at https://scikit-image.org/). This algorithm changes LKFs to an attenuated form with unitary width. Note that we use centerlines only to locate the center coordinates of the spatial filters but do not estimate the orientation of centerlines per se.

Estimating directional responses: Using dilated, directional filters (green rectangle in Fig. 2a) we estimate the local directional response for orientations between 0° and 180°, averaging over equally sized bins. Using small bins and the corresponding maximum response to determine the main orientation tends to produce noisy results. We therefore use relatively coarse 10° steps.

Determining the main orientation: Given that the filter response is computed only for specific rotation angles with relatively coarse steps, using the discrete angle with the maximum response (the modal value) directly would be inaccurate. To increase accuracy, we fit a parabolic function of the rotation angle to the directional filter response around the maximum response. The rotation angle corresponding to the maximum of the parabolic function is taken as the main local orientation.

### c. Intersection and line points classification

To enable the estimation of intersection angles between crossing LKFs we suggest a classification algorithm that separates the points on the centerline into intersection and line points. The conceptual framework of the classification is simple. Each point on the graph is associated with a histogram of directional response for a spatial filter (see Fig. 2b and section 3b). All histograms have at least one local maximum, and the histograms of points in the vicinity of corners or junctions have two or more local maxima. As an example, the location and directional response of two points P1 and P2 on an LKF centerline are shown in Figs. 2a and 2b. P1 is located on a line segment far from other LKFs and P2 at a junction. The histogram associated with P1 exhibits one maximum, the histogram for P2 two and three local maxima, dependent on the dimensions of the directional filter (see also section 4a). Our classification algorithm consists of the following five steps:

We assign all points on the centerline that have more than one local maximum in the directional response to the class of intersection points. The remaining points belong to the class of line points.

All points on the centerline form a graph. Using the graph, we compute the minimum distance

*d*along the centerline between each pair of points belonging to the class of intersection points (Dijkstra 1959).For each point in the class of intersection points, we compute the average of its minimum distance with all intersection points $d\xaf$.

We divide the class of intersection points into smaller groups of point that have approximately equal average minimum distance $d\xaf$ (red points in Figs. 2e and 2f). Using trial and error, we suggest that for the LKF dataset, the difference between the average minimum distances of points in such a group should not be more than five grid points.

In each of these groups of intersection points, we determine the centroid as the point closest to its actual centroid (pink and dark cyan points in Fig. 2e). This centroid point is considered to be the final intersection (and junction) point.

### d. Estimation of intersection angles

The simplest approach to estimate intersection angles is to calculate intersection angles as the difference between the corresponding orientations of the local maxima in the histograms of directional response. However, we find that this is not accurate because the LKFs are generally wider around intersections than elsewhere and as a consequence, the maxima in the directional response are broad. Instead, we calculate the difference between the main orientation of those points that are furthest away from each centroid but still within the same intersection-point cluster (the red points in Fig. 2e, step 3 of section 3c). This method is designed to work robustly as long as the filter size is smaller than the typical distance between adjacent features, intersection points, and junctions. Moreover, it is to be expected that the estimated intersection angles are more uncertain than the estimated orientations because the former inherit twice the uncertainty in orientation plus uncertainties related to the feature classification.

### e. Distance and similarity measures

A point-by-point comparison between LKFs of two different sets, for example, model results and observations, typically is not useful because even slight spatial offsets between observed and simulated structures will be counted as a double error (Mohammadi-Aragh et al. 2018). Instead, the comparison of frequency distributions of orientations and intersections angles of LKFs within the analyzed geographical (sub)domains is more meaningful for assessing the similarity of two sets of LKFs. Therefore, we complement the MDA by using measures that compare modeled versus observed frequency distributions of orientations and intersection angles. From a large number of existing similarity measures for histograms, many of which are reviewed in Cha (2007), the following three appear to be particularly useful for our purpose:

The simplest metric we use is the Euclidean distance (Cha 2007),

where

*s*is the number of bins in the discrete frequency distributions and*P*_{num}and*P*_{obs}are modeled and observed frequency distributions.The main advantage of the Bhattacharyya coefficient (Kailath 1967),

over the Euclidean distance is that it is scaled:

*ρ*_{BhC}∈ [0, 1];*ρ*_{BhC}= 1 indicates identical distributions, whereas*ρ*_{BhC}= 0 for a pair of nonoverlapping distributions.The previous two metrics consider each bin in isolation and are not sensitive to, for example, the distance between the frequency maxima (e.g., the orientation angle) in two distributions. The Earth mover’s distance (EMD) (Rubner et al. 2000; Pele and Werman 2008) overcomes this issue:

where

*g*_{kl}is a flow from bin*k*to bin*l*and*r*_{kl}is a ground distance (in our case in orientation or intersection angle) between the histogram bins. For intersection angles*y*_{k}and*y*_{l}at the center locations of two bins,*r*_{kl}= |*y*_{k}−*y*_{l}|, and for orientations*r*_{kl}= max(|*y*_{k}−*y*_{l}|, 180° − |*y*_{k}−*y*_{l}|). As an example, consider two very long, narrow linear LKFs with orientation angles of 30° and 150°, which are different by 120°. With our choice or*r*_{kl}, the circular nature of the orientation angle is taken into account and the ground distance in the histogram is 60°. Our implementation of EMD is based on the Python library PyEMD [Pele and Werman (2008, 2009); code available at https://github.com/garydoranjr/pyemd].

## 4. Application and evaluation of the MDA

In this section we illustrate the application of the MDA, first with data from an idealized synthetic field of LKFs, and then we compare data retrieved from satellite observations to numerical model simulation data.

### a. MDA of idealized linear features

To test the performance of the suggested methods, we first construct a field of synthetic linear features with clearly defined orientations and intersection angles (Fig. 2a). The features are designed such that they reflect some properties typically found in observed LKF patterns. The grid spacing is designed such that the width of linear features are about 3–5 grid boxes, similar to typical widths found in observations (see section 4b below). Note that the LKF width in gridded datasets may be unrealistic due to gridding of the originally Lagrangian linear features on a discrete raster (the grid). In the same way the discretization in our example introduces inaccuracy into the gridded representation of the idealized linear features.

We define three rectangular filters starting from a reference size of 37.5 km by 37.5 km (length and width). This reference size corresponds to three by three grid points of the RGPS data with grid scale of 12.5 km. The filters differ only by the dilation factors *S*_{lf} and *S*_{wf} that describe the ratio of length and width of the dilated filter to the reference length and width. We use *S*_{lf} ∈ {3, 6, 9} and *S*_{wf} ∈ {1, 3}, resulting in filters with lengths of 112.5, 225.0, and 337.5 km and widths of 37.5 and 112.5 km (e.g., the green rectangle in Fig. 2a). Histograms of the directional response are generated by convoluting the filters with the synthetic data along the centerlines (the cyan lines in Fig. 2a). The histograms for the two example points P1 and P2 marked in Fig. 2a are shown in Fig. 2b.

Grid point P1 represents a line point and P2 point belonging to an intersection point. The histogram of the directional response at point P1 (Fig. 2b, magenta lines) has only one peak because the feature is locally linear for the applied aspect ratios. The histogram of directional response at point P2 has more than one maximum. For point P2 the response histogram for the filter with (*S*_{lf}, *S*_{wf}) = (3, 1) is smooth indicating that the filter is too small in comparison to the width of the linear feature (Fig. 2b, thin cyan line); for (*S*_{lf}, *S*_{wf}) = (6, 1), the histogram has two easily distinguishable peaks (Fig. 2b, medium cyan line); for (*S*_{lf}, *S*_{lw}) = (9, 3), the histogram exhibits three maxima because the filter is too large, and therefore it is influenced by the adjacent zigzag features around approximately 135° as well (Fig. 2b, thick cyan line).

To illustrate how the results of the algorithm depend on the parameter choices and how they correspond to our visual impression, Figs. 2c and 2d show the estimated main local orientation of the linear features for two different dilation factors (see also Figs. 1 and 2 in Singh et al. 2016). For example, the zigzag feature on the right-hand side of the idealized pattern as a whole is clearly oriented along 90°, but consists of several segments in two directions of 45° and 135°. The filter with *S*_{lf} = 3 (i.e., a length scale of 112.5 km, Fig. 2c) is so short that it detects the local signal and neglects the large-scale orientation. As a result, the main local orientation only reflects the direction of segments, but not that of the entire feature. However, if we apply larger dilation factors (corresponding to a length scale of 337.5 km, Fig. 2e), the measured orientation of the zigzag feature considers all segments as a single feature with a uniform orientation angle of 90°. This illustrates that the multiscale orientation estimation is useful to analyze the LKFs for different scales without changing the binary map of the features.

The fact that the largest filter with (*S*_{lf}, *S*_{wf}) = (9, 3) does not detect the local orientation of 45° and 135° of the zigzag feature is also apparent in the frequency distributions of orientations for the whole synthetic domain (Fig. 3a), where peaks at these angles exist only for *S*_{lf} = 3 and *S*_{lf} = 6. The largest difference of intersection angles to the reference is found for angles above 75° with (*S*_{lf}, *S*_{wf}) = (9, 3) (Fig. 3b). Generally, primarily intersection angle estimates agree best with the reference for (*S*_{lf}, *S*_{lw}) = (6, 1).

### b. MDA of high-resolution sea ice deformation

To assess the potential of our MDA in a realistic context, we use average sea ice deformation fields obtained from RGPS data (Fig. 1a) on a particular day (13 February 2001) with many linear features in the deformation fields and from a simulation of the MITgcm (Fig. 1b) of the same day. As previously reported (Spreen et al. 2017; Mohammadi-Aragh et al. 2018), there are far fewer LKFs in the MITgcm data than in the RGPS data (Figs. 1e,f). For the specific case in Fig. 1, the ratio of observed LKF density to simulated LKF density is 2.90. From visual inspection, we estimate that the LKFs in the RGPS data are predominantly oriented between −10° and 0° and at 30° and 60° relative to the horizontal axis of the plot. The LKFs in the MITgcm field are mostly oriented in directions around 60° and 150°, with some LKF oriented between 0° and 10° to the North of the Canadian Arctic Archipelago. The model also fails to generate the short, scattered LKFs. Thus, the MDA and the distance metrics [Eqs. (1)–(3)] should reflect this difference of simulated orientation of intersection angles and RGPS derived data.

The estimated local orientation of the LKFs for two different dilation factors (*S*_{lf} = 3, *S*_{lf} = 6) and *S*_{wf} = 1 are shown in Figs. 4a–d. Although the LKF pattern of the RGPS data is complex and dense, the multidirectional estimation of orientation using the larger dilation factor of *S*_{lf} = 6 measures the main orientation with elevated frequency between 30° and 60° and peaks at 0° and 170° in good agreement with our visual impression (Figs. 4a,c and 5a). The probability distribution of LKF orientations extracted from RGPS data (Fig. 5a) exhibits the same maxima, but superimposed on a uniformly high frequency, in particular between 20° and 70°. The estimated frequency distribution of orientations in the MITgcm data (solid lines in Fig. 5a) shows several maxima near orientation angles of 10°, 60°, 90°, and 150°. We conclude that our MDA can estimate plausible patterns in frequency distributions of orientations, even if the observed and simulated fields differ considerably in terms of LKF density.

The difference in orientation distributions between the pair of deformation fields as quantified by the distance metrics [Eqs. (1)–(3), section 3e] depends slightly on the bin size of the distributions and the dilation factor (Table 1, top left). With increasing bin size, the Euclidean distance increases, the Bhattacharyya coefficient stays the same, and the EMD fluctuates slightly. With increasing dilation factor *S*_{lf}, the Euclidean distance increases, Bhattacharyya coefficient stays the same, and the EMD decreases.

Figures 4e and 4f show the grid points in intersections of LKFs and the centroids of intersection clusters for *S*_{lf} = 3. We have left out the intersection angles below 17.5°, shown as green dots, in our further analysis because these points appear mostly not to be real intersection points but rather an artifact associated with the filter size. For intersection angles above 17.5°, the location and number of intersection centroids is largely consistent with our visual impression for both the observed and the simulated data. The LKFs in the MITgcm field intersect often at almost right angles. In contrast, the RGPS data contain numerous acute-angled intersections. For the MITgcm data, the frequency distribution of the estimated intersection angles of LKFs with *S*_{lf} = 3 increases toward 90°, whereas intersection angles in the RGPS data are more frequent toward small angles (Fig. 5b). Note that the relative frequency of small acute intersection angles is higher when we estimate intersection angles with a smaller dilation factor of *S*_{lf} = 3.

The Euclidean distance between the probability distributions of intersection angles for 13 February 2001 tend to increase with increasing bin size and dilation factor (Table 1, top right). The EMD increases with increasing dilation factor. In this example, the Earth mover’s distance is particularly sensitive to the value of *S*_{lf} for intersection angles.

So far, we used only one pair of deformation fields corresponding to one particular date to illustrate the functioning of the MDA. To obtain more robust results, we assess the spatial characteristics of LKFs with our MDA for a dataset of the MITgcm and RGPS data at six selected dates, first individually and then aggregating, to reveal similarities between individual and multidate angle distributions. This approach provides some information on overall characteristics of the model simulations versus the observation fields. Finally, to address the operational needs of daily verification, we provide and discuss time series of relevant statistics over a four-month period.

The averaged orientation angle histograms (Fig. 6a) are much smoother compared to the single-date histograms (Fig. 5a), except for conspicuous maxima around 0° and 90°, in particular for a small dilation factor (*S*_{lf} = 3). These maxima are identified as artifacts caused by the alignment of the spatial filter with the computational grid directions, more specifically the small aspect ratio between the filter size and the grid spacing. The otherwise uniform histograms indicate that there are no generally preferred orientations of LKFs. This statement is valid only for the Arctic as a whole, whereas we expect that preferred orientations caused by land–sea geometry would arise if only particular coastal regions were analyzed.

On a pan-Arctic scale and averaged over numerous dates, the frequency distributions of LKF main orientations are expected (and measured) to be close to uniform, so that the comparison of corresponding distributions between observations and simulations is not very informative. In contrast, systematic intersection angles are expected to occur. These angles are closely linked to the sea ice rheology (Ringeisen et al. 2019) and thus provide insight into the flow properties of sea ice or their representation in a model.

In our simulations, the MITgcm underestimates the frequency of small intersection angles up to 45° compared to the RGPS data (Fig. 6b). At the same time there are too many intersections close to 90° compared to the satellite data. The preference toward small intersection angles revealed by our analysis for RGPS data is consistent with earlier results for the western Arctic Ocean where angles between 20° and 50° were reported to occur most frequently (Marko 1976; Walter and Overland 1993; Cunningham et al. 1994; Schulson 2004; Wang 2007) and also consistent with laboratory measurements (Schulson et al. 2006). Note that we have excluded intersection angles below 17.5° because the corresponding points are not real LKF intersections in many cases (cf. Fig. 4).

To quantify the accuracy of the multidate frequency distributions, we have again applied the set of metrics introduced in section 3e. As expected, the generally uniform orientation distributions are more similar between RGPS and MITgcm for all metrics, dilation factors and bin sizes compared to the single-date comparison (cf. Table 1 upper-left to lower-left statistics). Also, when it comes to intersection angles, all distance measures reveal a higher degree of similarity for the multidate distributions compared to the single-date distributions.

Finally, we compute the EMD and other relevant statistics for each date individually over a four-month period. Not surprisingly, the mean sea ice deformation of the RGPS and MITgcm data are positively correlated, but the MITgcm exhibits weaker deformation on average (Fig. 7a). The same holds for the LKF density (positive correlation and lower mean in MITgcm). In both datasets, the LKF density covaries with the mean deformation (Figs. 7a,b).

Turning to the corresponding EMD time series, we find several dates in which the angle distributions are insufficiently sampled in the MITgcm data. Therefore, we consider EMD for orientation and intersection angle solely for cases in which the ratio of observed versus modeled LKF density is smaller than 5:1. Furthermore, we compute EMD benchmark values for each date by computing the EMD of the single RGPS LKF pattern at that date relative to each MITgcm LKF pattern for all dates excluding the one of the RGPS LKF pattern, and then averaging. The degree to which the EMD is smaller than the benchmark EMD reveals how well the simulation is able to discriminate orientation and intersection angles, meaning that day-to-day variations of LKF spatial characteristics are reproduced by the simulation.

Regarding LKF orientation, the EMD is mostly considerably lower (average: EMD_{$(Slf=3)$} = 9.3, EMD_{$(Slf=6)$} = 9.6) than the benchmark EMD (average benchmark: EMD_{$(Slf=3)$} = 13.0, EMD_{$(Slf=6)$} = 13.5) (Fig. 7c), with a few exceptions. This indicates that the day-to-day spatial characteristics of LKF patterns in the RGPS data tend to be reproduced by the MITgcm. However, this is much less obvious for the intersection angles, where the EMD (average: EMD_{$(Slf=3)$} = 7.6, EMD_{$(Slf=6)$} = 8.7) exceeds the benchmark EMD (average benchmark: EMD_{$(Slf=3)$} = 7.8, EMD_{$(Slf=6)$} = 8.6) in a larger number of cases, sometimes by a factor of more than 2 (Fig. 7d). This suggests that, compared to orientation, the EMD for intersection angles more strongly reflects systematic differences between simulation and observations than day-to-day variations. Note that for both orientation and intersection angles, the discriminative power (EMD benchmark divided by EMD) is slightly larger for the shorter spatial filter. We conclude that modern sea ice models like the MITgcm setup used here start to be able to reproduce basic characteristics of LKFs, such as their orientation, and how these characteristics vary in time. However, higher-order characteristics such as intersection angles and their temporal evolution are not yet reproduced, calling for still higher resolutions and further sea ice model development.

## 5. Summary and conclusions

In this study, we introduce a multiscale directional analysis (MDA) as a new method for evaluating the simulated spatial patterns of linear kinematic features (LKF) in high-resolution sea ice model simulations. We envisage the application of our MDA to assess sea ice dynamics on scales ranging from weather to climate in uncoupled (ice–ocean) as well as coupled (ice–ocean–atmosphere) simulations. An MDA-based evaluation will complement established sea ice forecast verification methods by adding local and scale-dependent spatial information of geomorphological patterns of LKFs.

The MDA has a number of attractive features that distinguish it from previous methods. The approach offers an effective and simple way to evaluate spatial characteristics of LKFs. In particular, the MDA quantifies properties such as orientation for each LKF-classified grid point individually instead of assigning one value for a whole segment of LKF. Using different spatial filters enables one to assess sea ice dynamics on different spatial scales while properly accounting for the finite width of LKFs. Shearlet-based methods may also provide an interesting alternative if the width of the LKFs is not important for the comparison.

State-of-the-art algorithms are components of the MDA, but the general framework allows substituting them with alternative methods that are tailored to the specific application. The MDA consists of five components: (i) LKF detection, (ii) estimation of main orientation, (iii) feature classification, (iv) estimation of intersection angles, and (v) distance and similarities measures. In the context of the example analysis with RGPS satellite data and MITgcm simulation data, we now discuss the advantages and limitations of each component of the MDA separately.

Using EDIS (see section 3a), a unified definition of an LKF in both modeled and observed fields is obtained. An image separation algorithm identifies the LKFs as curvilinear features from a binary gridded field. With the algorithm it is possible to distinguish signal (LKFs) from noise (e.g., isolated grid points also classified as LKF by a thresholding method), which is an important prerequisite not only for obtaining accurate results for local properties in subsequent steps but also for meaningful statistics when local properties are aggregated into regional histograms. In our analysis, EDIS shows that the MITgcm simulation generates fewer than observed LKFs (Figs. 1e,f). This is linked to smoother sea ice deformation fields in the MITgcm simulation data compared to the RGPS data. Higher resolution of the simulation will also increase the number of LKFs.

The main local orientations of LKFs are estimated using rectangular filters of different length for different spatial scales. Existing (object based) methods treat the influence of neighboring features either by masking other LKFs, by connecting them based on some maximum orientation difference, or by clustering the linear features that are in close proximity and oriented approximately in the same direction (e.g., Linow and Dierking 2017). In contrast, the MDA is designed to work reliably even in places with complex patterns of intersecting LKFs and avoids the cost of feature segmentation and semantic postprocessing. At the center of our MDA implementation is the multiscale directional orientation method of Singh et al. (2016). We have, however, refined the method to obtain more accurate orientation estimates by fitting a parabolic function around the maximum of the directional response histogram.

Intersection angles of sets of LKFs are estimated based on an algorithm that classifies the LKFs as intersection and line points. Compared to alternative methods, the clustering approach has the advantage of minimizing unwanted effects of localized irregularities, thereby enabling more accurate intersection angle estimates. However, our analysis also reveals that some points are erroneously detected as intersection points, which leads to inflated frequencies of small angles (i.e., below ~17.5°). The quantified errors (the ratio of the number of excluded intersection points to the total diagnosed intersection points) for the modeled and observed fields are 0.167 and 0.063.

Finally, pairs of frequency distributions are compared using a set of distance measures including the Euclidean distance, the Bhattacharyya coefficient, and the Earth mover’s distance (EMD). The Euclidean distance and Bhattacharyya coefficient are bin-by-bin measures and take into account only the correspondence between bins with the same index. These measures are sensitive to the resolution of the binning: coarse binning entails low discriminative power, whereas fine binning can lead to spuriously low estimates of similarity where frequency distributions might be just slightly shifted (Komaty et al. 2014). According to our analysis, a meaningful compromise is a bin size of 5° for these two metrics. In the end, the more sophisticated EMD is found to be a more useful and robust distance metric because it captures histogram differences, including shifts, in a natural manner.

We test the performance of the MDA by estimating the spatial characteristics of a synthetic case with a geometrically simple network of linear features. The qualitative evaluation and the quantitative comparison with the clearly defined geometrical characteristics show agreement with human subjective judgement but reveal sources of uncertainties: the discretized nature of the defined reference linear features, lack of scale-dependent reference, imperfect centerline tracing in intersection regions, and uncertainties associated with the clustering. The results show that errors of intersection angles are higher than errors of orientation estimates.

Our analysis was then applied to a case study of pan-Arctic model simulation and observational data. We show that the sea ice component of the MITgcm in our 4.5 km configuration generates considerably fewer than observed LKFs, with a moderate level of agreement between observed and simulated LKF orientations and intersection angles. In particular, the frequency of small acute intersection angles is too low in the model simulations compared to the observations (see also Hutter and Losch 2020; Ringeisen et al. 2019).

LKF orientation and intersection angles and their dependence on different sea ice rheologies (in particular yield curves) and grid resolution has received growing attention recently (e.g., Hutchings et al. 2005; Wang and Wang 2009; Ringeisen et al. 2019). In this regard, the MDA will help to shed light on the dynamical behavior of sea ice and its representation in models. The MDA might also be useful to further extend our knowledge about the formation and evolution of the complex observed geometrical features in sea ice (Marko 1976; Cunningham et al. 1994) or to characterize the linear features in order to diagnose properties of the atmospheric forcing (Zhao and Liu 2007; Barry et al. 1989; Marko and Thomson 1975). The atmospheric forcing plays a pivotal role in the generation and predictability of LKFs (Mohammadi-Aragh et al. 2018). Since the MDA is an evaluation method sensitive to direction and spatial scale, it can be used effectively for such studies as the very first steps toward measuring the actual predictability. Although the focus of our study is on spatial characteristics of LKFs in sea ice cover, we anticipate that this method is applicable to other branch of science. In summary, we envisage that the MDA will be useful to understand sea ice dynamics and to provide additional guidance for the development of sea ice models and forecast systems.

## Acknowledgments

This study has been supported through the AWI strategy fund. M.M.A. acknowledges the financial support by the German Federal Ministry for Education and Research (BMBF) and forms part of the GROCE project (Greenland Ice Sheet/Ocean Interaction, Grant 03F0778C). All model simulations were carried out on the computer facilities of the ECMWF within under the special project “Potential sea ice predictability with a high-resolution Arctic sea ice-ocean model.” H.F.G. acknowledges the financial support by the Federal Ministry of Education and Research of Germany (BMBF) in the framework of SSIP (Grant 01LN1701A). We thank Daniel Hodyss, Barbara Casati, and the two anonymous reviewers for their comments that improved this manuscript. We are grateful for comments by Wolfgang Dierking and Sergey Danilov.

## REFERENCES

*J. Geophys. Res.*

*IEEE Trans. Geosci. Remote Sens.*

*Ann. Glaciol.*

*Cryosphere*

*Meteor. Appl.*

*Int. J. Math. Models Methods Appl. Sci.*

*J. Geophys. Res.*

*Proc. Int. Geoscience and Remote Sensing Symp.*, Pasadena, CA, IEEE, 1747–1749

*Quart. J. Roy. Meteor. Soc.*

*IEEE Trans. Geosci. Remote Sens.*

*Numer. Math.*

*J. Geophys. Res. Oceans*

*Appl. Comput. Harmon. Anal.*

*Meteor. Appl.*

*J. Phys. Oceanogr.*

*Annu. Rev. Fluid Mech.*

*Quart. J. Roy. Meteor. Soc.*

*Geophys. Res. Lett.*

*Mon. Wea. Rev.*

*Cryosphere*

*J. Geophys. Res. Oceans*

*Cryosphere*

*IEEE Trans. Commun. Technol.*

*Sci. Rep.*

*IEEE Trans. Instrum. Meas.*

*Int. Conf. on Curves and Surfaces*, Osnabrück, Germany, University of Osnabrück, 416–430

*Scaling Laws in Ice Mechanics and Ice Dynamics*, J. P. Dempsey and H. H. Shen, Eds., Kluwer Academic, 315–323

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*Remote Sens.*

*Ocean Modell.*

*IEEE Signal Process. Mag.*

*Cryosphere*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*Environ. Res. Lett.*

*Ocean Sci.*

*Mercator Ocean Quarterly Newsletter*, No. 31, Mercator-Ocean, Ramonville Saint-Agne, France, 13–21

*J. Geophys. Res.*

*Sci. Rep.*

*HPCMP Users Group Conf.*, Denver, CO, IEEE, 280–282

*J. Geophys. Res.*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*Computer Vision—ECCV 2008*, Lecture Notes in Computer Science, Vol. 5304, Springer, 495–508

*12th Int. Conf. on Computer Vision*, Kyoto, Japan, IEEE

*Exp. Fluids*

*Cryosphere*

*Int. J. Comput. Vis.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Neuroinformatics*

*Proc. Natl. Acad. Sci. USA*

*Image Process. On Line*

*Cryosphere*

*Ann. Glaciol.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Ann. Glaciol.*

*Eos, Trans. Amer. Geophys. Union*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Geophys. Res. Lett.*

*Drift, Deformation, and Fracture of Sea Ice: A Perspective Across Scales*. SpringerBriefs in Earth Sciences, Vol. 83, Springer Science and Business Media, 83 pp

*IEEE Trans. Image Process.*

*Commun. ACM*

*J. Oceanogr.*

## Footnotes

Denotes content that is immediately available upon publication as open access.