GPS-derived velocity and strain fields around Dome Argus, Antarctica

Abstract

Knowledge of the surface topography, velocity field and strain field at an ice-core site is critical to the accurate interpretation of ice-core records. At Dome Argus (Dome A), where a Chinese deep ice-core drilling project is being carried out, we have produced an accurate surface digital elevation model from GPS measurements in January 2013 at 47 sites. We identify two peaks at Dome A, with the northern peak ~7 cm higher than the southern peak. Repeat GPS measurements at 12 sites in 2008 and 2013 provide a surface velocity field around the dome. The surface velocity ranges from 3.1±2.6 to 29.4±1.2 cm a^–1, with a mean of 11.1 ~2.4 cm a^–1. The surface flow directions are near perpendicular to the surface elevation contours. Velocities from GPS are lower than derived from satellite radar interferometry (InSAR). From GPS velocities, the accuracy of velocity from the existing InSAR velocity field is determined, resulting in a standard deviation of 0.570 m a^–1 in speed and 117.5º in direction. This result is consistent with the reported accuracy of InSAR, showing the value of in situ GPS measurements for assessing and correcting remote-sensing results. A surface strain field for the drilling site over Dome A is calculated from 24 strain triangles, showing north–south extension, east– west compression and vertical layer thinning.

1. Introduction

The ongoing deep ice-core drilling at Dome Argus (DomeA) is vital to the understanding of climate change in Antarctica. The ice-core record may extend the paleo-climatic record to ~1 Ma BP (Xiao and others, 2008). Information about ice thickness, ice velocity, surface topography, snow accumulation history and strain field is needed in order to assess the age of the deep ice here. In 2008, the ice-sheet thickness and bed topography were measured by ground-based, ice-penetrating radar during the 24th Chinese National Antarctic Research Expedition (CHINARE; Sun and others, 2009), exhibiting a 3139 m thickness at the core location (Cui and others, 2010). The meteorological characteristics have been extracted from an automatic weather station (AWS) (Xiao and others, 2008; Ma and others, 2010). The surface mass balance (SMB)/snow accumulation rate was estimated using snow-pit and stake-array observations (Hou and others, 2007; Ding and others, 2011). Steep bedrock valley walls beneath Dome A were reported, and freezing from the base has important impacts on ice thickness and ice flow over Dome A (Bell and others, 2011). The surface topography of Dome A was GPS-surveyed in January 2004 (Zhang and others, 2007) and January 2007 (Cheng and others, 2009).

Despite the above efforts, there is little progress in obtaining a reliable velocity field near Dome A, making it difficult to model the mass balance and bed ice age here. Satellite radar interferometry (InSAR) has been used to estimate ice flows over the Antarctic ice sheet (Rignot and others, 2011), but the precision of InSAR here is affected by factors such as spatial–temporal correlations of SAR images over ice surfaces and snow accumulation effects. It has been reported (Mouginot and others, 2012) that InSAR-derived velocities may have errors in magnitude ranging from 1 m a^–1 in the interior regions to 17 m a^–1 in coastal sectors, and InSAR-derived flow direction may have errors ranging from <0.5⁰ in areas of fast flow to unconstrained direction in sectors of slow motion. The resolution is only ~450 m, and this can induce the challenge of geolocating the SAR images over remote regions of the ice sheet. Comparison of surface ice velocity and measured InSAR velocity over dome regions in Antarctica is rare. In general, ice velocities over the dome regions with smaller slopes (Vittuari and others, 2004; Wesche and others, 2007; Cheng and others, 2009) are low, requiring very precise field measurements to detect them. Differential GPS provides high-precision coordinates and is a useful tool for measuring surface topography in the Antarctic interior.

The key objective of this paper is to present GPS data collected during expeditions to Dome A in 2008 and 2013. The main results will be a digital elevation model (DEM), a velocity field and a strain-rate field. GPS results are compared with the InSAR result of Rignot and others (2011), and will be important for interpreting the deformation in the upper part of the ice sheet over Dome A.

2. Field Gps Survey and Data Processing

2.1. GPS data collection

