Structure-From-Motion Photogrammetry of Antarctic Historical Aerial Photographs in Conjunction with Ground Control Derived from Satellite Data

: A longer temporal scale of Antarctic observations is vital to better understanding glacier dynamics and improving ice sheet model projections. One underutilized data source that expands the temporal scale is aerial photography, specifically imagery collected prior to 1990. However, processing Antarctic historical aerial imagery using modern photogrammetry software is difﬁcult, as it requires precise information about the data collection process and extensive in situ ground control is required. Often, the necessary orientation metadata for older aerial imagery is lost and in situ data collection in regions like Antarctica is extremely difﬁcult to obtain, limiting the use of traditional photogrammetric methods. Here, we test an alternative methodology to generate elevations from historical Antarctic aerial imagery. Instead of relying on pre-existing ground control, we use structure-from-motion photogrammetry techniques to process the imagery with manually derived ground control from high-resolution satellite imagery. This case study is based on vertical aerial image sets collected over Byrd Glacier, East Antarctica in December 1978 and January 1979. Our results are the oldest, highest resolution digital elevation models (DEMs) ever generated for an Antarctic glacier. We use these DEMs to estimate glacier dynamics and show that surface elevation of Byrd Glacier has been constant for the past ∼ 40 years.


Introduction
Quantifying how the cryosphere responds to climate forcings is essential for predictions of ice sheet mass loss and sea level rise. Unfortunately, records of Antarctic glacier behavior only extend a few decades-the extent of the modern satellite era. As such, glacier characteristics such as ice velocity and surface elevation are only available from the early 1990s to the present. Based on our own data searches, infrequent data collection from older satellites (e.g., early generation Landsat and spy satellites) consists of coarse resolutions and/or significant cloud cover. As a result, little is known about Antarctic glacier behavior prior to 1990. In this study, we outline a methodology for using historical analog aerial imagery to quantify glacier behavior. By extending the temporal scale of glaciology data, we improve our understanding of ice sheet behavior in a changing climate, which is essential for physically-based sea level rise projections.
The 4th IPCC assessment report specifically identified the need for high-temporal resolution remote sensing data to improve understanding of glaciological processes [1]. The 5th IPCC assessment report (∼six years later) stated that much more effort was put into launching new polar-observing satellites as well as collecting additional airborne data [2]; however, the most current special IPCC report on the cryosphere states that longer observational records are still needed in order to better quantify and understand current ice sheet changes [3]. Aerial imagery has been collected since the 1930s in Antarctica [4,5], but only a small fraction of the data has been used in analytical research. These data were originally collected for topographic mapping, but have largely been neglected since; one major reason is that the imagery is inherently difficult to process.
Digitally-scanned historical images have been successfully used to render accurate glacial terrain surfaces elsewhere on the globe, but for Antarctica, this type of data processing has occurred only on the Peninsula e.g., [6][7][8][9][10]. Processing historical imagery in Antarctica is uniquely difficult because of the lack of exchangeable image file (EXIF) information and ground control. As a result, rendering elevations and orthoimages requires a great amount of manual effort. We build the necessary EXIF files from metadata extracted from the original flight records and use a new method of manually extracting ground control points (GCPs) from high-resolution digital elevation models (DEMs) to provide the missing ground control. This technique has been applied to small Arctic glaciers [11,12], but never to Antarctic glaciers, which are often wider than individual images (where images of just glacier ice have no suitable features for assigning GCPs). To assess the applicability of structure-from-motion (SfM) processing with remotely extracted GCPs for the purposes of calculating elevations, we apply this methodology to an extensive collection of 1970s aerial photographs of Byrd Glacier, East Antarctica to create high-resolution elevation and velocity products.
In December 1978 and January 1979, two sets of aerial images were collected over Byrd Glacier, as part of a large-scale research project [13]. The original image processing was carried out in the early 1980s, after which the film went missing for 36 years. In 2016, the original film negatives were located in unmarked film canisters in the basement of the USGS Library in Reston, Virginia. They were then sent to the EROS Center in Sioux Falls, South Dakota, where they were digitally scanned. We use SfM photogrammetry software to process the scanned images and compare the elevation and velocity results with the manually derived products from the 1980s [14]. We assign real-world geopositions to the images from GCP XYZ values that we extract from high-resolution satellite data. To the best of our knowledge, there are no other studies in Antarctica using SfM photogrammetry processing with manually extracted ground control to generate elevations from analog aerial imagery. Urbini et al. [15] do use some form of SfM processing with extracted GCPs, but only to orthorectify aerial imagery; no historical elevations are produced. Our workflow, which relies on modern processing techniques applied to historical images with little ground control, can also be applied to other historical aerial images. Over 300,000 aerial photographs were collected in Antarctica from 1930 [16,17] to 2000 [18], but they have remained untouched for decades, largely due to the arduous processing logistics. Here, we demonstrate the validity and accuracy of this elevation processing technique for historical Antarctic data, which can later be applied to this large dataset of historical images.

