57-Year Ice Velocity Dynamics in Byrd Glacier Based on Multisource Remote Sensing Data

Long time-series glacier ice velocity reflects the local climate changes and can be used to estimate mass balance (MB) changes, which is a critical parameter for understanding of glacier–climate interactions and prediction of sea level rise. However, due to the difficulty in image matching caused by the poor quality of historical satellite images, there was insufficient glacier ice velocity data before 1999. Here, we proposed a multiple-constraint dense image matching approach for mapping historical ice velocity based on the early poor-quality images from ARGON, Landsat-1, and Landsat-4/5. We successfully applied this method to Byrd Glacier to generate its historical ice velocity maps from 1963 to 1999. Additionally, ice velocity maps of Byrd Glacier from 2000 to 2014 were generated by IMCORR software using Landsat-7 and Landsat-8 images. Combining with the ice velocity maps from the Global Land Ice Velocity Extraction from Landsat-8 dataset since 2014, we obtained the ice velocity of Byrd Glacier for 57 years. Our results showed that the glacier experienced slight fluctuations in ice velocity, which may not be due to the calving events in the studied portion of the Ross Ice Shelf or the air temperature changes, but by the activity of subglacial drainage systems. Furthermore, Byrd Glacier showed a positive MB (average rate of 2.6 ± 2.0 Gt/year) from 1963 to 2020, indicating that global climate change may have a limited impact on it.


I. INTRODUCTION
O WING to global climate change, several large glaciers in Antarctica (e.g., Pine Island Glacier, Thwaites Glacier, Totten Glacier, and Cook Glacier) are changing rapidly by the complex interplay of external (e.g., atmospheric warming, intrusion of warm ocean water, and the presence of sea ice or ice mélange in front of floating ice shelves) and internal (e.g., bedrock and glacier topography, calving events, and subglacial lake drainage events) forces [1], [2], [3], [4]. The ice velocity of glaciers is not only an important parameter for estimating the mass balance (MB) of the Antarctic Ice Sheet (AIS) using the input-output method but can also be used to monitor the spatiotemporal ice dynamics of glaciers and increase the understanding of the driving forces behind ice velocity change [1], [5], [6], [7]. Unfortunately, several glaciers still lack sufficient velocity analyses, especially in East Antarctica [1], and most studies on glacier ice flow have focused only from the 1990s to the present [8], resulting in an incomplete understanding of long-term glacier dynamics. Ice velocity has been widely studied using optical or synthetic aperture radar (SAR) satellite remote sensing techniques [9], [10], [11], [12]. Owing to their active data acquisition during austral winter and high precision in ice flow detection, SAR-based techniques have played an important role in the understanding of interannual and seasonal ice velocity since the launch of ERS-1/2 satellites in the 1990s [1], [10]. The first Antarctic ice velocity measurements using optical remote sensing were acquired at discrete points in Byrd Glacier by manually tracking glacier surface features on Landsat-1 and Landsat-4 images [13]. With the improved spatial resolution and geolocation performance of available Landsat-7 and ASTER images, several reliable image matching methods and software, both in the frequency and space domains, have been developed to map the ice velocity field, for example, IMCORR and COSI-CORR. These methods have been widely applied in ice velocity studies of individual glaciers and small glacierized regions [7], [14]. With the launch of the Landsat-8 and Sentinel-2 satellites, it has become possible to map the ice velocities for the majority of the AIS [6], [9]. Thousands of Landsat-8 images have been collected since 2013, and Antarctic-wide annual or seasonal ice velocity datasets have been successfully generated [15], [16]. For example, the Global Land Ice Velocity Extraction from Landsat 8 (GoLIVE) dataset [17] has provided near-real-time ice velocity maps since 2013, and the Landsat 8 Ice Speed of Antarctica (LISA) dataset [18] has provided annual or multiyear composite ice velocity maps from 2013 to 2017. Additionally, the Inter-Mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) [19] is producing annual ice velocity maps since 1985; however, the coverage before 2013 was incomplete. Since 2000, ice velocities can be directly obtained from existing products or generated by common methods and software based on Landsat or Sentinel images. Historical satellite images (e.g., Landsat-1 and Landsat-4) with poor geolocation and internal geometric distortion were used for ice velocity studies in Antarctica before 1999. However, these images have now been discarded because of the need for extensive preprocessing and complex matching techniques [11]. In this article, we propose an image processing approach to improve the utilization of these historical satellite images to fill the data gap in historical ice velocity studies. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ As the largest ice shelf in Antarctica, the Ross Ice Shelf (RIS) is fed by large glaciers and ice streams flowing from both East Antarctic Ice Sheet (EAIS) and West Antarctic Ice Sheet (WAIS). Recent research shows that local, seasonal production of warm upper-ocean water near the ice front of the RIS drives rapid ice shelf melting east of Ross Island, where thinning would reduce buttressing of nearby EAIS glaciers as well as the more distant WAIS ice streams [20]. As the fastest ice stream entering the RIS in EAIS [21], prolonged changes in the flow dynamics of Byrd Glacier will impact the stability of the RIS and MB of EAIS [22].
Ice flow of Byrd Glacier is being studied for the past ∼60 years. Based on the earliest investigations carried out in 1960-1962, Swithinbank [23] reported ice velocities along a line of in situ survey markers across the glacier, and Alder [24] measured ice velocities at several points of the glacier center through aerial triangulation from sets of almost annual (348 days) aerial photographs. Brecher [25] obtained two sets of aerial photographs over a 56-day interval in 1978-1979, and measured 601 common moving points to generate the ice velocity of Byrd Glacier. However, due to the inaccessible or hazardous field environments in Antarctica, the logistical costs would be expensive even if in situ data collection is possible. Lucchitta et al. [13] first demonstrated that satellite observation is a valuable tool for ice velocity studies by manually measuring the displacements of surface features that remain preserved over long periods (several years). Subsequently, Stearns et al. [26] observed an acceleration of the glacier flow, whose onset coincided with the water discharge from two large subglacial lakes, suggesting that an active lake drainage system can cause large and rapid changes in glacier dynamics. Despite the abovementioned studies on ice flow velocities of Byrd Glacier, there was little continuous velocity data before 1999, resulting in a poor understanding of many key details of the long time-series dynamics of the glacier.
In this article, we proposed a multiple-constraint image dense matching approach for the production of ice velocity maps based on poor-quality historical remote sensing images collected before 1999. We successfully applied this approach to Byrd Glacier to produce ice velocity maps between 1963 and 1999. Finally, combined with the ice velocity maps generated after 2000 from the common ice velocity mapping software (e.g., IMCORR) and products (e.g., GoLIVE), we analyzed the long time-series of ice velocity changes of Byrd Glacier from 1963 to 2020 and compared them with the dynamics of the RIS (the RIS is fed by Byrd Glacier), the variability in environmental factors (e.g., air temperature near glacier surface), and the subglacial drainage activities.