In January 2008, 49 poles were set up to estimate SMB with 5 km × 5 km resolution, and to map the ice surface topography in a 30 km × 30 km area centered at Kunlun station (Fig. 1). Two Leica 1230 dual-frequency GPS receivers were used to collect data at 12 of the 49 poles, while the remaining poles were not occupied due to lack of time. A reference GPS station was established near Kunlun station, which also served as the base camp. The reference GPS receiver made continuous measurements from 13 to 26 January 2008 at 5 s intervals. Static GPS measurements were carried out at 12 poles for at least 20 min at 5 s intervals (Cheng and others, 2009) to estimate high-precision coordinates.

Fig. 1. Sketch map showing the locations of poles for GPS measurements over Dome A in 2013; the spatial resolution is 5 km.

In January 2013, GPS measurements were made at 49 poles installed in 2008 with the two Leica GS15 dual-frequency GPS receivers. Thus, there are repeat GPS measurements at the 12 poles occupied in both 2008 and 2013. The GPS data collected in 2013 enabled us to map the surface topography in a larger area than in 2008, and the repeat measurements at the 12 poles can be used to construct a surface velocity field. In the 2013 GPS survey, five reference GPS stations were set up with >~5 hours available for differential GPS (triangles in Fig. 1) with the other GPS measurements. These five reference GPS receivers collected continuous data over 7–9 and 12–14 January 2013 at 1 s intervals. The other GPS receivers were used as roving receivers. The measurements at each pole last at least 5 min. The slant height and vertical height of a pole were measured to determine the tilt angle, which was used to correct for height error. Two poles, P14 and P17, were missing, so in fact we collected GPS data at only 47 poles in 2013. Besides the five reference GPS stations, 30 of the 42 roving poles were occupied for >10 min. The coordinates at the 47 poles are listed in Table 1 and are used to construct an elevation model around Dome A (in comparison to that of 2008). The repeat GPS data at the 12 poles were used to construct a velocity field and a strain field around Dome A.

Table 1. GPS sites, the estimated coordinates and formal errors (68% confidence interval) at Dome A in 2013 and 2008, where σ and σ′ are for 2013 and 2008, respectively

Fig. 2. (a) Surface topography, and surface velocity field from GPS and InSAR over Dome A; large and small red rectangles show the 400 km² and 64 km² areas surveyed by Cheng and others (2009) and Zhang and others (2007), respectively. (b) Surface speeds from InSAR by Rignot and others (2011).

2.2. Data processing

2.2.1. Analysis of absolute motion

The minimum session length for the GPS data is 296 min for the five reference stations in 2013 and the reference station in 2008. The GPS data from these six reference stations were processed using the Version 10.50 GAMIT/GLOBK software (King, 2002), developed at the Massachusetts Institute of Technology, USA. Seven International GNSS (global navigation satellite systems) Service (IGS) stations around Antarctica, namely MAW1, CAS1, DAV1, DUM1, MCM4, OHI2 and SYOG, were used in the processing. During data processing, the following options were selected: (1) IGS SP3 precise ephemerides were used; (2) the IGS stations were tightly constrained (within 1 cm) to their International Terrestrial Reference Frame (ITRF) 2008 values, while the six reference stations at Dome A were loosely constrained (within 100 m); (3) an elevation cut-off angle of 15º was adopted; (4) antenna-phase center variation corrections were applied; (5) the ionospheric-free linear combination of the L1 and L2 frequencies was used; (6) corrections due the solid-Earth and frequency-dependent tides were applied;(7) the dry component of the zenith tropospheric delay was implemented using the Saastamoinen model with global mapping function; and (8) the wet component of the zenith tropospheric delay was estimated every 2 hours with global mapping function. The GAMIT solutions were then combined using the GLOBK software.

To account for the effect of tectonic motion, we corrected for the east–west component and north–south component of velocities at MAW1, CAS1, DAV1, DUM1, MCM4, OHI2 and SYOG using the result of Jiang and others (2009). The maximum absolute values of the east–west and north–south components of velocities at these seven stations are 1.5 and 1.3 cm a^–1, respectively. However, there is no rock at Dome A, and no long-term GPS data are collected near Dome A. Therefore, the plate motion cannot be obtained directly. To reduce the impact of plate motion, we fix the coordinates of seven IGS stations to the same epoch 2013.0 for both 2008 and 2013 for data processing with GAMIT.