Study Site
Byrd Glacier is located at ∼80.4 • S, where it flows through the Transantarctic Mountains and terminates in the Ross Ice Shelf (RIS) (Figure 1). The catchment basin of Byrd Glacier, the largest in the world, is 933,744 km 2 [19]. The glacier trunk is ∼100 km long and ∼25 km wide and is bounded by the Britannia Range to the north and the Churchill Mountains to the south. At the grounding line, the average ice thickness is ∼2000 m, and the ice is moving at a velocity of ∼850 m y −1 [20]. The surface of Byrd Glacier is heavily crevassed and surrounded by steeply sloped terrain. These characteristics make Byrd Glacier a good candidate for SfM processing because textured surfaces produce the best results from stereoparallax calculations. In this region of Antarctica, there is no vegetation, little snowfall, and no man-made structures, so the DEM construction of the topography is actually a digital terrain model (DTM), but is referred to as DEM throughout the rest of this paper.

Materials
We use a combination of optical imagery, GPS, and vector data to create and validate historical DEMs. To produce historical surface elevations, ground control must be applied to image geometry in order to calculate surface relief with absolute XYZ values. Our ground control is manually extracted using present-day high-resolution elevation data. The GPS and vector data are used to register the historical DEMs to absolute geolocations and elevations as well as assess the accuracy of the final results.
3.1. Aerial Photography 3.1.1. Imagery On 6 December 1978 and 31 January 1979 (56 days apart), 568 photographs were collected by the U.S.Navy. All images were acquired using a Wild-Heebrugg RC8 Aerial Camera [21] that captured vertical (0 • nadir) images on ∼23 × 23 cm or ∼9 × 9 in Kodak Safety Film (film type printed on each image). The along flight image overlap is ∼60% and the side overlap is ∼35% (Figure 1). The areal coverage of a photograph is ∼80 km 2 . The film was scanned by the USGS EROS Center at 1000 dpi [18], 25 µm, in 2016.
The plane collecting the imagery was a C-130H aircraft flown by the U.S. Navy VXE-6 at an altitude of ∼7000-8000 m according to the SfM photogrammetrical calculations. This altitude range with a lens focal length of 152.28 mm [21] means the images were captured at an approximate map scale of 1:48,000-1:50,000, or a ground sampling distance of 1.5-2 m. All flight lines were parallel with ice flow with the exception of one semiorthogonal flight line on 6 December 1978. A malfunction with the instrument lights at the edge of the frame occurred during the 1978 acquisition, resulting in white stripes in all images acquired on that day ( Figure 2).

Original Study
The aerial dataset was initially processed in the early 1980s, using manual photogrammetric techniques. A total of 322 images (152 from 1978 and 170 from 1979) were used with a stereo-comparator to estimate XYZ positions by a least-squares bundle block adjustment of the two blocks of photographs simultaneously [14]. The relative horizontal and vertical errors were estimated at 0.8 m and 1.8 m, respectively. In situ ground control was established from electronic traverse trilateration measured from JMR-1 TRANSIT Doppler receiver positioning equipment and ground survey [13]. A total of 13 GCPs were recorded over stable relief surrounding Byrd Glacier. Point elevations were generated from each set of images and surface ice velocities were calculated for the period of image collection. Velocities were derived from displacement values of manually identified 601 common features between the two time periods, with an approximated uncertainty of 5 m [14].

High-Resolution DEM
A high-resolution DEM is essential for remotely identifying accurate GCPs. The only high-resolution DEM currently available over Byrd Glacier is the Regional Elevation Model of Antarctica (REMA version 3.0 and higher; [22]). This dataset is produced by using panchromatic strips with a spatial resolution of 0.5 m collected from DigitalGlobe's World-View sensor over 2016-2017. REMA is publicly available at the Polar Geospatial Center (PGC) at a horizontal resolution of 2 m projected horizontally to WGS84 Antarctic Polar Stereographic and vertically to the WGS84 ellipsoid. The present-day elevation product is a series of REMA DEM strips which were generated from a surface extraction TIN-based search-space minimization (SETSM) algorithm after Noh and Howat [23]. These DEMs are created without ground control, which means that the geolocation is approximated from the sensor's rational polynomial coefficients (RPCs), and the subsequent error is entirely dependent on the accuracy of the RPCs [23]. Unprocessed WorldView imagery has a vertical accuracy of ∼5 m excluding terrain effects [24]. The online documentation provided by PGC reports the vertical error to be 3.5 m; however, Noh and Howat Noh and Howat [25] report accuracies of ∼2-3 m in all directions and ∼0.2 m for DEMs coregistered to dense ground control. An error estimate of ∼0.2 m in the Transantarctic Mountains is unlikely due to image collection geometry (e.g., high degrees off nadir) and sparse ground control. A comparison of an SETSM DEM with LiDAR data resulted in errors of 4-5 m in regions of steeply sloped terrain [26] and as a result, we adopted an uncertainty of 5 m.

