Estimating statistical errors in retrievals of ice velocity and deformation parameters from satellite images and buoy arrays

The objective of this note is to provide the background and basic tools to estimate the statistical error of deformation parameters that are calculated from displacement fields retrieved from synthetic aperture radar (SAR) imagery or from location changes of position sensors in an array. We focus here specifically on sea ice drift and deformation. In the most general case, the uncertainties of divergence/convergence, shear, vorticity, and total deformation are dependent on errors in coordinate measurements, the size of the area and the time interval over which these parameters are determined, as well as the velocity gradients within the boundary of the area. If displacements are calculated from sequences of SAR images, a tracking error also has to be considered. Timing errors in position readings are usually very small and can be neglected. We give examples for magnitudes of position and timing errors typical for buoys and SAR sensors, in the latter case supplemented by magnitudes of the tracking error, and apply the derived equations on geometric shapes frequently used for deriving deformation from SAR images and buoy arrays. Our case studies show that the size of the area and the time interval for calculating deformation parameters have to be chosen within certain limits to make sure that the uncertainties are smaller than the magnitude of deformation parameters.

1 Introduction

Sea ice drifts under the influence of wind and ocean currents. Spatial gradients in the sea ice motion lead to distortion of the sea ice cover, termed deformation. The retrieval of sea ice drift vectors and deformation parameters from pairs or sequences of satellite synthetic aperture radar (SAR) images has gained increased attention during recent years because of the growing availability of suitable data (e.g., Stern and Moritz, 2002; Karvonen, 2012; Berg and Eriksson, 2014; Komarov and Barber, 2014; Lehtiranta, 2015; Muckenhuber et al., 2016; Demchev et al., 2017; Korosov and Rampal, 2017). Sea ice kinematics is also studied based on data from arrays of buoys or GPS receivers (e.g., Lindsay, 2002; Hutchings et al., 2008; Hutchings et al., 2012; Itkin et al., 2017), which in addition can serve as reference in comparisons to motion vectors obtained from SAR images. The knowledge of spatially detailed motion and deformation fields is potentially useful in ice navigation to locate divergent or compressive ice areas, as complementary information for operational sea ice mapping, for validation of models for forecasting of ice conditions, and for assimilation into ice models (Karvonen, 2012). Such practical applications require that the errors of the retrieved drift and deformation parameters are known. For buoys, errors in drift measurements depend on the accuracy of position and time readings. The accuracy of deformation parameters is not only affected by errors in drift magnitude and direction but also by the size and shape of buoy arrays (e.g., Hutchings et al., 2012; Griebel and Dierking, 2018). Drift vectors derived from pairs of satellite images are the result of correlation techniques or object detections, while deformation parameters are calculated from spatial arrangements of adjacent drift vectors surrounding the area of interest, in a manner that is independent of the coordinate system. This means that drift and deformation errors do not only depend on the geolocation accuracy and spatial resolution of satellite images but also on the reliability and robustness of the drift retrieval algorithm. In this technical note we focus on the estimation of statistical errors for ice velocity and deformation. The issue of error estimation was repeatedly addressed in the past, scattered in a number of publications and restricted to single aspects related to the respective analysis (e.g., Lindsay and Stern, 2003; Hollands and Dierking, 2011; Bouillon and Rampal, 2015; Hollands et al., 2015; Linow et al., 2015; Griebel and Dierking, 2018), and is also addressed in a more recent analysis by Bouchat and Tremblay (2020). Our motivation is to provide the mathematical background, together with examples of applications and discussions of validity, in a broader context. We emphasize that here we deal with statistical errors but not with boundary definition errors as described, e.g., in Lindsay and Stern (2003), Bouillon and Rampal (2015), and Griebel and Dierking (2018). Although this note is specifically focused on retrievals of parameters characterizing sea ice kinematics, the mathematical framework is also applicable to movement and deformation of ice shelves and glaciers, or for model simulations of sea ice, glacier, and ice sheet dynamics.

In Sect. 2 we summarize the basics and provide equations for calculating errors of drift and deformation parameters: divergence, vorticity, shear, and total deformation. The equations are used in Sect. 3 to quantify the influence of different parameters such as geolocation and tracking errors, or shape and size of buoy arrays and grid cells. Conclusions are presented in Sect. 4.

2 Errors of drift and deformation parameters

In this section, we provide a short description of the estimation of errors and the computation of strain rates and then derive the statistical errors for drift velocity, polygon areas, divergence, shear, vorticity, and total deformation. The statistical errors quantify uncertainties that are introduced by random fluctuations in the measurements. If the random fluctuations are small, data are measured with a high degree of precision but not necessarily with high accuracy. The latter requires that the measured value is close to the true value, whereas precision refers to the reproducibility of a measurement (Bevington and Robinson, 2003, chap. 1).

2.1 Error propagation and calculation of deformation

The formula for error propagation is based on the splitting method, i.e., the decomposition of a measured variable x into its true value and the measurement error: $x=x_{\mathrm{true}}+x_{\mathrm{error}}$, where x_true is considered to be a constant and x_error is a random variable with expected value E(x_error)=0 and variance $E\left(x_{\mathrm{error}}^{\mathrm{2}}\right)={\mathit{\sigma }}^{\mathrm{2}}$. If a quantity Q is calculated from measured variables x_k, i.e., $Q=f(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_n)$, a Taylor series expansion can be applied to estimate the error of Q. Usually only the linear term is retained:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}Q=& f\left(x_{\mathrm{1},\mathrm{true}},x_{\mathrm{2},true},\mathrm{\dots },x_{n,\mathrm{true}}\right)\\ & +\sum _{i=\mathrm{1}}^n\left[{\displaystyle \frac{\partial f}{\partial x_i}}\left(x_{\mathrm{1},\mathrm{true}},x_{\mathrm{2},\mathrm{true}},\mathrm{\dots },x_{n,\mathrm{true}}\right)\right]\left[x_{i,\mathrm{error}}\right]\end{array}.\end{array}$$