The final geodetic coordinates and the estimated uncertainties of the five reference stations in 2013 are shown in Table 1. The maximum estimated uncertainties for the five stations in the X, Y and Z directions are 0.7, 0.8 and 2.0 cm, respectively. The mean uncertainties for the five reference stations in the X, Y and Z directions are 0.5, 0.6 and 1.8 cm, respectively. The WGS84 coordinates of the reference station in 2008 are 80º25^ʹ02.240332^˝ S, 77º0632.26741˝ E, 4091.2319 m, with the estimated uncertainties being 0.1,0. 1 and 0.4 cm, respectively.

2.2.2. Analysis of relative motion

The GPS data at the 12 poles (considered as roving stations) in 2008 were reprocessed to compute the site coordinates relative to the reference station. Version 2.50 Trimble Business Center (TBC) was used to obtain the baseline between a roving station and a base station. During data processing, IGS SP3 precise ephemerides were used, and an elevation cut-off angle of 15º was adopted. The 12 baseline solutions were all fixed. Table 1 shows the formal errors of the horizontal coordinates at the 12 poles, which were computed using the uncertainties of the baseline solutions and the uncertainty of the base station. The mean formal error at the 12 poles for both X and Y is 1.3 cm, the maxima reaching 3.0 and 3.5 cm.

Besides the five reference stations, static GPS data were collected for the other 42 poles in 2013. TBC was used to obtain the baseline solutions between the roving stations and the five base stations. The fixed solutions were obtained for all 42 baselines, and the formal errors of 42 poles (considering the uncertainty of tilting) are shown in Table 1. The mean formal errors of 42 poles in the X, Y and Z directions are 1.0, 1.0 and 3.5 cm, respectively, with the maxima 2.5, 3.0 and 8.6 cm, respectively. Generally, more accurate results are obtained with longer durations. As a rule (Leick, 2004), the vertical components from GPS have the least accuracy compared to other components, because all visible GPS satellites are above the horizon at a given receiver site. The use of a cut-off angle further decreases the vertical accuracy.

3. Results and Discussion

3.1. Surface topography around Dome A

The derived surface topography in the study area (Fig. 2a) was produced using the GPS-derived heights (Table 1) at the 47 poles in Figure 1. We interpolated the 47 heights on a grid using the minimum-curvature algorithm. Elevation contours around Dome A were plotted at a 2 m interval (Fig. 2a). The dome is saddle-shaped and flat. The range of heights within the 1156 km² area is <23 m, and most of the area is above 4082 m. There exist two peaks in the study region, which may be due to the existence of bedrock mountains and steep valley walls in the region shown in Cui and others (2010, Fig. 4b) and Bell and others (2011). Our result shows that the northern peak is ~7 cm (with a mean uncertainty of 3.5 cm) higher than the southern peak (4092.11 m vs 4092.04 m in Fig. 2a). The two red rectangles in Figure 2a show the areas surveyed by Zhang and others (2007) in early 2005 and by Cheng and others (2009) in early 2008.

Cheng and others (2009) showed that the southern peak was ~30 cm higher than the northern peak. However, the mean elevation formal error from Cheng and others (2009) is ~20 cm, much larger than our formal error of 3.5 cm. This smaller uncertainty is attributed to the improved surveying technique here over that used by Cheng and others (2009). Moreover, the SMB and the 20 cm elevation uncertainties of Cheng and others (2009) can both lead to the 30 cm elevation difference in 2008. To investigate this 30 cm difference, the SMB difference between the southern and northern peaks from 2008 to 2013 is computed, and is ~ 8 cm (43 cm vs 51 cm). With this 8 cm SMB (2008–13) and ~7 cm elevation difference between the northern and southern peaks in 2013, the elevation difference discrepancy between the southern and northern peaks in 2008 is estimated to be 1 cm, as compared with 30 cm by Cheng and others (2009). Hence, the 20 cm elevation formal error of Cheng and others (2009) is the main reason for the elevation difference between the southern and northern peaks in 2008.