A. Study Site
Byrd Glacier (80 • S, 160 • E), located in a 140-km long and 25-km wide fjord through the Transantarctic Mountains, has the largest catchment basin (ca. 1 101 725 km 2 ) in Antarctica and annually delivers a large ice mass into the RIS (ca. 22 Gt [5]), making it an influence factor to the MB of EAIS and resulting sea level rise [26]. It is the fastest ice stream entering the RIS with sustained rapid flow exceeding 850 m/year around its grounding line. The location of Byrd Glacier, including its catchment basin from the updated drainage inventory [5], the ice velocity map [27], and the grounding line derived from interferometric synthetic aperture radar [28], are shown in Fig. 1.

B. Data
Owing to the different time periods of satellite services, multisource satellite data were employed to study the changes in surface ice velocities of Byrd Glacier during 1963-2020 (Fig. 2). The earliest satellite images used in this article were collected from the CORONA program, the first-generation film-based U.S. reconnaissance satellite series (KH series, including three missions, namely, CORONA, ARGON, and LANYARD) from 1959 to 1972 [29], [30]. Among them, only the ARGON mission (KH-5) successfully captured images from Antarctica between 1961 and 1964 [31]. In this article, we used the orthorectified ARGON mosaic around the Antarctic coastline, produced by Kim [32], with least-squares-based bundle block adjustment. The geolocation error of this mosaic was within 1 pixel (140 m), providing a good opportunity for historical ice velocity mapping of Byrd Glacier. Other historical satellite images employed in this article include images from Landsat-1 Multi-Spectral Scanner (MSS) sensor band 4 (60-m resolution) taken in 1973, Landsat-4 MSS band 7 (60-m resolution) taken in 1983, and Thematic Mapper (TM) band 4 (30-m resolution) taken in 1989. After the 1990s, panchromatic images with a 15-m resolution acquired by the Landsat-7 carrying Enhanced Thematic Mapper Plus and Landsat-8 carrying Operational Land Imager sensors were carefully selected. All Landsat-1/4 images were processed at the systematic correction (L1GS) level, which was created when the location accuracy was not sufficient to apply terrain correction, and Landsat-7/8 images were processed at the systematic terrain correction (L1GT) level, which was created when location accuracy was sufficient to permit georegistration using a terrain model (www.usgs.gov/landsat-missions/landsatlevels-processing). After 2014, instead of calculating the ice velocity of Byrd Glacier, we adopted the existing maps in the GoLIVE dataset from Landsat-8 image pairs [17]. The Landsat-7 images collected after the failure of scan line corrector (SLC-off images) were destriped using the Landsat gap-fill plug-in in the ENVI software developed by the National Aeronautics and Space Administration [33]. The local linear histogram matching method was used to fill each SLC-off image based on a linear transformation of the near-time acquired images. The triangulation interpolation method was used in the absence of near-time images. In case of different resolutions of an image pair, the higher resolution image was downsampled to the resolution of the low-resolution image using a Gaussian kernel from the OpenCV library in Python.
To assess the influence of local climate on the ice velocity of Byrd Glacier, we obtained a time-series of climate data recorded by the Marilyn automatic weather station (AWS), located on the RIS at the closest distance of ∼140 km from the grounding line of the glacier [see Fig. 1(a)]. It is part of the Antarctic AWS  Project funded by the National Science Foundation Office of Polar Programs. Data collection began in 1991 with an update frequency of every 10 min or 3 h. Five meteorological data variables (temperature, pressure, wind direction, relative humidity, and vertical temperature difference) were available in tabular ASCII format via the University of Wisconsin's AWS Project website (http://uwamrc.ssec.wisc.edu/aws/).