Rock Mask
Coregistration of the present-day DEM and historical DEMs can only be accurately conducted over stable terrain like rock outcrops [27,28]. We use a publicly available Antarctic rock mask generated from a normalized difference snow index (NDSI) using the Landsat 8 OLI green and shortwave infrared (SWIR-1) bands [29]. The spatial resolution of those bands is limited to 30 m, so we verify all rock outcrop boundaries using the accompanying orthoimagery of the REMA DEM strips (posted at 2 m).

GPS Data
From 2010-2013, a network of NetR9 Trimble GPS receivers was deployed on Byrd Glacier, including three reference static sites on exposed bedrock. The data from these bedrock sites was used to evaluate the absolute accuracy of the SfM processed elevations. GPS data were processed using Litchen and Border's [30] GIPSY software package with high-precision kinematic data processing methods (e.g., Elósegui et al. [31,32]). By calculating the weighted root-mean-square scatter of the vertical position estimates over the entire deployment time, an uncertainty of 2-5 cm was determined. These are the only in situ XYZ measurements from the study site.

Ice Thickness
Ice thickness is calculated by subtracting radar-derived bed topography (assumed to be static over grounded ice) from surface elevations. Thickness values are a product of ground penetrating radar ice bed picks from data collected by the Center for Remote Sensing of Ice Sheets (CReSIS) during the austral summer 2011/12 [33]. We specifically use the radar depth sounder (MCoRDS V2) data that operates at 195 MHz with a bandwidth of 30 MHz [34]. We use approximately 25 of the flight paths that cross over the central trunk of Byrd Glacier in our analysis. The radar data is referenced vertically to the WGS84 ellipsoid and has a crossover error estimate of ∼30 m [34].

Methods
We use a combination of manual and automated techniques to generate historical DEMs. We begin by extracting ground control which we then apply to the SfM photogrammetric processing. Figure 3 is a diagram of the workflow for extracting GCPs and building DEMs. Table 1 lists the various processing parameters for each set of images. . Workflow for generating DEMs from historical imagery. On the left are the steps for manually extracting ground control points (the software used for coregistering is PyBob) and on the right, the SfM photogrammetry processing of the analog aerial images (the images processed using in MicMac and the point clouds are cleaned with CloudCompare). The orange shapes represent manual procedures; the teal shapes are automated steps.

Ground Control Points
Traditional ground control is typically established using in situ geodetic surveys. For most of Antarctica, this type of ground control is unavailable due to inaccessibility, highly sloped relief, and the harsh environment. Byrd Glacier is no exception; there is not enough in situ ground control to produce accurate digital elevations from the historical aerial images. Here, we describe a new methodology that utilizes manually extracted GCPs from a present-day DEM as ground control for the purposes of generating historical elevations. This method allows for GCP extraction of mountain peaks and other isolated terrain data that would otherwise be near-impossible to obtain. GCPs consist only of features of stable, unchanging terrain; we do not use surface features that may have undergone substantial motion over the past ∼40 years (e.g., ice, permafrost). While no ground control is physically applied to the glacier ice, the SfM processing is still able to calculate elevations based on surrounding GCPs from the exposed bedrock. We take into consideration that the Transantarctic Mountains are uplifting at a rate of 3 mm y −1 due to isostatic adjustment from the Last Glacial Maximum [35,36], but it is considered negligible as the cumulative uplift of 0.12 m over the ∼40 years of our study is well within our margin of vertical error.

Blunder/Cloud/Shadow Masking
Before GCPs can be identified, the REMA data is corrected for blunders. These blunders originate from effects of shadows, clouds, and the high degree of off-nadir with which the imagery was collected. Noh and Howat [25] acknowledge that the SETSM algorithm has trouble processing some surfaces like shadowed regions, but that this occurs when shadows are just the result of "low sun angles". Conversely, other studies have obtained more accurate elevation results from imagery collected at low sun angles [37], and we find miscalculated elevations from SETSM-processed imagery collected at both high and low sun angles. To remove potential blunders from shadows, we opt for a conservative approach of clipping all shadowed regions within the DEMs. Polygons are manually drawn over shadowed regions identified in orthoimagery, and then removed from the DEMs. We ensure that topographic displacement in the orthoimagery is considered by creating polygon boundaries that are larger than the actual shadowed region.
The SETSM algorithm does contain a cloud detection component [23], but it misses a significant percentage of cloud coverage. To ensure cloud-free data, and because it can be difficult to determine if the DEM contains clouds solely from viewing a DEM, we follow the same technique for shadow removal. Every orthoimage is inspected and manually delineated cloud cover is eliminated from the DEM data.
The final blunder type embedded within the present-day elevations is due to the high degree of off-nadir with which the imagery was collected. Large, off-nadir degrees (10 • -30 • ) means that terrain like the backside of mountains is outside the satellite sensor's field of view and not captured. However, the REMA data does contain derived elevation values outside of the field of view, which are erroneous, and were discarded. We identify the off-nadir and shadow zones from both the orthoimage and the DEMs. In the orthoimagery, these zones appear distorted with a wavelike pattern. In the DEMs, the regions contain either clusters of >4 pixels, all of the same elevation, or sharp rectangular boundaries. These "out of view" regions are also manually deleted from the DEM strips.