3.2. Surface velocity field around Dome A

With the repeated GPS measurements at the 12 poles, the surface velocity field and the corresponding uncertainties were calculated and are listed in Table 2. The mean velocity over Dome A is ~11.1 ±2.5 cm a^–1, with the maximum velocity reaching ~29.4 ±1.2 cm a^–1 at P01 in the northwest corner of the survey grid. The minimum surface velocity is 3.1 ±2.6 cm a^–1 at P26, the closest site to the northern peak (Fig. 2a). This radial ice flow over Dome A is different in character than that at the EDML drilling site (Wesche and others, 2007), which is parallel to the ice divide. The flow directions at these sites are consistent with the downslope motion of the ice sheet (perpendicular to the elevation contours). The velocity over Dome A, which is near the summit of the East Antarctic ice sheet, is higher than that over Dome C because the latter has much gentler slopes than the summit (Vittuari and others, 2004), 2.8% vs 12.0%.

Table 2. GPS velocity components, speed, direction and their uncertainties (obtained by applying a factor 5 to the formal errors), and InSAR velocity components, speed and direction

A further investigation was made to show the relationship between ice surface velocity and ice surface topography. The surface velocity field is listed in Table 2. The slope amplitude and direction at the same GPS sites are derived from the surface topography in Section 3.1 as follows. (1) The east– west component and north–south component grid of surface topography gradient are calculated. (2) The east–west component and north–south component of surface topography gradients at 12 sites are interpolated from the two grids, respectively. (3) The slope gradient amplitude and direction at the 12 sites are calculated. Velocity vs slope, and velocity direction vs slope direction are compared. The correlation coefficient between velocity and slope is 0.766, while a linear equation between them is shown in Figure 3a. A higher correlation coefficient, up to 0.926, between slope direction and velocity direction is found in Figure 3b, suggesting the direction of the ice movement correlates well with the direction of maximum surface gradient. The comparison indicates that the surface velocity is largely determined by the surface topography. However, the spatial scale may influence the correlation coefficients. Higher correlation coefficients may be obtained with a better spatial scale.

Fig. 3. The surface slope field over Dome A derived from the ice surface topography in Figure 2a. (a) The relationship between surface velocity field and surface topography through velocity vs slope gradient amplitude; (b) velocity direction vs slope direction.

As shown by Mouginot and others (2012, Fig. 3), the errors in surface velocity and direction from InSAR are up to5.5 m a^–1 and 1808, respectively, as determined from cross-calibration of InSAR satellite measurements. Because of a lack of ground data, there are few assessments of InSAR velocity over regions such as Dome A. We assessed the accuracy of the surface velocity field from InSAR over Dome A using our GPS-derived surface velocity field. The InSAR surface velocity field in the experiment region was from Rignot and others (2011) with 450 m resolution. The velocity field from InSAR over Dome A is shown in Figure 2b. The mean velocity from InSAR is 1.425 m a^–1. The standarddeviation of velocity from the InSAR results is 1.438 m a^–1, indicating that the velocity over Dome A is spatially unevenly distributed. The surface velocities at 12 sites were interpolated from InSAR velocities, and are shown in Figure 2a.

Significant differences exist between the GPS and the InSAR results. The InSAR velocities always exceed those from GPS, with a mean difference ~1.075 m a^–1 and a standard deviation of 0.570 m a^–1. They are 2.0–31.3 times larger than the GPS velocities, and on average 14.3 times larger. At 10 of the 12 poles, they are nine times larger. The directions from GPS and InSAR are not consistent, and the mean and the standard deviation of the direction differences are 33.88º and 117.5º, respectively. We found consistent directions from GPS and InSAR at four poles, while the two directions are opposite at three poles. We attribute the discrepancies between InSAR and GPS to error in InSAR, DEM error, geolocation of InSAR, and interpolation errors. Our GPS assessment suggests that one must be cautious about the large uncertainties in surface velocity from InSAR, especially near major ice divides.