A. Ice Flow Velocity
For optical satellite images, the ice velocity is derived by matching the features on two coregistered images to measure their displacements over time. The ice velocity can then be calculated using the following equation: where v is the annual ice velocity (m/year), v x and v y are the annual ice velocities in two perpendicular directions of the image (m/year), d i and d j are the displacements in the two directions in pixels, r is the image resolution (m), and Δt is the time interval between the two images (year). The main step in ice velocity mapping is the image matching of the corresponding surface ice features on two images of the same area collected at different times. Image quality is an important factor in this process for a favorable result. The low spatial and radiometric resolutions and sparse temporal coverage of the early satellites (before 1999) (ARGON, Landsat-1, and Landsat-4) resulted in significant difficulties for surface texture recognition and feature tracking using common software, such as IMCORR and COSI-CORR (Fig. 7). We proposed a multipleconstraint dense image matching approach for these poor-quality images by integrating feature-and grid-based matching. In the following parts of this section, we introduce the image preprocessing and approach of the proposed method for poor-quality image matching.
For the satellite images collected after 2000 (e.g., Landsat-7/8), an experiment was conducted (Section IV-A-2) to compare the performance of the proposed method with that of IMCORR.
Both methods were successful in producing ice velocity maps, indicating that the quality of Landsat-7/8 images was sufficient to produce massive and reliable matches using publicly available software. However, because of the time-consuming seed point selection and mismatch elimination of the proposed method, IMCORR software [34] was employed owing to its high quality and data processing efficiency.
1) Image Preprocessing: Image preprocessing is a prerequisite for image matching due to the weak contrast and low geolocation accuracy of historical satellite images, and includes image enhancement and geolocation.
Most of the Antarctic surface is covered by snow and ice, resulting in weak-textured remote sensing images. We used a contrast-limited adaptive histogram equalization algorithm to extend the original grayscale range and suppress noise for image enhancement. Geolocation is necessary to remove the positioning error of the images to ensure accurate georeferencing. For the poor-quality Landsat images before 1999, which were usually L1GS products without terrain correction and had large positioning errors of up to dozens of kilometers [35], we manually selected control points on the stationary features (e.g., bare rocks) from the Landsat Image Mosaic of Antarctica product [36] for horizontal control and RAMP DEM Version 2 [37] for vertical control, with the assistance of the orthorectification module in the PCI Geomatic software for georectification [38], [39]. We manually checked the geolocation accuracy after georegistration by randomly selecting five inspection points. If the geolocation error of one of the inspection points was larger than one pixel, the control points for the image geolocation were manually adjusted. After orthorectification, the average geolocation error of Landsat images before 1999 was 0.74 pixels. For the high-quality images after 2000, which were terrain-corrected L1GT products with high location accuracy, all Landsat-7/8 images in each pair were coregistered by manually selected tie points over stationary regions.
2) Multiple-Constraint Dense Image Matching for Poor-Quality Imagery Before 1999: The long study period (∼40 years, from 1963 to 1999) and the scarcity of early historical satellite images (only six in total) resulted in a very long-time interval of the two images to be matched for ice velocity estimates (mostly up to 10 years), e.g., 1963-1974, 1974-1983, and 1989-1999. During this long time period, the surface flow features traveled varying distances depending on the different ice velocity of the regions, causing significant changes in their structures, or even disappearances in case of large flow velocity. This led to difficulties in setting the search window for image matching. Additionally, the poor quality and low contrast of the surface textures on the optical images further reduce the feasibility of existing image matching algorithms. To acquire reliable matching, we modified the method of multipleconstraint-based robust matching [40], which was proposed for landslide monitoring using poor-texture images. One of the most important considerations of this method is to reduce the search range of feature tracking to improve matching reliability. Dense image matching was performed using geometric and similarity constraints. For geometric constraints, we constructed Delaunay triangulated irregular networks (TINs) using seed points to initially describe the spatial patterns of ice flow. Then, featurebased image matching was performed by assuming consistent ice flow within the TINs, where we could restrict the search area for a matching point within a reasonable range using affine transformation. Based on this, grid-based image matching was employed to ensure high spatial resolution of the matches. In these two steps, a similarity constraint based on normalized cross-correlation (NCC) was used to determine the best match. Fig. 3 shows a flowchart of the multiple-constraint dense image-matching method. First, sparse seed points with strong textures (e.g., bare rocks and crevasses) were manually selected and matched to generate the initial TINs in the first image (master image) and the homologous TINs in the second image (slave image) as the initial geometric constraint. The ice velocity in a triangle should be approximately consistent. Therefore, the spatial patterns of the ice velocity should be considered during seed point selection. We generally selected fewer seed points around bare rock devoid of ice flows, while more seed points were selected in the fjord of Byrd Glacier, where ice flow increases significantly.
Second, feature-based image matching was performed. To extract uniformly distributed feature points from both bare rocks and poor-textured areas, the Shi-Tomasi corner detector [41] was used in this article after comparison with other methods. Corresponding matches of the feature points in each triangle of the first image were then predicted in the homologous triangle of the second image, based on an affine transformation between the two triangles. The best match was then determined using the similarity constraint (NCC) from a range of possible locations around the predicted matches. To improve matching reliability, we switched the first and second images and reperformed feature matching as a bilateral matching. Only the feature points that could be matched both ways were retained. Subsequently, a mismatch elimination strategy was applied to remove mismatches according to a predefined NCC threshold and local regional consistency (LRC) principle (described below, similar to the principles in [6] and [11]). Consequently, the reserved feature points were used to generate expanded TINs as new geometric constraints in the next step.
Third, to obtain more matches and ensure a high spatial resolution of the ice velocity map, we performed grid-based matching under the updated geometric constraints. The gridbased matching approach iterates through each pixel of the first image at equal spacing to obtain subscenes centered on each pixel (e.g., 24 × 24 pixels). Each subscene was then compared to a range of possible matching locations in the predefined search range (estimated based on the resolution and time interval of the two images) of the second image to obtain the displacement in the central pixel of the subscene in the first image. The best match at the subpixel level was determined by mathematical interpolation of the primary peak in the surface generated by the NCCs of all pixels in the second range. Similarly, the mismatch elimination strategy described in the previous step was employed to remove mismatches.
In the mismatch elimination strategy, we considered the complete spatial pattern of glacier ice velocity reflected by the seed points. NCC coefficients and flow directions of the seed points were calculated as references for the following steps. First, low-correlated matches were removed by an appropriate NCC coefficient threshold (0.1-0.2), which was set depending on the image pair quality (refer to the lowest NCC coefficient of seed points) and was obtained by trial and error. According to the LRC principle, ice flow should follow glaciological processes, and the magnitude and direction of the velocity vectors derived from the displacements of matches generally remain consistent within a close neighborhood. Therefore, first, if a glacier did not change dramatically for a long time (i.e., the difference between ice velocities of seed points and the existing ice velocity map on the glacier fjord was <20% in magnitude and 8°in direction), a few mismatches were removed because of the conspicuous difference between their velocities with the reference map. Here, an ice velocity map with high spatial resolution and complete coverage could be selected as a reference map, e.g., the comprehensive MEaSUREs ice velocity map (version 2) (derived from both SAR and Landsat-8 optical images [27]) and the ITS_LIVE multiyear composited map (generated by Landsat optical images [19]). We manually checked the ice velocity of Byrd Glacier on these ice velocity maps and found no outliers. Furthermore, the ice velocity directions were consistent with the glacier streamline on the optical images. However, this step should be skipped if this stability assumption was not satisfied, or the existing ice velocity map was unavailable. The mismatches would then be removed in the next steps based on the LRC principle of the ice velocity and the manually selected seed points. Second, we calculated the mean and standard deviation within a neighborhood (e.g., 3 km). When the standard deviation was large (generally ∼30 m/year), indicating large different velocities in a small region, we eliminated mismatches by manually inspecting the image texture and considering the glacier velocity patterns reflected by the seed points. For the remaining matches, velocities exceeding two standard deviations of the mean were deleted. The second part of the work can be reiterated to maintain the regional consistency of the ice flow field.
For the images from 2000 to 2014, the commonly used IMCORR software [11], [34] was employed to generate reliable results in much less time compared with our proposed method (Fig. 8). Considering the approximately 1-year interval of the image pairs, we used a reference window size of 32 pixels and a search window of 256 pixels to track the displacement of the surface features. Finally, all ice velocity maps were produced at a resolution of 500 m using the natural neighbor spatial interpolation method [42] implemented in ArcGIS.
3) Adjustment of Geolocation and Error Analysis: During image coregistration, it is difficult to manually identify subpixel information, which could introduce an artificial offset of within one pixel between images in a single velocity estimate. Similar to ice velocity estimation based on optical satellite images (e.g., Fahnestock et al. [11] and Gardner et al. [6]), we removed the remaining error based on ice velocities over exposed rock regions using an Antarctic exposed rock map derived from Landsat-8 images [43]. We calculated the mean and standard deviation of the ice velocities over the bare rock regions (v x , σ v x in x-direction and v y , σ v y in y-direction). The mean was used to correct the overall ice velocity so that its averages in the x and y directions over bare rock were zero. The standard deviations were considered to be the velocity measurement errors in the x and y directions. According to the error propagation, the velocity error over the bare rock regions was calculated using (4), which was considered for the entire ice velocity map