Coregistering/Mosaicking
To ensure accurate and uniform elevation values throughout the mosaic, each DEM strip is first coregistered as a way of standardizing the collective dataset. Registering DEMs is necessary because while the DEMs all have the same ground resolution, their grids are not aligned and even subpixel misalignment causes considerable bias when estimating surface elevation change [27,28,38]. Other errors in an elevation dataset, outside of processing-atmospheric effects, steep relief slopes, different sensors used to collect data, time of year-can also cause erroneous elevation values [28]. Coregistering one DEM to another minimizes these errors [27].
All present-day DEM strips are resampled from 2 m to 5 m spatial resolution using cubic convolution and projected to WGS84 Antarctic Polar Stereographic. The resampling helps to smooth any missed erroneous elevations and fills some gaps from blunder extraction. With the exception of the two furthest up-flow DEMs, each strip is then coregistered to its immediate upflow newly registered neighbor using the Python module PyBob, ref [39]. This coregistration method (after Nuth and Kääb [27]) is a robust analytical solution that uses the pairs' elevation difference residuals and the terrain's aspect and slope to correct for bias and errors between DEM pairs. The coregistered DEMs strips are then mosaicked together to create one DEM of the entire glacier. A list of all REMA DEM files in the mosaic is included in the Supplementary Materials. Throughout the rest of this paper, the mosaicked DEM is referred to as "present-day elevations".

Identifying GCPs
We evaluate potential GCP locations of stable terrain from WorldView imagery that was orthorectified using REMA data by PGC. The imagery can not be used to determine the absolute locations of GCPs because the high degree of nadir with which the data were collected caused the planimetric accuracy to vary throughout the entire ortho dataset. Data collected at high (>10 • ) degrees off-nadir contain the strongest effects from topographic displacement and occlusion. Therefore, the orthoimagery is only used as a guide. To avoid feature misplacement, GCPs are chosen exclusively at local elevation maxima which are easily identifiable within the DEMs. This approach ensures that GCP placement error is a product of the DEM and not the orthoimages.
SfM processing only requires three GCPs for stereo processing, but increasing the number of GCPs used ensures improved bundle adjustment results [11]. A larger number of correctly placed GCPs means the bundle adjustment is more likely to correct misplacements. The exact number of GCPs necessary for accurate elevations varies from project to project, but research shows that more GCPs can improve accuracy to twice that of the image's pixel ground sample distance (GSD) [40,41]. In classic photogrammetric processing, accurate results are achieved with evenly distributed GCPs in a triangular grid [41]. This arrangement reduces the horizontal distance to any GCP. In our case study, exposed bedrock is not evenly distributed throughout our study site making a perfectly distributed grid of triangulated ground control unrealistic.
We use 141 and 162 GCPs for the 1978 and 1979 image sets, respectively. Other studies using SfM processing for glacier elevations also find that implementing a large number of GCPs increases the overall accuracy of the final DEM [42,43]. The historical DEM's absolute orientation is established by applying at least one GCP, though as many as ten GCPs, for every image containing exposed bedrock. This means that there is one GCP added every 17 km 2 . A map of GCP distribution is available in the Supplementary Materials.

Elevation Processing
The historical surface reliefs are constructed with an open-source SfM photogrammetry software called MicMac. MicMac is a project hosted at the National Institute of Geographic and Forestry Information and the National School of Geographic Sciences [44]. The software provides a strategic framework that includes manual and automated processing steps which is necessary when calculating elevations from historical analog imagery.

Interior Orientation
Before performing stereo calculations with the imagery, a uniform geometry is first established for every image. Analog film geometry is in a gridded coordinate system where the columns are X and the rows Y in millimeter units. The coordinates of the principal point from one digitally scanned aerial image is not guaranteed to be identical to the next due to asymmetrical scanning and possible physical distortion of the film, likely in this case caused by nonideal storage conditions. We determine the image coordinate system by measuring the distance between the N-S and W-E fiducial marks (vertical and horizontal directions, respectively, Figure 4a) as well as the fiducial mark XY positions in Photoshop. These parameters are then applied to each image's four protruding fiducial marks, manually identified in MicMac, to establish the imagery's interior boundary. Following this, the images are resampled to a resolution of 25 µm which creates a homogeneous geometry throughout the historical dataset. Figure 4b is a schematic showing the result of normalizing a skewed digitized photograph. Correcting these misalignments assured that all of the images contain the same pixel size and grid alignment, without which automatic tie point identification would be impossible.

Relative Orientation
Using the normalized images, MicMac identifies common features (e.g., tie points) within each image using the Scale-Invariant Feature Transformation (SIFT) algorithm [45]. All of the images contain black borders and some flight metadata text. Unfortunately, the software treats the border and text as ground features and includes them in the initial tie point extraction. A mask is created to exclude the text and borders for each image, and the white stripes visible in the 1978 dataset. After masking, there is a total of 3.7 million tie points for the 1978 data and 2.7 million tie points for the 1979 data. The additional million tie points in the 1978 set are attributed primarily to the one semiorthogonal flight flown in 1978. Relative locations of the tie points are estimated using Fraser's [46] calibrated distortion model. Camera distortion parameters are calculated using the least squares method which iterates until the calculated residuals reach an optimal value. The resulting average residuals for the two sets of data are 0. 91 (1978) and 0.95 (1979) of a pixel (∼1.2 m).