3.3. Surface strain field around Dome A

A surface strain field over Dome A was determined using the velocities at the 12 poles. Using this surface velocity field, with x as the northward component and y as the eastward component, we numerically determined the strain rate of two neighborhood points as (Cuffey and Paterson, 2010)

(1)

where ∆v_x and ∆v_y are the differences of the x and y velocity components in the pair of GPS survey points, and ∆x and ∆y are the corresponding distances in the x and y components. Considering the practical GPS positioning errors, we multiplied the formal errors by a factor 5 when estimating the uncertainties of the strain rates in Tables 3 and 4.

Table 3. Velocity differences, surface strain rates and their uncertainties for pairs of survey points

Table 4. Directions θ, surface strain rates έ and their uncertainties for the strain triangles

To compute the surface strain rate, we divide the area into six rectangles shown in Figure 4. We formed four triangles in each rectangle. Therefore, there are 24 triangles for the six rectangles. From Eqn (1), the strain rates were calculated for all possible pairs of neighboring points, yielding 29 values (Table 3). The distances among pairs range from 10 087.78 to 18 536.04 m, while the distances in x and y components range from 200.88 to 15 777.40 m, and from 332.27 to 10 934.63 m, respectively. In the 29 pairs, 14 of the resulting strains have errors less than 5.0 × 10^–6 a^–1 in both x and y components. Most strain rates in x and y components for the EDML drilling site (Wesche and others, 2007) exceed1.0 × 10^–5 a^–1, and the larger distances over Dome A may lead to these smaller surface strain rates.

Fig. 4. The investigation region over Dome A is divided into six rectangles, numbered 1–6. For each rectangle, four triangles are formed with different combinations of points. For example, the four triangles in rectangle 1 are P01/P22/P24, P03/P22/P24, P01/P03/P22 and P01/P03/P24. Surface strain rate and mean strain vectors of each triangle are calculated. Strain vectors s are plotted for all 24 triangles. The elevation contour is 2 m.

We calculated the average strain of each triangle, yielding 24 values (Table 4). The positive strain rates in x

Using the average strain rates of each triangle, we calculated the rotation of the principal components as

(2)

where the angle θ corresponds to the axis of maximum strain rate , and the angle θ + 90⁰ corresponds to the axis of minimum strain . The calculation was repeated for each strain triangle, and the result shows that the directions of maximum strain rate vary from 60^º to 180^º (Table 4).

With the average strain rates of each triangle, the strain-rate magnitudes for the two directions are calculated from

(3)

The calculation is repeated for every strain triangle. The strain-rate magnitudes in the direction of the maximum strain rate vary from (1.6 ± 2.8) × 10^–5 to (10.4 ± 1.4) × 10^–5 a^–1. The strain-rate magnitude in the direction of the minimum strain rate can be up to half that in the direction of the maximum strain.

Using the incompressibility condition, we calculated the flow-induced vertical strain rate for the 24 strain triangles using (Cuffey and Paterson, 2010)

(4)

The result shows that the vertical strain rate ranges from (–5.91 ± 1.44) × 10^–5 to (6.89 ± 3.38) × 10^–5 a^–1. We obtained six positive triangles in the southeast part of Dome A.

Among the 24 triangles, the mean elevation of three triangles, 5 (P05/P24/P26), 13 (P24/P43/P45) and 18 (P24/P26/P45), exceeds 4089 m. These three triangles lie in an area with least terrain variation in Dome A. The directions of the maximum strain rate from 13 and 18 are both close to 08, hence the maximum strain rate is mostly equal to strain rate in the x component. The average strain rates of the three triangles in the x, y and combined directions are(6.37 ± 1.67) × 10^–5, (–1.93 ± 1.34) × 10^–5 and (–2.48 ± 1.39) × 10^–6 a^–1, respectively. The average flow-induced vertical strain rate έ _z of the three triangles is (–4.45 ± 2.16) × 10^–5 a^–1. Furthermore, the direction of the maximum strain rate is 176.95⁰, with the maximum strain rate(6.43 × 2.31) × 10^–5 a^–1.