The variance is obtained by moving the first term to the left-hand side, squaring both sides, and applying the expected value operator E() (Bevington and Robinson, 2003). This operation results in

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\sigma }}_Q^{\mathrm{2}}=& \sum _i{\left({\displaystyle \frac{\partial f}{\partial x_i}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_i^{\mathrm{2}}+\sum _{i\ne j}\sum \left({\displaystyle \frac{\partial f}{\partial x_i}}\right)\left({\displaystyle \frac{\partial f}{\partial x_j}}\right){\mathit{\sigma }}_{ij}\\ & i=\mathrm{1},n,j=\mathrm{1},n,\end{array}\end{array}$$

where ${\mathit{\sigma }}_i^{\mathrm{2}}$ is the variance of x_i and σ_ij the covariance of x_i and x_j. If we can assume that the errors are uncorrelated, the second term on the right side of Eq. (2) is zero. We will use the notation “uncertainty” synonymously with “standard deviation of the absolute error”.

Deformation parameters are calculated from different combinations of the components of the velocity gradient tensor ($\partial u/\partial x,\partial v/\partial x,\partial u/\partial y,\partial v/\partial y)=(u_x,v_x,u_y,v_y)$ (Leppäranta, 2011), here given in a Cartesian coordinate system, where u(x,y) and v(x,y) are the velocity components in the x and y direction at position (x,y). We have

$$\begin{array}{}\text{(3a)}& {\displaystyle }\text{divergence}{\dot{\mathit{\epsilon }}}_{\mathrm{div}}=u_x+v_y,\text{(3b)}& {\displaystyle }\text{vorticity}{\dot{\mathit{\epsilon }}}_{\mathrm{vrt}}=v_x-u_y,\text{(3c)}& {\displaystyle }\text{shear}{\dot{\mathit{\epsilon }}}_{\mathrm{shr}}=\sqrt[\mathrm{2}]{{\left(u_y+v_x\right)}^{\mathrm{2}}+{\left(u_x-v_y\right)}^{\mathrm{2}}},\text{(3d)}& {\displaystyle }\text{and total deformation}{\dot{\mathit{\epsilon }}}_{\mathrm{tot}}=\sqrt[\mathrm{2}]{{\dot{\mathit{\epsilon }}}_{\mathrm{div}}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{\mathrm{shr}}^{\mathrm{2}}}.\end{array}$$

Divergence and shear are the two invariants of the symmetric deformation tensor. The dimension of $\dot{\mathit{\epsilon }}$ is velocity change per length unit, hence [time]⁻¹. For ease of reference, we briefly repeat the physical meaning of different velocity gradient combinations (after Cuffey and Paterson, 2010; Leppäranta, 2011): imagine a rectangle with its sides L_x and L_y parallel to the x and y axes of a 2D Cartesian coordinate system. In this case the gradients u_x and v_y are normal strain rates, leading to an extension or contraction of the rectangle in the respective direction. The normal strain along the x axis, e.g., is $\mathrm{\Delta }{L}_x\left(t\right)/L_x=u_x\mathrm{\Delta }T$. Here ΔT is the time interval $\mathrm{\Delta }T=t-t_{\mathrm{0}}$ during which the effect of deformation is analyzed, and L_x+ΔL_x is the side length at time t₀+ΔT. The sum u_x+v_y is the divergence or convergence, dependent on the sign. The expression u_y+v_x is linked to the change in shape of the rectangle (pure shear). The normal shear, u_x−v_y, quantifies the change in length difference between the sides of the rectangle. The vorticity (v_x−u_y), which is twice the rotation rate, describes the rotation about an axis vertical to the x–y plane (positive counterclockwise) without change in shape. Let the rectangle be located in a temporally constant velocity field with, e.g., u_x=0.1 d⁻¹, v_y=0.05 d⁻¹, u_y=0, v_x=0; then the divergence is ${\dot{\mathit{\epsilon }}}_{\mathrm{div}}=\mathrm{0.15}$ d${}^{-\mathrm{1}}=\mathrm{15}$ % d⁻¹. Assuming that the sides of the rectangle are L_x and L_y at time t₀, its area A₀=L_xL_y increases to (L_x+u_xL_xΔT) ($L_y+v_y{L}_y\mathrm{\Delta }T)=A_{\mathrm{0}}(\mathrm{1}+u_x\mathrm{\Delta }T\left)\right(\mathrm{1}+v_y\mathrm{\Delta }T)=\mathrm{1.155}{A}_{\mathrm{0}}$ for ΔT=1 d. Since only the difference u_x–v_y contributes to the square root (Eq. 3c), ${\dot{\mathit{\epsilon }}}_{\mathrm{shr}}=\mathrm{0.05}$ d${}^{-\mathrm{1}}=\mathrm{5}$ % d⁻¹ is the normal shear (Hutchings et al., 2012).

The deformation of a region R (covered by the buoy array or grid cell) with area A is calculated from the spatial averages of the velocity gradient components over the region R, in Eq. (4) indicated by an overbar. For the u_x component, for example, the expression is (Thorndike, 1986)

$$\begin{array}{}\text{(4)}& \overline{u_x}={\displaystyle \frac{\mathrm{1}}{A}}\int \underset{R}{\int }{\displaystyle \frac{\partial u}{\partial x}}\mathrm{d}a={\displaystyle \frac{\mathrm{1}}{A}}\underset{C}{\oint }u\mathit{n}{\mathit{e}}_x\mathrm{d}l.\end{array}$$

Here da and dl are the differentials for area and length, n is the outward normal vector to the perimeter C of R, and e_x is the unit vector in the x direction. This is Green's theorem, which relates a line integral along a closed curve C to the area integral over a plane region R bounded by C. The application of the theorem requires that the velocity components u and v have continuous first-order partial derivatives on R. In a Cartesian coordinate system, the calculation of the velocity gradient in the x direction is carried out using

$$\begin{array}{}\text{(5)}& u_x={\displaystyle \frac{\mathrm{1}}{A}}\underset{C}{\oint }u\mathrm{d}y\cong {\displaystyle \frac{\mathrm{1}}{\mathrm{2}A}}\sum _{i=\mathrm{1}}^n\left(u_{i+\mathrm{1}}+u_i\right)\left(y_{i+\mathrm{1}}-y_i\right)\end{array}$$

and the other components of the velocity gradient tensor accordingly. In Eq. (5) we have omitted the overbar above u_x. The sum comes from the trapezoid rule for integration, taking n points around the perimeter of R, where ($u_{i+\mathrm{1}}+u_i)/\mathrm{2}$ is the estimate of u on the ith segment, ($y_{i+\mathrm{1}}-y_i$) is dy, i is the summation index which traces the boundary in a counterclockwise sense, n is the number of vertices for the grid cell (or number of buoys), and A is the area of the grid cell (or of the polygon spanned by the buoy array). Here, $u_{n+\mathrm{1}}\equiv u_{\mathrm{1}}$ and $y_{n+\mathrm{1}}\equiv y_{\mathrm{1}}$ (closed polygon). The velocity gradients are implicitly averages over R. This will also be the case for our estimates of the deformation parameters, Eqs. (3a)–(3d).

The velocity vectors may be obtained from an array of buoys, where the buoys' positions are regarded as the vertices of a polygon. The displacement of a buoy is usually calculated from the distance between distinct positions, and the velocity is determined as the displacement divided by the time period between position fixes. When using pairs of satellite images, sea ice deformation is obtained from the displacements of recognizable structures or patterns in these images. These are referred to as ice structures from here on. In the reference image, a grid can be constructed by connecting the center positions of adjacent ice structures by lines. If movements of single ice structures differ between acquisitions of image 1 and image 2, the shapes and sizes of grid cells have changed in the second image. It is the presence of velocity gradients due to locally varying physical forces that causes the deformation. In practice the movement of sea ice is obtained using different methods (e.g., Holt et al., 1992; Stern and Moritz, 2002; Karvonen, 2012; Muckenhuber et al., 2016; Korosov and Rampal, 2017), which determine the spatial distribution and density of the displacement vectors. The vectors can be regularly spaced on the crossing points of horizontal and vertical grid lines as a result of pattern matching algorithms in an Eulerian approach, or they can be irregularly distributed, which is typical for the Lagrangian approach applied in feature or buoy tracking (see Fig. 1).

Figure 1 Eulerian grids (a, b) are reinitialized at every time step to a regular configuration. Lagrangian grids (c) evolve over time without being reinitialized.

The errors discussed in the following subsections can be traced back to errors in the position of reference points (i.e., vertices of a grid, or buoys). Lindsay and Stern (2003) denote this error type as geolocation error. On a horizontal plane two coordinates (e.g., x, y or latitude, longitude) determine the positions of the start and end points of the displacement, respectively. The distance $d=\sqrt[\mathrm{2}]{{\left(x^{\prime }-x\right)}^{\mathrm{2}}+{\left(y^{\prime }-y\right)}^{\mathrm{2}}}$ is prone to the errors of the coordinate readings. Its uncertainty is ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$, assuming ${\mathit{\sigma }}_{\mathrm{coord}}={\mathit{\sigma }}_x={\mathit{\sigma }}_y={\mathit{\sigma }}_{x^{\prime }}={\mathit{\sigma }}_{y^{\prime }}$ and no correlation between coordinate measurements at the end points (see Eq. 2). When displacements are retrieved from a pair of SAR images, one needs to consider position and tracking uncertainties, i.e., ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$, respectively. The latter arises from the fact that in a satellite image details of structures on a pixel scale may be difficult to match between images 1 and 2. In this case the uncertainty in displacement (which here is the distance between positions of a fixed point on an ice structure in images 1 and 2) is ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}=$ 2${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$. For buoy arrays, ${\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$ is zero, since a buoy remains fixed relative to the ice floe on which it was deployed.

In a SAR image, the geolocation (position) error is caused by the inaccuracies of the parameters describing the satellite orbit as a function of space and time. In general, the error caused by these inaccuracies is uniform across the image with only small local variations. Hence the assumption of independent geolocation errors is not valid if distances between moving objects are small. Holt et al. (1992) give a correlation length of 10 km for the uncertainty of the geolocation error, σ_coord, but correlation lengths of up to 100 km may be possible (Ronald Kwok, personal communication, 2020). Deformation parameters from SAR image pairs are usually calculated over regions that are on the order of 10 km or less across. With correlation lengths of ≥10 km, geolocation errors at all pixels in the region are almost equal, which means that geolocation error variances σ_coord are small (as is discussed in Sect. 3.4.1). It is hence reasonable in many cases to regard the geolocation errors in image 1 and image 2 as constant biases and to assume that σ_coord=0 (Sect. 3.4.2). When calculating the distance between two points with identical geolocation errors, we obtain hence ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$. Differences between the biases in image 1 and 2 affect the retrieval of ice drift. Deformation, on the other hand, is related to the relative change in size and shape of a given area between acquisitions of image 1 and image 2. The relative area change is independent of the regionally constant difference between the biases and depends only on the error variances. Therefore, deformation can be estimated with sufficient accuracy even if geolocation errors are large.

2.2 Uncertainty of drift velocity

The deformation is calculated from components of the velocity gradients according to Eq. (5). Hence, we have to consider the uncertainty in the measurements of velocity components u_i and v_i. The components are calculated from $u=d_x/\mathrm{\Delta }T$ and $v=d_y/\mathrm{\Delta }T$, where $d_x=(x^{\prime }-x)$ and $d_y=(y^{\prime }-y)$ are the displacements in the x and y direction, respectively, and ΔT is the time interval needed for the position change from (x, y) to (x^′, y^′). Considering that errors in measuring time and positions are not correlated, we obtain from Eq. (2), taking into account a possible tracking error:

$$\begin{array}{}\text{(6a)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_u^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{T}^{\mathrm{2}}}}{\mathit{\sigma }}_{{\mathrm{d}}_x}^{\mathrm{2}}+{\left({\displaystyle \frac{-d_x}{\mathrm{\Delta }{T}^{\mathrm{2}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{T}^{\mathrm{2}}}}\left(\mathrm{2}{\mathit{\sigma }}_x^{\mathrm{2}}+{\mathit{\sigma }}_{{\mathrm{tr}}_x}^{\mathrm{2}}+u^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}\right)\end{array},\text{(6b)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_v^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{T}^{\mathrm{2}}}}{\mathit{\sigma }}_{{\mathrm{d}}_y}^{\mathrm{2}}+{\left({\displaystyle \frac{-d_y}{\mathrm{\Delta }{T}^{\mathrm{2}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{T}^{\mathrm{2}}}}\left(\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}+{\mathit{\sigma }}_{{\mathrm{tr}}_y}^{\mathrm{2}}+v^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}\right)\end{array},\end{array}$$

where ${\mathit{\sigma }}_{{\mathrm{d}}_x}$ and ${\mathit{\sigma }}_{{\mathrm{d}}_y}$ are the uncertainties of the displacements (distances) in the x and y direction and ${\mathit{\sigma }}_{{\mathrm{tr}}_x}$ and ${\mathit{\sigma }}_{{\mathrm{tr}}_y}$ are the corresponding components of the tracking error. If the uncertainty in timing, ${\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}$, is not zero, the assumption that ${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}$ is only valid if u²=v². The uncertainty in speed U (i.e., the magnitude of velocity vector U) can be computed using Eq. (6), replacing ${\mathit{\sigma }}_u^{\mathrm{2}}$ with ${\mathit{\sigma }}_U^{\mathrm{2}}$, ${\mathit{\sigma }}_{{\mathrm{d}}_x}^{\mathrm{2}}$ with ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$,${\mathit{\sigma }}_{{\mathrm{tr}}_x}^{\mathrm{2}}$ with ${\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$, and u with U, considering that $U=d/\mathrm{\Delta }T=\sqrt[\mathrm{2}]{u^{\mathrm{2}}+v^{\mathrm{2}}}$ and $d=\sqrt[\mathrm{2}]{{\left(x^{\prime }-x\right)}^{\mathrm{2}}+{\left(y^{\prime }-y\right)}^{\mathrm{2}}}$. When calculating the relative error variance σ_U∕U, one obtains Eq. (A1) in Hutchings et al. (2012).

If, on the other hand, both components of the vector U are determined separately (hence considering magnitude and direction), the result is different:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_U^{\mathrm{2}}={\left({\displaystyle \frac{\partial U}{\partial u}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_u^{\mathrm{2}}+{\left({\displaystyle \frac{\partial U}{\partial v}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_v^{\mathrm{2}}={\left({\displaystyle \frac{u}{U}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_u^{\mathrm{2}}+{\left({\displaystyle \frac{v}{U}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_v^{\mathrm{2}}.\end{array}$$

If one substitutes Eq. (6a) and (6b) into Eq. (7) and sets ${\mathit{\sigma }}_{\mathrm{d}x}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{d}y}^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$, one obtains

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_U^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}}}\left({\displaystyle \frac{u^{\mathrm{4}}+v^{\mathrm{4}}}{u^{\mathrm{2}}+v^{\mathrm{2}}}}\right).\end{array}$$

If σ_ΔT cannot be neglected, and if u=0 and v=U or v=0 and u=U, the second term of Eq. (8) yields $U^{\mathrm{2}}({\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}$/ΔT²), which is the uncertainty in speed given above. If, on the other hand, u=v and hence U²=2u², the second term is $\mathrm{0.5}{U}^{\mathrm{2}}({\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}/\mathrm{\Delta }{T}^{\mathrm{2}})$. This result may be viewed as if independent measurements of the two components u and v reduce the uncertainty contribution of ${\mathit{\sigma }}_{{\mathrm{\Delta }}_T}^{\mathrm{2}}$.

2.3 Uncertainty of polygon area

The uncertainty of an area measurement is needed for application of Eq. (5) and equations presented in the following sections. The starting point for calculating the variance of error for the measurement of an area is the surveyor's area formula valid for a polygon with an outline consisting of n segments in a plane spanned by the x and y axis:

$$\begin{array}{}\text{(9)}& A={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^n(x_i{y}_{i+\mathrm{1}}-x_{i+\mathrm{1}}{y}_i).\end{array}$$

Here $x_{n+\mathrm{1}}\equiv x_{\mathrm{1}}$ and $y_{n+\mathrm{1}}\equiv y_{\mathrm{1}}$ (closed polygon), i is the summation index, and the boundary is traced in a counterclockwise sense. We have to consider that each coordinate appears twice in the sum of Eq. (9). When i=k we have, e.g., for x: x_ky_k+1, and when $i=k-\mathrm{1}$ we have $-x_k{y}_{k-\mathrm{1}}$. For the law of error propagation, we need the derivatives:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\partial A}{\partial x_k}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(y_{k+\mathrm{1}}-y_{k-\mathrm{1}}\right)\text{ and }\\ & {\displaystyle \frac{\partial A}{\partial y_k}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(x_{k+\mathrm{1}}-x_{k-\mathrm{1}}\right)\end{array},\end{array}$$

where k is the index of the derivative. Hence, we obtain

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_A^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _{i=\mathrm{1}}^n\left[{\mathit{\sigma }}_{i\mathit{\_}x}^{\mathrm{2}}{\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{i\mathit{\_}y}^{\mathrm{2}}{\left(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right].\end{array}$$

We can assume that coordinate uncertainties ${\mathit{\sigma }}_{i\mathit{\_}x}^{\mathrm{2}}={\mathit{\sigma }}_{i\mathit{\_}y}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ are equal and the same for all measured positions. The uncertainty of the area is then

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_A^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{4}}}\sum _{i=\mathrm{1}}^n\left[{\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}+{\left(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right].\end{array}$$

Examples of applying Eq. (12) on basic polygons are shown in Fig. 2. Arbitrarily shaped triangles and quadrangles, which are basic patterns for arrays of three or four buoys and for grid cells in satellite images when applying the Lagrangian approach, are shown at the bottom. The x–y coordinate system is here oriented such that the calculation of the uncertainty is easy. For any orientation of the triangle or quadrangle, side lengths and distances can be derived from the coordinates (x,y) of the edge points. For squares and equal-sided right-angled triangles, which are typical grid cells when retrieving ice drift from satellite images in a Eulerian approach, the uncertainty is directly proportional to the area. If a square grid cell is split into two triangles (as in Fig. 1), the uncertainty in area of each triangle is half that of the square.

Figure 2 Application of Eq. (12) to different geometrical figures: rectangle, equal-sided right triangle, rhombus, regular hexagon, triangle, and quadrangle.

For an assessment on how the polygon shape affects the magnitude of uncertainty we require that the enclosed area remains constant. The areas of a square with side length L and a right-angled triangle with two sides of length L_T are equal if $L_T=\sqrt{\mathrm{2}}L$. In this case we get ${\mathit{\sigma }}_A^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}{L}^{\mathrm{2}}$ for both square and triangle, which means that in this particular case the increase in number of vertices does not result in a decrease in σ_A. For a hexagon with A=L², on the other hand, one obtains $s^{\mathrm{2}}=\mathrm{2}{L}^{\mathrm{2}}/\mathrm{3}\sqrt{\mathrm{3}}$ and ${\mathit{\sigma }}_A^{\mathrm{2}}=\mathrm{1.44}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}{L}^{\mathrm{2}}$ (where s is the length of a line segment on the boundary of the hexagon; see Fig. 2). The issue of adding more vertices while keeping the shape of the polygon is addressed in Sect. 3.6.

The question arises of how large the smallest detectable area change is in a SAR image. To address this question, we assume a square grid cell with its vertices on the positions of adjacent displacement vectors and its sides parallel to the x and y axes of a Cartesian coordinate system. The cell covers m×m square-shaped pixels of side length Δx. The minimum possible change is to move one edge point by the side length of 1 pixel, either in the x or y direction. This adds the area of a right triangle with legs Δx and mΔx(Δy=Δx), and the change in the area is $\mathrm{\Delta }A=\left(m\mathrm{\Delta }{x}^{\mathrm{2}}\right)/\mathrm{2}$, i.e., 100∕(2m) percent of the original area (mΔx)². Hence the larger the number of pixels in the area, the smaller the detectable relative area change. However, until now we assumed that the position error is zero, but we have to consider the uncertainty of the area estimate, which is ${\mathit{\sigma }}_A^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}{m}^{\mathrm{2}}\mathrm{\Delta }{x}^{\mathrm{2}}$ for a square with L=mΔx. To be sure that a detected area change is real, ΔA needs to be larger than σ_A or ${\mathit{\sigma }}_{\mathrm{coord}}<\frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathrm{2}}}\mathrm{\Delta }x$.

2.4 Uncertainties for divergence, shear, vorticity, and total deformation in fixed grids

We consider a grid with displacement or drift velocity vectors on the vertices. For calculating the deformation parameters, we need the velocity gradients $u_x,u_y,v_x,v_y$, obtained from Eq. (5). Formally, the gradients depend on the area A, positions (x_i, y_i), and velocities (u_i, v_i); see Sect. 2.5. Here we assume that the georeferencing of the satellite images is accurate. In this case, the positions (x_i, y_i) of vertices and the area of each grid cell are known precisely, which means that σ_coord=0 and σ_A=0. The displacement or velocity vectors, however, have an uncertainty related to the tracking error. With $\partial u_x/\partial u_k=(y_{k+\mathrm{1}}-y_{k-\mathrm{1}})/\mathrm{2}A$ and again considering that two terms in the sum Eq. (5) include u_i, the uncertainty of the velocity gradient in the x direction is (Griebel and Dierking, 2018)

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_{ux}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_u^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\sum _{i=\mathrm{1}}^n{\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}},\end{array}$$

with analogous equations for the other gradient components. The divergence is ${\dot{\mathit{\epsilon }}}_{\mathrm{div}}=u_x+v_y$, Eq. (3a), and the corresponding uncertainty is ${\mathit{\sigma }}_{\mathrm{div}}=\sqrt[\mathrm{2}]{{\mathit{\sigma }}_{ux}^{\mathrm{2}}+{\mathit{\sigma }}_{vy}^{\mathrm{2}}}$, if u_x and v_y are independent. Throughout this section we assume that ${\mathit{\sigma }}_U={\mathit{\sigma }}_u={\mathit{\sigma }}_v$ and σ_ΔT= 0; hence the error variance for the divergence is

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_U^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\sum _{i=\mathrm{1}}^n\left[{\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}+{\left(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right]\\ & ={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}\mathrm{\Delta }{T}^{\mathrm{2}}}}\sum _{i=\mathrm{1}}^n\left[{\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}+{\left(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right]\end{array}.\end{array}$$

Equation (14) resembles the uncertainty for a polygon, Eq. (11). Since the position uncertainty σ_coord is set to zero, the uncertainty of velocity U is only a function of the tracking uncertainty σ_tr; see Eq. (8) (assuming σ_ΔT=0). For the vorticity Eq. (3b) one obtains ${\mathit{\sigma }}_{\mathrm{vrt}}=\sqrt[\mathrm{2}]{{\mathit{\sigma }}_{vx}^{\mathrm{2}}+{\mathit{\sigma }}_{uy}^{\mathrm{2}}}$ and thus the same expression as for the divergence. The shear rate is given by Eq. (3c). Calculating the derivatives with respect to the velocity gradient components and applying the law of error propagation yields

$$\begin{array}{}\text{(15a)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{shr}}^{\mathrm{2}}=& {\displaystyle \frac{{\left(u_x-v_y\right)}^{\mathrm{2}}}{{\dot{\mathit{\epsilon }}}_{\mathrm{shr}}^{\mathrm{2}}}}\left({\mathit{\sigma }}_{ux}^{\mathrm{2}}+{\mathit{\sigma }}_{vy}^{\mathrm{2}}\right)\\ & +{\displaystyle \frac{{\left(u_y+v_x\right)}^{\mathrm{2}}}{{\dot{\mathit{\epsilon }}}_{\mathrm{shr}}^{\mathrm{2}}}}\left({\mathit{\sigma }}_{uy}^{\mathrm{2}}+{\mathit{\sigma }}_{vx}^{\mathrm{2}}\right).\end{array}\end{array}$$

With $\mathit{\phi }=\mathrm{1}/\mathrm{2}\arctan \left(\right(u_y+v_x)/(u_x-v_y\left)\right)$, which gives the principal direction of shear, and using Eq. (3c) and relations ${\cos }^{\mathrm{2}}\left(\arctan \left(x\right)\right)=\mathrm{1}/\left(\mathrm{1}+x^{\mathrm{2}}\right)$ and ${\sin }^{\mathrm{2}}\left(\arctan \left(x\right)\right)=x^{\mathrm{2}}/\left(\mathrm{1}+x^{\mathrm{2}}\right)$, Eq. (15a) can be expressed as

$$\begin{array}{}\text{(15b)}& {\mathit{\sigma }}_{\mathrm{shr}}^{\mathrm{2}}={\cos }^{\mathrm{2}}\left(\mathrm{2}\mathit{\varphi }\right){\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}+{\sin }^{\mathrm{2}}\left(\mathrm{2}\mathit{\varphi }\right){\mathit{\sigma }}_{\mathrm{vrt}}^{\mathrm{2}}.\end{array}$$

Since ${\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{vrt}}^{\mathrm{2}}$ and cos²(2φ)+ sin²(2φ)=1, the error variances are equal for divergence, vorticity, and shear. For the total deformation, Eq. (3d), we need the derivatives $\partial \left({\dot{\mathit{\epsilon }}}_{\mathrm{tot}}\right)/\partial \left({\dot{\mathit{\epsilon }}}_{\mathrm{shr}}\right)$ and $\partial \left({\dot{\mathit{\epsilon }}}_{\mathrm{tot}}\right)/\partial \left({\dot{\mathit{\epsilon }}}_{\mathrm{div}}\right)$, with which we obtain

$$\begin{array}{}\text{(16a)}& {\mathit{\sigma }}_{\mathrm{tot}}^{\mathrm{2}}={\displaystyle \frac{{\dot{\mathit{\epsilon }}}_{\mathrm{shr}}^{\mathrm{2}}}{{\dot{\mathit{\epsilon }}}_{\mathrm{tot}}^{\mathrm{2}}}}{\mathit{\sigma }}_{\mathrm{shr}}^{\mathrm{2}}+{\displaystyle \frac{{\dot{\mathit{\epsilon }}}_{\mathrm{div}}^{\mathrm{2}}}{{\dot{\mathit{\epsilon }}}_{\mathrm{tot}}^{\mathrm{2}}}}{\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}.\end{array}$$

If we define θ= arctan (ε_shr∕ε_div) (Stern et al., 1995), Eq. (16a) can be rewritten as

$$\begin{array}{}\text{(16b)}& {\mathit{\sigma }}_{\mathrm{tot}}^{\mathrm{2}}={\sin }^{\mathrm{2}}\left(\mathit{\theta }\right){\mathit{\sigma }}_{\mathrm{shr}}^{\mathrm{2}}+{\cos }^{\mathrm{2}}\left(\mathit{\theta }\right){\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}.\end{array}$$

The angle θ gives the relative contributions of divergence and shear: pure divergence is θ=0^∘, uniaxial extension is θ=45^∘, pure shear is θ=90^∘, uniaxial compression is θ=135^∘, and pure convergence is θ=180^∘. Since the uncertainties for shear and divergence are of equal magnitude, it follows that

$$\begin{array}{}\text{(16c)}& {\mathit{\sigma }}_{\mathrm{tot}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{shr}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{vrt}}^{\mathrm{2}}.\end{array}$$

In the following, we assume that σ_ΔT can be neglected and that the standard deviations for the velocity components u and v are equal. Using Eq. (14) for a square cell, we obtain the following for the uncertainty of the divergence:

$$\begin{array}{}\text{(17)}& {\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_U^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\left(\mathrm{4}{L}^{\mathrm{2}}+\mathrm{4}{L}^{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}}{L^{\mathrm{2}}\mathrm{\Delta }{T}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{L^{\mathrm{2}}\mathrm{\Delta }{T}^{\mathrm{2}}}},\end{array}$$

with A=L², ${\mathit{\sigma }}_U^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}/\mathrm{\Delta }{T}^{\mathrm{2}}$, and $\mathrm{\Delta }T=t-t_{\mathrm{0}}$ as above. Since the position uncertainty is zero in the case investigated here, ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ (which equals $\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}$; see Sect. 2.1) depends only on the tracking error (compare to Eq. 17 in Lindsay and Stern, 2003).

2.5 Uncertainties of deformation parameters, general case

For an array of buoys, we have to consider errors of the area, the buoy velocity components u and v, and the coordinates (x, y) of each buoy position. The general case does also apply to SAR images if geolocation error variances cannot be neglected. A buoy array consists of single buoys arbitrarily positioned over a plane. When connecting all buoy positions with lines, a polygon of area A is formed in which distances between adjacent buoys are usually different. The starting point is Eq. (5). In the following equations summation bounds from i=1 to n are omitted. We note that the equations in this section have been independently derived by Bouchat and Tremblay (2020) as well.

For the uncertainty in u_x we obtain

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\mathit{\sigma }}_{ux}^{\mathrm{2}}=& {\mathit{\sigma }}_A^{\mathrm{2}}{\left({\displaystyle \frac{\partial u_x}{\partial A}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_u^{\mathrm{2}}\sum {\left({\displaystyle \frac{\partial u_x}{\partial u_i}}\right)}^{\mathrm{2}}\\ & +{\mathit{\sigma }}_y^{\mathrm{2}}\sum {\left({\displaystyle \frac{\partial u_x}{\partial y_i}}\right)}^{\mathrm{2}}.\end{array}\end{array}$$

With

$$\begin{array}{rl}& {\displaystyle \frac{\partial u_x}{\partial A}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{A}^{\mathrm{2}}}}\sum \left(u_{i+\mathrm{1}}+u_i\right)\left(y_{i+\mathrm{1}}-y_i\right),\\ & {\displaystyle \frac{\partial u_x}{\partial u_k}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}A}}\left(y_{k+\mathrm{1}}-y_{k-\mathrm{1}}\right),\text{and}{\displaystyle \frac{\partial u_x}{\partial y_k}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}A}}\left(u_{k+\mathrm{1}}-u_{k-\mathrm{1}}\right),\end{array}$$

Eq. (19) reads

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\mathit{\sigma }}_{ux}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_A^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{4}}}}{\left[\sum \left(u_{i+\mathrm{1}}+u_i\right)\left(y_{i+\mathrm{1}}-y_i\right)\right]}^{\mathrm{2}}\\ & +{\displaystyle \frac{{\mathit{\sigma }}_u^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\sum {\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}\\ & +{\displaystyle \frac{{\mathit{\sigma }}_y^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\sum {\left(u_{i+\mathrm{1}}-u_{i-\mathrm{1}}\right)}^{\mathrm{2}}.\end{array}\end{array}$$

The first term on the right side is calculated from line segments connecting adjacent vertices ($i+\mathrm{1},j+\mathrm{1}$) and (i,j), and the second and third terms are calculated from chords between (i+1, j+1) and (i−1, j−1). Assuming ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}={\mathit{\sigma }}_x^{\mathrm{2}}={\mathit{\sigma }}_y^{\mathrm{2}}$; ${\mathit{\sigma }}_U^{\mathrm{2}}={\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}$ (the latter follows from ${\mathit{\sigma }}_T^{\mathrm{2}}/\mathrm{\Delta }{T}^{\mathrm{2}}\approx \mathrm{0}$) one obtains for the divergence:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}={\mathit{\sigma }}_{ux}^{\mathrm{2}}+{\mathit{\sigma }}_{vy}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_A^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{4}}}}\left\{{\left[\sum \left(u_{i+\mathrm{1}}+u_i\right)\left(y_{i+\mathrm{1}}-y_i\right)\right]}^{\mathrm{2}}\right.\\ & \left.+{\left[\sum \left(v_{i+\mathrm{1}}+v_i\right)\left(x_{i+\mathrm{1}}-x_i\right)\right]}^{\mathrm{2}}\right\}\\ & +{\displaystyle \frac{{\mathit{\sigma }}_U^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\left[\sum {\left(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}\right)}^{\mathrm{2}}+\sum {\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right]\\ & +{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\left[\sum {\left(u_{i+\mathrm{1}}-u_{i-\mathrm{1}}\right)}^{\mathrm{2}}+\sum {\left(v_{i+\mathrm{1}}-v_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right],\end{array}\end{array}$$

where the first term can be written as $\frac{{\mathit{\sigma }}_A^{\mathrm{2}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)}{A^{\mathrm{2}}}$, considering Eq. (5). For the vorticity, only the first term is different:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{vrt}}^{\mathrm{2}}={\mathit{\sigma }}_{uy}^{\mathrm{2}}+{\mathit{\sigma }}_{vx}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_A^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{4}}}}\left\{{\left[\sum \left(u_{i+\mathrm{1}}+u_i\right)\left(x_{i+\mathrm{1}}-x_i\right)\right]}^{\mathrm{2}}\right.\\ & \left.+{\left[\sum \left(v_{i+\mathrm{1}}+v_i\right)\left(y_{i+\mathrm{1}}-y_i\right)\right]}^{\mathrm{2}}\right\}\\ & +{\displaystyle \frac{{\mathit{\sigma }}_U^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\left[\sum {\left(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}\right)}^{\mathrm{2}}+\sum {\left(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right]\\ & +{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{4}{A}^{\mathrm{2}}}}\left[\sum {\left(u_{i+\mathrm{1}}-u_{i-\mathrm{1}}\right)}^{\mathrm{2}}+\sum {\left(v_{i+\mathrm{1}}-v_{i-\mathrm{1}}\right)}^{\mathrm{2}}\right].\end{array}\end{array}$$

Here, the first term can be written as $\frac{{\mathit{\sigma }}_A^{\mathrm{2}}\left({\mathit{u}}_y^{\mathrm{2}}+{\mathit{v}}_x^{\mathrm{2}}\right)}{A^{\mathrm{2}}}$. The first terms in Eqs. (20) and (21), right side, consider that the relative error variance of the area affects the magnitude of the average velocity gradients. The second term is the variance of divergence/vorticity of the velocity field in a fixed grid where positions of vertices are known precisely, Eq. (14). The last term takes into account the effect of uncertainties in the positions of buoys in the field of velocity vectors. The velocity is usually determined from buoy positions separated by a time interval $\mathrm{\Delta }T=T_{i+\mathrm{1}}-T_i$. However, within ΔT also the buoy array changes its area and shape. Hence an alternative approach would be to determine the average velocity from positions at T_i−1, T_i, and T_i+1 and link it with the geometric properties of the buoy array at time T_i. For the shear and total deformation, the results are formally equal to Eqs. (15a) and (15b) as well as (16a) and (16b), where now σ_ux, σ_uy, σ_vx, and σ_vy are calculated using Eq. (20) and analogous expressions. Note that in this case the uncertainties of divergence, vorticity, shear, and total deformation differ from one another, unless ${\mathit{\sigma }}_{ux}^{\mathrm{2}}={\mathit{\sigma }}_{uy}^{\mathrm{2}}={\mathit{\sigma }}_{vx}^{\mathrm{2}}={\mathit{\sigma }}_{vy}^{\mathrm{2}}$. In practical applications, they can be evaluated numerically. This requires the knowledge of uncertainties σ_coord for buoys and σ_coord, σ_tr for satellite images.

3 Discussion

Equations (20) and (21) together with Eqs. (15) and (16) provided above indicate that statistical uncertainties are not only influenced by geolocation and tracking errors but also depend on the shape and size of grid cells and buoy arrays. In the following discussion we consider magnitudes of geolocation and tracking errors reported in the literature and selected squares and triangles as examples for grid cells in SAR images (Lindsay, 2002; Bouillon and Rampal, 2015) and for splitting large buoy arrays into smaller units (Hutchings et al., 2012; Itkin et al., 2017). The effect of combining several cells is investigated. Finally, we focus on the range of validity of the equations derived in Sect. 2 and alternative methods of analysis.

3.1 Typical magnitudes of deformation parameters

The statistical uncertainties have to be related to the typical magnitudes of the deformation parameters. According to Leppäranta (2011, p. 70) the total deformation of drifting ice typically varies between around 90 % d⁻¹ in the marginal ice zone to 0.9 % d⁻¹ in the central Arctic. For the vorticity, magnitudes up to 9 % d${}^{-\mathrm{1}}=\mathrm{1}/\mathrm{2}$ (0.09) revolutions d${}^{-\mathrm{1}}=\mathrm{16.2}$^∘ d⁻¹ were observed. Hutchings et al. (2012, Figs. 4 and 7) analyzed displacements of an array of 24 buoys deployed in the Weddell Sea on first- and second-year ice with concentrations above 90 %. For divergence, they found most values between −90 % d⁻¹ and 90 % d⁻¹ at a spatial scale of 10 km, as well as at a 60 km scale mainly between ±25 % d⁻¹ and up to 35 % d⁻¹ for the shear. Note that spatial scales are mentioned here since they affect the observed magnitudes of deformation (e.g., Marsan et al., 2004). Itkin et al. (2017) observed exceptional events of strong divergence and shear of up to 200 % d⁻¹ from buoys in an area north of Svalbard (their Fig. 4), but over most of the measurement period, magnitudes were lower. At scales of 15 km or less, values for divergence covered the range ±20 % d⁻¹ over several days to weeks, but also variations of about ±100 % d⁻¹ occurred for 3 weeks. Shear was close to zero for a few days but varied mainly from 20 % d⁻¹ to 70 % d⁻¹ for 3 weeks. At measurement scales larger than 60 km, the magnitudes of divergence and shear were lower than at ≤15 km scale, with the exception of very short periods during which the opposite was the case. Magnitudes for divergence were roughly at ±10 % d⁻¹ with occasional minima and maxima in the range of ±100 % d⁻¹, and for shear most values were ≤10 % d⁻¹ with a few peaks at about 100 % d⁻¹. Based on merged velocity measurements from buoys and different satellite sensors, Lindsay (2002) provided a table for monthly averaged values of divergence (−0.6 % d⁻¹ to 0.5 % d⁻¹), shear (0.9 % d⁻¹ to 4 % d⁻¹), and vorticity (−2.3 % d⁻¹ to 3.2 % d⁻¹) from the Beaufort Sea at a scale of 100 km. Stern and Moritz (2002, Fig. 4) used SAR images and found decreasing magnitudes for the divergence for increasing spatial scales from 50 km×50 km to 200 km×200 km in the Beaufort Sea. Magnitudes were largest between August and February with minima/maxima between −5 % d⁻¹ and 5 % d⁻¹ at a scale of 50 km, decreasing at larger scales. Note that the uncertainties resulting from the equations given in the subsections below have to be multiplied by 100 to obtain a value in percent per time unit.

3.2 Uncertainties for areas of simple geometric shape

In general, the uncertainty of the deformation parameters depends on the ratio ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/A^{\mathrm{2}}$ (since σ_A and σ_U are functions of σ_coord); hence for given geolocation and tracking errors it decreases with increasing area. The first term in Eqs. (20) and (21) is smallest if, for given area and velocity gradients, σ_A is at a minimum. For an arbitrary triangle with sides a, b, c, the uncertainty ${\mathit{\sigma }}_A^{\mathrm{2}}$ is $\mathrm{0.25}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}(a^{\mathrm{2}}+b^{\mathrm{2}}+c^{\mathrm{2}})$ (see Fig. 2). Of all triangles with the same base and the same area A, the equal-sided triangle with $a=b=c$ has the smallest perimeter and hence the lowest uncertainty, which is ${\mathit{\sigma }}_A^{\mathrm{2}}=\sqrt{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}=\mathrm{1.73}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ for a unit area. (This follows from the equations for the area of the equal-sided triangle which is $A=\frac{\sqrt{\mathrm{3}}}{\mathrm{4}}{a}^{\mathrm{2}}$ and for the uncertainty ${\mathit{\sigma }}_A^{\mathrm{2}}=\frac{\mathrm{3}{a}^{\mathrm{2}}}{\mathrm{4}}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ if A=1). In case of rectangles and rhombi, squares have the smallest perimeter (see Fig. 2). In both cases the uncertainty is ${\mathit{\sigma }}_A^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ for a unit area, and hence larger than for the equal-sided triangle. For the regular hexagon, which is composed of six equal-sided triangles, one obtains ${\mathit{\sigma }}_A^{\mathrm{2}}=\frac{\mathrm{5}}{\mathrm{2}\sqrt{\mathrm{3}}}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}=\mathrm{1.44}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ (Fig. 2). So the progression of ${\mathit{\sigma }}_A^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ from triangles to squares to hexagons goes from 1.73 A to 2.00 A to 1.44 A.

3.3 Uncertainties in time

The accuracy of time readings for the acquisitions of satellite images is on the order of subseconds. The product of sea ice drift velocity and uncertainty of time reading appears on the right-hand side of Eq. (6): $\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}+u^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}$. Average sea ice drift velocities range mostly from 0 to 0.35 m s⁻¹ (Rampal et al., 2009). Kræmer et al. (2015) determined instantaneous line-of-sight ice drift velocities, using Doppler frequency measurements from SAR, and found values as large as 0.4–0.6 m s⁻¹. If we assume a maximum value of u=1 m s⁻¹ and a maximum uncertainty of time readings of 1 ms, the term $u^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}$ on the right side of Eq. (6) is 10⁻⁶ m² at the most. It can be neglected compared to the typical values of terms ${\mathit{\sigma }}_x^{\mathrm{2}}$, ${\mathit{\sigma }}_y^{\mathrm{2}}$ and ${\mathit{\sigma }}_{{\mathrm{tr}}_x}^{\mathrm{2}}$, ${\mathit{\sigma }}_{{\mathrm{tr}}_y}^{\mathrm{2}}$ in Eq. (6) (see Sect. 3.4 for a discussion of the effect of position and tracking errors). The uncertainty σ_ΔT of the GPS time (used both for buoys and satellites such as Sentinel-1) is given as better than 1 ms (see, e.g., https://www.atomic-clock.galleon.eu.com, last access: 4 September 2020, and https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-1-sar/pod/egp, last access: 4 September 2020). Similar considerations apply to Eq. (8). Hence, in Eqs. (20) and (21) we have ${\mathit{\sigma }}_U^{\mathrm{2}}=(\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}})/\mathrm{\Delta }{T}^{\mathrm{2}}$ both for velocity retrievals from satellite image pairs and buoy arrays. For given position and tracking errors, the second term in Eqs. (20) and (21) decreases with increasing time interval ΔT and area A. The third term involving the coordinate uncertainty σ_coord also decreases with increasing area A.

Another issue that has to be considered is the time synchronization between individual buoys in an array. Differences of a few seconds may be possible in practice. In the following discussion we assume that position data of all buoys are exactly synchronized but also discuss an example for which this was not the case in Sect. 3.5.

3.4 Deformation retrievals from square grid cells

Here we first focus on the retrieval of deformation parameters calculated from square grid cells in SAR images or from square-shaped buoy arrays. For SAR images, we consider the case in which geolocation errors may have slight variations, hence σ_coord≠0. If a square of side length L, with sides parallel to the x and y axes, is positioned in a spatially varying velocity field as shown in Fig. 3, the uncertainty of the divergence is

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{3}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{L^{\mathrm{2}}}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{L^{\mathrm{2}}}}\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right)\\ & +{\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{L}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{L}^{\mathrm{2}}}}.\end{array}\end{array}$$

This follows from Eq. (20) with the velocities given in Fig. 3 at the edges 1–4 of the square. The uncertainty of the vorticity is from Eq. (21):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{vrt}}^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{3}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{L^{\mathrm{2}}}}\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right)+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{L^{\mathrm{2}}}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)\\ & +{\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{L}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{L}^{\mathrm{2}}}}.\end{array}\end{array}$$