B. Ice Discharge and Mass Balance
MB of the ice sheet can be estimated by the input-output method, which compares the surface mass balance (SMB) over the interior basin with the ice discharge (D) by the glacier across the grounding line [5], and can be expressed as [44] Here, the SMB generally provides mass input to the surface of the ice sheet, including precipitation, condensation, and deposition [44]. D provides the majority of the mass output of the ice sheet. In this article, we adopted the ice discharge calculation approach introduced by Gardner et al. [6]. The flux gate (FG) (Fig. 1) was selected to move inland from major shear zones, avoid glacier shear zones with poorly constrained velocities, and follow nearby radio-echo sounding (RES) flight lines from which valid ice thickness data can be extracted. In ice discharge estimation, the distances between sampling nodes must be smaller than or equal to the sampling density of the ice thickness data and the resolution of the ice velocity data to avoid introducing additional integration errors. Here, each FG segment was defined between two gate nodes that were separated every ∼280 m, which is smaller than the 500-m resolution of our ice velocity maps. The unmeasured flux in the area between the FG and the grounding line was considered to be approximately equal to the SMB in this area [6]. The ice thickness of the FG on Byrd Glacier was the RES measurements data in December 2011 from the Center for Remote Sensing of Ice Sheets [37] and a small amount of BedMap-2 data [45]. Similar to Rignot et al. [5], we used the ice thickness values in the calculation of multiyear ice discharge. Owing to the incomplete ice velocity data at the FG in 1963-1974, we compared the ice velocity with the 1974-1983 result near the grounding line of Byrd Glacier to estimate a scaling factor for the calculation of ice discharge for 1963-1974 (based on an ice velocity of 839.8 m/year in 1974-1983 and 819.8 m/year in 1974-1983, a linear scaling of 1.02 was obtained); this scaling method was also adopted by [5]. For the SMB over the interior basin, we use the updated output products from the Regional Atmospheric and Climate Model version 2.3 (RACMO2.3p2; van Wessem et al. [46]), which incorporates upper-air relaxation, a revised topography, and tuned parameters in the cloud scheme to generate more precipitation toward the AIS interior and modified snow properties to reduce drifting snow sublimation and increase surface snowmelt. RACMO2.3p2 data of 1979-2016 were applied at a spatial resolution of 27 km in the AIS [46]. For this duration, the annual SMB of Byrd Glacier changed very little and had a negligible increasing trend (a slope of 0.06 Gt/year), except an outlier in 2011. Thus, we considered the SMB in the glacier basin before 2011 to be more stable, and rapidly changing after 2011. The SMB beyond the duration covered by RACMO2.3p2 (1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978) was estimated by taking the decadal average of 1979-1988, and the SMB after 2016 was regarded as the 5-year average of 2012-2016.
Error in the MB evaluation arises from the SMB calculation and ice discharge estimation. First, the SMB uncertainty was estimated according to the elevation-SMB relative bias approach [46], [47]. This was because the modeled RACMO2.3p2 SMB was validated with 3234 in situ SMB observations distributed across the entire AIS [46]. After validation, the SMB bias was calculated and binned into 500-m surface elevation intervals (0-250, 250-750, etc.). The annual SMB related uncertainty in each RACMO2.3p2 grid cell was considered to be in 3%-16%, depending on the surface elevation bin where the grid cell was located. Considering the insignificant surface elevation change around Byrd Glacier basin [48], a multiyear REMA DEM (200-m resolution) mosaic was used. Second, uncertainties in ice discharge were calculated by error propagation, considering the errors in ice thickness and ice velocity data. The ice thickness errors depend on different data sources, among which the errors of the Bedmap-2 product [45] were mostly >200 m. For the RES data, we used the errors collected by Gardner et al. [6], which were 10-110 m, similar to the errors from the crossover analysis in Gogineni et al. [21]. Consequently, uncertainty of the ice discharge of Byrd Glacier was estimated as ∼1.5 Gt/year, similar to the uncertainty of the SMB, and the uncertainty of the MB was ∼2.0 Gt/year.