Absolute Orientation
GCPs identified from the present-day DEMs are applied to the relative coordinate system of the tie points to relate the tie points to real-world (absolute) ground positioning.
The GCPs are manually positioned in the aerial imagery in MicMac. The interior image information is then transferred to absolute XYZ values of the GCPs. Following this, a bundle adjustment is executed to minimize error between light-ray bundles and GCP locations by adjusting the camera's orientation and position [47]. The residuals calculated by the bundle adjustment represent the distance the GCPs are estimated as "offset" from the manually placed GCPs. For the 1978 and 1979 image set, the average residuals are 1.00 and 0.96 pixels, respectively (see Table 2). These pixel residual values are comparable to errors quoted by other studies using historical aerial imagery [48]. The bulk application residuals provide one resource for the overall error of the GCPs and the distance to its corresponding point in the point cloud. The final product produced from the SfM photogrammetry processing is a point cloud of the elevation values from each time stamp. We filter each point cloud for spurious points using the software CloudCompare (Version 2.11) [49]. Spurious points are XYZ values outside of the known local minimum and maximum elevations, as well as points that are clearly erroneous elevations. Potential causes of such points are bad radiometric calibration, homogeneous features, film distortion, and inaccurately placed GCPs. After manually removing the inaccurate points, we convert the filtered point clouds into a gridded raster format. Grid cell values are an average of the points that fall within a the spatial extent of given cell [49].
The initial point clouds generated by MicMac have a 1 m spatial resolution, but after editing for spurious points, the final DEMs are gridded to 10 m. The coarser resolution is still an accurate representation of the terrain, while also acting to fill in data gaps that would otherwise remain at the 1 m resolution.

DEM Registration
Before any validation is carried out, and the point clouds are converted to grid format, both point clouds are registered to the GPS data in CloudCompare. Aerial images are draped over the point cloud in CloudCompare to increase the accuracy of GPS placement to the DEMs ( Figure 5). We recognize that three in situ data points is less than the optimal number of GCPs needed to accurately register an image [50], but these points are the only absolute elevation data of stable terrain available for Byrd Glacier.

Velocity
The white stripes in the 1978 imagery make producing an orthorectified mosaic impossible. This means that traditional methods of estimating velocities from optical imagery (e.g., feature tracking [54]) cannot be employed. Our solution instead is to use a similar method to Ryan et al. [55], but rather than use the surface elevations, we use feature tracking of surface slope patterns, a byproduct of the DEMs. Velocity determined from feature tracking techniques requires that elements within a dataset remain the same from one time period to the next in order to track their displacement. While elevations will change between epochs, over a short time period (56 days) the relative shape of the terrain, such as the slope, remains fairly constant. We estimate glacier velocities using the Optically Sensed Images and Correlation (COSI-Corr) software [56][57][58]. The tracking is done using patches of specified number of pixels that depict patterns of surface slope [59]. We use a patch size of 64 × 32 pixels. The same patterns of slope are located on both DEMs, which means the amount and direction of slope horizontal displacement is quantifiable. A signal-to-noise ratio (SNR) raster is also included with the velocity results. The SNR raster is a grid of correlation percentage which allows us to assess the quality of the estimated feature tracking. Like other studies using COSI-Corr to estimate feature displacement (e.g., [60,61]), we discard any velocity results with a corresponding SNR value <0.9. We also eliminate any bearings not congruent with ice flow which we determine from previous studies of Byrd Glacier's velocities (e.g., [62][63][64]) (∼5% of the results). The remaining velocity values are interpolated with ordinary kriging to a resolution of 500 m.

Data Validation
To determine how well the DEM generation methodology performed, the vertical error of the historical elevations is estimated. We assess the error by comparing the stable terrain of the historical DEMs with present-day in situ elevations, present-day DEMs, and one historical DEM to another. A map of these differences is available in the Supplementary Materials. We also compare the SfM elevation results with those calculated by Brecher [14].