Uncertainties of shear and total deformation can be calculated using Eqs. (15b) and (16b) as weighted averages of the error variances of (i) divergence and vorticity and of (ii) shear and divergence, respectively. The second term in Eq. (22) and first term in Eq. (23) indicates that the uncertainties of divergence and vorticity are affected by contributions from pure shear. The third and fourth term of Eqs. (22) and (23) are independent of the velocity gradients and are only a function of position and tracking error, time interval between position measurements, and size of the square. The fourth term is equal to Eq. (17) in Lindsay and Stern (2003). In general, it is more realistic to assume that arrays of four buoys are arbitrarily shaped quadrangles. As mentioned in Sect. 1, drift vectors from SAR image pairs are irregularly spaced if calculated using feature tracking (e.g., Komarov and Barber, 2014; Muckenhuber et al., 2016; Demchev et al., 2017). While σ_tr, σ_coord, and ΔT are constant, σ_A and A depend on the size and shape of the quadrangle that changes from grid cell to grid cell (Figs. 1c and 2). In this case the most convenient approach for calculating deformation parameters is the application of Eqs. (20) and (21) together with Eqs. (15) and (16). We emphasize, however, that the heterogeneous spatial distribution of drift vectors is regarded as a disadvantage for evaluating and analyzing sea ice deformation, since the latter is a scale-dependent process (Korosov and Rampal, 2017).