C. Ice Front Position
Frontal retreat would reduce the buttressing that ice shelves provide to grounded upstream ice, which may be a possible reason for interannual acceleration [3]. As Byrd Glacier entered the RIS, we manually mapped the ice front positions of a part of the RIS [Figs. 1(a) and 10(a)] using ARGON mosaic and Landsat images from 1963 to 2020. Each Landsat image was georegistered to the Landsat-8 image in January 2019 using stable control points over the exposed rocks. After mapping these ice front positions, changes in the ice shelf front length were quantified using the curvilinear box method [49], [50]. We generated a series of reference boxes whose lateral sides were two parallel curves approximating the ice flow line [51]. The upflow side was fixed at an arbitrary position perpendicular to the parallel curves, and the downflow edge was the ice front position of the ice shelf. Then, change in the average width of the ice front was quantified by comparing the areas of the reference boxes and dividing them by the distance of the parallel curves. Here, the ice front position of part of the RIS [ Fig. 1(a)] was ∼420 km away from the grounding line of Byrd Glacier.

A. Ice Flow Velocity
We used 20 multisource satellite images to generate 19 ice velocity maps from 1963 to 2014. The temporal resolution was mainly determined by the availability of satellite images, which After preprocessing, the glacier surface features (e.g., crevasses and depressions) were enhanced (Fig. 4) and the geolocation errors of two MSS images in 1974 and 1983 were reduced from ∼17 km and ∼300 m to 47.1 and 55.7 m, respectively. Fig. 5 shows an example of the multiple-constraint dense image matching process for the Landsat image pair of 1974 and 1983. First, 157 seed points were manually selected and matched on the 2 images to generate 298 triangles as initial TIN constraints, with an average area of 14 976 km 2 for each triangle. Second, 64 462 and 64 988 feature points were extracted from the two images using the Shi-Tomasi detector. Among them, 36 411 were matched from the first image to the second image, 36 345 were reverse matched using the NCC constraint, and 28 823 bimatched feature points were obtained. During mismatch elimination, considering the lowest NCC coefficient of 0.12 of seed points, we first set a conservative threshold of 0.11 to ensure that the correct matches were not removed. Here, 646 points were first removed owing to this empirical NCC coefficient threshold. Then, because the stability assumption of Byrd Glacier was satisfied by comparing the maximum difference between estimated ice velocities in the seed points and the reference map (here, the MEaSUREs ice velocity map v2 or the ITS_LIVE ice velocity map), 1039 points on the glacier fjord (velocity greater than 200 m/year) were removed according to the LRC principle because of the apparent difference in velocity (20% in magnitude and 8°in direction) from the reference map, and 4523 due to the velocity magnitude difference (676 from manual inspection, 3484 due to larger than two standard deviations within a 3-km neighborhood). Thus Finally, five ice velocity maps from 1963 to 1999, with a resolution of 500 m, were generated using the natural neighbor spatial interpolation method [42] implemented in ArcGIS (Fig. 6). As one of the major outlet glaciers of AIS flowing into the RIS through the Transantarctic Mountains, the ice flow at Byrd Glacier gradually converges upstream, passes through the narrow outlet with dramatic acceleration from ∼200 to ∼850 m/year at the grounding line, and gradually slows down when reaching the RIS. From the three-dimensional view [ Fig. 6(f)], mountains (up to 3 km high) exist on both sides of Byrd Glacier with a large number of bare rocks. A continuous elevation drop from the inland ice sheet (100 km upstream of the grounding line) to the floating part of the glacier was observed (Fig. 6), together with a corresponding dramatic ice flow increase from <200 to >800 m/year.