Elevation Comparisons
When the historical images were originally processed in the 1980s, the results were projected horizontally and vertically in the International 1924 datum [14]. The SfM DEMs produced in this study are projected in reference to the WGS84 ellipsoid datum. The only way to compare the vertical data from these two elevation sets is to transform the 1980s elevations to the SfM DEM's ellipsoid. The 1980s XYZ values are transformed to the WGS84 ellipsoid using a local geodetic datum with the Standard Molodensky three-parameter (∆Φ, ∆λ, ∆h) formula [65]: where Φ Local , λ Local , and h Local are local geodetic latitude, longitude, and elevation. The specific reference ellipsoid and shift parameters are from the local geodetic datum 1987 Camp McMurdo Area, Antarctica [66].
Brecher [14]'s elevation data was generated from 472 irregularly distributed discrete points for the trunk of Byrd Glacier. Interpolating the point data to a 500 m grid produces "bulls-eyes" in the final raster due to the sparseness of the points. The coarse resolution (500 m) makes it difficult to directly compare with the high-resolution (10 m) SfM DEMs. We therefore instead calculate the elevation differences at the discrete point data locations rather than each interpolated grid cell.
The elevations derived in the 1980s do not have spatial coverage beyond the glacier boundaries [14], which means that the original published dataset does not have any XYZ values on stable terrain. This issue was resolved in Schenk et al. [67] by reprocessing small sections of the imagery along four regions on the edge of the glacier from diapositives, but that is beyond the scope of our study. Consequently, we are unable to register Brecher [14] results to the SfM results, and we cannot perform an absolute direct comparison. Instead, a relative comparison is implemented from the mean and standard deviation values of the vertical differences between the Brecher [14] and the SfM elevations.

Velocity Comparisons
We compare the SfM velocity results with the speed contours published in Brecher [14]. The vector speed isolines are converted to a raster grid with 500 m cell size using ordinary kriging interpolation. Like the Brecher [14] elevation dataset, the Brecher [14] velocity spatial coverage also does not extend to any stable terrain.

Grounding Zone
We use the new historical DEMs to assess the stability of Byrd Glacier by estimating the grounding line migration from 1979 to 2017. A grounding line is the location where the course of the glacier's orward movement shifts from ice flow over bedrock to ice flow floating on the ocean. The grounding zone is the region where this transition occurs. Byrd Glacier flows over a retrograde bed at the grounding line [68] which has shown in models to be more influential to grounding line migration due to changes in ice dynamics [69]. Migration of the grounding line and zone are indicators of a glacier losing or gaining mass [70] and thus in or out of steady-state. The grounding zone is bounded by the point of flexure (upper limit: the hinge position of tidal influence) and the hydrostatic equilibrium boundary (lower limit: complete flotation) [52,71]; their locations can differ by as much as 10 km [53].
We estimate total migration of the grounding zone by calculating the hydrostatic equilibrium boundary using the historical 1979 and present-day 2017 elevation datasets. Determining migration of the hydrostatic equilibrium boundary is a valid substitution for assessing actual grounding line migration because the parameters that regulate grounding line shift (e.g., ice thickness change, basal melt) are the same as those that influence hydrostatic equilibrium boundary change. Boundaries are estimated using the 1979 SfM and present-day DEMs with the ice thickness data derived from the manually picked CReSIS echograms. Hydrostatic equilibrium is modeled by: where, h is the orthometric surface elevation, H the ice thickness, ρ w is the density of ocean water (1028 kg m −3 ) and ρ i the density of ice (917 kg m −3 ). We convert the presentday and historical DEMs to orthometric heights using Pavlis et al. [72] EGM2008 model. This portion of Byrd Glacier is entirely blue ice [73], so a firn correction is not necessary in our calculations [74]. A glacier is in hydrostatic equilibrium where Equation (2) is approximately zero.
To estimate the total mean migration of the hydrostatic equilibrium boundary, we use a similar method to Moon and Joughin [75] where change is quantified from differencing the areas of two reference boxes and dividing by the glacier's width. The reference box is an open-ended rectangle whose upflow side is an arbitrary position above the grounding line and the sides define the glacier's width. The open-ended edge is the hydrostatic equilibrium boundary.

Basal Evolution
Highly detailed surface elevation data from different time periods provides the opportunity to evaluate the basal dynamics of the floating portion of a glacier. We conduct a modest study of the subglacial conditions by investigating localized hypsometry on the glacier's surface. In particular, we focus on the evolution of a basal crevasse based on the changes we observe from a depression located directly overhead. A surface depression forms over large basal crevasses as a response to hydrostatic adjustment. We observe from Landsat data spanning the 1970s to the present several depressions that remain visible for decades within the trunk and continue to persist as features in the ice surface as ice flows from Byrd Glacier to the Ross Ice Shelf's terminus. We know the Byrd Glacier depressions directly overlay basal crevasses from radar echogram data and observations of other basal crevasses/surface depressions that have manifested post-1979 data collection. Over time, the geometry of these depressions evolves (e.g., changes in concavity) and we assume this to be a result of the changing geometry of their underlying basal crevasses. The basal crevasse we focus on was 165 m tall in 2011 during the time of the CReSIS radar data collection. To estimate the amount of geometry change the crevasse has experienced from 1979 to 2017, we use a similar method to Le Brocq et al. [76] where the basal melt (b) is a product from the continuity equation: u is surface velocity, which we obtain from Scambos et al. [77], x is in reference to the distance the basal crevasse traveled over 38 years, and a, the accumulation rate, 0.27 m y −1 [78,79]. We do not have exact H values for these two time periods, so we assume that the column of ice overhead of the basal crevasse is in hydrostatic equilibrium and use Equation (2) to estimate ice thickness using the two period DEMs as h.
As with Equation (2), this region of Byrd Glacier is also blue ice [73], so we disregard the firn correction.