Figure 3 Uncertainty of divergence and vorticity for a square in a spatially varying velocity field with gradients u_x, u_y, v_x, v_y.

3.4.1 Geolocation error and uncertainties in SAR images

When ice drift is retrieved from images of modern SAR systems, the contribution of those terms that depend on σ_coord∕L can usually be neglected, as we will show below. For Envisat ASAR stripmap and wide-swath mode images (IM and WSM), e.g., Small et al. (2005) reported differences between measured positions of reflectors and their positions in the SAR image of 1.63±0.82 m in azimuth (considering bistatic correction) and 2.02±0.51 m in slant range for normal imaging mode in single-look complex format. Ground-range products require the transformation from slant to ground range as an additional step. When judging the effect of position errors on the uncertainty of divergence and vorticity, the systematic bias (mean error) of positions affects all vertices of a grid cell in the same way; hence only the standard deviation σ has to be considered as geolocation uncertainty. Considering the σ values of position errors given above, we use a value of 1 m as a conservative estimate of the azimuth and ground-range position uncertainty for IM. For ground-range WSM images, the accuracy of positioning was better than 1 pixel. If we assume that the ratio σ[m]∕σ[pixel] is approximately the same for IM and WSM, the uncertainty for the latter is about 7 m at maximum. In the study of Hollands and Dierking (2011), e.g., resolution pyramids and cascades are used for retrieving sea ice displacements from Envisat ASAR IM and WSM data. For the level of highest spatial resolution, the side length of the grid cells (distance between adjacent displacement vectors) was 300 m for IM and 1200 m for WSM. Hence, the corresponding ratios ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/L^{\mathrm{2}}$ are on the order of ${\mathrm{1}}^{\mathrm{2}}/{\mathrm{300}}^{\mathrm{2}}\approx {\mathrm{10}}^{-\mathrm{5}}$ and ${\mathrm{7}}^{\mathrm{2}}/{\mathrm{1200}}^{\mathrm{2}}\approx \mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{5}}$, respectively. For modern SAR systems such as TerraSAR-X and Sentinel-1, the positioning accuracy is even better than for Envisat (e.g., Schubert et al., 2008; Schubert et al., 2017). The geolocation error of older SAR systems, however, is larger. In their analysis of drift and deformation products from the RADARSAT Geophysical Processor System (RGPS), Lindsay and Stern (2003) report geolocation errors of 225 and 277 m for RADARSAT ScanSAR images (to be treated as bias; see above). For the combined geolocation and tracking uncertainty ${\mathit{\epsilon }}_{\mathrm{RGPS}}=\sqrt{\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}$, they found a value of 286 m. With a tracking uncertainty of 100 m, the geolocation uncertainty is hence 190 m. The initial grid cells used for the RGPS are squares of 10 km side length, but they change their shape in successive time steps since the RGPS drift and deformation products are based on the Lagrangian approach. The ratio ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/L^{\mathrm{2}}$ is approximately 200²∕10 000${}^{\mathrm{2}}=\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{4}}$. The third and fourth term in Eqs. (22) and (23) can be directly computed from position and tracking error, the time interval ΔT between image acquisitions, and the grid cell size. The ratio between the fourth and the third term is ${\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}/\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$. In the following section, the relative contributions of single terms in Eqs. (22) and (23) are discussed.