2) Comparison Experiments of the Proposed Method With Common Software or Existing Productions:
To assess the reliability and effectiveness of the proposed method, we conducted two comparison experiments based on both poor- (Fig. 7) and high-quality (Fig. 8) images.
For poor-quality images, we compared the proposed method with two publicly available software packages (IMCORR and COSI-CORR) based on an image pair between January 1974 and November 1983 (Fig. 7). It was observed that our method obtained massive matching points with reasonable spatial patterns compared with IMCORR and COSI-CORR, indicating that the glacier experienced high flow velocity and changed smoothly along its fjord, while the mountain region moved relatively slowly. We compared the histogram of the matches within a 16 × 8 km window near the grounding line obtained using the three approaches. Our method matched 463 For high-quality images, we compared the ice velocity maps produced using the proposed method and IMCORR with one of the maps from the GoLIVE dataset (Fig. 8). The three ice velocity maps were based on a pair of Landsat-8 images from Here, 65 591 feature points and 532 854 grid points were matched by the proposed method under 95 manually selected seed points, and 334 009 matching points were generated using IMCORR software. Based on these matching points, two ice velocity maps were produced with a resolution of 300 m, similar to that of the GoLIVE map. From the ice velocities over bare rock, the errors of the three velocity maps were calculated to be 4.34 m/year (proposed method), 5.23 m/year (IMCORR), and 4.86 m/year (GoLIVE). We then generated two ice velocity difference maps by subtracting the GoLIVE map from the two produced maps. The histograms of the two difference maps were obtained after removing 1% of outliers. Finally, we estimated an average difference of 0.35 m/year between the velocities of the proposed method and the GoLIVE map, and 1.68 m/year between IMCORR and the GoLIVE maps. These differences were within the error of the three maps. The proposed method obtained more matching points than IMCORR software and produced an ice velocity map with a smaller difference with the GoLIVE map. However, in terms of time consumption, the proposed method required approximately 1 day for computation and 1 day for manual inspection, which was significantly longer than the 3 h duration of IMCORR.
Overall, the proposed method successfully produced historical ice velocity maps, when commonly used methods failed (Fig. 7). However, the proposed method was relatively less efficient when processing modern high-quality images (after 2000) because of the requirement for manual effort, which was originally designed for poor-quality images (Fig. 8).
3) Long-Time Behavior of Byrd Glacier: To evaluate the long time-series of changes in ice velocity of Byrd Glacier, we extracted the mean velocities [ Fig. 9(a)] of three boxes [Box 1, Box 2, and Box 3 in Fig. 1(b)] each with a side length of 3 km over the glacier centerline, showing an increase-decrease ice flow change trend along the centerline [ Fig. 1(b)]. All three boxes were distributed along the approximate centerline of the glacier (profile AA'). Box 2 was on the grounding line (the critical boundary between grounded ice and the ocean

B. Mass Balance
Based on a long time-series of ice velocity measurements, we obtained the MB of Byrd Glacier from to 1963 to 2020 [ Fig. 9(b)]. Over the last 57 years, the ice discharge across the grounding line from the catchment of Byrd Glacier to the RIS showed slight change, with an average rate of 18.6 ± 1.5 Gt/year and a standard deviation of ∼0.3 Gt/year. Owing to the large range of SMB change (maximum of 32.7 ± 2.1 Gt/year and minimum of 17.7 ± 1.1 Gt/year) and the smaller range of discharge (∼1.2 Gt/year), the MB of Byrd Glacier was almost dominated by the SMB, changing from −0.9 ± 2.1 to 14.5 ± 2.5 Gt/year. Overall, Byrd Glacier gained a total mass of 159.4 ± 15.3 Gt at an average rate of 2.6 ± 2.0 Gt/year over the last 57 years.

C. Ice Front Position
We observed the ice front position in a portion of the RIS [white box in Fig. 1(a) Fig. 10(b)]. After calving, the advancing speed of the ice front increased ∼16.2% to become ∼720 m/year from 2003 to 2020.

D. AWS Temperature
Based on a time-series of climate data during 1991-2020, recorded by the Marilyn AWS, we generated monthly and austral summer (December to February) air temperatures at Byrd Glacier. The glacier experienced a slight fluctuation with an average of approximately −23.2°C (Fig. 11). The summer air temperature had an average of approximately −9.14°C, increased from 1993 to 2011 and decreased later. Additionally, the monthly temperature was almost <0°C, with a maximum of approximately −3.6°C.

V. DISCUSSION
Several reliable in situ measurements exist for the ice velocity of Byrd Glacier. We compared these data with our historical satellite-based ice velocity maps to discuss the differences between them. Although a few studies introduced stake-based observations in the 1960s [23], [24], the exact locations of these in situ observations were not released to the public. On December 6, 1978 andJanuary 31, 1979, two sets of aerial photographs were obtained by Brecher [25] for photogrammetric processing, and the total displacements of 601 points (470 on the glacier) were measured to estimate the ice velocity of the glacier. Child et al. [8] regenerated two DEMs using aerial datasets to produce an ice velocity map by tracking the surface slope patterns and found that the velocity values along the main track agreed well with those obtained by Brecher [25]. Hughes and Fastook [53] obtained similar results in their field survey. We compared our results (from January 1974 to November 1983 using Landsat images) with the velocities released by Brecher [25] (Fig. 12). We transformed the coordinates of these points from the International 1924 ellipsoid datum into the WGS 1984 datum using an abridged Molodensky method [54] based on the parameters of the local geodetic datum, the Camp McMurdo Area, Antarctica [55]. Notably, 80 points were removed from the comparison because there were no corresponding matched points within 500 m. The scatter diagram of the two velocities on the remaining points [ Fig. 12(c)] showed good consistency, with a near-unitary slope (slope = 1.02) and coefficient of determination (R 2 = 0.95) for the fitted line. Fig. 12(a) shows the velocity difference map at each point, with the color representing the difference values. Notably, some outliers are mainly in the glacier shear zone with largely changed ice velocity owing to the low resolution of the historical images. The average Landsat ice velocity (January 1974-November 1983) was ∼26.5 m/year less than the average aerial velocity (December 1978-January 1979) [ Fig. 12(d)]. We further compared these two ice flow velocities along the profile of centerline AA' [ Fig. 12(a)], where each of the velocity values was calculated at a sampling distance of 1 km along the profile by averaging the valid velocity measurements within a 3-km length-side rectangle [ Fig. 12(b)]. It was observed that the two measurements were generally consistent upstream, and the difference gradually appeared from inland to the coast, with a maximum of ∼40 m/year near the grounding line. This difference was probably owing to the disparate time interval of the two measurements (56 days for aerial surveys and ∼10 years for satellite image coverage), during which the ice velocity might experience a short-time fluctuation.
The proposed method has the potential to improve the spatial resolution of the ice velocity map obtained from historical low-quality images. Based on feature-and grid-based matching strategies, we produced massive reliable matching points to generate historical ice velocity maps at a high resolution of 500 m (∼8 times the 60-m resolution of MSS images or ∼17 times the 30-m resolution of TM images), which performed better compared to the existing products, e.g., 50 times for LISA products (750-m spatial resolution of ice velocity maps vs. 15-m source image resolution [18]) and 20 times for the GoLIVE products (300-m spatial resolution of ice velocity maps vs. 15-m source image resolution [17]). However, because of the requirement for manual efforts in the selection of seed points and mismatch elimination, the efficiency of our proposed method is relatively  lower than that of the widely used software. For example, our method took approximately 2 days to generate ice velocity maps based on a pair of Landsat-8 images, while IMCORR software only took ∼3 h (Fig. 8). Thus, one of the limitations of our proposed method is that it is not suitable for large-scale ice velocity mapping over the entire ice sheet or for processing high-quality Landsat-7/8 or Sentinel-2 images. The method can produce historical ice velocity maps, where commonly used methods (e.g., IMCORR and COSI-CORR) fail, and fill the data gap because of poor image quality (Fig. 7). Owing to the infrequent data collection of early satellite and cloud coverage, available image pairs generally have a long-time interval, and the glacier surface features may change dramatically or disappear due to glaciological processes over this period (e.g., surface crevasse evolution, ice front calving, etc.). Another limitation is that the method may fail when processing image pairs with a long-time interval in very fast-moving glaciers, such as Thwaites Glacier in West Antarctica. Finally, owing to the initial constraints from the manually selected seed points, the method was tolerant to the illumination conditions of the images.
The spatial pattern of the ice velocity of Byrd Glacier was controlled by the basal terrain and amount of ice funneling through the trunk [56]. We observed a rapid basal elevation decrease from ∼0 to ∼−3000 m at a distance of 110-65 km upstream of the grounding line (Fig. 13), after which the bed elevation continuously increased. This type of retrograde bed was also found in other glaciers and was considered potentially unstable if the grounding line retreated inland [57]. Along the profile, the ice velocity of Byrd Glacier gradually accelerated from ∼230 to ∼830 m/year, except for a constant part of ∼430 m/year from 90 to 60 km upstream of the grounding line. Finally, after the ice flow reached the grounding line, its velocity gradually decreased. We observed that the acceleration stopped ∼30 km before the basal valley bottom (a depth of ∼ 3000 m and a distance of    Fig. 10(b)], they only caused the loss of the apparent "passive shelf ice" (Fig. 1), which can be removed owing to little or no influence on the buttressing and ice shelf dynamics [58]. Additionally, from the diagnostic modeling experiments, the change in buttressing from the regions near the calving fronts of the RIS showed insignificant influence (∼1%) on the grounding-line ice velocity of Byrd Glacier [59]. According to satellite observations and hydrostatic equilibrium calculations based on DEMs, little to no change in the grounding zone of Byrd Glacier has occurred over the last 40 years, indicating its steady state [8], [60]. Temperature in the Marilyn AWS was just above freezing point within a few hours during the austral summer, suggesting a low chance of sustained surface melting in summer for Byrd Glacier. Furthermore, the highest summer air temperature was recorded in 2012 when Byrd Glacier experienced a relatively low ice velocity.  Fig. 1, and the blue dashed line indicates its surface elevation from a REMA DEM [63]. Brown dashed line indicates bed topography from BedMachine [57]. All elevation was referenced to the WGS84 ellipsoid.
These phenomena suggest no remarkable correlation between air temperature changes and glacier ice velocity. A subglacial drainage system, including two large active subglacial lakes (Byrd 1 and Byrd 2 ) and 14 small lakes in Byrd Glacier catchment ( Fig. 1), was discovered through surface deformation analysis using ICESat altimetry data [61]. Byrd 2 was filled gradually from 2003 to November 2005, and drained sometime between March and November 2006, then recharged from mid-2007; Byrd 1 shared a similar pattern of surface displacement, but with an approximately 6-month delay [26], [61]. These sudden and short-lived drainage events were considered to have caused an acceleration of ∼80 m/year (10%) in Byrd Glacier between November 2006 and February 2007 [26]. Instead of finding this large and rapid change in ice velocity, we observed a relatively high average ice velocity (∼6%) between February 2005 and December 2006 [Box 2 in Fig. 9(a)] as a possible result of the sudden and short-lived drainage event. The changes in ice dynamics in Byrd Glacier might be related to the activity of the subglacial drainage systems. Other factors, such as grounding line migration, decreased backstress from the ice shelf, and ice thickness changes due to oceanic melt, may alter the velocities over a much longer timeline than those estimated during the acceleration. Moreover, any changes in dynamics that Byrd Glacier is experiencing external forces would also be likely to affect the other outlet glaciers draining the Transantarctic Mountains. Table AI shows the comparison of MB estimates between our results and the previous articles (e.g., Stearns [22] and Rignot et al. [5]). In this article, nearly all the estimated discharge rates were below the SMB in the corresponding time, resulting in a state of slightly positive MB of 2.7 ± 2.0 Gt/year during 1963-2020, compared with the nearly zero or negative values from Rignot et al. [5] 29 Gt/year, which was attributed to an overvalued accumulation of 44.27 ± 3.32 Gt/year from the previous RACMO2/ANT version [62]. Here, usage of a different version of the RACMO models was a reason for the different MB results. The ice discharge from this article (average of 18.6 ± 1.5 Gt/year and a standard deviation of ∼0. 3 Gt/year during 1963-2020) was overall slightly smaller than those from Rignot et al. [5] (average of 22.27 Gt/year with a standard deviation of ∼0. 5 Gt/year in 1979-2017). Ice discharge is related to the ice velocity and ice thickness at the grounding line. By comparing our ice velocity maps with those of Rignot et al. [27], the difference between ice velocity at the grounding line and FG were found to be small (∼4%). Stearns [22] calculated ice thickness from hydrostatic equilibrium, Rignot et al. [5] employed the BedMachine dataset, and this article used ice thickness of FG almost entirely from flight data [6]. Thus, the difference in discharge between these studies is probably due to the ice thickness at the grounding line, where narrow submarine valleys have been challenging to resolve for radar sounding for decades [57]. In general, little change in discharge and a slightly positive MB from this article indicate that Byrd Glacier has been relatively stable for the past 57 years, and may not have been significantly impacted by the global climate change.

VI. CONCLUSION
In this article, we proposed a multiple-constraint dense image matching approach to estimate the ice velocity of glaciers, based on poor-quality historical satellite images with low spatial resolution that were acquired before the 1990s. By integrating multisource satellite remote sensing images and existing ice velocity products, reliable and long time-series of ice velocity maps of Byrd Glacier from 1963 to 2020 were generated. This method could be used for ice velocity studies of other glaciers in the cryosphere.
Our results indicate that the spatial pattern of the ice velocity of Byrd Glacier is mainly affected by bed terrain. Ice velocity of Byrd Glacier slightly fluctuated with time, which was not affected by the calving events in the portion of the RIS or the air temperature changes; however, activity of the subglacial drainage systems may have contributed. Besides, the MB of Byrd Glacier was almost dominated by the SMB and was slightly positive at an average rate of 2.6 ± 2.0 Gt/year from 1963 to 2020, indicating that global climate change may have a limited impact on Byrd Glacier. APPENDIX   TABLE AI  LIST OF SURFACE MASS BALANCE, ICE DISCHARGE, AND MASS BALANCE ESTIMATED IN DIFFERENT STUDIES