DEM Products
The point cloud generated from the 1978 historical aerial photographs contains ∼52 million points, and the point cloud from the 1979 mosaic contains ∼33 million points. Both DEMs cover ∼2400 km 2 over the lower trunk of Byrd Glacier (see Figure 6 for results). The extents of the two mosaics differ slightly due to differences in image overlap and flight path course. The SfM processing performs best with 60% image overlap [80,81] and terrain coverage consisting of exposed bedrock and crevassing. A significant portion of Byrd Glacier is covered in surface crevasses, so results for the majority of the main trunk are generated. However, for large expanses of the surrounding solely snow-covered terrain, the processed elevations were of incorrect values. These values were elevations that were out of range for the known surface heights and were eliminated manually in point cloud format within the CloudCompare software.
The resulting interpolated point data may be seen in Figure 7a Figure 8). One possible cause for the large mean elevation difference is the ellipsoidal transformation. The XYZ shift calculated from the standard Molodensky method has a total accuracy of 5-10 m [82,83]. Due to the coarseness of the point dataset's areal distribution and that it has no spatial coverage of stable terrain, we were unable to coregister the Brecher [14] discrete elevations to the SfM generated ones. Studies have shown that coregistering elevation data results in a possible alignment correction of 5-70% [27], making it likely that registering the 1980s processed data with the SfM elevations would exhibit a significantly smaller mean difference.   Figure 9a. The normalized median absolute deviation (NMAD)-a robust estimate of the mean's standard deviation [84,85]-is 7.04 m. Given that the distribution of elevation differences is quite linear, it makes since that this NMAD value is so close to that of the standard deviation [86]. Both SfM DEMs are also coregistered and compared to the present-day DEM over stable terrain. The mean difference for 1978 and 1979 is 2.06 m and 1.71 m with a standard deviation of 12.52 m and 14.33 m, respectively. The NMAD values are 6.09 m (1978) and 3.85 m (1979). It is likely that the NMAD values being closer to the mean elevation differences than to the standard deviations is an indicator that the outliers in the data are causing the increased standard deviation values. Outliers in Figure 9b,c are due to the present-day DEM along cliff edges, which is where REMA's SETSM algorithm has difficulty rendering accurate elevation values [23].

Velocity
Velocities derived from feature tracking of surface elevation slopes (Figure 10a) produces values along the main trunk that range from ∼600-875 (m y −1 ). These values agree well with Brecher [14] and the field study conducted at the same time as the aerial image collection [13].

Velocity Differences
We calculate the velocity differences by subtracting the slope velocity from the original velocity data. The mean and standard deviation of the velocity differences (N = 6974) are −21.93 (m y −1 ) and 61.63 (m y −1 ), respectively. The majority of the large difference values (±100 (m y −1 )) are located on the edge of the glacier, which is where several of the COSI-Corr results contain NoData gaps (Figure 10c). The outer edges of the glacier consist of slow-moving smooth surfaces where it is difficult to extract velocity data when using cross-correlation techniques [28]. These large difference values are likely the result of interpolating the slope velocities, but may also be due to the Brecher [14] velocities which were estimated without registering one dataset to the other. We remove all differences outside two standard deviations (N = 6447) making the new mean −9.68 (m y −1 ) and the standard deviation 37.64 (m y −1 ) (see Figure 11). Unfortunately, there are no available satellite images or other remotely sensed data of Byrd Glacier between December 1978 and January 1979 with which to validate both sets of velocity results.

Grounding Zone
We find the present-day hydrostatic equilibrium boundary to have migrated slightly down-flow of the 1979 one. Both present-day and historical equilibrium boundaries are mapped in Figure 12. We calculate an overall mean difference of 390 m in down-flow migration. The stability of Byrd Glacier's grounding line is of particular interest because it was believed to have been receding since collection of the air photos [87]. Upon further investigation, it has been determined that the basis of these conclusions come from estimating grounding line migration by comparing the grounding zone's flexure point from 1978/79 field results with the hydrostatic equilibrium boundary calculated by Rignot and Jacobs [88] using 1997 satellite altimetry data. In assuming these two separate constraints of the grounding zone were the grounding line, Reusch and Hughes [87] concluded that the grounding line had migrated 20 km upstream. Were this the case, this much upflow migration would be an indicator that Byrd Glacier has been losing mass and is potentially unstable. Our results show an average migration of 390 m for the hydrostatic equilibrium boundary over the last ∼40 years. In Konrad et al. [70] calculation of actual grounding line migration, they estimate an even smaller migration rate of a 0.4 m yr −1 downflow shift. This amount of migration is minor and indicative of a glacier in mass balance. One possible cause for the altering boundary locations is due to changes in the ice bottom surface from new basal crevassing, evolving subglacial tunnels, and basal melt, which will alter an ice column's thickness. The ice bottom surface data were collected in the 2011-2012 austral summer, so it is expected that any changes in basal crevasses or melt would mean the basal surface from 1979 would be altered 33 years later. The basal features present at the grounding line in 1979 are now ∼30 km downflow and if fewer features have formed since then, the ice thickness 38 years later has increased, which will cause the boundary to migrate downflow. Having fewer deformational features along the basal surface-especially those resulting from periodic changes in ice dynamics (e.g., basal crevasses)-is another indicator that a glacier is in steady-state. It seems highly probable that the conclusion drawn by Reusch and Hughes [87] was incorrect and likely due to comparing two separate components of the grounding zone. Our results demonstrate that little to no change of Byrd Glacier's grounding zone has taken place over the last ∼40 years and the glacier is exhibiting stable condition.