3.4.2 Examples: uncertainties versus true magnitudes of deformation

According to Sect. 3.1, a value of ±1 d⁻¹ can be regarded as a large divergence rate which is rarely exceeded in reality. Large values of shear were at about 0.7 d⁻¹. Considering the numbers for divergence and shear given in Sect. 3.1 we can deduce that the terms $\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)$ and $\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right)$ in Eqs. (22) and (23) are <1 d⁻² in most cases and at larger length scales and weak deformation more likely on the order of 10⁻¹ or 10⁻² d⁻². This means that ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/L^{\mathrm{2}}$ and $\mathrm{3}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/L^{\mathrm{2}}$ can be used as upper bounds for the first and second term in Eqs. (22) and (23) (see Table 1).

Table 1 Magnitudes of terms 1 to 4 in Eqs. (23) and (24).

Hollands and Dierking (2011) found tracking errors between 0.8 and 1.6 pixels (their Tables 3 and 4, standard deviations), which corresponded to 20–40 m for IM (pixel size 25 m) and 120–240 m for WSM (pixel size 150 m). With σ_coord=1 m for IM and 7 m for WSM, the ratios between the fourth and third term in Eqs. (22) and (23) are hence 200–800 for IM and 147–588 for WSM. In this case the first three terms can be neglected compared to the fourth (see Table 1, columns 2 and 3, in which the range from minimum to maximum values for the fourth term is estimated using corresponding combinations of ΔT and σ_tr). With a grid cell size of L=300 m (IM) and 1200 m (WS), and time differences ranging from 1.2 to 5.8 d for IM image pairs and from 2 to 6 d for WSM image pairs, the uncertainties σ_div and σ_vrt were between 2.4 % d⁻¹ and 14 % d⁻¹ for IM and 3.5 % d⁻¹ and 12.7 % d⁻¹ for WSM (calculated for each image pair listed in Table 1 of Hollands and Dierking, 2011, with the corresponding tracking errors from their Tables 3 and 4). Comparing these values to the typical magnitudes of divergence and vorticity in Sect. 3.1, the respective uncertainties are too large in areas of weaker deformation.

Lindsay and Stern (2003) calculated deformation parameters for the RGPS initial velocity grid (L=10 km) and a time interval ΔT of 3 d. They use a tracking error of 100 m for RADARSAT ScanSAR images (pixel size 100 m) and assumed that the geolocation error can be regarded as bias with zero uncertainty. Hence, only the fourth term of Eqs. (22) and (23) is used (their Eq. 17), and uncertainties for divergence and vorticity are 0.5 % d⁻¹ (Table 1, column 4). However, when considering the uncertainty of the geolocation error mentioned in Sect. 3.4.1, the fourth term contributes less than the other three terms (Table 1, column 5). Only if terms $\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)$ and $\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right)$ are of magnitudes <0.001 d⁻², the first and second term can be neglected compared to the third term.

At first sight, larger time intervals and grid cells seem to be advantageous to keep the uncertainties of deformation parameters at a low level. However, larger time intervals may cause problems in the retrieval of the ice drift field, since ice structures, which serve as reference for the retrieval, may change or even vanish with time. Larger grid cells may smooth out local variations of deformation.

If the first and second term in Eqs. (22) and (23) can be neglected, i.e., when magnitudes of deformation parameters are low (which is most likely for measurements over larger spatial scales and for weak deformation events), we can determine the minimum grid cell size that is required to keep the uncertainties of divergence and vorticity below a given threshold. If we assume an uncertainty threshold of 1 % d⁻¹, then the third and fourth term of Eqs. (22) and (23) tells us that the ratio between combined position and tracking uncertainty and grid cell size should satisfy $\sqrt{\mathrm{4}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}/L\le \mathrm{0.01}\left[{\mathrm{d}}^{-\mathrm{1}}\right]\times \mathrm{\Delta }T\left[\mathrm{d}\right]$. If σ_coord≪σ_tr we obtain ${\mathit{\sigma }}_{\mathrm{tr}}/L\le \mathrm{0.01}\left[{\mathrm{d}}^{-\mathrm{1}}\right]\times \mathrm{\Delta }T\left[\mathrm{d}\right]/\sqrt{\mathrm{2}}\cong \mathrm{0.007}\left[{\mathrm{d}}^{-\mathrm{1}}\right]\times \mathrm{\Delta }T\left[\mathrm{d}\right]$. For ΔT=1 d, this means a grid cell length of roughly 150×σ_tr (uncertainty 1 % d⁻¹) or larger (uncertainty <1 % d⁻¹).

3.5 Deformation retrievals from triangular grid cells or buoy arrays

Also triangles are used for calculations of deformation parameters in SAR images (e.g., Bouillon and Rampal, 2015; Griebel and Dierking, 2018), and they form the smallest units of buoy arrays (e.g., Hutchings et al., 2011; Hutchings et al., 2012). Using the same approach as for the square above, we obtain the following for a triangle with its base a parallel to the x axis (Fig. 4):

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}\left(a^{\mathrm{2}}+b^{\mathrm{2}}+c^{\mathrm{2}}\right)}{h_a^{\mathrm{2}}{a}^{\mathrm{2}}}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)\\ & +{\displaystyle \frac{\left(\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}\right)\left(a^{\mathrm{2}}+b^{\mathrm{2}}+c^{\mathrm{2}}\right)}{\mathrm{\Delta }{T}^{\mathrm{2}}{h}_a^{\mathrm{2}}{a}^{\mathrm{2}}}}\\ & +{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{h_a^{\mathrm{2}}{a}^{\mathrm{2}}}}\left[\left(u_x^{\mathrm{2}}+v_x^{\mathrm{2}}\right)\left(a^{\mathrm{2}}+a_{\mathrm{1}}^{\mathrm{2}}-a{a}_{\mathrm{1}}\right)\right.\\ & \left.+\left(u_y^{\mathrm{2}}+v_y^{\mathrm{2}}\right)h_a^{\mathrm{2}}+\left(u_x{u}_y+v_x{v}_y\right)\left(\mathrm{2}{a}_{\mathrm{1}}-a\right)h_a\right].\end{array}\end{array}$$