Kunlun station, where the drilling site is located, is within triangles 18 and 19 (P24/P26/P47). To estimate a strain rate representative of the drilling site, we calculated the average strain rate of triangles 18 and 19. The estimated strain rates in the x, y and combined directions are (5.40 ± 1.51) × 10^–5, (–3.16 ± 1.99) × 10^–5 and (–0.34 ± 1.36) × 10^–5 a^–1, respectively. The estimated flow-induced vertical strain rate έ _z is(–2.24 ± 1.36) × 10^–5 a^–1. The direction of the maximum strain rate is 177.878, with the maximum strain rate(5.44 ± 1.68) × 10^–5 a^–1.

4. Discussion and Conclusion

The surface velocity field and surface strain field around an ice-core site are critical to the interpretation of the ice-core records. Although remote-sensing tools such as InSAR may estimate the velocity field around a desired ice-core site, their precision is affected by factors such as snow accumulation. They require validation by in situ observations. As the surface velocity around Dome A is low and is hard to detect by InSAR, our GPS result is indispensable in assessing the InSAR result.

With the GPS measurements carried out at 47 sites in 2013, a DEM around Dome A was obtained. Since our GPS-derived elevations have cm-level accuracies that are far better than the 20 cm accuracy of Cheng and others (2009), the DEM presented here is more convincing.

The surface velocity field over Dome A was derived using GPS measurements at 12 sites in both 2008 and 2013. The directions of ice flow are radial and conform to the downslope motion of the ice sheet. A comparison is made between surface velocity field and surface slope. The result shows that the correlation coefficient can be up to 0.93, which indicated that the surface velocity field is mainly determined by surface topography. When assuming surface velocity as balance velocity, this analysis provides constraints for ice thickness, surface slope and snow accumulation data, and provides an important insight into the probable flow pattern of the ice sheet.

With the GPS surface velocity field, the estimated accuracy of InSAR surface velocity in speed and direction is about 0.57 m a^–1 and 117.5^º, respectively. The labeled accuracy of the InSAR surface velocity fields is usually assumed to exceed 1 m a^–1 (Rignot and others, 2011), so the actual accuracy is better than the labeled accuracy. As the maximum surface speed is ~ 0.29 m a^–1, the errors in speeds are comparable to the speeds near ice divides. In this case, the error in flow direction is unconstrained, as the error increases as the velocity decreases (Rignot and others, 2011). Therefore, the labeled accuracy in the InSAR surface velocity field is reasonable.

With the surface velocity field from GPS measurements, the surface strain fields are calculated from 24 triangles. The overall direction of the maximum strain rate varies from 60^ºto 180^º. The flow-induced vertical strain rate in the southeast of Dome A is positive, while the negative s are shown in the other parts of Dome A.

To estimate a strain rate representative of the drilling site, we calculate the average strain rate of two triangles where the drilling site is located. Similar strain rates are shown for the two triangles. The estimated strain rate in x, y and combined directions is (5.40 ± 1.51) × 10^–5, (–3.16 ± 1.99) × 10^–5 and (–0.34 ± 1.36) ×10 ^–5 a^–1, respectively. Furthermore, the direction of the maximum strain rate is 177.878, with the maximum strain rate (5.44 ± 1.68) × 10^–5 a^–1. Hence, the strain rate at the drilling site is similar to that in flat regions of Dome A. The dilatational force and the compressing force acts in the direction of 177.87^º and 87.87º as the maximum and minimum components, respectively. The estimated flow-induced vertical strain rate of the drilling site is (–2.24 + 1.36) × 10^–5 a^–1, which also indicates layer thinning is required in the vertical component to achieve balance.

The GPS-derived surface velocity field at the 12 sites provides a critical quality control for InSAR results in Antarctica, especially over dome regions. The surface strain field is characterized in layer thickening, which has to be accounted for to yield a correct interpretation of ice-core data. In order to capture correct signatures of ice dynamics and mass balance in an area such as Dome A, it is highly important to combine various data types such as surface velocity, surface topography, ice thickness and snow accumulation. With surface velocity, the accumulation, ice thickness and density profile, it is possible to give an indication of the dome’s long-term stability.