NASA's RADARSAT Geophysical Processor System (RGPS) uses sequential synthetic aperture radar (SAR) images to track the trajectories of some 30 000 points on the Arctic sea ice for periods of up to 6 months. Much of the Arctic basin is imaged and tracked every 3 days. The result is a highly detailed picture of how the sea ice moves and deforms. The points are initially spaced 10 km apart and are organized into four-cornered cells. The area and the strain rates are calculated for each cell for each new observation of its corners. The accuracy of the RGPS ice tracking, area changes, and deformation estimates is needed to make the dataset useful for analysis, model validation, and data assimilation. Two comparisons are made to assess the accuracy. The first compares the tracking performed at two different facilities (the Jet Propulsion Laboratory in Pasadena, California, and the Alaska SAR Facility in Fairbanks, Alaska), between which the primary difference is the operator intervention. The error standard deviation of the tracking, not including geolocation errors, is 100 m, which is the pixel size of the SAR images. The second comparison is made with buoy trajectories from the International Arctic Buoy Program. The squared correlation coefficient for RGPS and buoy displacements is 0.996. The median magnitude of the displacement differences is 323 m. The tracking errors give rise to error standard deviations of 0.5% day−1 in the divergence, shear, and vorticity. The uncertainty in the area change of a cell is 1.4% due to tracking errors and 3.2% due to resolving the cell boundary with only four points. The uncertainties in the area change and deformation invariants can be reduced substantially by averaging over a number of cells, at the expense of spatial resolution.
The RADARSAT Geophysical Processor System (RGPS) is a sophisticated computer system that takes RADARSAT synthetic aperture radar (SAR) images of the Arctic Ocean and produces estimates of the sea ice motion and deformation. The goal of the RGPS is to provide datasets to improve our understanding of the impact of sea ice on climate. Sea ice motion from historical sources such as drifting ships and Russian ice stations revealed the general circulation of the Arctic Ocean, including the Beaufort Gyre and the Transpolar Drift Stream. Since 1979 the International Arctic Buoy Program (IABP) has provided ice motion data from a network of buoys, enabling details of the mean circulation and changes in the circulation pattern to be detected (Colony and Thorndike 1984; Rigor et al. 2002). Sea ice deformation drives changes in the ice thickness distribution: divergence of the ice opens cracks where new ice can grow, and convergence causes ice to pile up in ridges. Thus ice motion and deformation are important for understanding the momentum, mass, and energy balance of the Arctic Ocean. Current methods for measuring ice motion and deformation (e.g., buoys) do not offer the spatial detail available from SAR imagery, and current methods for measuring ice thickness are limited to a few point values (upward-looking sonar buoys) or values along occasional transects (submarine sonar). The RGPS offers the unique opportunity to obtain fields of ice motion and deformation, independent of cloud cover and with broad spatial coverage, yet it is still a new and experimental satellite data processing system. Here the accuracy and limitations of its products are documented in order that they may be useful to a broad range of scientists.
The RGPS products include ice motion, ice deformation, ice age and thickness distributions, radar backscatter distributions, winter multiyear ice concentration, summer open water concentration, and the dates of melt onset and freeze-up (Kwok et al. 1995; Kwok 1998). The products build on each other. The ice motion is computed by tracking a large number of points in sequential images with the maximum correlation technique (Fily and Rothrock 1987; Kwok et al. 1990). Individual points are tracked at 3-day intervals (typically) for an entire season to determine their trajectories. The points define the boundaries of 10-km Lagrangian “cells” whose areas are then computed. Area changes of the cells lead to inferences about the creation or destruction of young ice. Finally, the thickness of new and young ice is computed using the local surface air temperature history and a simple freezing–degree–day parameterization.
The RGPS is funded by the National Aeronautics and Space Administration (NASA) and was developed and built at the Jet Propulsion Laboratory (JPL) in Pasadena, California. The Alaska SAR Facility (ASF) in Fairbanks operates a second RGPS and processes data in a complementary effort to that of JPL. Input to the RGPS consists of ScanSAR images from the Canadian RADARSAT satellite, which carries a C-band (5.3 GHz) SAR. The active microwave radar penetrates clouds and darkness to image a swath on the earth 460 km wide with a pixel size of 100 m. The 24-day exact-repeat polar orbit (inclination 98.6°) has subcycles of 3 and 7 days, with more frequent coverage at high latitudes where swaths overlap substantially. The data acquisition plan for RGPS calls for complete Arctic coverage every 3 days, and this is generally achieved in the western Arctic. The raw SAR data are received, processed, and geolocated at ASF to create the ScanSAR images. The eastern Arctic is completely imaged every 6 days. The raw data are collected at the Tromso Satellite Station and the tapes are sent to ASF. The centers have processed almost 3 yr of images from the western Arctic Ocean (Alaska station mask) and 6 months from the eastern (Tromso station mask). Some minor data gaps exist but spatial coverage is generally good (Fig. 1). (The RGPS data are currently available online at www-radar.jpl.nasa.gov/rgps/radarsat.html.)
RGPS products have already been used to study a number of aspects of Arctic sea ice. Kwok and Cunningham (2002) used data products from the first winter (1996/97) to compute the area and volume of new ice production, and the areal coverage of multiyear ice. Kwok (2002) then compared the area and volume production from the second winter (1997/98) with that of the first winter, finding substantial differences. Lindsay (2001) used RGPS ice thickness estimates as the basis for computing the albedo of the ice cover. The ice motion products computed by RGPS in the vicinity of the drifting ice station SHEBA (1997/98) (Perovich et al. 1999) have been used by Lindsay (2002) and by Stern and Moritz (2002) to compute the ice deformation there, which in turn has been used in a Lagrangian model of the ice cover (Lindsay 2003). Richter-Menge et al. (2002) compared in situ stress measurements with RGPS-derived ice deformation near the SHEBA site. Hopkins (1999) used the ice motion products to drive a granular model of sea ice. Other researchers are using RGPS products to validate both isotropic and anisotropic sea ice models (Hibler 2001; Lindsay et al. 2003).
The RGPS products constitute a large and complex dataset, with mixed spatial and temporal information (see the product specification document, Kwok et al. 1999). The sampling is irregular in both time and space, so that any estimates of the ice conditions at a particular day and location must account for the local sampling characteristics. There are five seasonal datasets currently available for regions covered by the Alaska station mask and the one season available for the Tromso station mask. The Alaskan mask covers the Beaufort, Chukchi, and East Siberian Seas north to within 200 km of the Pole, and the Tromso mask covers the Laptev, Kara, and Barents Seas and Fram Strait north also to within 200 km of the Pole (Table 1). The satellite configuration and orbit precludes imaging close to the Pole.
Examples of the spatial coverage for the Alaska station mask are shown in Fig. 1. The first panel shows the coverage near the beginning of the SHEBA experiment (Perovich et al. 1999), as it appeared in May 1998. Trajectories initialized both in November 1997 and in January 1998 are included. These later trajectories filled in coastal areas that were not ice covered in the fall. The data give excellent coverage in the central Beaufort Sea. Within about 200 km of the coast, the trajectories have an initial 25-km spacing instead of the 10-km spacing farther out. This was done to save processing time, because the large deformation rates of the coastal regions make them difficult to track. This difficulty is also reflected in the loss of most of the coastal points north of Alaska. Gaps in the spatial coverage caused by divergent motion at leads have been filled with new trajectories, so not all trajectories date back to the time of initialization. However, a number of spatial gaps remain unfilled where some trajectories failed to obtain the minimum 15-day sampling interval. Figure 1b shows the trajectories near the end of August 1998. These trajectories were initialized in a new grid that covered the entire region in May 1998, but difficulty in tracking the ice during the summer, when the backscatter contrast is low, has caused the loss of over half of the trajectories by late summer. Liquid surface water from melting snow and melt ponds dominates the radar backscatter, obscures the signal from the ice underneath, and lowers the contrast in the images.
There are four types of RGPS products: Lagrangian ice motion, ice deformation, backscatter histograms, and ice age and thickness distributions. Each is produced for monthly intervals. In addition, the date of melt onset is also provided. The first two are cumulative in that all the trajectory or deformation information is included from the beginning of the season in each monthly product. The two distribution products include only distributions determined during the month. The ice age and thickness are not provided for the summer, because there is no direct relationship between ice deformation and new-ice growth in the summer months. This paper reviews the accuracy of the more fundamental products, the Lagrangian ice motion product (trajectories) and the deformation product (cells). A comparison of the thickness distributions to thin-ice thickness estimates derived from AVHRR imagery is the subject of a companion paper (Yu and Lindsay 2003, manuscript submitted to J. Geophys. Res.).
2. Lagrangian ice motion
The ice trajectories for individual points are provided in the Lagrangian ice motion (LIM) product. This product provides the time, location, and a tracking quality flag for each observation for each trajectory, as well as information about the images used in the ice tracking. It is the most fundamental product, from which all other products are derived. The intervals between observations range from 1.5 h (one orbit) to 15 days. Any trajectory that fails to obtain an observation for 15 days is no longer tracked. The locations are given in kilometers on a polar stereographic grid centered at the Pole and aligned with the x axis along the 45°E meridian.
a. General discussion of errors
The errors in the displacements computed for the trajectories can arise from two sources: errors in determining the location of the ice in the second image that corresponds to the point in the first image, and errors in the geolocation of either of the two images. Tracking errors, the first type, arise from ambiguities in determining the location at which the correlation between the backscatter from the first image and from the second image is maximum. This ambiguity can arise from deformation in the ice, from rotation of the ice, and from low contrast (variance) in the SAR images. Related to the tracking errors are discretization constraints: the SAR images used as input to the RGPS consist of regular arrays of square pixels of finite size (100 m), so that horizontal and vertical positions are constrained to multiples of 100 m. This must be kept in mind when considering the size of the tracking errors, and when considering small relative motions between points (small deformation).
where a1 and a2 are the initial and final locations, ega1,2 are the geolocation errors, and efa is the tracking error. The error variance in the displacement estimate is expressed as the sum of the error variances for tracking and geolocation
One of our goals is to determine this error variance for the RGPS data. The magnitude of the displacement difference for two locations made with different measurement systems is
and the error variance in the difference is
where the errors are all assumed to be independent. The uncertainty in the RGPS estimated displacement is then obtained by evaluating ɛ2Δab and making assumptions about the tracking and geolocation errors of the position vectors.
One approach is to estimate the error variance in the displacement difference ɛ2Δab by
where the sum is over the N pairs of displacement measurements. The two components are combined by including both in the summation. Assuming the mean component difference is zero, then
where rmsd is the root-mean-square difference of the magnitude of the vector difference in the two displacements. Because it is derived from the squares of the differences, it is very sensitive to a few large outliers.
A second approach is to use a statistic that is less dependent on outliers, such as the median or 95th percentile order statistics. In this approach the order statistic of Δab is evaluated from the comparison datasets and compared to Monte Carlo simulations of Δab to give a robust estimate of the sum of the error terms in (4).
Two comparison datasets are analyzed that address the two sources of error. One is a comparison between two RGPS tracking efforts with the same input SAR images. In this comparison the geolocation errors are not independent and indeed cancel so that the tracking error can be estimated from
The second is between buoy displacement measurements and RGPS displacement estimates. This comparison includes geolocation errors for both datasets, but the tracking error for the buoys ɛfb is zero (the buoy remains fixed relative to the ice it sits on). The combined tracking and geolocation error ɛRGPS = 2ɛ2gr + ɛ2fr is then obtained from (4) and (6):
where ɛgr is the geolocation error of the RGPS, ɛfr is the tracking error, and ɛbuoys = ɛgb is the positioning (geolocation) error of the buoys.
The difference in displacement found from two estimates (either two RGPS estimates or an RGPS estimate and a buoy estimate) occasionally exhibits large values. These large displacement differences may be due to error or may be due to real differences in the ice motion found from two closely spaced points that are located on opposite sides of an active lead at a time of substantial shear or that experience large rotation. These few large difference values lead to much larger root-mean-square (rms) differences than mean or median differences.
Each observation of a trajectory is given a quality indicator Q_flag between 0 and 222. The flag indicates how good the correlation is between the backscatter fields from the two images in the vicinity of the trajectory point and if manual intervention has been used to track the point. Values of Q_flag between 1 and 6 imply no operator intervention, with 1 indicating the highest correlation and 6 the lowest. Thus points with Q_flag = 1 have the highest reliability or confidence in their computed positions. Values of Q_flag above 6 indicate operator intervention, that is, manual help in selection of corresponding points in the second image. In the May 1998 product, 52% of the observations have Q_flag = 1, 42% have Q_flag between 2 and 6, and 6% have Q_flag > 6.
b. Comparison of two independent tracking efforts
The original RGPS facility at JPL is mirrored by the facility at ASF that began operation in 1999. The ASF facility processed the winter data for the Alaska station mask from November 1997 to mid-May 1998. The first large-scale summer tracking exercise was performed at JPL for data beginning at the start of May 1998 and also covering the entire Alaska station mask. Thus there is a 12-day period during which the data from the two facilities overlap in time. The trajectories from the ASF data were initialized in November and had drifted and deformed for 6 months at the time the summer tracking began. Those for the JPL data were initialized with a new 10-km grid at the beginning of May. These two datasets offer an opportunity to determine the errors that arise from independent tracking efforts. The calibrated and geolocated SAR images (input data) were identical for the two facilities during the overlap period, the software for the tracking was the same, but the operators and the points on the ice that were tracked differed.
All pairs of displacements were found in which the observation time intervals were identical and the distance between the ASF and JPL points was initially less than 10 km. A total of 126 064 pairs were found. Table 2 shows that for all points less than 2 km apart the mean displacement difference (vector magnitude) was 0.090 km, smaller than the pixel size of the SAR images. The displacement difference is related to the quality flag. The mean displacement difference for cases in which both facilities assigned the same quality flag and for which the initial separation distance was 2 km or less increased from 0.032 to 0.332 km as the quality flag increased from 1 to 6. Flags 6 or smaller indicate that no manual intervention occurred and include 96.2% of the ASF dataset and 98.8% of the JPL dataset for the matched pairs. The quality flag is related to the level of correlation found in the SAR images in the vicinity of a point after tracking. In regions of deformation, this correlation would be less because of the changes occurring in the backscatter patterns, so that the increase in the mean difference with increasing quality flag values may also be related to real deformation or rotation events.
The dependence of the differences on the quality flag shows that the differences are very small when no operator intervention is required, but that when it is required there were significant differences of opinion between the operators on how to intervene. The mean difference in these cases was 2.586 km. Operator intervention entails the operator indicating manually where a point from the first image should be located in the second image in regions of high deformation (or in the summer, low contrast). These are places where the system was not able to automatically determine the location of maximum correlation between the two images. It is not surprising that there should be differences between operators in these cases because deformation can greatly confuse the scene interpretation. Note that there were just 26 cases (out of 7580) in this comparison dataset for which both systems required operator intervention.
The rms difference is sensitive to the few large outliers. For a maximum initial separation of 2 km (all quality flags), only 6.1% of the displacement differences exceed 0.100 km (one pixel) and the 95th percentile is just 0.224 km, yet the rms difference (0.477 km) is over 5 times as great as the mean difference. The large displacement differences in the maximums may represent the occurrence of shear cracks between the two points as well as difficulty deciding on which side of a crack a point should be placed when a crack occurs very near the initial location of a trajectory.
The locations of all the comparison points are shown in Fig. 2a, and the dependence of the displacement differences on the separation distance is shown in Fig. 2b. The mean difference falls below 1 pixel for separations less than 1 km. The median difference is zero for separations of 1.5 km or less. The rms difference is large for all separations, but for only the best quality flag (Q_flag = 1) the rms difference is much lower and is comparable to the mean difference for all quality flags. To extrapolate the rms difference to zero separation, we determined a linear least squares fit of the squares of the differences for points with separations less than 5 km (N = 47 968). The square root of the intercept of this fit yields an rms difference of 0.366 km at zero separation.
c. Comparisons with buoys
A second quality measure of the trajectories may be found in comparisons to buoy motion. We compared the displacements measured by RGPS to buoy positions obtained from the IABP. The position error for the buoys reported by the IABP is 0.3 km (Thorndike and Colony 1980). The buoy locations were estimated for the RGPS measurement times from the 3-hourly IABP data with linear interpolation. This linear interpolation adds just a few percent to the nominal buoy position error. As we will see, the RGPS data suggest that the position errors for the buoys are smaller than this nominal value. We assembled all pairs of buoy and RGPS displacement measurements within 30 km of each other (initial locations) and with the time intervals for the displacements up to 12 days.
In one comparison of a single buoy trajectory and the nearest RGPS trajectory, the buoy appears to remain on the same rigid plate, or floe, as the nearby RGPS trajectories (Figs. 3a and 3b). The distance to the nearest trajectory remains nearly constant and the mean difference in the displacement measured by the buoy and the nearest trajectory is small, 0.7 km. The other nearby trajectories (not shown) also move with the buoy and show little relative motion, although some solid-body rotation is apparent and accounts for much of the displacement differences. A second example is of a buoy located near a strong shear line (Figs. 3c and 3d). The distance to the nearest RGPS trajectory changes little, but for three observations the difference in displacement between the two is quite large, about 6 km. The relative motion of nearby buoys shows a shear line, and it is apparent that the buoy and the nearest RGPS trajectory were located on different rigid plates during a single shearing event.
The RGPS and buoy trajectories are not as closely matched as in the comparison of the JPL and ASF May 1998 products. Table 3 shows the mean, rms, median, 95th percentile, and maximum differences in displacement for all RGPS trajectories within either 2 or 5 km of a buoy for each of the five RGPS datasets considered. The improvement for the closer trajectories is due to the decreasing likelihood of a deformation event occurring between the buoy and the RGPS trajectory as well as a smaller solid-body rotational component for the displacement differences. For initial separation distances less than 5 km, the distribution of the displacement difference magnitudes is lognormal except for a heavy tail involving about 5% of the observations. The 95th percentile values show that although the mean differences are small, there are a few differences that are quite large for both separation distances. These few large differences account for the large rms differences compared to the means. The difference between the buoys and the RGPS trajectories is larger in regions of high deformation, but the old problem of separating out real error from differences caused by deformation between the buoy and a trajectory is also more severe. The summer dataset exhibits the largest differences. This is due, at least in part, to the greater deformation in the ice pack at small scales in the summer and the consequent increased likelihood of real deformation occurring between a buoy and a nearby RGPS point. It bodes well that the statistics of the displacement differences for the Tromso dataset are not significantly worse than for the others, in spite of the large deformation rates in the European sector of the Arctic.
The dependence of the displacement difference on the quality flag is seen in Table 4 for initial distance ranges of less than 2 km. Here all five datasets are combined. The mean and median differences are again substantially smaller than the rms differences but not to the same degree. The displacement differences increase from a mean of 0.464 km for the best quality flag, to 2.231 km for quality flags greater than 6 (implying manual intervention). The difference has little dependence on the quality flag for flags 6 or less and there is not as large an increase in the displacement differences when operator intervention occurs compared to the previous case. The larger quality flags, including those indicating manual intervention, are more likely at times of large deformation, so some of the observed displacement differences are likely real, not error. The comparisons between the RGPS displacement estimates and those of the buoys are illustrated in Fig. 4. The smaller number of pairs does not allow for a robust determination of an intercept of the squares of the differences versus separation distance.
The RGPS dataset initialized for November 1997 consisted of mostly multiyear ice. A large region of open water in the Chukchi and Beaufort Seas could not be seeded with points to track at that time. Therefore another dataset was initialized for these seas in January 1998 once this region had frozen and it consisted of mostly first-year ice. It has the same fraction of observations with Q_flag = 1 or Q_flag < 7 as in the November 1997 dataset so we conclude that the tracking algorithm worked equally well whether tracking multiyear or first-year ice.
Here analysis is focused on the displacement differences, but the correlation in the displacement vectors (expressed as complex numbers) between the buoys and the RGPS estimates is quite high, with a squared correlation coefficient of 0.996 for all cases and 0.998 for those not requiring manual intervention. Such a high correlation is possible because the mean magnitude of the displacement, 18.2 km, is large compared to the displacement differences. The error in the displacement directions is usually small because of this large mean displacement.
In a separate ice tracking effort using RADARSAT images, Wilson et al. (2001) compared sea ice displacement in northern Baffin Bay derived from two sources: sequential RADARSAT ScanSAR images, and beacons deployed on the sea ice. Features in the ScanSAR images were tracked using the Canadian Ice Service “Tracker” algorithm, which uses a correlation technique similar to that of the RGPS. The beacons transmitted their positions hourly via satellite using the same system as the IABP buoys (Service Argos). The rms Tracker errors in this study were 3.8 km in winter/spring and 6.8 km in summer/fall. The large errors probably reflect the dynamic ice conditions of northern Baffin Bay relative to the central Arctic Ocean.
d. Displacement error variance
For the comparisons between the two May 1998 RGPS datasets, the initial locations of the displacement vectors are known exactly (or by definition) and the uncertainty in the final locations are assumed to be the same. The tracking error is given by (7). For all quality flags and separations extrapolated to 0 km, the rms difference is 0.366 km and ɛf = 0.184 km, but for Q_flag = 1, about half of the points, the rms difference is 0.071 km (Table 2) and ɛf = 0.035 km, much less than one pixel. For all those points that did not require manual intervention, the rms difference is 0.266 km and ɛf = 0.133 km.
The rms difference is very sensitive to a few outliers, but the order statistics, such as the median or the 95th percentile value, are less sensitive. We performed a Monte Carlo simulation of the displacement difference with 50 000 trials assuming normally distributed tracking errors for each component and with the components discretized to 0.1 km. The QQ plots show that the component differences in the comparison datasets are nearly normally distributed, allowing for the discretization, if the smallest and largest 2% are trimmed. The observed median value for all comparisons less than 2 km apart is 0 because of the discretization, so the less robust 95th percentile value is used for comparison with the simulations. The simulation gives a 95th percentile value the same as the observed, 0.224 km, for the range of values between ɛf = 0.087 and 0.102 km. This is consistent with the estimate of Kwok and Cunningham (2002) of 0.100 km based on the standard deviation of the observed area changes, and their value seems to be a good choice for the tracking error.
For the buoy comparisons the combined tracking and geolocation error is given by (8) with ɛbuoys = 0.3 km. For all quality flags and for separation distances less than 2 km, rmsd = 1.431 km and ɛRGPS = 0.919 km. For Q_flag ≤ 6, rmsd = 1.002 km and ɛRGPS = 0.567 km.
As with the May 1998 comparison, a more robust estimate is made with the order statistics. This time the median is used because the observed value, 0.323 km, is not smaller than the discretization. The error in the displacement in (4) for the RGPS–buoy comparison is
We performed another Monte Carlo simulation of the displacement difference with 50 000 samples assuming normally distributed errors for the components of both the buoys and the RGPS with ɛ2RGPS = 2ɛ2gr + ɛ2fr = 0.1 km and ɛbuoys = 0.3 km. The QQ plots show that the component differences in the comparison datasets are nearly normally distributed if the smallest and largest 5% are trimmed. Since the median is independent of the very large and small values, the nature of the distribution tails is not important for this analysis. The simulation gives a median value of the right-hand-side displacement difference of 0.629 km. Because this is greater than the value found for all observations (0.323 km, the left-hand side) the assumed error variance for the RGPS or for the buoy (or both) is too large. Even for ɛRGPS = 0, the simulated median is 0.502 km, greater than the observed value. This suggests that the error variance for the buoys is less than what we have assumed and we can only determine an upper bound for ɛRGPS from the median statistic. If the buoy error were 0.100 km (smaller than the nominal value but perhaps realistic for the smoothed 3-hourly values used here; I. Rigor 2002, personal communication), the uncertainty in the displacement that provides the best match with the observed median difference is ɛRGPS = 0.286 km.
These data can also be used to put an upper bound on the buoy errors. Even if the RGPS tracking and geolocation errors were zero, the buoy positioning errors must be less than 0.200 km. A realistic but still optimistic estimate of the combined RGPS geolocation and tracking errors is 0.200 km, implying a maximum buoy positioning error of 0.158 km, almost half the nominal value listed by Thorndike and Colony (1980). The improvement may be due to better technology found on the Argos satellites nearly 20 yr after the Thorndike and Colony report, but it is also likely due to the smoothing used by the IABP in constructing the 3-hourly dataset.
Geolocation errors in the SAR image products contribute to the buoy-derived displacement differences. These image products (input data) were generated at ASF and were common for the tracking at both the JPL and ASF facilities. Geolocation errors in the input images would not be detected in comparisons of the displacement vectors calculated by the two facilities. The absolute geolocation error of the ScanSAR Wide B images determined by ASF was 0.277 km (Martyn et al. 1999). This is comparable to the estimate of the tracking and geolocation errors combined, so that the comparison datasets suggest that the geolocation error of the images is less than previously thought. Wilson et al. (2001) found similar geolocation errors of 0.225 km in the ScanSAR images they worked with.
Because the geolocation errors cancel in the ASF/JPL comparisons, our best estimate of the error standard deviation of individual displacements, including geolocation errors, are those determined from the buoys. However, the geolocation errors of the SAR images do not enter the deformation calculations, because the deformation rates are calculated from relative displacement vectors determined from a common set of images.
3. Ice deformation
a. The deformation products
The RGPS ice deformation products are also cumulative so that the product for the latest month contains all the information from the previous months back to the time the analysis was initialized. The initial areas of the cells are either roughly 100 or 625 km2, depending on the spacing of the trajectories. Additional trajectories are added when an edge becomes more than twice its initial size. Although trajectories are added throughout the 7 months, no new cells are added after the initialization of the tracking grids. The deformation of a cell is determined whenever the constituent trajectories all have common observation times.
The monthly deformation products include the current observation time, the time interval since the last observation for the cell, the current area, the change in the area over the interval, the center position of the cell, the displacement of the center, and the four components of the strain tensor (not the strain rate tensor). The number of trajectories that make up the cell and the trajectory identification numbers are not included. Figure 5 is a map of all of the cells at the end of May 1998. The 10% of the cells with the largest total deformation rate are shown in black. This threshold is 4.3% day−1. Note the lines of deformation, also called linear kinematic features, that delineate approximately rigid plates in some regions. The large area with no data in the Beaufort Sea is caused by advection of the cells that covered the sea at the time of initialization in January 1998. Cells more than about 500 km from the coast of Alaska were initialized in November 1997.
Summary plots for all of the cells in the May 1998 (ASF) deformation product that had an initial area of near 100 km2 are shown in Fig. 6. While 65% of the cell observation intervals are 3 days or less, 5% are 6 days or more. A few of the sample intervals for cells are large, up to 50 days, although the maximum sample interval for trajectories is 15 days. This is because it is occasionally not possible to find common observation times for all corners for an extended period. Only 0.15% of the observations have intervals greater than 15 days. The cell areas have been corrected for map distortion and are consequently 2%–6% greater than the area determined from the positions reported for the corner points.
Also shown in Fig. 6 are histograms of the area change, the divergence, the shear, and the vorticity of the velocity vectors. The distribution of total area change, as a fraction of the initial area is slightly asymmetric, indicating more increased area than decreased. This might be expected given the resistance to convergence supplied by the ice strength. Some of the areas are more than twice as large as the initial areas and some show negative areas.
Negative areas arise from the formula for cell area (14) when the vertices of the cell move in such a way that the segments connecting them cross one another. We do not consider these negative areas to be errors so much as flags to indicate highly distorted cells. Area changes and area fractions for bins of age, ice thickness, and backscatter for these cells are meaningless. In the May 1997 and May 1998 products, 1.9% and 1.5%, respectively, of the observations show negative area. It is also important to note that cells that have begun to turn into a butterfly have a computed area made up of a positive portion and a negative portion. These cells, which cannot be identified unless the area goes negative, also have area-change values that are not meaningful. To ensure more accurate deformation estimates, we often restrict our attention to those cells with areas within 50% of their initial areas.
The spiky nature of the area change and deformation histograms arises from the discrete 100-m pixels in the SAR images. For the initial very narrow distribution of cell areas, the allowable deformation values for small deformations are limited. As the distribution of cell areas and shapes is broadened over time, the allowable deformation values for small deformations broadens and a smooth underlying distribution results; however, the inherent uncertainty in the deformation due to the discrete nature of the grid remains. For example, if the initial cell is 10 km × 10 km, the smallest area change detectable is 0.5%, or about 0.17% day−1 over a 3-day interval (the largest spikes with near-zero divergence). This discretization uncertainty is a substantial fraction of the standard deviation of the divergence itself. Because of the broad tails, the standard deviations of the invariants are greater than what might be guessed from the plotted distributions. For divergence the standard deviation is about 1% day−1 and for shear and vorticity about 3% day−1. The magnitudes of the shear and vorticity are larger so the discretization uncertainty is less consequential. About 5% of the deformation observations indicate very near zero deformation (within ±0.05% day−1). The spikes are not evident in the final area change because the distributions of cell shapes and areas are broad by the end of 6 months and 1-pixel changes in the shape of the cell can imply different area changes.
b. Deformation errors
The strain rate of the ice arises from spatial variability of the ice velocity, where the velocity components u and υ are computed from the displacement components and the time interval between the SAR images. The spatial gradient in the velocity field has three invariants (invariant with respect to coordinate system): divergence, shear, and vorticity. These are computed from the components of the strain rate tensor:
The strain rate components are computed by the RGPS from approximations of the line integral around the boundary of each of the cells. For example,
where n is the number of vertices for the cell and A is the area of the cell. The area is computed as
where xi and yi are the position coordinates of the vertices, which trace the boundary in a counterclockwise sense. These equations are derived from the trapezoid rule for integration. If the velocity is nonlinear along any side, an error arises in Eq. (13). In particular, if the velocity is discontinuous, such as along a shearing lead, an error will occur. Errors in the deformation rates then arise from either tracking errors in the estimated displacements (velocities) or in boundary-definition errors originating in the approximation of the line integrals.
c. Uncertainty from propagation of errors
The theory of propagation of errors (Bevington 1969) may be used to estimate the impact of displacement errors on the deformation estimates. The uncertainty in the area of a cell depends on the uncertainty of the positions of each of the n vertices and the size and shape of the cell. The error variance of the area (and the area change) is
where ɛi is the error standard deviation of the position of each of the vertices of the polygon. For a square cell with n points defining the boundary and identical error variances for the positions, the expression simplifies to
where ɛf is the tracking error and L is the length of any one of the n segments. If n = 4, L = 10 km, and ɛf = 0.1 km (from the winter–summer comparison), then ɛA = 1.4 km2. This value is confirmed by numerical simulations. Kwok and Cunningham (2002) also performed numerical simulations and obtained similar results. In the May 1998 winter product, 65% of the area changes did not exceed this value. This indicates that significant area changes are confined to just some of the cells at any one time. This uncertainty in the estimated cell area, even for cells experiencing no deformation, can lead to the spurious production of ice. However, Kwok and Cunningham (2002) found in simulations that this source of error was less than 1% of the total ice production.
The errors in the invariants depend on the tracking errors and the elapsed time between the images Δt (the timing error of the images is negligible). The error in the vorticity, shear, and divergence are all the same and proportional to the area change error
If L = 10 km, Δt = 3 days, and ɛf = 0.1 km, then ɛshr = ɛvrt = ɛdiv = 0.5% day−1.
d. Boundary-definition errors
The cells are typically defined by the four corner points, or occasionally additional points if an edge grows to more than double its initial length. With just four points defining the boundary of a cell, significant errors may result in deformation determinations because the four points may not adequately define a material boundary in the presence of localized lines of deformation (leads). Thorndike (1986) examined the error in estimating the deformation on scales of hundreds of kilometers. He imagined computing ∂u/∂x (and the other partial derivatives) for a circular region, using the line integral formulation, in which the velocity (u, υ) is known at three or more points on the circle. His results are based on the large-scale velocity correlation functions for a homogeneous and isotropic field. He shows that there is not much additional improvement (reduction in error) in using more than six measurement points around the circle when the radius varies from 100 to 500 km. But because the correlation functions are not well known for separations less than 100 km, his analysis and methods do not apply in the present case of estimating the error in small-scale deformation.
Numerical experiments show that if a shear crack passes through a field of square cells with no convergence or divergence, spurious area changes are calculated. Some cells exhibit opening (positive area change) and others closing (negative area change) and the amount depends on the relative motion across the crack and the angle the crack makes with the cell edges. With randomly oriented cracks, the rms error per unit crack length is 15.0% of the sliding distance. For example, if the sliding distance is 1 km, spurious openings averaging 0.15 km2 km−1 of crack will be calculated as well as spurious closings of an equal magnitude. If the cells are defined with eight points instead of four, the error drops almost in half to 8.4%. The bias toward additional opening and closing created by the limited number of points used to define a cell would cause the rate of ice production to be overestimated. This overestimation may be offset by the 3-day sample intervals, which miss some opening and closing events.
The limitations of using just four points to determine the area changes of a cell can be further investigated using a special 5-km Eulerian dataset produced at JPL for the Surface Heat Budget of the Arctic Ocean (SHEBA) project. This dataset consists of 5-km cells within a 200-km box centered on the drifting icebreaker Des Groseilliers (Stern and Moritz 2002). The ship was used as a base of operations for the SHEBA ice camp (Perovich et al. 1999). The special SHEBA dataset was tracked in an Eulerian fashion: for each pair of SAR images, the initial points lie on the same regular 5-km grid. This allows the area change of a 10 km × 10 km square to be computed in either of two ways: (i) using only the four corners, initially spaced 10 km apart, or (ii) using the corners and the midpoints, initially spaced 5 km apart. In 10 623 comparisons of method (i) versus method (ii), the standard deviation of the difference in the computed area change is 3.4 km2 (3.4% of the initial area) and the squared correlation is only R2 = 0.46. The difference in the area change can be partitioned into two factors: an error due to tracking ɛA,f, and a boundary-definition error ɛA,bd. The factor ɛA,f can be computed from (15), and it equals 1.2 km2. If the variance of the difference in the areas is assumed to be the sum of the error variances from the two sources, ɛ2ΔA = ɛ2A,f + ɛ2A,bd, then the standard deviation of the boundary-definition error is 3.42 − 1.22 = 3.2 km2. This can be compared with the error arising only from tracking the corners of the cell, 1.4 km2. The SHEBA 5-km dataset clearly shows that the boundary-definition error is substantially larger than the tracking error.
e. Comparison of May 1998 ASF and JPL deformations
The coincident tracking efforts from May 1998 can also be used to compare the computed deformation rates. The comparisons are not expected to be exact because even for cells with centers collocated exactly, the positions of the constitutive trajectories of the two datasets are different. The cell shapes in the winter product, which have evolved for a 6-month period, will often be greatly distorted compared to those of the summer product, initialized just at the beginning of May. Nevertheless, the comparison of the two products gives a good indication of the magnitude of the difference in the deformation obtained for a single location but with different boundary points. For cells with center locations less than 2 km apart, one-fifth the nominal width of the cells, and for which the cell areas are within 20% of the initial 100 km2 values, the shear is correlated with an R2 value of 0.72 (N = 2177). For just a 1-km separation of the cell centers and cell areas within 10% of the initial value, the squared correlation increases to 0.81 (N = 282). While these correlation values are significant, the correlations of divergence values are not (R2 = 0.16 and 0.04, respectively), leading to the conclusion that the computed divergence for a small region is highly dependent on the specific locations of the points used in the computation.
The analysis of the trajectory errors is formulated in terms of the difference in the computed displacement of the ice as determined from two comparison datasets. In the first comparison, the datasets were two RGPS products generated at JPL and ASF for May 1998. The input SAR images and software were identical for the two systems. The only differences were in the initial locations on the images of the displacement vectors and the personnel that provided guidance to the systems when manual intervention was required. If no manual intervention was required and the initial locations were identical, the computed displacements were identical from the two systems. With small separations of the initial locations, the two systems produced only slightly different results as long as no manual intervention was required. Manual intervention is needed when there is substantial local deformation or rotation. In these cases there was often a substantial difference in opinion between the operators of the two systems. The median difference in these cases was 1.030 km (Table 2). These differences reflect the difficulty in tracking ice in highly deformed areas. They also highlight the fact that the errors are not independent of deformation rates. We believe the errors are higher where the ice is deforming rapidly. The best estimate of the tracking error, 0.100 km (one pixel), is consistent with the observed 95th percentile of the displacement differences.
The second comparison is with buoy-determined displacements. In this comparison the vector correlation R2 value is 0.996 (N = 402). The rms difference is 1.431 km. If this value is used and an uncertainty in the buoy positions of 0.3 km is assumed, the uncertainty in the RGPS displacement is 0.919 km. However, the distribution of the displacement differences has a much heavier tail than would be expected for normally distributed position errors for the components, a tail attributable to occasional real ice deformation near the location of the buoys. If the median statistic, which is insensitive to outliers, is used instead of the rms difference, the displacement differences are small enough to infer that the buoy errors are less than 0.3 km; and if the buoy errors were assumed to be 0.100 km, the RGPS uncertainty would be 0.286 km. This is an upper bound for the RGPS displacement error including both geolocation and tracking errors because the buoy errors are likely larger than 0.100 km. The tracking errors for an individual trajectory accumulate over time, but the error in the displacements used for individual deformation calculations do not.
The displacement measurements obtained by the RGPS may be compared to those of two other satellite-based systems. The passive microwave Special Sensor Microwave Imager (SSM/I) has been used by a number of investigators to measure the ice motion (Agnew et al. 1997; Emery et al. 1997; Kwok et al. 1998; Kwok and Rothrock 1999; Lindsay 2002). Either the 85-GHz or the 37-GHz channels may be used with either 1- or 2-day intervals between images. The low spatial resolution of the images (12 and 25 km, respectively) and uncertainties in the timing of the image data contribute to errors of 5–12 km in the estimated velocity vectors (Kwok et al. 1998). The spatial coverage is usually quite good in the Arctic Ocean during winter, but tracking is not possible in the summer months due to loss of contrast of the wet ice surface. The Advanced Very High Resolution Radiometer (AVHRR) is also used to track sea ice (Emery et al. 1991). Either visible-band or thermal-band images may be used; however, both bands are subject to loss of surface tracking due to clouds. The resolution of the images is 1–4 km and a 1-day interval is commonly used. The accuracy of the tracking is about 1.6 km when compared to buoy motions (Lindsay 2002). The AVHRR tracking errors are more than 5 times as large as those of the RGPS and the SSM/I errors are more than 40 times as large.
Errors in the computed deformation, including area change, arise from relative tracking errors and from the approximation of the line integral around a region with a limited number of measurements of the displacement, what we call the boundary-definition error. The theory of the propagation of errors allows for an estimate of the error that arises just from tracking. If we assume the tracking error ɛf = 0.100 km, then ɛdiv = ɛshr = ɛvrt = 0.5% day−1. This is about half the standard deviation of the divergence computed from the cell strain rates (about 1% day−1, depending on the dataset and the screening for outliers) implying a low signal-to-noise ratio. However, the deformation in the ice is concentrated in quasi-linear features (leads or slip lines) and most of the cells experience very little deformation at any one time [Fig. 5 and also Stern et al. (1995) show pictures]. Those that do experience strong deformation might be thought of as having a larger signal-to-noise ratio, but the large deformation events are also subject to larger boundary-definition errors. The error in the area change due to tracking for the 10-km cells is 1.4%.
The number of points in the cell used to define the boundary and to determine the deformation also contributes to errors in the deformation estimates. The special SHEBA dataset with 5-km cells allows an estimate of this error by comparing the area changes for a 10-km cell found with either 4 or 8 points. The rms difference in the area changes is 3.4 km2. Monte Carlo simulations of the tracking error indicate that most of this difference arises from the boundary definition errors, which amount to 3.2 km2. This source of error is more than twice as large as the error from tracking alone (1.4 km2). The error variances from the two sources, tracking and boundary definition, are additive, so that the total error standard deviation is 3.5 km2 or 3.5%. This suggests that the Kwok and Cunningham (2002) estimate that the spurious ice production from area change errors is less than 1% of the total, based on tracking errors alone, is too low. The error in ice production is still quite small because of the large net divergence in the regional ice movement.
There are a very few outliers in the dataset that exhibit unrealistic behavior. For example, 50 observations (out of 2.7 × 106) have speeds over 50 km day−1 for 1 to 3 days, much greater than the maximum buoy speeds of about 30 km day−1. While it may be possible that some of these observations are valid, they should be screened on a case-by-case basis and compared to their neighboring trajectories. Similarly, a few of the deformation and area change values are unrealistic, often arising from highly distorted cells. Calculations of statistics that are sensitive to outliers should be performed only after screening the dataset. The only measure of a cell's distortion available in the dataset is the area, so a conservative screening for the area at the beginning and end of each observation interval is suggested.
The errors in the deformation and area change can be reduced by grouping cells into larger areas (Kwok and Cunningham 2002; Lindsay et al. 2003). The errors in the strain-rate invariants are reduced because ɛA [from (16)] grows more slowly than A in (17) as more cells are aggregated, and the boundary-definition error is reduced as well because of the larger number of points used to determine the boundary. The deformation from a set of cells is properly averaged by computing an area-weighted sum of the deformation-tensor components and then computing the invariants. If the cells are approximately the same size and form a large square, the error of the aggregate will be reduced by a factor of 1/N3/4, where N is the number of cells averaged together. The major difficulty in expanding the area analyzed is to ensure that the observation times of all the constituent cells are identical. For a few cells this is not too hard, but larger areas are increasingly difficult because the sizes of the regions analyzed from individual satellite passes are not always large.
The RGPS datasets are a remarkable accomplishment and offer a new and unique view into the kinematics of Arctic sea ice. Our validation efforts have shown that the ice tracking is in general excellent, at least as good as the tracking errors for buoys. Accurately determining the ice deformation at small scales is extremely demanding, but areal averaging improves the accuracy at the expense of spatial resolution. The deformation determined by the RGPS offers a new opportunity to investigate the nature of the ice deformation in response to wind and ocean forcing and to directly test the performance of models in simulating the deformation.
The buoy position data were obtained from I. Rigor and M. Ortmeyer of the International Arctic Buoy Program. The RGPS data were obtained from the Polar Remote Sensing Group of the Jet Propulsion Laboratory, Pasadena, California, and from the Alaska SAR Facility, Fairbanks, Alaska. We also acknowledge the reviewers for their helpful comments. This work has been supported by the NASA Cryospheric Sciences Program.
Corresponding author address: Ronald W. Lindsay, Polar Science Center, Applied Physics Lab, University of Washington, 1013 NE 40 St., Seattle, WA 98105. Email: email@example.com