Sides b and c, height h_a, and segments a₁ and a₂ are shown in Fig. 4. For the vorticity, the sum $(u_x^{\mathrm{2}}+v_y^{\mathrm{2}})$ in the first term has to be replaced by $(u_y^{\mathrm{2}}+v_x^{\mathrm{2}})$. Equation (24), which is shown here for an acute triangle (all internal angles <90^∘), is also valid for an obtuse triangle (one internal angle >90^∘) setting a₁ negative and a₂ to zero. For a right triangle with b=a, $c=\sqrt{\mathrm{2}}a$, h_a=a, and a₁=a, Eq. (24) yields

$$\begin{array}{}\text{(25a)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{6}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{a^{\mathrm{2}}}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{a^{\mathrm{2}}}}\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right.\\ & \left.+\left(u_x{u}_y+v_x{v}_y\right)\right)+{\displaystyle \frac{\mathrm{8}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{a}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{a}^{\mathrm{2}}}}.\end{array}\end{array}$$

However, if the right angle is placed at the left side of the triangle, i.e., c=a, $b=\sqrt{\mathrm{2}}a$, h_a=a, and a₁=0, the resulting equation changes to

$$\begin{array}{}\text{(25b)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{6}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{a^{\mathrm{2}}}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{a^{\mathrm{2}}}}\left[u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right.\\ & \left.-\left(u_x{u}_y+v_x{v}_y\right)\right]+{\displaystyle \frac{\mathrm{8}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{a}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{a}^{\mathrm{2}}}}.\end{array}\end{array}$$

Similarly as for the grid of squares, the contributions of terms 1–3 of Eqs. (25a) and (25b) can be neglected when geolocation uncertainties are much smaller than tracking uncertainties. When comparing the third and fourth terms of Eqs. (25) and (22) one finds that the squared uncertainty of a right triangle is 2 times the squared uncertainty of a square for a=L and identical σ_coord, σ_tr, and ΔT, which can be attributed to the reduced coverage of the triangle over the varying velocity field. For an uncertainty of 1 % d⁻¹, we obtain a value of ≤0.005 [d⁻¹] ×ΔT [d] for the ratio σ_tr∕a, corresponding to a base length a of 200×σ_tr if ΔT=1 d.

Figure 4 Uncertainty of divergence for a triangle in a spatially varying velocity field with gradients u_x, u_y, v_x, v_y. The height h_a is 2A∕a (A can be calculated from Heron's formula), and $a_{\mathrm{1}}=c^{\mathrm{2}}-h_a^{\mathrm{2}}$. Side a is the base of the triangle.

The uncertainty of divergence for the equal-sided triangle ($c=b=a$, $h_a^{\mathrm{2}}=\mathrm{3}{a}^{\mathrm{2}}/\mathrm{4}$, and $a_{\mathrm{1}}=a/\mathrm{2}$) is

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{div}}^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{6}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{a^{\mathrm{2}}}}\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)+{\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{a^{\mathrm{2}}}}\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right)\\ & +{\displaystyle \frac{\mathrm{8}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{a}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}}{\mathrm{\Delta }{T}^{\mathrm{2}}{a}^{\mathrm{2}}}}.\end{array}\end{array}$$

Note that compared to a square of length L, the area of an equal-sided triangle with base L is 0.433 A_square. The area of an arbitrary triangle with constant base increases when changing its shape from the equal-sided to the right triangle.

3.5.1 Uncertainties in position and temporal sampling

For buoys, the tracking error is zero. Itkin et al. (2017) quoted 25 m as geolocation accuracy for stationary buoys but used 50 m to account for effects of buoy drift. One of us (Hutchings) analyzed the position errors of GPS receivers in the Fairbanks (Alaska) region. The errors were normally distributed for position data collected at the same location for several days. The relative position error between pairs of GPS receivers, which has to be used for deformation calculations, was 2 m over distances of 1–10 km. Reported time intervals between acquisitions of buoy positions range from 10 s to 3 h (Hutchings, 2012; Itkin et al., 2017) with uncertainties in time less than milliseconds (see above). Hutchings at al. (2012), however, mention also a time error of 30 s, which was due to the acquisition times of the buoys not being exactly time coincident. In such an exceptional case, the second term on the right-hand side of Eq. (8) may have to be considered. If the ice drifts in the x direction (i.e., v=0), the right-hand side of Eq. (8) reads ($\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}+u^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}})/\mathrm{\Delta }{T}^{\mathrm{2}}$ (σ_tr=0 for buoys). Here we ask the following: what is the maximum value of the drift velocity for which the term $u^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}$ can still be neglected? Our criterion for neglecting it is that its value is 1 % of $\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ or less. Then the velocity u must be equal to or smaller than 424 m h⁻¹ if we assume that σ_coord=25 m and σ_ΔT=30 s. If σ_coord=2 m, the result is u=34 m h⁻¹. The speed of sea ice drift ranges mainly between 0 and 1.3 km h⁻¹, with possible extreme values around 3.6 km h⁻¹ (see Sect. 3.3), which means that the term $u^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Delta }T}^{\mathrm{2}}$ has to be taken into account in most cases. Conversely, we may ask how large the acceptable maximum temporal sampling error is so that the second term is negligible (i.e., <1 % compared to the first term). With $u_{max}=\mathrm{3.6}$ km h${}^{-\mathrm{1}}=\mathrm{1}$ m s⁻¹ and σ_coord=2 m one obtains σ_ΔT=0.3 s, and for σ_coord=25 m it is σ_ΔT=3.5 s.

3.5.2 Optimal sizes of buoy arrays

In this section we ask how large the area of a triangle-shaped buoy array has to be chosen to keep the uncertainty for deformation below a given threshold. We assume that the temporal sampling error can be neglected. The time interval ΔT is set to the temporal sampling rate of buoy positions. For buoy arrays, the tracking error is zero. With a given threshold for the uncertainty of divergence, e.g., one can use Eqs. (25) and (26) to calculate base a of right-angle or equal-sided triangles. Solutions of these simple cases can serve for approximately fixing the optimal area size for triangles of arbitrary shapes. For such triangles, the corresponding Eq. (24) cannot be directly solved since they need to be described by additional geometric parameters besides base a.

The first two terms of Eqs. (25) and (26) require the knowledge of the sea ice velocity field and its gradients. We will here focus on cases for which these terms can be neglected. This requires that $\mathrm{8}\mathrm{\Delta }{T}^{-\mathrm{2}}\gg \mathrm{6}(u_x^{\mathrm{2}}+v_y^{\mathrm{2}})$ or $\mathrm{8}\mathrm{\Delta }{T}^{-\mathrm{2}}\gg \mathrm{6}(u_y^{\mathrm{2}}+v_x^{\mathrm{2}})$. Itkin et al. (2017) analyzed deformation for constellations of three buoys using temporal sampling intervals of ΔT₁=1 h and ΔT₂=3 h, which results in $\mathrm{\Delta }{T}_{\mathrm{1}}^{-\mathrm{2}}=\mathrm{576}$ d⁻² and $\mathrm{\Delta }{T}_{\mathrm{2}}^{-\mathrm{2}}=\mathrm{64}$ d⁻². For a large fraction of measured divergence and shear data we can assume that $\left(u_x^{\mathrm{2}}+v_y^{\mathrm{2}}\right)$ and $\left(u_y^{\mathrm{2}}+v_x^{\mathrm{2}}\right)$ are smaller than 1 d⁻² (see Sect. 3.4.2) and neglect the first two terms in Eqs. (25) and (26). At low magnitudes of deformation this is also justified for ΔT₃=24 h, which gives $\mathrm{\Delta }{T}_{\mathrm{3}}^{-\mathrm{2}}=\mathrm{1}$ d⁻².

Using only the third term $(\mathrm{8}/\mathrm{\Delta }{T}^{\mathrm{2}})\times ({\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/a^{\mathrm{2}})$ the uncertainty of the divergence can be expressed as ${\mathit{\sigma }}_{\mathrm{div}}=\mathrm{71}/a$ h⁻¹ for ΔT₁=1 h and 24∕a h⁻¹ for ΔT₂=3 h, where σ_coord=25 m and the value for base a has to be given in meters. In the following we accept an uncertainty of 5 % or less relative to the majority of the magnitudes of divergence derived in Itkin et al. (2017), which are ≤0.4 d${}^{-\mathrm{1}}=\mathrm{0.017}$ h⁻¹. Hence the uncertainty is σ_div=0.00085 h⁻¹, which means that base a of the triangle has to be larger than 83.5 km for ΔT₁=1 h and 28 km for ΔT₂=3 h. If one calculates the divergence using only the position change after 24 h, the required base is 2.95∕a h⁻¹, and for σ_div=0.00085 h⁻¹ one obtains a=3.5 km. Hence, by choosing a larger time interval, acceptable uncertainties can be obtained over smaller spatial scales. If positions acquired at shorter time intervals are available, they can be used for controlling the temporal evolution of the ice drift. Using σ_coord=2 m instead of 25 m in the example given above, we obtain ${\mathit{\sigma }}_{\mathrm{div}}=\mathrm{5.7}/a$ h⁻¹ for ΔT=1 h, i.e., a base length of 6.7 km for a single measurement with σ_div=0.00085 h⁻¹ and 2.2 km for one measurement with ΔT=3 h. Itkin et al. (2017) used triangle arrays with the smallest distance between two buoys of 2 km and the largest of 70 km.

Since the area and shape of the triangle change under the action of continuous stress, the uncertainty does not simply decrease by a factor of $\mathrm{1}/\sqrt{n}$, i.e., with the number n of buoy position readings. If we assume that the three-buoy array keeps the shape of an equal-sided triangle for 24 h, with an increase in side length from a₀ to 1.1a₀ (i.e., the area of the triangle increases by a factor of 1.05), the uncertainty of the last single measurement at the end of the 24 h period is lower by a factor of $\mathrm{1}/\mathrm{1.1}=\mathrm{0.91}$ (Eq. 26, third term). Here it is assumed that the divergence is constant, the ratio u_x∕v_y is fixed by the ratio between base a and height h_a of the triangle, and the vorticity is zero, i.e., $u_y=v_x=\mathrm{0}$. As mentioned above, the position error may be as small as 2 m.

3.6 Combination of grid cells or buoys

The combination of grid cells or several buoys is one possibility to lower the uncertainty of the area σ_A. In general, the uncertainty of deformation rates is reduced when they are evaluated over a larger area, as can be deduced from the equations provided in Sect. 2.4 and 2.5. However, the uncertainty of the area, σ_A, appears explicitly only in the equations derived for the general case, Sect. 2.5. In Sect. 3.4 and 3.5 we showed that the terms including σ_A can be neglected since velocity gradients observed for sea ice are usually small. Since, on the other hand, the change in the area inside a buoy array or of a grid cell can also be used to quantify deformation (Lindsay and Stern, 2003), it is worthwhile to have a closer look at the effect of combining several grid cells or buoys.

Because buoy arrays rarely reveal simple shapes such as squares or right triangles, the uncertainties in area have to be calculated numerically using Eqs. (11) or (12). Hutchings et al. (2012), e.g., used 22 buoys, which they split into arrays of approximately equilateral triangles but also into arrays of 6, 9, and 12 buoys. Here we discuss combinations of squares and triangles.

First we investigate the effect of splitting a square or a right triangle into smaller units. We start with a square window covering N×N cells, i.e., we have 4N displacement vectors around it. Let L^′ be the length of each side of a square covering several grid cells (Fig. 5). We divide each side of the square into N segments of equal length. If N=2 then each side of the square is 2 segments of length $L^{\prime }/\mathrm{2}$, and correspondingly for N>2 it is $L^{\prime }/N$. The term $\mathrm{\Sigma }(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}{)}^{\mathrm{2}}$ is zero if both x_i+1 and x_i−1 are located on the vertical sides of the square. On the top and bottom sides parallel to the x axis, N−1 terms in the summation contribute $(\mathrm{2}{L}^{\prime }/N{)}^{\mathrm{2}}$ for each side (indicated by green bars in Fig. 3). In addition, each corner contributes $(L^{\prime }/N{)}^{\mathrm{2}}$ (blue bars). The total contribution is hence $\mathrm{2}(N-\mathrm{1})(\mathrm{2}{L}^{\prime }/N{)}^{\mathrm{2}}+\mathrm{4}(L^{\prime }/N{)}^{\mathrm{2}}=\mathrm{4}(\mathrm{2}N-\mathrm{1})(L^{\prime }/N{)}^{\mathrm{2}}$. The term $\mathrm{\Sigma }(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}{)}^{\mathrm{2}}$ contributes the same amount. Hence application of Eq. (12) yields

