On producing sea ice deformation data sets from SAR-derived sea ice motion

We propose a method to reduce the error generated when computing sea ice deformation fields from synthetic aperture radar (SAR)-derived sea ice motion. The method is based on two steps. The first step consists of using a triangulation of the positions taken from the sea ice trajectories to define a mesh on which a first estimate of sea ice deformation is computed. The second step consists of applying a specific smoother to the deformation field to reduce the artificial noise that arises along discontinuities in the sea ice motion field. This method is here applied to RADARSAT Geophysical Processor System (RGPS) sea ice trajectories having a temporal and spatial resolution of about 3 days and 10 km, respectively. From the comparison between unfiltered and filtered fields, we estimate that the artificial noise causes an overestimation of about 60 % of opening and closing. The artificial noise also has a strong impact on the statistical distribution of the deformation and on the scaling exponents estimated with multifractal analysis. We also show that a similar noise is present in the deformation fields provided in the widely used four-point deformation RGPS data set. These findings may have serious implications for previous studies as the constant overestimation of the opening and closing could lead to a large overestimation of freezing in leads, salt rejection and sea ice ridging.


Introduction
Sea ice motion can be retrieved from satellite synthetic aperture radar (SAR) images using cross-correlation techniques and feature tracking algorithms (Kwok et al., 1990;Fily et al., 1990;Hollands and Dierking, 2011). Sea ice deforma-tion is then estimated by computing the spatial derivatives of the sea ice motion. The most popular data set providing both sea ice motion and deformation is the RADARSAT Geophysical Processor System (RGPS) data set (Kwok, 1998). It covers the western Arctic for the period 1996-2008 at temporal and spatial resolution of about 3 days and 10 km, respectively.
As deformation determines sea ice opening (i.e., positive divergence) and closing (i.e., negative divergence), it may be used to estimate important global quantities, such as the ice production in leads, with some assumptions on sea ice growth and redistribution (Kwok et al., 1995). Using the RGPS data set, Kwok (2006) estimated that deformationrelated ice production is about 25-40 % of the winter ice production in both the perennial and seasonal ice zones. Kwok et al. (2008) also showed that the deformation-related ice production derived from the RGPS data set is up to 2 times higher than the one estimated by numerical models, implying a potential underestimate of the associated sea ice-ocean feedbacks.
In addition to essential information about sea ice opening and closing, the analysis of sea ice motion and deformation also gives a particular insight into the underlying physics controlling the sea ice dynamics and provides precious information with which to validate sea ice models. Marsan et al. (2004) described how the statistics of sea ice deformation vary as a function of spatial scale, while Rampal et al. (2008) generalized these scaling properties to both the spatial and temporal domains. Stern and Lindsay (2009) and Herman and Glowacki (2012) documented the seasonal and interan-664 S. Bouillon and P. Rampal: On producing sea ice deformation data (2009,2011) showed that classical sea ice models do not capture these statistical properties.
The estimation of these global quantities (e.g., total opening/closing) and statistical properties (e.g., spatial scaling exponents) may be impacted by errors in sea ice deformation data. Uncertainty on deformation is usually seen as a consequence of motion tracking errors that depend on the algorithm and parameters used. Lindsay and Stern (2003) estimated the standard deviation of the error in area change to be about 1.4 km 2 for a 10 by 10 km cell when the tracking error (i.e., tie point) is about 100 m. This error estimate is equivalent to the level of significance of 0.005 per day for 3-day intervals estimated by Kwok and Cunningham (2002), and used to determine the error on ice production as being less than 1 % of the total.
However, two other sources of error can be identified. Both are linked to the definition of the boundary of the cell (usually quadrangle) over which deformation is computed. Lindsay and Stern (2003) showed that unrealistic deformation is often obtained when this boundary is too irregular. Also, spurious openings and closings (that we will refer to as artificial noise hereafter) are caused by unfavorable orientation of the cell boundary relative to the discontinuities in the sea ice motion field, also called dynamic discontinuities, slip lines or linear kinematic features (LKFs). Lindsay and Stern (2003) evaluated the standard deviation of the error in area change due to the boundary definition to be about 3.2 km 2 for a 10 by 10 km square cell, which is more than twice the error from tracking mentioned above. Kwok (2006) stated that this artificial noise would lead to an overestimation of the ice volume production, although no precise number was given.  proposed to reduce this error by combining cells together, but this solution reduces the benefits of having highresolution data and reduces the spatial range over which one could perform scaling analysis. The error in estimating deformation has also been studied by Thorndike (1978) for large spatial scales (from 100 to 500 km), but such an analysis is only valid for homogeneous and isotropic fields and does not apply to small-scale deformation as explained in .
This paper proposes a method to avoid unrealistic values and to significantly reduce the noise obtained when computing sea ice deformation from SAR-derived motion and presents an example of its application to sea ice trajectories coming from the RGPS data set. The complete method is described in Sect. 2. In Sect. 3, we discuss the quality of the obtained deformation fields and analyze the impacts of removing the artificial noise on the estimated global opening/closing and on the spatial scaling of the deformation. Section 4 concludes the paper with a discussion on potential improvements of the method and on implications of our findings for the existing literature. Figure 1. Example of the divergence (a, c, e) and shear (b, d, f) obtained for the single-crack test case at a normalized resolution of 0.1. The relative displacement parallel and normal to the crack (black line) are set to 0.01 and 0, respectively. (a) and (b) correspond to the unfiltered deformation fields, (c) and (d) to the deformation fields filtered with a classical smoothing kernel and (e) and (f) to the deformation fields filtered with our smoother. Triangles in white show the kernel defined for the triangle in green. With both smoothers the kernel corresponds to cells that can be reached by crossing a maximum of n edges (here n = 7). The classical smoother takes all the cells into account, whereas our smoother only takes into account the cells whose deformation is above a given threshold.