Basal Evolution
The surface depression in Figure 13 migrated ∼30.5 km over the 38-year period between data collections. The ice column thickness calculations estimate 735.0 m in 1979 and 610.5 m in 2017. The depression depth in 1979 is ∼37 m and it shallows to ∼24 m in 2017-assuming uniform change, that is a depth decrease of ∼0.34 m y −1 . The estimated melt rate is 5.4 m y −1 , which is a total of 203.5 m for the timeline of this study. The basal crevasse, seen in the echogram in Figure 13b, was ∼165 m in 2011; using the calculated melt rate, the crevasse has undergone a ∼40% reduction in height in 38 years.
Besides depressions, there are other distinct surface features on floating ice that are also indicators of basal features and processes. For instance, elongated, longitudinal surface striations on ice shelves are suggested to represent subglacial channels [76,89], which direct the course of fresh basal water outflow and increase the amount of basal surface exposed to ocean-induced melt [90]. By studying changes in geometry of these features, our results provide insight into fluctuating ice dynamics, subglacial hydrology, and ocean temperatures [91]. The surface depression in Figure 13 appears to be connecting two of these channels, which means that the underlying basal crevasse could be acting as a junction point for lateral freshwater flow. This is a first ever glance at evolving basal conditions for this long ofa temporal scale. Further analysis of these changes will yield a better understanding of the external forces influencing Byrd Glacier's dynamics and the broader context in which they fit in an increasingly warming climate.

Conclusions
We find that using manually extracted GCPs to apply absolute locations to historical aerial imagery in SfM photogrammetry yields accurate results without requiring in situ, evenly spaced ground control, or other traditional photogrammetric parameters like flight altitude and photo center XY. The stable terrain elevation comparisons of both historical DEMs with the present-day DEM have an average vertical difference of ∼1.89 m. The gridded DEM spatial resolution and error estimate are of comparable quality to digitally processed modern-day data so that the SfM historical results may be used in further analysis of glacier dynamics.
The process requires a large manual component that is nontrivial, but cost-effective, and the results are DEMs of spatial resolutions unavailable to researchers at the time of original analog processing. Reprocessing historical data with modern SfM methodologies has produced DEMs of glacier surfaces that have never before been observed in Antarctica. Such high resolutions allow for better analysis of glacier dynamic evolution by better constraining and validating time-dependent elevation modeling. We show this by estimating the hydrostatic equilibrium boundary of Byrd Glacier using both the historical and present-day DEMs. Our results contrast previous studies that concluded the grounding line had migrated significantly inland, which inferred Byrd Glacier that was a risk for being unstable. The high-resolution SfM results show that the grounding zone of Byrd Glacier has instead exhibited stability for the past ∼40 years. Unfortunately, for the rest of Antarctica, high-density aerial image coverage from a single flight campaign is rare; however, over 300,000 trimetrogon aerial images were collected from 1946-2000, offering at least partial coverage for most Antarctic glaciers, many of which have never before been investigated. By using long-neglected historical aerial imagery, scientists will finally be able to expand the temporal scale of elevation change for the Antarctic Ice Sheet from an average of ∼25 years [92] to a possible ∼75 years.
Supplementary Materials: The historical DEM datasets are available online at http://dx.doi.org/1 0.17632/j8bjhmfg4j.1. An additional document is available containing every file used in this study with a www.mdpi.com/xxx/s1 to access it.
Author Contributions: This study was conceptualized by S.F.C. and L.A.S. The digital processing was conducted by S.F.C. and L.G. The data validation was performed by S.F.C. H.H.B. did the 1980 processing and provided the original point elevation data. All authors contributed to the writing and the editing of this paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by NSF grants ANT-1255488 and ANT-1543530. The original air photo collection was funded by NSF grant DPP-7722204.

Acknowledgments:
A special recognition is paid to Tony Meunier of the USGS who found the lost aerial photos after investigating 40 film canisters. The aerial photography and Landsat 8 OLI data are freely available from the USGS/EROS's Earth Explorer database. The REMA dataset was funded by NSF-OPP awards 1543501, 1810976, 1542736, 1559691, 1043681, 1541332, 0753663, 1548562, 1238993, and NASA award NNX10AN61G and the 2 m strips may be downloaded from PGC's website. The CReSIS data were funded by NSF grant ANT-0424589, NASA grant NNX10AT68G, and by the University of Kansas. The raw data products are obtainable from their website. All data products produced in this research are available by request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.