$$\begin{array}{}\text{(27a)}& {\mathit{\sigma }}_A^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}(\mathrm{4}N-\mathrm{2})(L^{\prime }/N{)}^{\mathrm{2}}).\end{array}$$

Since each side of the square is divided into N segments, the total number of points defining the boundary is n=4N. With $L=L^{\prime }/N$, we can rewrite Eq. (27a) as ${\mathit{\sigma }}_A^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}(n-$2)L². However, the notation in Eq. (27a), using N and L^′ instead of n and L, is preferable because it explicitly shows that ${\mathit{\sigma }}_A^{\mathrm{2}}$ decreases as N increases for a fixed L^′. Note that Eq. (27a) is valid for buoy arrays. In case of SAR images, the tracking error has to be considered as well. When ${\mathit{\sigma }}_A^{\mathrm{2}}$ is estimated for an area that deforms between acquisitions of SAR image 1 and 2, and the variance of the position error ${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ can be set to zero (see Sect.  2.1), ${\mathit{\sigma }}_A^{\mathrm{2}}=$ 0 for image 1. In image 2, however, it is ${\mathit{\sigma }}_A^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{tr}}^{\mathrm{2}}(n-\mathrm{2})L^{\mathrm{2}}$ (Lindsay and Stern, 2003).

Figure 5 Derivation of equation 30 in the x direction for N=3. Green and blue bars indicate terms to be considered in the derivation of Eqs. (28a) and (28b).

For a right triangle, we have only two contributions from the corners instead of four as for the square (Fig. 5). In the x direction, e.g., the term $x_{i+\mathrm{4}}-x_{i+\mathrm{2}}$ is zero. Hence the total contribution in the x and y direction is $\mathrm{4}(N-\mathrm{1})(\mathrm{2}{L}^{\prime }/N{)}^{\mathrm{2}}+\mathrm{4}(L^{\prime }/N{)}^{\mathrm{2}}$ or

$$\begin{array}{}\text{(27b)}& {\mathit{\sigma }}_A^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}(\mathrm{4}N-\mathrm{3})(L^{\prime }/N{)}^{\mathrm{2}}.\end{array}$$

This can be written as ${\mathit{\sigma }}_A^{\mathrm{2}}/\left({\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}{L^{\prime }}^{\mathrm{2}}\right)=(\mathrm{4}N-\mathrm{3})/N^{\mathrm{2}}$, which takes the values 1, 5∕4, 1 for $N=\mathrm{1},\mathrm{2},\mathrm{3}$, and then decreases as N increases. Note that the uncertainty initially increases from N=1 to N=2, and an improvement over N=1 is not reached until N=4.

In SAR applications, the question is whether it is preferable to use, e.g., the smallest possible (“elementary”) square cell (determined by the resolution of the ice drift field) with four drift vectors at the edges or to combine adjacent cells. Formally, the uncertainty in area for the elementary cell is $\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}{L}^{\mathrm{2}}$, and for a cell with side length $L^{\prime }=N\times L$, covering N×4 drift vectors, Eq. (27a) yields ${\mathit{\sigma }}_A^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}(\mathrm{4}N-\mathrm{2})L^{\mathrm{2}}$. Hence the uncertainty of the area increases when elementary cells are combined. However, since also the cell area increases by a factor of N², the single terms in Eqs. (13)–(21) that include the factor A⁻² decrease. When all terms except the fourth in Eqs. (22) and (23) can be neglected, it is immediately clear that the uncertainties of divergence and vorticity decrease if several elementary cells are combined into a larger square. The effect of local variations of the drift field on the deformation rate, however, can be considered in more detail when elementary cells (or smaller units of buoys arrays) are used for the calculations.

For buoy arrays it may be of advantage to use a larger number of buoys along the outline of a polygon. Here we study the example of an isosceles triangle with two sides of equal length (Fig. 6), which, e.g., comes closest to the array/subarrays used by Hutchings et al. (2012). The term $\mathrm{\Sigma }(x_{i+\mathrm{1}}-x_{i-\mathrm{1}}{)}^{\mathrm{2}}$ of Eq. (12) results in $(\mathrm{6}N-\mathrm{4.5})(L^{\prime }/N{)}^{\mathrm{2}}$; for the term $\mathrm{\Sigma }(y_{i+\mathrm{1}}-y_{i-\mathrm{1}}{)}^{\mathrm{2}}$ we obtain $(\mathrm{8}N-\mathrm{6})(h/N{)}^{\mathrm{2}}$. The areal uncertainty is hence

$$\begin{array}{}\text{(28)}& {\mathit{\sigma }}_A^{\mathrm{2}}=({\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}/\mathrm{2})\left[\right(\mathrm{3}N-\mathrm{2.25}\left)\right(L^{\prime }/N{)}^{\mathrm{2}}+(\mathrm{4}N-\mathrm{3})(h/N{)}^{\mathrm{2}}].\end{array}$$

Compared to an array consisting of three buoys at the edges of the triangle, the uncertainty can be reduced for N≥4, i.e., at least 12 buoys are required along the outline of the triangle. This also applies to the use of SAR images, when drift fields are retrieved from triangular cells.

Figure 6 Application of Eq. (12) on a triangle with two equal sides for N=3.

If the shape of an array with many buoys approximately approaches the shape of a circle with radius r, and if the sum of two line segments s connecting vertices with summation index i+1 and i differs only slightly compared to the chord length s_c between vertices i+1 to i−1, the uncertainty of the area can be estimated as follows. We require that $s_{\mathrm{c}}^{\mathrm{2}}\approx (\mathrm{2}s{)}^{\mathrm{2}}$. According to Eq. (12) the uncertainty in area is

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{\mathit{\sigma }}_A^{\mathrm{2}}=& {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{4}}}\sum _{i=\mathrm{1}}^n{s}_{\mathrm{c}i}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{4}}}n{s}_{\mathrm{c}}^{\mathrm{2}}\approx {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{\mathrm{4}}}n\mathrm{4}{s}^{\mathrm{2}}\\ & \approx {\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}n{\left(r{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{n}}\right)}^{\mathrm{2}}={\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}}{n}}{\mathit{\pi }}^{\mathrm{2}}{r}^{\mathrm{2}},\end{array}\end{array}$$

in agreement with Jansson and Persson (2014, Eq. 29). Here we use the relationship $s=\mathrm{2}\mathit{\pi }r/n$ and take into account that both the even chords (i−1, i+1) with $i=\mathrm{1},\mathrm{3},\mathrm{5},\mathrm{\dots }$ and the odd chords with $i=\mathrm{2},\mathrm{4},\mathrm{6},\mathrm{\dots }$ each approximate the total perimeter of the circle (see, e.g., Fig. 2, hexagon). To calculate the number of chords required to fulfill Eq. (29), we demand that $n^{\prime }{s}_{\mathrm{c}}(\mathrm{1}+e)=\mathrm{2}\mathit{\pi }r$, with $n^{\prime }=n/\mathrm{2}$, where e is the error between the perimeter of a regular polygon and a circle. With $s_{\mathrm{c}}/r=\mathrm{2}\sin (\mathit{\pi }/n^{\prime })$ the condition is sin($\mathit{\pi }/n^{\prime }\left)\right(\mathrm{1}+e)=\mathit{\pi }/n^{\prime }$. If $n^{\prime }=\mathrm{10}$ (i.e., a circled-shaped array with 20 buoys), e is <0.017.

3.7 Validity

It has to be kept in mind that the fundamental Eqs. (1), (2), (4), and (5) that we used for estimating the statistical uncertainties in the retrieval of deformation parameters are based on simplifying assumptions. Hence it is necessary to consider their range of validity when applying them.

3.7.1 Truncation error

The right-hand side of Eq. (5) for estimating u_x is based on the trapezoid rule applied to the contour integral on the left side. The trapezoid rule is exact if u is linear in x and y; otherwise, the nonlinear part of u gives rise to a truncation error. Define segment k of the contour integral to be the straight line from (x_k, y_k) to (x_k+1, y_k+1), and define $\mathrm{\Delta }{x}_k=x_{k+\mathrm{1}}-x_k$ and $\mathrm{\Delta }{y}_k=y_{k+\mathrm{1}}-y_k$. Then segment k of the contour integral ∮udy is estimated by ($u_{k+\mathrm{1}}+u_k)\mathrm{\Delta }{y}_k/\mathrm{2}$, as in Eq. (5), and the associated error is