Method
The method we developed is based on two steps. The first step consists of defining a mesh by doing a triangulation of a set of tracked points. For each individual triangular cell, the deformation is calculated using the motion of its three nodes estimated from the tracking procedure. The second step consists of applying a specific smoother to the obtained deformation fields to reduce the artificial noise.

Application to simple test cases
In order to present the method, we first define a simple setup on a square domain having a normalized area equal to 1. In this domain, tracked points are distributed uniformly with a mean distance d between them (see for example Fig. 1 with d = 0.1). d is hereafter called the normalized resolution.
In the first test case, a single crack is defined (black line in Fig. 1). This crack passes by the center of the domain and makes an angle θ with the horizontal x axis. We want to simulate a discontinuous displacement field that is induced by the presence of that crack. To do so, we keep the points located below the crack (lower part of the domain in Fig. 1) fixed, and we require the points above the crack (upper part of the domain in Fig. 1) to move with the same displacement. The two components of the imposed displacement, u p and u n , correspond to the displacement parallel and normal to the crack, respectively.
The first step of the method is to perform a Delaunay triangulation of these points to generate a mesh on which deformation is computed. The spatial derivatives of the displace-ment are obtained by calculating the following contour integrals as in Kwok et al. (2008) around the boundary of each triangle: where A is the cell area. For example, u x is approximated by where m = 3 and subscript m + 1 = 1. The shear shear and divergence div deformations are computed as In the case of a slip line, u n is set to 0. No opening or closing should occur and shear should have the same value along the crack. Figure 1a and b show the divergence and shear computed for such a case, with u p = 0.01 and u n = 0. The divergence field exhibits spurious positive (opening) and negative (closing) values along the slip line. The shear field also exhibits some noise, but that is hardly visible in the figure. This artificial noise generates an overestimation of the total opening (and closing). In our test cases, the opening (and closing) error is defined as the absolute difference between the total opening (and closing) computed from the geometry of the problem and the total opening (and closing) computed from the integration of the positive (and negative) divergence given by the method. Repeating the slip line experiment 100 times, with θ varying from − arctan(0.2) to + arctan(0.2) and with different meshes, we find that the root mean square (rms) error per unit crack is about 20 % of the sliding distance u p for both the opening and closing. In other words, with a 100 km long crack and a sliding distance of 1 km, the artificial opening (and closing) would be about 20 km 2 . It is particularly interesting to note that this error does not depend on the normalized resolution d (we tested with d equal to 0.1, 0.01 and 0.001).
When repeating the same test case with quadrangles instead of triangles, we found a rms error of about 18 % of the sliding distance u p . For comparison, Lindsay and Stern (2003) found an error per unit crack of about 15 % of the sliding distance, for a similar test case on a mesh made of square cells. This analysis shows that using triangles generates an increase of only about 10 % of the opening (and closing) error compared to using quadrangles. This increase of the error is minor compared to the advantages of using triangles. Triangulation methods are more flexible. It roughly doubles the number of deformation estimates and increases the resolution at which deformation is defined. For the rest of this paper, we thus only present results on triangular meshes, but the method could also be applied to other type of meshes.
In order to remove the artificial noise in the deformation fields, one could apply an isotropic smoother over all the cells of the mesh. We here denote C the list of all the cells, and for each cell c ∈ C we define the kernel K c ⊂ C as the subset of cells that can be reached by crossing a maximum of n edges. An example of kernel with n = 7 is shown in Fig. 1c and d. The size of the kernel is noted |K c |. For the example shown in Fig. 1c and d, |K c | is equal to 87. The components of the filtered deformation are then defined as an area-weighted average over the cells of the kernel. For example, the filtered value for u x on the cell c is defined as This method reduces part of the artificial noise but is not appropriate since it ruins the localization of the shear and adds unreal deformation to non-deforming cells. It also modifies the spatial scale at which the deformation is defined, resulting in a modification of the value of the shear along the crack.
With the single-crack case, the area-weighted average of the shear for the cells cut by the crack is found to be inversely proportional to n. We propose a better method based on the fact that the deformation is by nature constant along a linear kinematic feature. Averaging motion derivatives along these features could then filter out the noise without spoiling the information on the real deformation. Contrary to the isotropic smoother presented here above, the mean value of the shear obtained with the second method in the single-crack case does not vary as a function of n. In other words, the scale at which the deformation is defined remains constant with the second method, which can be seen as an anisotropic smoother. In the case of a regular mesh and a single crack aligned with the x axis, we verified that the area-weighted average of the shear along the crack is strictly constant whatever the value of n. For the rest of this study, we consider that our method does not modify the spatial scale at which the deformation is defined. In some cases this assumption may not be verified, for example, when two cracks intersect or when the treated cells do not belong to a linear kinematic feature. Two metrics will be used to analyze the error caused by intersecting cracks and to evaluate the proportion of treated cells that can actually be considered as cracks.
To detect the cells that are involved in the mapping of each linear kinematic feature, we define a threshold for total deformation ( 2 div + 2 shear ). Only the cells whose total deformation is above the threshold are selected to build the smoothing kernels (see tering is applied on the cells where deformation is below the threshold. We denote S the list of all the selected cells. For each cell s ∈ S, we define the kernel K s ⊂ S as the subset of cells that can be reached by crossing only selected contiguous cells and a maximum of n edges. |K s | is the size of the kernel. In the case of the single crack the size of the kernel is alway equal to 2n + 1, except for the kernels whose center is close to the boundary of the domain. The kernel size may then be as low as n + 1. Our method preserves the localization of the deformation by avoiding mixing the deformation between LKFs (i.e., cells where the deformation is intense) and the surrounding rigid plates (i.e., cells where deformation is almost 0). Moreover, the way the smoothing kernels are built using contiguous cells ensures that deformation between LKFs that are not connected will not be averaged together. The proposed method relies on two parameters: the deformation threshold that determines which cells are selected and parameter n that determines how far we extend the kernel. In our test cases, the threshold value is chosen to be small enough to select all the deforming cells. For application to real data, the choice of this parameter is critical and is detailed in Sect. 2.2. The impact of parameter n on the total error, defined as the sum of the opening and closing errors, is shown in Fig. 2 (line with disk symbols). This error, when normalized by the sliding distance u p , decreases from about 40 % to a residual error that depends on the normalized resolution. For a resolution of 0.1, the residual error (i.e., the error remaining for n > d −1 ) is about 5 % as shown in Fig Points below the principal crack are fixed. Points above the principal crack experience the same displacement u n , perpendicular to the principal crack (here u n = −0.0025) but have distinct tangent components: u p for the block on the left (here u p = 0.01) and u p for the block on the right of the secondary crack (here u p = 0.0125). Triangles in white show the kernel defined for the triangle in green. In this example, the parameter n is equal to 3. malized resolution, whereas the initial error does not depend on the normalized resolution.
The two other curves in Fig. 2 (with square and triangle symbols) correspond to the normalized errors found for experiments considering a secondary crack as shown in Fig. 3. The domain is now divided into three blocks. Points below the principal crack are still fixed. Points above the principal crack experience the same displacement u n , perpendicular to the principal crack, but have distinct tangent components u p for the block on the left or u p for the block on the right of the secondary crack.
To get one crack opening while the other is closing, u p is defined as u p − u n . The example in Fig. 3 is given for u p = 0.01 and u n = −0.0025, so the principal crack should be closing whereas the secondary crack should be opening. Before filtering, the computed divergence field is highly polluted by the noise. Once the deformation is filtered (here with n = 3), the divergence field better matches the expected opening and closing. At the intersection of the two cracks, however, the solution may be incorrect, as the method does not distinguish cracks when they intersect and thus averages deformation over cells belonging to different cracks. It should be also noted that at the intersections of two cracks the size of the kernel |K s | may be as high as 3n + 1 (for threebranch intersections as in Fig. 3) or 4n + 1 (for four-branch intersections).
This mixing of intersecting cracks explains why the normalized error (triangle and square symbols in Fig. 2), after having rapidly decreased for small n as in the single-crack case, starts to increase for larger n. This simple test case shows that the shape of this function depends on the ratio u n u p and that the optimal value for n would be 4 for u n u p = 1 8 and 2 for u n u p = 1 4 . From this analysis, we identify n = 3 as an optimal value as it is the only value for which the error is reduced by at least a factor of 3 in any of the test cases presented here. In real cases, to define an optimal value for n is more difficult as it would depend on the number of intersecting cracks and on the local ratio between divergence and shear. For this study, we chose to use a constant parameter n, and its reference value is fixed at n = 3. To validate the choice of the method's parameters (i.e., n and the deformation threshold), we present in Sect. 3 another metric based on a multifractal scaling analysis of the deformation fields.

Application to RGPS sea ice trajectories
The RGPS Lagrangian displacement product provides trajectories of sea ice "points" initially located on a 10 km regular grid (http://rkwok.jpl.nasa.gov/radarsat/lagrangian.html). The positions of these points are updated when two successive images are available and treated by the tracking algorithm. The time interval between two updates is typically 3 days. Spatial coordinates are given in the SSM/I polar stereographic projection, with the origin of the Cartesian grid located at the North Pole and the negative y axis aligned to the 45 • W meridian.
The RGPS Lagrangian deformation product provides the deformation of each cell (which is quadrangle) of the original grid. The deformation of a cell is updated each time the position of all its nodes are updated. This method has a serious problem because cells may become so distorted that spatial derivatives are ill-defined. As the RGPS deformation data set does not provide for each cell the position of its node, it is not possible to filter the data to avoid this problem. This problem is specific to the RGPS deformation product and would not appear if each pair of images was treated separately with its own grid as in the GlobICE Image Pair product (http: //www.globice.info) and in the ENVISAT Geophysical Processor System (EGPS) (http://rkwok.jpl.nasa.gov/envisat/).
To tackle these problems, we reprocessed the RGPS Lagrangian displacement product to build a new deformation data set called the RGPS Image Pair Product. We first identify the tracked points corresponding to each pair of images (i.e., the set of points whose position has been updated at the exact same date and with the same time interval). We generate a Delaunay triangulation of these points. Then we compute the deformation over what we consider as being wellshaped cells, i.e., only for triangles having an area between 5 and 400 km 2 , their angles higher than 5 • or all their edges shorter than 25 km. We also only keep meshes if they have at least 200 nodes, and we discard single and pairs of triangles that are not connected to other cells. Figure 4 shows an example of a mesh and a sea ice divergence field after the processing of the data corresponding to one pair of images. Artificial noise, characterized by a succession of highly negative and positive values, is clearly visible and, as expected, is mainly located along lines.
Using triangles instead of quadrangles roughly doubles the number of deformation estimates and increases the resolution of the deformation product up to 7 km. Another advantage is that triangulations can be made on any set of points (if they are not all aligned), which is not the case with quadrangulation (Bremner et al., 2001). If the tracked points are given on a regular grid, quadrangulation could be easily performed and could be preferred. However for most of the available data sets (for example GlobICE and EGPS), the data are not given on a grid but as a list of points. The method presented here based on triangles is then very flexible and can be applied to many different sources of data. In the next section, the unfiltered and filtered deformation fields obtained on triangular meshes are compared to the RGPS deformation fields. The smoothing procedure is not applied to the RGPS deformation fields because it requires knowing the neighbors of each cell and this information is not present in the RGPS Lagrangian deformation data set.
To apply the smoother, we first need to detect the cells that are suspected to map the location of linear kinematic features. Thomas et al. (2008) proposed using a threshold for the shear based on the level of noise resulting from the motion tracking error. Instead, here we use a fixed threshold based on the total deformation (as in the simple test case presented above) to give more weight to the cracks suffering from artificial divergence. Cells showing total deformation greater than the threshold are thus selected and others simply not taken into account for the filtering procedure. Figure 5 shows the unfiltered total deformation rate and the selected cells (those with their edges in black) for a threshold equal to 0.02 per day.
Decreasing the threshold increases the number of selected cells and finally leads to excessive smoothing. Increasing the threshold reduces the number of selected cells and finally splits linear features into disconnected pieces for which the smoother is not efficient anymore. Indeed if a kernel only contains one cell, the smoother does not modify the value of the deformation over that cell.
To quantify the effect of this threshold on the quality of the selection, we define an index based on the size of the smoothing kernels |K s | (i.e., the number of cells that can be reached by crossing only selected cells and a maximum of n edges). As explained in Sect. 2.1, the size of the kernel |K s | could be as low as n + 1 for single cracks when the center of the kernel is at the boundary of the mesh and as high as 4n + 1 for two intersecting cracks. We then define the quality index as the percentage of treated cells having a kernel size between n + 1 and 4n + 1. For the example of Fig. 5, the quality index is equal to 89 %.
We explored the sensitivity of this quality index to the threshold value for the entire winter season 2006-2007 and with the parameter n equal to 3, which is the reference value defined in Sect. 2.1. For deformation thresholds equal to 0, 0.01, 0.02, 0.03, 0.04 and 0.05 per day, median quality indices are equal to 33, 78, 78, 76, 74 and 72 %, respectively. Based on this quality index, the threshold values 0.01 and 0.02 per days are the best. The value of 0.02 per day is chosen as the reference value for the deformation threshold. To quantify the range of the quality index obtained with this reference value, we look at the percentage of pairs of images for the entire winter 2006-2007 for which the quality index is lower than 50 %; our finding is that only 14 % of the pairs of images have a quality index lower than 50 %. To further validate the choice of the model parameters, a metric on a multifractal scaling analysis of the deformation fields is proposed in Sect. 3.  Figure 6 shows the sea ice divergence field after the application of the smoother with the parameter n equal to 3. Compared to the unfiltered divergence field shown in Fig. 4, the filtered field exhibits much less artificial noise and its interpretation is now much easier.

Results and discussion
In this section, we compare the original RGPS deformation data to the unfiltered and filtered versions of our RGPS Image Pair data set. A consistency check based on spatial scaling analysis is proposed and the differences between the three data sets in terms of spatial scaling and total opening/closing are discussed.
To compare the original RGPS deformation data with the unfiltered and filtered deformation data produced by our method, we generate composite pictures of the deformation rates for specific periods. The periods have to be long enough to ensure a good spatial coverage but not too long so as to avoid mixing incoherent information. For this study, we select the data for which the time of the first and second images -noted t k−1 and t k , respectively -are within a period of 8 days centered on a target date and for which the time interval, t = t k − t k−1 , is between 1 and 6 days. For the RGPS data set, we add a criterion to reject cells larger than 400 km 2 .
Selected cells may cover the same area but correspond to different dates and time intervals. This redundancy may impact statistical distribution and scaling analysis, so we apply a second selection step. We first define a regular grid at a resolution of 20 km. For each box of this grid, we find the cells whose center is in the box and we keep only those whose date, defined as (t k + t k−1 )/2, is the closest to the target date. This selection step creates some gaps in the coverage but is necessary to ensure a minimum consistency of the composite fields. Note that no averaging or interpolation is done during the generation of the composite deformation fields. Figure 7 shows the divergence rate for the period 2-10 February 2007 given by the RGPS Lagrangian deformation data set. Some features are so polluted by a succession of highly negative and positive values that it is very difficult to identify where cracks are opening, closing or sliding. Figure 8 shows the unfiltered divergence rate for the same period obtained after the first step of our method. As in the RGPS data set, the artificial noise is important and mainly located along linear kinematic features. Figure 9 shows the filtered divergence rate obtained with a deformation threshold of 0.02 per day and n = 3. The reduction of the noise makes much easier the identification of the opening and closing cracks.

Consistency check
To evaluate our method in a more quantitative way, we propose a metric based on a spatial scaling analysis. Scaling analysis is a powerful tool to characterize sea ice dynamical behavior and has been successfully used in previous studies to reveal the power-law scaling of sea ice deformations Rampal et al., 2008).
The spurious noise in the deformation fields corresponds to high values of deformation and is potentially present for any active linear kinematic features. This noise may then impact the distributions of shear and absolute divergence and modify their mean (first-order moment) but even more their standard deviation (second-order moment) and skew- ness (third-order moment). Moreover this noise is the highest at the resolution of the data but rapidly decreases for larger spatial scales. We thus expect that the presence of noise in the deformation fields will have a strong impact on the result of the scaling analysis, especially for the smallest scales and the highest-order moments of the distribution. Indeed, we assume that the power-law model for the spatial scaling of sea ice deformations has no physical reason to not hold over several orders of magnitude. This assumption is based on Weiss and Marsan (2004), who showed that the powerlaw model for the spatial scaling of the open-water density, which can be directly related to sea ice divergence, is valid down to 0.2 km. Therefore, any significant departure from the power-law model when approaching the spatial resolution of the data could be seen as an indicator of the remaining noise in the deformation field.
To perform the scaling analysis of sea ice deformation, we implemented a coarse graining method similar to the one proposed by Marsan et al. (2004) and applied it to the unfiltered and filtered versions of our RGPS Image Pair data set. Sea ice shear and absolute divergence rates are computed at different spatial scales ranging from 7 to 700 km. For the lowest scale, which is also the scale of the triangular cells, all the cells are taken into account. For the other scales, the coarse graining procedure covers the domain with boxes of different sizes (14, 28, 56, 112, 224, 448 and 896 km). The boxes actually overlap since a distance equal to half the box size separates their respective centers. For each box, we select the cells that have their center in the box. When the sum of the area of those cells is greater than half the box area, the deformation over the box is defined by averaging the spatial derivatives of the displacement weighted by the surface of each cell. The spatial scale for this new estimate of the deformation is the square root of the sum of the cell areas. The shear and absolute divergence for each box are then reported as a function of the spatial scale on a log-log plot (see Fig. 10 for the absolute divergence rate). The mean value ˙ L (where˙ L is either the shear rates or the absolute divergence rates, computed at scale L) relates to scale L following a power law. The powerlaw exponent is evaluated by applying a linear regression of the logarithm of ˙ L vs. the logarithm of L. Due to the finite size of the domain the power-law model is not expected to hold for the largest scales. For this reason we restrict the power-law regression of the data to the spatial scales 7, 14, 25, 50, 100 and 200 km.
The filtered shear and absolute divergence closely follow the power-law model for the spatial scaling as their firstorder moments are well aligned with the power-law fit for the spatial scales ranging from 7 to 200 km (see right panel of Fig. 10 for the absolute divergence). This is not the case for the unfiltered deformation fields (see left panel of Fig. 10), and we explain this strong departure from the power-law model by the presence of artificial noise. If the power-law fits were computed only from 50 to 200 km, spatial scaling exponents would be similar for the filtered and unfiltered data. Furthermore, the other moments ˙ q L of the distributions (see Fig. 11 for the absolute divergence rate, with q, the moment order, ranging from 0.5 to 3) computed from the unfiltered deformation fields also exhibit a strong departure from the power law, whereas the moments computed from the filtered deformation fields are well aligned with the power-law fits.
By definition, if a scaling holds for a given range of scales, it should be respected for any pair of scales within this range. To evaluate the deviation from the power-law scaling, we compute the power-law exponents for each pair of successive spatial scales (i.e., from 7 to 14 km, from 14 to 25 km, and so on) and we take the minimum and maximum values of these exponents. Those values as well as the exponents previously obtained with the whole range from 7 to 200 km are reported as a function of the moment order q in Fig. 12. The relationship between the power-law exponents and the moment order q is called the structure function β(q) and is defined by ˙ q L ∼ L −β(q) . The minimum and maximum exponents define the bars around β(q).
To check that the reference values for the model parameters are well chosen, we look at the deviation from the power law. This deviation is evaluated by the min-max error. For each moment order, the min-max error is defined as the difference between the minimum and maximum exponents obtained with any pairs of spatial scales within the defined range. It other words, it is the length of the bar drawn in Fig. 12. Applied to the composite fields used as an example here, we find that using a threshold for the total deformation Structure function β(q) corresponding to the exponents of the power-law relationship between the absolute divergence rate and the spatial scale: |˙ div | q ∼ L −β(q) . The bars on the graph indicate the deviation from the power law as they correspond to the minimum and maximum power-law exponents obtained for two successive spatial scales. of 0.02 per day gives the lowest min-max error for the highest order. For the parameter n, the lowest min-max error for the highest moment order is obtained with n = 2, but n = 3 is better for the other moment orders. The application of our metric to this single example tends to indicate that the reference values are well chosen. However, this metric should be applied to a larger number of examples to really identify the best values for the parameters.

Discussion
Comparing the original RGPS deformation to the unfiltered deformation allows us to evaluate the impact of using a triangulation to define well-shaped triangular cells. As in Lindsay and Stern (2003), we observe unrealistic values for the shear and divergence rates retrieved from the RGPS deformation data set. For the period 2-10 February 2007, the composite picture made from RGPS has maximum opening, closing and shear rates equal to 1.73, −6.73 and 66.47 per day, respectively. These extreme values arise from highly distorted cells. A very small fraction of the data set is polluted by these unrealistic values; however it has a high impact on the multifractal scaling analysis, particularly when looking at the highest moment orders of the distributions. In Marsan et al. (2004) and Stern and Lindsay (2009), additional constraints on the initial and current size of the cells were applied and the cells with total deformation higher than 1 per day were not taken into account. In many other studies based on the RGPS deformation data set, the presence and impacts of these unrealistic extreme values are simply not discussed.
The simple fact of redefining a new mesh from the actual position of the RGPS nodes allows us to avoid badly shaped cells and then to significantly reduce the number and magnitude of extreme values. For the same period, the composite picture obtained from the unfiltered version of our RGPS Image Pair data set has maximum opening, closing and shear rates equal to 0.63, −1.17 and 1.97 per day, respectively. The smoother also logically decreases the extreme values. For this example, the filtered composite picture has maximum opening, closing and shear rates equal to 0.13, −0.20 and 0.73 per day, respectively.
Comparing the filtered and unfiltered deformation allows us to analyze the impact of the artificial noise. From the scaling analysis for the total deformation, shear and absolute divergence, we found that the scaling exponents estimated from the unfiltered fields are systematically larger in magnitude by about 100 % for the absolute divergence and by about 50 % for the shear and total deformation. In the example corresponding to Fig. 10, the power-law exponent for the absolute divergence is −0.38 for the unfiltered field, instead of −0.20 for filtered data (for the shear: −0.17 instead of −0.1; for the total deformation: −0.19 instead of −0.12). For each moment, we observe that using unfiltered data leads to a systematic overestimation of the scaling exponents of about 100 % for the absolute divergence and 50 % for the shear and total deformation. We also performed the multifractal scaling analysis on the original RGPS deformation data set with the same constraints on the data as in Stern and Lindsay (2009), and we found that the departure from the power law is similar to the one observed for the unfiltered deformation data set.
The impact of the artificial noise is also seen on the structure function β(q) (see Fig. 12). As a consequence of the systematic overestimation of the scaling exponents, its curvature, which indicates the degree of multifractality of the deformation fields, is found to be twice as high for the unfiltered divergence field (0.20) than for the filtered one (0.11). For the shear and total deformation, the overestimation of the curvature is about 50 %, with a value of 0.16 instead of 0.1.
Differences are also seen in the cumulative distribution of the closing and opening rates (see Fig. 13 for the period 2-10 February 2007). Differences between the RGPS and the unfiltered deformation may be due to differences in the coverage and the selection of the data, but they also come from the difference in resolution (10 km for the RGPS instead of 7 km) and from the impact of distorted cells included in the RGPS data set. Differences between the filtered and unfiltered deformation induce a modification of the shape of the distribution. The distribution of the filtered divergence field is closer to an exponential distribution (linear in the semi-log plot), while the distribution of the unfiltered divergence field is clearly a stretched exponential.
Finally, we compare the three data sets by computing the total area that has been opened and closed. For the original RGPS deformation data, 40 000 km 2 has been opened during the period 2-10 February 2007, whereas 39 000 km 2 have been closed. For our unfiltered data, we find lower values of 30 000 and 38 000 km 2 , respectively. For the filtered data these numbers drastically drop down to 15 000 and 24 000 km 2 , respectively. In this example the artificial noise is thus responsible of an overestimation of the opening and closing of about 100 and 60 %, respectively. Over the entire winter season 2006-2007, the cumulative opening and closing are both 60 % higher in the unfiltered data than in the filtered data.

Conclusions
A method is proposed to derive accurate sea ice deformation fields from SAR-derived motion products. The first step of the method consists of a triangulation of the tracked points to generate a mesh of triangular cells on which a first estimate of deformation is computed. The second step consists of applying a smoother to the deformation fields. The method relies on two parameters: a deformation threshold and the size of the smoothing kernel.
By applying the method to idealized test cases, we show that using triangles instead of quadrangles induces an increase of only about 10 % of the opening and closing error, whereas our smoothing method reduces the opening and closing error by at least a factor of 3 in any of the test cases presented here. The sensitivity to the value of the threshold used to detect deformation features is analyzed with a quality index, and the efficiency of our method is assessed using a metric based on a spatial scaling analysis and comparison between the unfiltered and filtered deformation fields.
The proposed method is used to produce a new deformation data set called RGPS Image Pair Product. Compared to the RGPS deformation data set, the RGPS Image Pair data set does not exhibit unrealistic large values caused by badly shaped cells. Moreover, our method drastically reduces the artificial noise arising along dynamic discontinuities.
By comparing the unfiltered and filtered deformation fields for winter 2006-2007, we estimate that this artificial noise may cause an overestimation of the opening and closing of about 60 %. We also estimate that the spatial scaling exponents as computed in Marsan et al. (2004) and Stern and Lindsay (2009) could have been overestimated by about 100 % for the absolute divergence and by about 60 % for the shear and total deformation.
The findings of the present study indicate that errors in sea ice deformation fields retrieved from SAR-derived motion may have been strongly underestimated, leading to po-tential significant biases on the estimates of sea ice production, salt rejection and sea ice ridging that one may find in the literature.
The method proposed here is applicable to other sea ice drift data sets, as provided, for example, by the GlobICE project. The method can handle Lagrangian trajectories or displacement between pairs of images. The same method could be applied to buoy trajectories when their spatial resolution is high enough, as with nested arrays of buoys (Hutchings et al., 2011(Hutchings et al., , 2012. The method proposed here could be modified to better manage intersecting cracks and to adapt its parameters depending on the local fields. However, substantial improvements may also come by combining, within tracking algorithms, the detection of dynamic discontinuities and the computation of sea ice deformation as proposed by Thomas et al. (2008). A complete validation using independent data sets should also be carried out.