$$\begin{array}{}\text{(30)}& e_k=-{\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\left(u_{xx}\mathrm{\Delta }{x}_k^{\mathrm{2}}+\mathrm{2}{u}_{xy}\mathrm{\Delta }{x}_k\mathrm{\Delta }{y}_k+u_{yy}\mathrm{\Delta }{y}_k^{\mathrm{2}}\right)\mathrm{\Delta }{y}_k,\end{array}$$

where the partial derivatives are evaluated at some point on segment k (Atkinson, 1989). As can be seen, if u is linear in x and y on segment k then e_k=0. Similar error expressions apply to the estimates of the other velocity derivatives.

Higher-order estimates for u_x could be derived, but they would not necessarily be more accurate because the ice motion may not be continuously differentiable to higher order; e.g., u_xxx and higher derivatives may not exist. Higher-order estimates would only be more accurate for sufficiently differentiable fields.

3.7.2 Spatial resolution

Equation (5) provides an area-averaged estimate of u_x. The question arises as to whether the spatial resolution (i.e., the area) is small enough to capture the spatial variability in u(x, y). One way to answer this question is to subdivide the region into smaller pieces and repeat the calculation of u_x for each piece. If the variability of u_x from piece to piece is large then the subdivision of the original area was necessary; otherwise, it was not. In practice, subdividing a region means adding new data points, which is not always possible, unless the original region is purposely chosen to consist of the union of several smaller pieces. An alternative method for determining whether the spatial resolution is adequate is given at the end of Sect. 3.8 below.

3.7.3 Temporal sampling

What temporal sampling is necessary to resolve changes in the sea ice velocity field? The velocity may be decomposed into a mean field and a fluctuating part (Thorndike, 1986). Rampal et al. (2009) showed that the variance of the fluctuating part has two regimes separated by a timescale of ∼1.5 d. Since buoys deployed on sea ice report their positions every few hours or less, their sampling frequency is sufficient to resolve the velocity and its fluctuations. The revisit time of modern satellite constellations such as Sentinel-1 is less than a day at the high latitudes of the poles, but older systems with 3 d sampling may have missed some of the deformation caused by spatial variations in those fluctuations.

3.7.4 Correlation of errors

We have assumed that different error sources are uncorrelated, and hence we have ignored the second term on the right-hand side of Eq. (2). While it is often true that spatial errors are uncorrelated with temporal errors, it may not always be the case that spatial errors are uncorrelated with each other. For example, for the distance $d=\sqrt[\mathrm{2}]{{\left(x^{\prime }-x\right)}^{\mathrm{2}}+{\left(y^{\prime }-y\right)}^{\mathrm{2}}}$ between two points (x^′, y^′) and (x, y), the full error variance of d is given by

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}={\left({\displaystyle \frac{\partial d}{\partial x^{\prime }}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{x^{\prime }}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial d}{\partial x}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}+{\left({\displaystyle \frac{\partial d}{\partial y^{\prime }}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{y^{\prime }}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{\partial d}{\partial y}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_y^{\mathrm{2}}+\mathrm{2}\left({\displaystyle \frac{\partial d}{\partial x^{\prime }}}\right)\left({\displaystyle \frac{\partial d}{\partial x}}\right){\mathit{\sigma }}_{x^{\prime }x}\\ & +\mathrm{2}\left({\displaystyle \frac{\partial d}{\partial y^{\prime }}}\right)\left({\displaystyle \frac{\partial d}{\partial y}}\right){\mathit{\sigma }}_{y^{\prime }y}+\mathrm{2}\left({\displaystyle \frac{\partial d}{\partial x^{\prime }}}\right)\left({\displaystyle \frac{\partial d}{\partial y^{\prime }}}\right){\mathit{\sigma }}_{x^{\prime }{y}^{\prime }}\\ & +\mathrm{2}\left({\displaystyle \frac{\partial d}{\partial x}}\right)\left({\displaystyle \frac{\partial d}{\partial y}}\right){\mathit{\sigma }}_{xy}+\mathrm{2}\left({\displaystyle \frac{\partial d}{\partial x^{\prime }}}\right)\left({\displaystyle \frac{\partial d}{\partial y}}\right){\mathit{\sigma }}_{x^{\prime }y}\\ & +\mathrm{2}\left({\displaystyle \frac{\partial d}{\partial x}}\right)\left({\displaystyle \frac{\partial d}{\partial y^{\prime }}}\right){\mathit{\sigma }}_{x{y}^{\prime }}.\end{array}\end{array}$$

If the coordinate uncertainties are all equal (${\mathit{\sigma }}_x={\mathit{\sigma }}_y={\mathit{\sigma }}_{x^{\prime }}={\mathit{\sigma }}_{y^{\prime }}={\mathit{\sigma }}_{\mathrm{coord}})$ and the covariances are all equal (σ_xy= ${\mathit{\sigma }}_{x^{\prime }{y}^{\prime }}={\mathit{\sigma }}_{x^{\prime }y}={\mathit{\sigma }}_{x{y}^{\prime }}={\mathit{\sigma }}_{y^{\prime }y}={\mathit{\sigma }}_{x^{\prime }x}=c$), then we obtain ${\mathit{\sigma }}_d^{\mathrm{2}}=$ 2${\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}-\mathrm{2}c$. Since the correlation between, e.g., x and y (and correspondingly for all combinations above) is $\mathit{\rho }={\mathit{\sigma }}_{xy/}\left({\mathit{\sigma }}_x{\mathit{\sigma }}_y\right)=c/{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$, we obtain ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}=\mathrm{2}{\mathit{\sigma }}_{\mathrm{coord}}^{\mathrm{2}}$ (1−ρ). In this case, a positive correlation serves to reduce ${\mathit{\sigma }}_{\mathrm{d}}^{\mathrm{2}}$ while a negative correlation serves to increase it. Since position errors are more likely to be positively correlated (due to systemic bias), ignoring the correlation terms is actually a conservative approach to error estimation.

3.7.5 Velocity discontinuities

When calculating uncertainties of deformation parameters, it is implicitly assumed that the sea ice velocity does not have discontinuities within the polygon in which the deformation is being estimated. This is because we use Eq. (5), which is based on Green's theorem. Numerous observations of the sea ice velocity field show narrow shear zones or “linear kinematic features” (e.g., Kwok, 2003; Marsan et al., 2004; Kwok, 2006) across which the velocity jumps abruptly, as a result of stresses in the ice that create leads and ridges. Some researchers, e.g., Griebel and Dierking (2017) have proposed methods to detect and isolate these discontinuities in the velocity field to avoid smoothing effects when averaging adjacent velocity vectors (e.g., for replacing outliers).

When applying Eq. (5) over an area with a discontinuity in the velocity field, a step-like function occurring between two positions r_i+1 and r_i with r_i=(x_i, y_i) is instead represented by a linear gradient. As the interval Δr is decreased, the gradient increases. Hence, there is a numerical scaling effect: e.g., divergence and shear increase when calculated on grids of velocity vectors with higher spatial resolution. A discontinuity can be defined by a threshold for the difference of the velocities on both sides of it. The threshold depends on realistic values of velocity gradients in sea ice and on the spatial resolution of the grid. The detection of possible discontinuities in a discrete field of velocity vectors, e.g., retrieved from SAR images, is helpful for the interpretation of the magnitudes of deformation.

3.8 Alternative methods of analysis

In Sect. 2, the area-averaged velocity derivatives in a region are obtained by estimating contour integrals of the velocity around the boundary of the region. Two alternatives to this boundary integral (BI) method are briefly discussed here: the least squares (LS) method and the finite difference (FD) method.

In the LS method, the velocity components u and v are modeled as linear functions of x and y, plus error. Suppose velocities (u_k, v_k) are given at locations (x_k, y_k) for k=1 to n. The linear model is

$$\begin{array}{}\text{(32a)}& {\displaystyle }{u}_k=A+B{x}_k+C{y}_k+{\mathit{\epsilon }}_k,\text{(32b)}& {\displaystyle }{v}_k=D+E{x}_k+F{y}_k+{\mathit{\delta }}_k,\end{array}$$

where the constants A, B, C, D, E, F are chosen to minimize the variance of the errors ε_k and δ_k. The velocity derivatives u_x and u_y are then B and C, while v_x and v_y are E and F. The next step is to check whether the linear model accounts for a reasonable fraction of the variance in u_k and v_k by computing the squared correlation and then whether the linear model does in fact provide a good fit to the data (by examining the spatial pattern of the errors ε_k and δ_k) or whether a quadratic or other nonlinear model is more appropriate.

The FD method provides an estimate of u_x (and the other velocity derivatives) at a single point, based on Taylor series expansions of u and v about that point. For example, suppose we have velocities u_k+1 and u_k−1 at locations (x_k+1, y) and (x_k−1, y), where $x_k=x_{\mathrm{0}}+k\mathrm{\Delta }x$. Then an estimate of u_x at (x_k, y) is

$$\begin{array}{}\text{(33)}& \begin{array}{rl}& u_x^{\mathrm{FD}}={\displaystyle \frac{u_{k+\mathrm{1}}-u_{k-\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }x}}={\displaystyle \frac{u\left(x_k+\mathrm{\Delta }x,y\right)-u(x_k-\mathrm{\Delta }x,y)}{\mathrm{2}\mathrm{\Delta }x}}\\ & =u_x+{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}{u}_{xxx}\mathrm{\Delta }{x}^{\mathrm{2}}+O\left(\mathrm{\Delta }{x}^{\mathrm{4}}\right),\end{array}\end{array}$$

where the derivatives on the right-hand side are evaluated at (x_k, y). The first term on the right-hand side is the true value of u_x at (x_k, y); the rest of the terms are the truncation error, i.e., error $=u_x^{\mathrm{FD}}-u_x=(\mathrm{1}/\mathrm{6})u_{xxx}\mathrm{\Delta }{x}^{\mathrm{2}}+$ higher-order terms.

In summary, the BI method provides area-averaged estimates of u_x, u_y, v_x, v_y; the LS method provides the best linear models of u and v, from which u_x, u_y, v_x, v_y follow; and the FD method provides point estimates of u_x, u_y, v_x, v_y.

For a rectangular region with velocities given only at the four corners, it turns out that all three methods give the same estimates of u_x, u_y, v_x, v_y, assuming the FD estimate is made at the center of the rectangle. For a general configuration of points, the three methods give different estimates. Note that in the BI method, velocities inside the boundary of the region are ignored. In the LS method, velocities farther from the mean location $\left(\overline{x},\overline{y}\right)$ have greater weight in determining the slope of the linear model. The FD method is most appropriate for regularly spaced square grids, whereas the BI and LS methods are equally applicable to irregular grids.

The LS method can be used as a diagnostic tool to determine whether the spatial resolution of the velocity data adequately captures the variability of the velocity field. Analysis of the spatial pattern of the LS residuals (errors) by standard methods (autocorrelation) reveals whether the linear velocity model is in fact a good fit to the velocity data or not. If it is a good fit, then the spatial resolution is adequate and the truncation error in the BI method is small. If it is not a good fit and sufficient data are available, the region should be divided into smaller pieces and the calculation repeated for each piece. The BI method should be used to calculate the actual (area-averaged) velocity derivatives, since it does not depend on a model that needs to be checked for goodness of fit.

4 Conclusions

In this study we derived equations for calculating the uncertainty of different deformation parameters for a given area, using displacement vectors retrieved from SAR images or buoy arrays. In the most general case, presented in Sect. 2.5, errors in measurements of position (geolocation error), velocity (determined from displacement), and area size have to be considered. Uncertainties in velocity and area size can be related to uncertainties in position measurements and (for velocity) time readings (Sect. 2.2 and 2.3). When retrieving displacements from pairs of SAR images, a tracking error has to be considered additionally.

In Sect. 3, uncertainties of divergence and vorticity are derived based on the general equations introduced in Sect. 2, assuming squares and triangles as outlines for the area over which deformation is calculated. We chose these geometric shapes since they have been frequently used in past and recent studies of deformation in sea ice. The major findings are as follows.

• The equations reveal that the uncertainties in divergence and vorticity increase with the magnitudes of the velocity gradients and with the geolocation and tracking errors. They decrease with increasing size of the area and the time interval ΔT used for calculating the velocity gradients (Sect. 3.4 and 3.5). These results agree with the recent work of Bouchat and Tremblay (2020). Since uncertainties of shear and total deformation are weighted averages of divergence and vorticity (Sect. 2.5), the conclusions drawn for the latter are also valid for the former.

• Since geolocation errors in SAR images are usually correlated over scales of ≥10 km they can be treated as a constant bias. In this case, position uncertainties are relatively small and may even be set to zero (Sects. 2.1, 2.4, 3.7.4).

• Geolocation errors in imaging modes of modern SAR systems are smaller than their spatial resolution (see Sect. 3.4.1). Errors in time readings of buoy positions and SAR image acquisitions are negligible in most cases. For buoy arrays, the magnitude of the position error may not be negligible. Here, the reader is advised to check the manual for the position sensor and pay attention to whether the error is given as standard deviation or in another format.

• The tracking error that needs to be considered for displacement fields retrieved from SAR images is on the order of the length of 1 pixel, as several studies have shown. If the geolocation error can be neglected relative to the tracking error, a good approximation for the uncertainty of divergence and vorticity valid for a square with side L or a triangle with base L is $\mathit{\sigma }=a\times {\mathit{\sigma }}_{\mathrm{tr}}/(\mathrm{\Delta }T\times L)$, where $a=\sqrt{\mathrm{2}}$ for the square and a=2 for the triangle. If squares or triangles are small, the ratio σ_tr∕L and hence the uncertainty is large.

• For a given threshold of acceptable uncertainty we estimated the necessary size of rectangular grid cells in SAR images and triangular buoy arrays, focusing on divergence and vorticity as examples (Sect. 3.4.2 and 3.5.2.). At larger temporal sampling rates, the areas can be made smaller.

• The area uncertainty of the smallest possible (elementary) cell, determined by the position of three or four adjacent displacement vectors at the edges of a triangle or square, is smaller than for a group of adjacent elementary cells with more displacement vectors on the perimeter around the group (Sect. 3.6). If, on the other hand, for an area of fixed size a variable number N of displacement vectors can be selected, the area uncertainty normally decreases with increasing N. For triangles, however, we found that the area uncertainty with six displacement vectors is larger than the one with three (see Sect. 3.6 for details).

• In Sect. 3.7 and 3.8 we provided thoughts concerning the validity of the derived equations, which assume that the velocity field inside elementary cells is continuous and can be approximated by a two-dimensional linear function. By including second-order terms or carrying out least-square fits over subregions of the velocity field, the validity of linearity can be judged. In the former case the second-order terms need to remain below a certain threshold; in the latter, the correlation coefficient should be large. Discontinuities in the velocity field should be detected before deformation is calculated to allow their impact to be assessed and to consider appropriate strategies to alleviate their impact.