In a comparative analysis between DOS and FLAASH corrected images for bathymetry mapping, DOS delivered superior performance with the most useful bands being coastal, green, yellow, and NIR 2 [149]. Malinowski et al. [150] analyzed patterns of localized flooding using 12 different decision tree (DT) methods comprising of different PBIA and GEOBIA approaches on WV-2 imagery. DOS delivered more accurate reflectance for water than ATCOR. The authors highlight that NIR 2 is the most water sensitive band in the WV-2 arsenal. This could have potential for mapping saturated snow facies, but in the current study, we observe that all facies can be mapped by BGRN1 with better accuracies than CYRN2. Collin et al. [149] compared Hyperion and WV-3 data using spectral combinations of BGR (3 wavelengths), BGRNIR (4 wavelengths), CBGYR (5 wavelengths), BGRNIRs (8 wavelengths), and full 16 band spectrums in an Artificial Neural Network (ANN) structure for mapping coastal systems. The authors suggest that most of the over-reflecting discrepancies between the sensors are a result of the coarser resolution of Hyperion when compared to WV-3. Between the spectral subsets, however, the 16-band full spectrum subset delivered maximum accuracy for WV-3. In the current study, however, the VNIR

(8 wavelength) subset delivered highest accuracy across all processing schemes and classification methods. One reason may be the little to no target information contained in the SWIR bands from SWIR 2 to SWIR 8. Only SWIR 1 could be used for delineating dirty ice, melting snow, and melting glacier ice. However, this was performed to utilize the SWIR band in the combination and not because it could not be alternated by a VNIR alternative. When matching field spectra to WV-2 bands to map asphalt roads, blue, green, red, and NIR 2 reported maximum usability [151]. Variations in roofing materials were distinguishable from the yellow band onwards using pansharpened WV-2 imagery [152].

Karimi et al. [31] mapped clean ice, periglacial debris cover, and supraglacial debris cover using 4 band WV-2 (CBGR) and Landsat-TM imagery. The authors performed an intersensor pansharpening by enhancing Landsat-TM using PAN WV-2 data, and QUAC was used to derive reflectance. The 4 band WV-2 imagery is reported to be effective in classifying ice due to its high reflectivity. The authors also differentiated the temperatures of shadowed areas, periglacial, and supraglacial debris using Landsat-TM. The goal of the current study is to only utilize the attributes of the optical imagery to determine the variations in output according to different processing schemes and classification methods. This allows for an independent robust analysis of the image itself in varying scenarios. Moreover, manual digitization of the glacier boundary in the current study did not necessitate distinction between peri- and supraglacial debris. Tiwari et al. [153] utilized WV-2 imagery solely for data interpretation purpose to identify features in the ablation zone of the Bara Shigri glacier to facilitate supraglacial debris mapping using ASTER thermal, GDEM, and TerraSAR-X imagery. Buhler et al. [154] identified wet snow and wind transported snow using WV-2 imagery. However, their analysis was performed on non-atmospherically corrected data. The NIR 2 band which was suggested to enable greater differentiation between snow surfaces did not individually enhance classification in the current study but aided in the differentiation between dry and wet snow (whenever available). Jawak et al. [155] extracted snow cover using a comparative analysis of 14 spectral index ratios on principle component sharpened and GS sharpened WV-2 imagery. Their spectral indices using the additional bands yellow, red edge, and NIR 2 were most significant in extracting snow cover. However, in the current study, only NIR 2 was moderately effective. Yellow and red edge did not add any significant improvement and could be easily replaced by green or red band. Gray et al. [156] highlighted an interesting application of WV-2/3 imagery for mapping green and red snow algae. The imagery was converted to reflectance using the 6S radiative transfer model. Algal distributions documented on field were used to train SAM classifier to distinguish between green, red, dirty, and clean snow. The authors highlighted that processing routines can induce biases in the final algal distributions. This is true even in the case of the current study. Modulations in spectral reflectance due to changing atmospheric corrections, resolution improvement by pansharpening, selection of classification algorithm (PBIA), and even operator skill (GEOBIA) play a large role in the final thematic classifications.

#### *4.5. Final Inferences and Recommendations*

The statistical evaluation tabulated and described in Sections 4.1 and 4.2. is summarized and illustrated to outline our findings in Supplementary Figures S1–S4. Among atmospheric corrections (Supplementary Figure S1), we find that DOS (OA = 0.64) is the most reliable correction across all image subsets. When analyzed between PBIA and GEO-BIA, DOS is preferred for PBIA, whereas QUAC is most suitable for GEOBIA. Among the PBIA classification methods, MXL and WTA work best with DOS. In GEOBIA, rule sets 2 and 3 work best with QUAC, whereas rule set 1 performs well with FLAASH. Between the spectral band combinations classified by PBIA, VNIR\_DOS is the most accurate, whereas in GEOBIA, BGRN1\_QUAC, VNIR\_DOS, and VNIR\_QUAC yielded maximum OA of 0.83. Among pansharpening methods (Supplementary Figure S2), HCS (0.64) delivers an overall better performance than GS. However, this is due to the aggregated extremely poor performance of GS in PBIA. In GEOBIA, both GS and HCS improve accuracy, with GS

delivering the maximum OA of 0.79. Among the GEOBIA rule sets, rule set 3 shows the maximum improvement in accuracy with GS scoring an OA of 0.86. Among the PBIA classifiers, MHD yields the highest OA with HCS. Between the spectral band combinations, the VNIR\_HCS combination is the best, but with a very low OA of 0.29. The VNIR\_GS and BGRN1\_GS combination via GEOBIA classification yields the maximum OA of (0.79). A comparison between GEOBIA rule sets and PBIA classifiers (Supplementary Figure S3) outlines that GEOBIA is the most accurate with rule set 3 delivering the highest OA of 0.86. Within the PBIA algorithms, MXL delivers the maximum OA of 0.61. Overall findings of the spectral band combinations (Supplementary Figure S4) suggest that VNIR is the best performing subset with an aggregated OA of 0.59. The same trend is maintained among the PBIA and GEOBIA classifications with the VNIR subsets achieving OA of 0.37 and 0.81, respectively.

In Supplementary Figure S5, we provide processing methodologies which may be followed when using specific spectral band combinations with the tested image processing and classification methods. The most accurate method is the VNIR\_VHR subset atmospherically corrected by QUAC, pansharpened by GS, and classified in GEOBIA by a combination of spatial and spectral object features using multiresolution segmentation and rule set 3 (Paper 2). The most efficient method is found to be the VNIR\_VHR subset atmospherically corrected by DOS and classified by MXL.

#### *4.6. Significances and Limitations*

In the current study, we find that the DOS atmospheric correction performs best for VNIR imagery using PBIA. When utilizing GEOBIA, the QUAC correction may yield satisfactory results, but the DOS correction outperforms all in the overall analysis. Pansharpening does not improve PBIA but causes a significant improvement in GEOBIA classifications. Rule set definition must include both spatial and spectral attributes. SWIR data were the most useful for mapping dirty ice, but ultimately did not improve the classification process. Scene-based adjustments will always be profound when trying to extract maximum information from satellite data (Paper 2 [60]). Moreover, shadowed snow has continued to remain a problem (Paper 2 [9]) even with the addition of SWIR bands. Mixed pixel classifications in PBIA are a persistent phenomenon in an image with highly diverse spectral [152]. Overlap of spectral signatures can also cause misclassifications [140]. In glacier facies mapping, the main driver is the accurate characterization of spectral features into consequent classes. This places an onus on the operator for effective identification and comparison with literature (in absence of field data) and with field spectra (whenever available). PBIA can also result in fractured nature of classification due to spectral complexity on VHR imagery [157]. This can be reduced by GEOBIA [150]. The key for effective GEOBIA is segmentation. The multiresolution algorithm parameters selected here were applicable across all image subsets and processing schemes. This was constant even when the thresholds and features used for mapping facies were not consistent. The 'optimal' segmentation parameters selected here may be transferable and need to be tested further. The current study acknowledges the lack of field data for verification. However, the equalized random sampling approach utilized here [5] ensures a robust analysis. The classification of facies and their manifestation in the region was compared with previous studies to qualitatively assess the mapping methods. Furthermore, the focus of the study is the compounding effect of preprocessing routines on the resultant classification of glacier facies using GEOBIA and PBIA, not the most accurate map of glacier facies. The major limitation is the time needed to manually design rule sets. Moreover, the cognitive skill and performance of an analyst may play a role in the accuracy of rule set-based classification. While an existing rule set may be tested exactly for its transferability, development of a rule set for the same facies may be determined according to the differing knowledge, skills, and psychological variability of the analyst [31]. To the best of our knowledge, this is the first study to test the effects of GEOBIA and PBIA using various image processing techniques for mapping glacier surface facies using VHR satellite data. While the current results are tested

on VHR VNIR\_SWIR data, the results of the study must be tested on medium resolution satellite data for asserting the scalability of each rule set and classifier. We have not utilized an ancillary data as in-scene methods utilizing the properties of the available image to the maximum are necessary for understanding how facies can be mapped in different scenarios as well as reducing the reliance on ancillary data for accurate results. Moreover, it is well known that scene-based adjustments are necessary for implementing the same mapping method on different satellite images over different areas. Our results highlight that, in addition to sensor- and glacier-based adjustments, image processing-based adjustments are also necessary due to the impact of different methods of atmospheric correction and pansharpening on the properties of the pixels/objects of the same image.

#### Future Directions

The SVM classifier, although not tested here, is shown to have excellent capacity for detailed information extraction [136]. In this study, all indices for GEOBIA were created through a trialand-error approach, while this is effective for an exhaustive study testing imagery and processing routines without field data. Shorter applications benefiting from field spectra can utilize an optimum index factor (OIF), which filters spectral bands based on standard deviation and correlation coefficients to determine optimal wavelengths [158,159]. Samsudin et al. [160] used the OIF to select spectral bands for creating indices for mapping degraded roof materials through WV-3 imagery. Alternatively, Shahi et al. [151] used a stepwise discriminant analysis (DA) to select bands for creating spectral indices for identifying asphalt roads. Resolution enhancement in the current study performed established and tested methods of pansharpening. Vivone and Chanussot [161] outline the potential of hypersharpening, which fuses a multispectral or hyperspectral image with a higher resolution multispectral/hyperspectral image. This could be potentially explored for multilevel resolution enhancements to characterize supraglacial features. To produce facies maps similar in detail to USGS mineral maps, we need a combination of hyperspectral shape matching [162], field spectra, laboratory spectra [163], and an expert mechanism of multispectral and hyperspectral wavelength feature characterization [140]. Papers 1, 2, and the current work account for the image processing and classification approaches which may be most suitable for facies mapping. Next, we intend to expand the methodology and examine classifications of SVM, random forest (RF), and neural network classifiers (NNC) in continuation of the current work to deliver a comprehensive machine learning based assessment of glacier facies mapping for the study sites—beyond which the current work will involve field investigations to produce detailed maps and rigorous comparison against all previous attempts.

#### **5. Conclusions**

This study focused on analyzing the impact of pre-processing routines and information extraction methods for mapping glacier surface facies across two different study areas. Three atmospheric corrections and two pansharpening algorithms were tested on VHR VNIR WV-2/3 for Ny-Ålesund, and Chandra-Bhaga Basin, Himalaya. The atmospheric correction methods included DOS, QUAC, and FLAASH. The pansharpening methods included GS and HCS. Following pansharpening, all processing schemes were divided into four spectral band subsets. Subset 1 is BGRN1, subset 2 is CYRN2, subset 3 is the VNIR range, whereas subset 4 is the entire VNIR\_SWIR range. These subsets were then subjected to PBIA and GEOBIA classification methods. The PBIA process involved the utilization of conventional classifiers such as MHD, MXL, MD, SAM, and WTA; and advanced classifiers such as ACE, CEM, MF, MTMF, MTTCIMF, OSP, and TCIMF. The GEOBIA process involved multiresolution segmentation followed by development of three rule sets. Rule set 1 was developed using only object spectral information, rule set 2 was developed using only object contextual and spatial features, and rule set 3 was developed using a combination of both spatial and spectral features. The parameters used for segmentation consist of a scale parameter set at 5, shape set at 0.9, and compactness at 0.4. The features identified for classification are mean, quantile, standard deviation, min. pixel value, max. pixel value, edge contrast of neighbor pixels, number of overlapping thematic objects, relative

border to, and customized ratios/arithmetic features. Among the PBIA classifications, the MXL classification on VNIR subsets achieved the highest F1 score of 0.50. MTTCIMF is the worst classifier yielding the lowest accuracy. However, MTTCIMF increased in performance when pan sharpened using the HCS method, whereas all other classifiers show a decrease in accuracy when using pan sharpened data. For the PBIA classifications, the best performing image subset is the VNIR. The trend of reliable performance among the atmospheric corrections for PBIA is DOS (OA = 0.48) > QUAC (OA = 0.40) > FLAASH (OA = 0.38). Pansharpening appears to be detrimental to VHR PBIA of surface facies and is not recommended. In the case of the GEOBIA classifications, the reliability trend among atmospheric corrections is QUAC (OA = 0.82) > DOS (OA = 0.81) > FLAASH (OA = 80). Among the pansharpening methods, GS (OA = 0.79) is more suitable than HCS (OA = 0.76) for the GEOBIA classifications. Rule set 3 delivers the best classification result, whereas the subsets with the highest accuracy (common OA = 0.79) are the VNIR and BGRN1 subset. This highlights the effectivity of the conventional spectral bands for characterizing glacier facies through both PBIA and GEOBIA. While GEOBIA delivers an overall greater performance than PBIA, the latter is more efficient. GEOBIA is limited by the human intervention for developing the rule sets. However, the current study highlights that GEOBIA may be better suited to map glacier surface facies using VHR VNIR data through a combination of spatial and spectral attributes. The segmentation parameters were consistent across all processing schemes and subsets and may be transferable to other VHR VNIR based facies mapping applications. The mapping procedures outlined in the current study may be applied on medium resolution satellite data. However, sub-pixel classification may be necessary to enhance the characterization of facies in coarser resolutions [70]. The direct alternative would be to utilize fine spatial resolution imagery [82] such as that used in the current study combined with the proposed methodology to map facies.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14246311/s1, Table S1: The selected glaciers of the study, their areal extents, and GLIMS reference IDs. The extents were calculated from the delineated shapefiles using the geometry calculator in ArcGIS, Table S2: Abbreviations/nomenclature of processing schemes used in the current study. AC: Advanced Classifiers; CC: Conventional Classifiers. Table S3: Measures of accuracy used in the current study. TP: Samples are those that are in the positive class and are correctly classified, TN: Samples that are correctly classified as negative, FP: Samples that are not truly of the positive class but are incorrectly mapped as positive, FN: Samples that are mapped as negative when they are positive. Figure S1: Highlights the impact of atmospheric corrections progressively at each level of image processing, Figure S2: Highlights the impact of pansharpening progressively at each level of image processing, Figure S3: Displays the variations in overall accuracy of GEOBIA rule sets and PBIA algorithms. Figure S4: Recommendations of the most practical processing strategies if using any of the image processing routines and mapping methods.

**Author Contributions:** Conceptualization and experiment setup, S.D.J. and S.F.W.; methodology and manuscript design, S.D.J. and S.F.W.; software and processing, S.D.J. and S.F.W.; validation, S.D.J. and S.F.W.; result analysis and discussion, S.D.J. and S.F.W.; resources, A.J.L. and K.B.; writing—original draft preparation, S.F.W. and S.D.J.; writing—review and editing, S.D.J., A.J.L. and K.B.; visualization, S.F.W.; supervision, S.D.J., A.J.L. and K.B.; project administration, S.D.J., A.J.L. and K.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Freely available data used in the current study-(1) ASTER GDEM v2. Downloaded from: Gdex.cr.usgs.gov/gdex/ (accessed on 2 February 2017). The data are now moved to GDEM v3: (reviewed on 12 March 2022) ASTER GDEM is a product of Japan's Ministry of Economy, Trade, and Industry (METI) and NASA. (2) Arctic DEM. Available online: Pgc.umn.edu/ data/arcticdem/ (accessed on 21 January 2019). Sentinel 2A (10 m)-Himalaya: Copernicus Sentinel data 2021, processed by ESA; Svalbard: Copernicus Sentinel data 2022, processed by ESA). Natural Earth: Free vector and raster map data @ naturalearthdata.com.

**Acknowledgments:** The authors would like to thank DigitalGlobe/MAXAR for providing highresolution satellite data. Chandra–Bhaga Basin WorldView-2 © 2014 Maxar; Ny-Ålesund WorldView-3 © 2016 Maxar. The authors also thank HOD, Civil Engineering, and Director, MIT, MAHE for their support. S.D.J and A.J.L would like to thank M. Ravichandran, Secretary, MoES and Thamban Meloth, Director NCPOR for their encouragement and support. The authors thank the anonymous reviewers for their constructive criticism which helped develop and improve the manuscript.

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

#### **References**

### *Article* **Combined Use of Aerial Photogrammetry and Terrestrial Laser Scanning for Detecting Geomorphological Changes in Hornsund, Svalbard**

**Małgorzata Błaszczyk 1,\*, Michał Laska 1, Agnar Sivertsen <sup>2</sup> and Shridhar D. Jawak <sup>3</sup>**

**Abstract:** The Arctic is a region undergoing continuous and significant changes in land relief due to different glaciological, geomorphological and hydrogeological processes. To study those phenomena, digital elevation models (DEMs) and highly accurate maps with high spatial resolution are of prime importance. In this work, we assess the accuracy of high-resolution photogrammetric DEMs and orthomosaics derived from aerial images captured in 2020 over Hornsund, Svalbard. Further, we demonstrate the accuracy of DEMs generated using point clouds acquired in 2021 with a Riegl VZ®-6000 terrestrial laser scanner (TLS). Aerial and terrestrial data were georeferenced and registered based on very reliable ground control points measured in the field. Both DEMs, however, had some data gaps due to insufficient overlaps in aerial images and limited sensing range of the TLS. Therefore, we compared and integrated the two techniques to create a continuous and gapless DEM for the scientific community in Svalbard. This approach also made it possible to identify geomorphological activity over a one-year period, such as the melting of ice cores at the periglacial zone, changes along the shoreline or snow thickness in gullies. The study highlights the potential for combining other techniques to represent the active processes in this region.

**Keywords:** structure-from-motion; terrestrial laser scanning; digital elevation model; Svalbard; SIOS

#### **1. Introduction**

The Hornsund area in southern Spitsbergen (Svalbard) is the focus of wide glaciological [1–4], hydrological [5–9], snow [9], permafrost [10,11], geomorphometric [12–14] and biological [15] studies performed due to intensive alterations in the surveyed terrain. Digital mapping of such an environment and its related phenomena has become a fundamental requirement to keep track of all such alterations [16]. The most up-to-date, accurate and precise digital elevation models (DEM) for the study area are generated based on high-resolution satellites and aerial photos [17]. From these, the DEM with the highest spatial resolution (2 m) and accuracy (standard deviation = 0.6 m) is derived from aerial photographs taken in 2011. However, to monitor climate-induced change in the territory, there is a need for more frequent and accurate data on terrain elevation [17].

The aim of the paper is to examine and combine two DEMs produced from aerial image sets and terrestrial laser scanning (TLS). All data created within the project are further available to the scientific community through the Svalbard Integrated Arctic Earth Observing System (SIOS). Aerial images for the studied area were provided by SIOS through a dedicated call of proposals (https://sios-svalbard.org/AirborneRS (accessed on 23 December 2021)). One of SIOS's missions is to reduce the environmental footprint of scientific observations in Svalbard [18]. To achieve this mission, SIOS supports and

**Citation:** Błaszczyk, M.; Laska, M.; Sivertsen, A.; Jawak, S.D. Combined Use of Aerial Photogrammetry and Terrestrial Laser Scanning for Detecting Geomorphological Changes in Hornsund, Svalbard. *Remote Sens.* **2022**, *14*, 601. https:// doi.org/10.3390/rs14030601

Academic Editor: Michael Lim

Received: 23 December 2021 Accepted: 25 January 2022 Published: 26 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

coordinates the usage of unmanned aerial vehicles (UAVs) and a crewed Dornier aircraft to acquire imagery and hyperspectral data for the Svalbard research community to support scientific projects (https://sios-svalbard.org/AirborneRemoteSensing (accessed on 23 December 2021)). The project started in 2020 and so far, no publications have presented the accuracy of the data collected within the framework of the project. In this study, we examine and present the accuracy of products such as orthomosaics and DEMs obtained during the mission. The aerial data over two Hornsund zones, the Fuglebergsletta and Werenskioldbreen areas, were collected in June 2020 (Figure 1). The mission took place at the beginning of the ablation season, when some snow cover was present in the numerous gullies. Additionally, the flight level limited the sidelap of images, causing data gaps in the DEM. To address these limitations, a long-range TLS campaign was performed in August 2021. Long-range terrestrial laser scanning is an emerging method for the monitoring of complex and rough terrain such as mountain slopes, outcrops and deformations [16,19–23]. Here, we used the ultra-long laser scanner Riegl VZ®-6000, which has been successfully used in cryosphere studies such as glacier mass balance measurements in China [24,25]; the mass balance of very small glaciers in the Swiss Alps [26]; snow distribution at a glacier located in the Ötztal Alps, Austria [27]; glacier snowline determination [28]; or relationships between different climate forcing and flows of Helheim Glacier, Greenland [29]. However, at this point there are not many studies on the accuracy of ultra-long laser scans [23,26,30], especially over 2 km [27]. Further, to our knowledge there has not yet been a study published concerning the quality control of TLS-collected data with the Riegl VZ®-6000 instrument in the polar region. To fill this gap, we analyzed the accuracy of the relative and absolute registration of four scans collected over the complex Fuglebergsletta area, the most studied area near the Polish Polar Station.

**Figure 1.** Location of the study area. The blue outlines represent the extent of aerial imagery acquired during the crewed aircraft campaign in 2020; the red polygon presents the range of data taken by terrestrial scanner Riegl VZ-6000 in 2021; the red asterisks illustrate the positions of the laser scanner; the orange dots represent ground control points and the yellow dots represent checkpoints measured by GPS in 2021; the blue rectangle presents the location of the Polish Polar Station, Hornsund.

Terrestrial laser scanning (TLS) was next applied to complement the aerial photogrammetry. The integration of 3D modeling techniques is advantageous for obtaining the most complete and useful object coverage for many application areas [19]. Integration of aircraftbased, helicopter-based and UAV photogrammetry with TLS is popular, especially at large or inaccessible sites [21,23,30–32]. In this paper, digital aerial photogrammetry and TLS were combined to produce a continuous DEM with the highest possible resolution and accuracy. Thanks to that, this DEM, based on aerial photos and scanning (TLS), can be successfully applied to many environmental applications, such as hydrology modeling, glacier change detection, quantifying depositional and erosional processes in dynamic and complex fluvial systems, the evolution of the landscape, geomorphology and the monitoring of landslide displacement [13,33].

#### **2. Study Area**

The study area is located in the southern part of Spitsbergen, Svalbard (Figure 1). The area is characterized by diverse surface coverage such as relatively flat topography, mountains, tidewater glaciers and land-based glaciers. Narrow coastal plains with raised marine terraces surround the fjord shores [34]. Further inland, mountain massifs range in elevation up to 763 m above sea level (a.s.l.). The steep rock walls are cut by a system of deep chutes and gullies [13]. Two regions are analyzed in this paper. The eastern part is Fuglebergsletta, containing the Fugleberget massif and Hansbreen tidewater glacier; the retreat of the Hansbreen caused the exposure of wide lateral moraine ridges with buried glacial ice [35]. The western studied part is the Werenskioldbreen area; the landbased Werenskioldbreen glacier's foreground is flat with active glaciofluvial landforms and moraine. The vicinity of the Polish Polar Station, located on the shore, makes these two areas the focus of numerous environmental studies.

#### **3. Materials and Methods**

#### *3.1. Aerial Imagery*

#### 3.1.1. Data Preprocessing

Imagery acquisition over the research area was carried out during a SIOS crewed aircraft campaign on 22 June 2020. The aircraft, a Dornier DO228, is fitted with an RGB camera (Phasone IXU-150, Schneider LS 55 mm f/2.8) and a hyperspectral imager (VNIR-1800, Norsk elektrooptikk; [18]). The RGB camera and the hyperspectral imager can acquire images with a ground resolution of 0.1 m and 0.3 m from a flight altitude of 1000 m, respectively. During the Hornsund campaign, 622 high-resolution RGB photos were acquired. Out of these, 326 were related to the Fuglebergsletta area (covering Polish Polar Station and Hansbreen surroundings) and 296 to the foreground and frontal part of Werenskioldbreen (Figure 1). The images were all acquired at a flight altitude of 1000 m a.s.l., or approximately 935 m above ground level, resulting in a ground resolution of approximately 0.085 m per pixel. Due to the variable lighting conditions during the mission, preliminary exposure compensation was performed (Figure 2). The entire dataset was converted from the RAW format (48 bit IIQ) to JPEG (24-bit).

**Figure 2.** Processing workflow for airborne data and terrestrial laser scanning (TLS).

#### 3.1.2. Ground Control Points and Checkpoints

There was no opportunity to place dedicated ground control points (GCPs) and checkpoints before the flight campaign. Therefore, natural points were measured during the fieldwork in 2021, between 13–18 August (Figure 1). The vast majority concerned centroids of the characteristic large boulders which could have been correctly defined in the previously taken aerial photos. Fifty-nine points were measured in WGS 84 UTM 33X projection by Leica GX1230GG receivers. Measurements were performed at GNSS post-processing mode (30 min sessions, Phase Fixed solution) and RTK-GNSS mode, with corrections from the reference station at the Polish Polar Station (PPS). The GCPs at Fuglebergsletta were within the 4 km range from the PPS. The GCPs at the Werenskioldbreen area were 11–12.5 km from the PPS. The average horizontal and vertical accuracy related to the postprocessing measurements was 0.0003 m and 0.0008 m, respectively. The average horizontal and vertical accuracy related to the RTK mode was 0.01 m and 0.024 m, respectively.

#### 3.1.3. Data Processing and Quality

The attitude and position of the camera were recorded for each image, using an onboard navigation system (Applanix POS AV 410). The flight data, together with the GNSS correction measurements from the Polish Polar Station, were postprocessed using Applanix PosPac MMS 8 software. The total accuracy of the postprocessed camera position and orientation were within 0.025 m and 0.01 degrees, respectively. The aerial images were processed in Agisoft Metashape 1.7.5 software (https://www.agisoft.com/ (accessed on 23 December 2021)), using the postprocessed camera position and orientation for each image. The structure-from-motion method was used to obtain high-resolution DEMs and orthomosaics (Figure 2). Structure from motion (SfM) is a photogrammetric technique for estimating three-dimensional models from two-dimensional images, collected in proper overlap and coupled with local motion signals [36]. It is widely used in modern photogrammetry, especially low-level, where UAVs are applied [37].

At the following stages of data processing, the images in full resolution (8280 × 6208 pixels) were applied. Next, 30 GCPs were used in the aerotriangulation process (Figure 1). In the

Fuglebergsletta area, the selected 27 GCPs were characterized by an appropriate distribution, presented in Figure 1. The GCP density was 1.17 GCP/km2. In the Werenskioldbreen area, only three GCPs were measured, mainly located in the southern part of the study. The GCP density was only 0.11 GCP/km2. The aggressive depth filtering was used to generate dense point clouds, which produced the best results among the available levels. The remaining artefacts were removed manually. Due to the negligible number of vertical objects and the use of very high-quality images, estimated by the software at above 0.80, the total number of removed points was insignificant, below 0.1%.

The remaining 29 points measured by GPS, which did not participate in the aerotriangulation process, were used as checkpoints (Figure 1) to validate the vertical accuracy of DEM and horizontal accuracy of the orthomosaic generated for Fuglebergsletta. The point cloud data were exported into the grid format (\*.tif) and compared against the GCPs in ArcGIS software.

#### *3.2. Terrestrial Laser Scanning (TLS)*

#### 3.2.1. Long-Range Terrestrial Laser Scanning

The TLS survey was carried out on 15 August 2021, using the Riegl VZ®-6000 longrange terrestrial laser scanner (Figure 1). The scanner has an effective maximum range of around 6 km, operates at a wavelength of 1064 nm and uses class 3B laser beams, so can be used in snow and glacier studies thanks to high rates of reflection (80%) from snowand ice-covered terrain [23,24,29]. The result of a single laser scan is a large quantity of 3D data, usually of the order of several million point measurements, each with an x, y, z and intensity value. This dataset is termed a point cloud and is the raw 'product' of scanning [38]. The laser footprint size estimated from the beam divergence (0.12 mrad for Riegl VZ®-6000) gives spot sizes of 15 mm at 60 m, 120 mm at 1000 m, 240 mm at 2000 m and up to 720 mm at distances of 6000 m [39].

Although the range of 6 km theoretically covers most of the area of interest, in order to increase the density of scans in the far field and reduce shadowing problems [23], a total of four scans were acquired from various positions to cover the area as uniformly as possible. The scanning frequency was set to 30 kHz and column and line resolutions to 0.002◦. Riegl instrument uses laser light, so is relatively slow compared to the phase difference scanners [40]. Therefore, the number of the positions and scan resolutions were chosen to take advantage of good atmospheric conditions, which is not common in harsh Arctic conditions. Associated phenomena affecting long-distance estimations are air temperature, pressure and humidity. Therefore, we used meteorological data from a nearby meteorological station at the Polish Polar Station during measurements.

#### 3.2.2. Point Cloud Registration

Relative and absolute positioning of scans [38] (Figure 2) was performed using RiSCAN Pro 2.12.1 (http://www.riegl.com/products/software-packages/riscan-pro/ (accessed on 23 December 2021)). The maximum measured distance by TLS within the study area was about 7.8 km. However, to minimize the effect of the laser footprint size for long distances, during the registration process we used only parts of point clouds within the range of 3 km from the laser scanner position. The data were also filtered to remove noises, and the unstable area of the glacier was omitted during the registration of the scan.

Fine cloud-to-cloud registration was performed by a Multi-Station Adjustment (MSA) tool in RiSCAN Pro [32]. The software identifies common plane patches from different scan positions, links them together and minimizes the errors between all these linked plane patches by using an iterative matching algorithm [41]. The standard procedure in MSA uses a plane patch filter [32]; however, our study used triangulated polydata to increase the number of planes in the registration process. As a result, all point clouds from the four stations were merged into a single point cloud.

An absolute registration was performed in the last step to fit the merged point clouds into the reflectors scanned during the field campaign. Five Riegl flat reflectors (50 mm diameters) were distributed spatially as evenly as possible around each of the four scanner positions and measured using RTK-GNSS in the UTM 33X zone. Leica GX1230GG receivers operating in the RTK mode used corrections from the reference station at the Polish Polar Station. The average horizontal and vertical accuracy of RTK GNSS measurements of all 20 reflectors was 0.01 m and 0.02 m, respectively.

#### 3.2.3. Validation of TLS

Absolute positioning allows laser scans to be integrated with the other registered data [38]. However, before such integration, validation of the data should be performed. The densities of TLS point clouds strongly depend on the distance to the scanning system and the chosen scan frequency [23]. Therefore, to uniform the data, point clouds were first filtered with the Octree algorithm in RiSCAN Pro 2.12.1 and then exported to the LAS format. Next, the data were imported to ArcGIS, converted into the raster format and compared against the 19 GCPs and checkpoints measured in 2021 (see Chapter 3.1.2) which were captured during scanning.

#### *3.3. Integration of Aerial and TLS Based Data*

To assess the difference between aerial-based and TLS-based DEMs, we used the freely available Multiscale Model-to-Model Cloud Comparison (M3C2) plugin for opensource CloudCompare software [42]. M3C2 is an algorithm used for multitemporal point cloud distance calculations. It estimates local positions in two input point clouds by using the surface normal vectors to determine the median point within a cylinder of defined radius [42,43]. For the point-cloud-based strategy, the M3C2 algorithm is chosen as the bestestablished method, especially in earth sciences, when dealing with irregular surfaces [44]. Here, we used the algorithm to compare the data and find areas where differences between both point clouds were significantly larger than the estimated accuracy of the DEMs generated from those point clouds. Buildings and noises were removed from the point clouds. Normal calculations were performed at a fixed scale (D = 2 m) using the core point file.

Both DEMs contain data gaps; therefore, we combined both DEMs to create a continuous gapless product with the best possible accuracy and resolution. Missing data in aerial-based DEM stem from too low a sidelap of images. TLS data gaps arise from two primary sources: a line-of-sight obstacle resulting in occlusion (holes in the data as some foreground interfered with the scanner's line of sight) and a dropout resulting from a specular reflective or absorbent surface, preventing the energy from a given laser pulse from returning to the TLS instrument [38,45]. Furthermore, TLS data has inhomogeneous point density throughout the area [16,22,30,42,46]. Thus, we used the aerial-based DEM with relatively homogeneous point distribution as the reference dataset. Next, the data gaps and areas covered by snow were replaced by the data from TLS. In a few small zones without any data, TLS data were interpolated using the Inverse Distance Weighting (IDW) algorithm [16].

#### **4. Results**

### *4.1. Digital Elevation Model and Orthomosaics Based on Aerial Imageries*

#### 4.1.1. Fuglebergsletta

As a result of aligning the images, 325 from 326 were correctly aligned. Camera locations and image overlaps are presented in Figure 3A. Details on data processing are shown in Table 1. The final resolution of the DEM equals 0.169 m and the point density was 35.2 points/m2. The root mean square error for all 27 GCP locations calculated by the software was 0.0018 m. Based on the improved dense clouds, digital elevation models were generated. The DEM was exported to raster format with 0.169 m resolution, the best possible resolution with this software, to increase the possibility of mapping even the smallest geomorphological features [16]. The vertical quality of the DEM was assessed based on 29 GCPs that were not used in the aerotriangulation process (Figure 4A). The standard deviation equaled 0.14 m, with a mean value of −0.22 m and a median value of −0.19 m. The maximal vertical error of the DEM equaled 0.54 m. There is a small systematic error in the DEM. In general, the elevation of the DEM generated in Agisoft is slightly higher than the measured elevation of the checkpoints.

**Figure 3.** Camera locations (black dots) and image numbers and overlaps over Fuglebergsletta (**A**) and Werenskioldbreen area (**B**).

**Table 1.** Processing settings and final resolution of the products.

After the final colour calibration of the entire dataset, an orthomosaic was produced (Figure 4B). The resolution of the orthomosaic was two times higher than the DEMs, which were the basis for its creation, and equalled 0.0843 m. The horizontal accuracy of the orthomosaic was estimated based on the same GCP points as those used for the DEM accuracy assessment. The standard deviation of the horizontal accuracy equalled 0.10 m, with a mean value of −0.12 m and a median value of −0.10 m.

In the low-image-overlap areas, the insufficient color blending over merged images may be noticeable. Even using color and white balance calibration in postprocessing, this effect was difficult to eliminate. For flat areas located below c. 140 m a.s.l., the applied flight parameters were sufficient to assemble the models. However, the mountain slopes and peaks on both study sites were characterized by either insufficient or a complete lack of coverage (Figure 3A). This resulted in broad blank areas (Figure 4) in the DEM and orthomosaic. Parts of the model composed only of a pair of images occasionally presented coarse graining or a deep generalization of elevation data (Figure 5A). Similar effects occur even with multiple overlaps on watercourses, surface water (Figure 5A) and the sea. For this reason, the DEM shows negative point values, reaching –15 m. Artefacts were also generated by moving objects, such as animals (Figure 5B).

**Figure 4.** Digital Elevation Model (**A**) and orthomosaic (**B**) for the Fuglebergsletta area. Orange points represent GCPs used during processing in Agisoft, while yellow points served as checkpoints to estimate the quality of the products.

**Figure 5.** Artefacts described in the text: coarse graining or a deep generalization of elevation data (**A**), artefacts generated by reindeers (**B**) and noises over water bodies (**A**,**C**).

#### 4.1.2. Werenskioldbreen Area

As a result of aligning the images, 285 out of 296 photos for Werenskioldbreen were aligned correctly. Both the automatic and manual alignment functions failed to merge ten images representing the southern part of the glacier. It was caused by low overlap over the Angellfjellet ridge (594 m a.s.l.) and partial cloud cover over this area, making it difficult to detect tie points correctly. Camera locations and image overlaps for the second area are presented in Figure 3B. Details on data processing are described in Table 1. The final resolution of the DEM equals 0.174 m and the point density is 33 points/m2. The root mean square error for three points used in aerotriangulation was 0.0011 m. The resolution of the DEM in raster format and orthomosaics generated for the zone was 0.174 and 0.087 m, respectively (Figure 6). The resolution difference between orthomosaics for Werenskioldbreen and Fuglebergsletta was caused by the variation between mean ground level and flight altitude. Similar to the DEM for Fuglebergsletta, noises over water bodies (e.g., Figure 5C) are also present for the Werenskioldbreen area.

No checkpoints were applied to assess the final absolute vertical accuracy of the DEM and horizontal accuracy of the orthomosaic. However, we compared the DEM and orthomosaic generated using three points in aerotriangulation and without using any GCP. Results show significant improvement of the former products. The DEM accuracy measured over the three GCPs increased from 1.77 m to 0.22 m. The horizontal accuracy of the orthomosaic increased from 1.29 m to around 0.1 m (the orthomosaic was shifted towards the southeast).

**Figure 6.** Digital Elevation Model (**A**) and orthomosaic (**B**) for the Werenskioldbreen area. Orange points represent GCPs used in processing in Agisoft.

#### *4.2. TLS*

In this study, we only used point clouds within 3 km from the laser scanner position (Figure 7; see Chapter 3.2.2). The four scans were registered using nearly 540,000 plane patches detected throughout the area (Figure 7). The final point cloud had over 117 million points, with a density of around 500–600 thousand points per m<sup>2</sup> in the near-field region of the scans and just a few points per m2 at the distance of 3 km from the scanner. The RMS error obtained during the relative registration of all four point clouds was 0.09 m.

**Figure 7.** Terrestrial laser scanning data over the study area. The yellow polygon in (**A**) indicates the area covered by all scans, while the brown polygon presents the data range used in the registration procedure. Yellow asterisks in (**B**) illustrate the positions of the laser scanner. Orange and yellow dots represent GCPs and checkpoints that served to estimate the quality of the raster DEM generated from TLS data. Red dots in (**C**) present the example positioning of the RIEGL reflectors regarding the scanner position.

The final point cloud was then georeferenced utilizing 20 RIEGL reflectors (Figure 7C), with an RMS error of 0.13 m. The lowest accuracy in the northerly direction (0.2 to 0.3 m) was noticed for the two lateral scans, numbers 2 and 4. The lowest accuracy in the easterly direction was for scan number 1 (0.21 to 0.25 m). Height accuracy was high for all scans (of order 0.01–0.05 m).

The georeferenced cloud point was next filtered using an Octree algorithm in RiSCAN Pro 2.12.1 with a cell length of 0.16 m and converted into the raster format in ArcGIS (Esri, California) with 0.16 m resolution. As the areas of DEM interpolated in the shadow are characterized by lower accuracy, NoData values were assigned to the raster where holes in the data were present. Therefore, the DEM contains data gaps and keeps information only at spots measured directly by the laser scanner. The vertical quality of the TLS-based DEM was assessed based on the 19 independent points (Figure 7B). The standard deviation equalled 0.31 m and median value −0.19 m. The maximal vertical error of the DEM was 0.93 m.

#### *4.3. Integration of Aerial DEM with TLS*

The resulting M3C2-based point cloud distances are presented in Figure 8. Generally, the differences between both point clouds over the surveyed area were in the range of the vertical accuracy of the product (SD equalled 0.14 m and 0.31 m for aerial and TLS data, respectively). The map of surface differences shows local changes ranging from −4 m to 2 m (Figure 8A). The largest recorded differences, up to −4 m, are noted over the gullies and depressions filled by snow in 2020. Therefore, before integrating both DEMs, all areas covered by snow were removed from the aerial-based point cloud (Figure 8B,C). We also documented the erosion of lateral moraine of Hansbreen (around −0.7 m), which consists of an ice core covered with debris, as well as erosion along the shoreline. There are elevation changes of about 1 m concentrated over the mountain slope in the western part of the analyzed area. Numerous large boulders probably affect the quality and accuracy of the TLS-based DEM in that zone. As the scanner was located on marine terraces, around 300 m below, the laser beam was directed significantly upwards, and reflected off the wall of these boulders, rather than the top. Regardless of these local differences, we combined both DEMs into one continuous product. The final dataset (Figure 8D) is composed of a mosaic of two DEMs in raster format. All data gaps and snow areas in the aerial-based DEM were filled with the TLS-based DEM, with a blending option over a two-meter seamline. Moreover, we interpolated DEM values for a few remaining holes, where no data from both campaigns were present and interpolation of the DEM was made under the building and infrastructure of PPS (Figure 8D).

**Figure 8.** Comparison of the multitemporal point clouds: aerial-based point cloud from images taken in 2020 and TLS-based point cloud collected in 2021. (**A**) M3C2-calculated distance between aerial-based and TLS dataset; (**B**) orthomosaic presenting data gaps and snow cover over the land in 2020; (**C**) vertical difference of the point clouds with snow cover area eliminated from further data integration; (**D**) the final DEM integrated from both DEMs. Black polygons present the areas where DEM was interpolated.

#### **5. Discussion**

The two approaches to acquiring the terrain data discussed in this paper are helpful in studying the landform topography and different environmental processes in the Hornsund area. However, both technologies have their limitations. The best vertical accuracy was noted for the aerial-based DEM. However, the limited sidelap of the mission caused data gaps in DEM over mountain slopes. Such settings should be considered when planning subsequent missions in the future, especially when the area of interest is land relief of the mountain regions.

Further, as GCPs and checkpoints were not established prior to the photograph surveys, we chose large, outsize boulders to serve as calibration and validation points [47,48]. Using such existing 'stable' features is a limitation of studies that use photogrammetric methods to produce DEMs [48]. For the first studied zone, Fuglebergsletta, we used enough GCPs to generate the DEM and orthomosaic and assess their accuracy. For the Werenskioldbreen area, only three points would have been used to create the DEM and orthomosaic. Thus, more significant distortions in the northern part of the study must be expected due to the lack of GCPs in that part of the model. Notwithstanding this, we generated the best products with the best resolution and accuracy for the Hornsund area. Further, it is possible to supplement them with better-representing GCPs and recalculate the dataset in the future. We would also like to underline here that regardless of the vertical accuracy of DEMs, the very high resolution of the DEMs gives the possibility to map even the smallest geomorphological features, compute geomorphometric indices and carefully estimate geomorphometric parameters [16].

The accuracy of the TLS-based DEM is slightly lower than the aerial-based product. One of the properties that strongly influenced the long-range TLS was an inhomogeneous point density throughout the area of interest and positional uncertainty due to the laser beam width [22,30,42,46]. Laser footprint size at long 3 km distances equals around 0.36 m. Furthermore, on inclined surfaces the laser footprint becomes deformed to an ellipse [22]. The uncertainty caused by these large footprints are not yet quantified, but is expected to be in the order of decimeters and strongly depends on the distance from the scanner to the surface [27]. Laser scanner data are susceptible to data gaps in locations not in the direct line-of-sight of the scanner, resulting in 'range shadows' which inherently added uncertainty to the derived product [22,38,42,45,49]. Further, the configuration of the scan position and Riegl flat reflectors [30] was limited by terrain conditions, which would increase the error in relative alignment and absolute registration of scans. Considering all these limitations, the final vertical and horizontal accuracy of the TLS-based DEM of order 0.30 m is not unexpected to be in the lower-decimeter range [27] and is the best possible for the complex natural terrain. The DEM at interpolated areas has a quality lower than estimated in the study. In the future, viewshed analysis before TLS fieldwork could remove the shadowing effect [50,51] and improve the accuracy of the data.

The snow cover present in the aerial imagery taken in 2020 and data gaps in both point clouds led us to integrate both datasets [32,52]. Usually, within point-cloud-based strategies, the data acquired in at least two epochs are used to determine the geometric changes between them [30,43,44,52]. Here, apart from the accuracy assessment, we compared and combined the data acquired within one year. The distances between both point clouds were distributed unevenly over the whole area. The differences lower than the vertical accuracy of DEM (SD = 0.31 m) indicated the stable area. The final DEM for Fuglebergsletta was composed of the aerial-based DEM, with data gaps and the snow-covered regions filled by TLS data. One needs to keep in mind that the expected accuracy of such a product is lower than singular point clouds acquired over two years, as the combined DEM comes from the datasets collected in 2020 and 2021. Merging two DEMs from two different years is challenging, as the study area is undergoing continuous changes in land relief. However, temporal changes in most of the area are below the accuracy of both DEMs, and therefore the DEM can be successfully applied to hydrological, glaciological, biological and geomorphometric studies, when observed changes are larger than the accuracy of the

DEM. In other words, the merged DEM that we have generated represents the topography of the region's independent annual surface elevation changes, since the incorporated data were collected during the following two years. In the process of merging two datasets, the annual elevation changes may not reflect in the resulting DEM because of suppression of high-density elevation data recorded from TLS and photogrammetric methods [53,54].

To conclude, the final point clouds, DEMs and orthomosaic are the best terrain products with the highest spatial resolution for the Hornsund area and can be used for further environmental studies. They allowed us to identify local geomorphic processes, such as melting of ice cores at the periglacial area or shoreline erosion. However, for studies aiming to identify changes within a specified time range or slow geomorphometric processes on slopes, we recommend using data from the same types of sources (aerial/TLS) or stemming from a more extended period, such as DEMs generated from aerial imagery taken in 2011 [17]. In addition, we would like to underline here that all the DEMs should be first co-registered [55,56] and the final error must be defined to differentiate between real changes and data noise [22]. This is especially important for the DEM of the Werenskioldbreen area, as only three GCPs located along an almost-straight line were used in the data processing. This suggests that despite the high relative accuracy of the DEM, the absolute accuracy could be lower and decrease with the distance from measured GCPs; therefore, co-registration is needed when comparing with the other data.

#### **6. Conclusions**

The Hornsund area is the site in the focus of several studies due to the vicinity of the Polish Polar Station. It is a complex terrain under constant changes, although it lacks repeated datasets with spatial resolution high enough to map even the smallest geomorphological features. In the present study we derived aerial-based DEMs, TLSbased DEMs and orthomosaics, and assessed their accuracy and usefulness for further environmental studies.

Aerial data provided by SIOS were used to generate the DEMs and orthomosaics for the Fuglebergsletta and Werenskioldbreen areas. The products lack gaps over the mountain slopes due to the low sidelap during the flight. Nevertheless, the derived DEM has sufficient quality to study different geomorphometric features. Long-range TLS, despite limitations such as data gaps, complex terrain which limits scan positions and large laser footprint sizes at long distances, can also serve for future analysis of the different geomorphological processes. In particular, the area close to the scanner position with very dense point clouds can be used for future studies of 3D displacements and minor geomorphological features.

This work demonstrates that combining other techniques could be an option, especially in remote polar areas, when data acquisition depends on many factors. In order to remove the aforementioned holes in the datasets from aerial and TLS surveys, we integrated both DEMs over the area close to the Polish Polar Station. Comparison of the point clouds revealed a few small zones under local geomorphic processes, such as the melting of ice cores at the periglacial zone and changes along the shoreline. The differences in both products over the rest of the area were generally lower than the accuracy of both DEMs. This allowed us to combine the data and create the final continuous DEM without gaps for the Fuglebergsletta area. After the proper co-registration, all point clouds and DEMs can be successfully applied in further studies of landform evolution and hydrological, geomorphometric and other processes in the region.

**Author Contributions:** Conceptualization, M.B.; methodology, M.B., M.L., A.S. and S.D.J.; flight planning, mission execution and aerial data acquisition A.S., S.D.J., M.L. and M.B.; data curation, M.B. and M.L.; writing—original draft preparation, M.B. and M.L.; writing—review and editing, S.D.J. and A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The fieldwork has been financed by the Centre for Polar Studies of the University of Silesia, Poland. The SIOS KC would like to acknowledge the support from the Research Council of Norway (RCN) through projects No. 291644 & No. 269927 and the Norwegian Space Agency (NoSA) through

contract No. NIT 12.20.5. The Hyperspectral sensor onboard the Dornier aircraft is funded by the SIOS-InfraNor project (No. 269927).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. This data can be found here: DEM for Fuglebergsletta area from Dornier images: https://ppdb.us.edu.pl/ geonetwork/srv/eng/catalog.search#/metadata/69987fe3-f0ed-4fdc-b597-41dd5ae7d21a (accessed on 23 December 2021), Orthomosaic for Fuglebergsletta area from Dornier images: https://ppdb.us. edu.pl/geonetwork/srv/eng/catalog.search#/metadata/9b980572-2bcb-401f-a2e6-2a7a16428dd7 (accessed on 23 December 2021), DEM for Werenskioldbreen area from Dornier images: https://ppdb. us.edu.pl/geonetwork/srv/eng/catalog.search#/metadata/87a556ed-6b8f-4d06-9c20-1a5b6b23d280 (accessed on 23 December 2021), Orthomosaic for Werenskioldbreen area from Dornier images: https: //ppdb.us.edu.pl/geonetwork/srv/eng/catalog.search#/metadata/c99085cb-7e4c-43dd-866b-f2f9 09b44937 (accessed on 23 December 2021), Terrestrial laser scanner (TLS) data for Fuglebergsletta area: https://ppdb.us.edu.pl/geonetwork/srv/eng/catalog.search#/metadata/8bb03c41-35de-46 b5-87c1-bd6c82fa2914 (accessed on 23 December 2021), Combined DEM for Fuglebergsletta area from Dornier images and terrestrial laser scanner (TLS): https://ppdb.us.edu.pl/geonetwork/srv/eng/ catalog.search#/metadata/eec5dd69-20c0-46c6-aa16-8bb9cece3e85 (accessed on 23 December 2021).

**Acknowledgments:** The authors would like to acknowledge the Institute of Geophysics PAS and project POLARPOL for enabling the use of terrestrial laser scanner Riegl VZ®-6000. The studies were carried out as part of the scientific activity of the Centre for Polar Studies (University of Silesia in Katowice) with the use of research and logistic equipment (GNSS receivers and supporting equipment) of the Polar Laboratory of the University of Silesia in Katowice. The authors would like to thank Jacek Jania and participants of the polar expedition of the University of Silesia in Katowice 2021 for assistance with the fieldwork. Thanks are directed to the staff of the Institute of Geophysics PAS and the Polish Polar Station, Hornsund for assistance during the field campaign. All the authors acknowledge the SIOS Knowledge Centre (SIOS-KC) for funding and conducting flight campaigns over the study area.

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

#### **References**

#### *Article*

## **Airborne Validation of ICESat-2 ATLAS Data over Crevassed Surfaces and Other Complex Glacial Environments: Results from Experiments of Laser Altimeter and Kinematic GPS Data Collection from a Helicopter over a Surging Arctic Glacier (Negribreen, Svalbard)**

**Ute C. Herzfeld 1,\*, Matthew Lawson 1, Thomas Trantow <sup>1</sup> and Thomas Nylen <sup>2</sup>**

**Citation:** Herzfeld, U.C.; Lawson, M.; Trantow, T.; Nylen, T. Airborne Validation of ICESat-2 ATLAS Data over Crevassed Surfaces and Other Complex Glacial Environments: Results from Experiments of Laser Altimeter and Kinematic GPS Data Collection from a Helicopter over a Surging Arctic Glacier (Negribreen, Svalbard). *Remote Sens.* **2022**, *14*, 1185. https://doi.org/10.3390/rs14051185

Academic Editors: Shridhar Jawak, Andreas Kääb, Veijo Pohjola, Hiroyuki Enomoto, Geir Moholdt, Kjell Arild Høgda, Małgorzata Błaszczyk, Bo N. Andersen, Ann Mari Fjæraa, Bartłomiej Luks, Roberto Salzano and Frode Dinessen

Received: 4 October 2021 Accepted: 30 January 2022 Published: 27 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** The topic of this paper is the airborne evaluation of ICESat-2 Advanced Topographic Laser Altimeter System (ATLAS) measurement capabilities and surface-height-determination over crevassed glacial terrain, with a focus on the geodetical accuracy of geophysical data collected from a helicopter. To obtain surface heights over crevassed and otherwise complex ice surface, ICESat-2 data are analyzed using the density-dimension algorithm for ice surfaces (DDA-ice), which yields surface heights at the nominal 0.7 m along-track spacing of ATLAS data. As the result of an ongoing surge, Negribreen, Svalbard, provided an ideal situation for the validation objectives in 2018 and 2019, because many different crevasse types and morphologically complex ice surfaces existed in close proximity. Airborne geophysical data, including laser altimeter data (profilometer data at 905 nm frequency), differential Global Positioning System (GPS), Inertial Measurement Unit (IMU) data, on-board-time-lapse imagery and photographs, were collected during two campaigns in summers of 2018 and 2019. Airborne experiment setup, geodetical correction and data processing steps are described here. To date, there is relatively little knowledge of the geodetical accuracy that can be obtained from kinematic data collection from a helicopter. Our study finds that (1) Kinematic GPS data collection with correction in post-processing yields higher accuracies than Real-Time-Kinematic (RTK) data collection. (2) Processing of only the rover data using the Natural Resources Canada Spatial Reference System Precise Point Positioning (CSRS-PPP) software is sufficiently accurate for the sub-satellite validation purpose. (3) Distances between ICESat-2 ground tracks and airborne ground tracks were generally better than 25 m, while distance between predicted and actual ICESat-2 ground track was on the order of 9 m, which allows direct comparison of ice-surface heights and spatial statistical characteristics of crevasses from the satellite and airborne measurements. (4) The Lasertech Universal Laser System (ULS), operated at up to 300 m above ground level, yields full return frequency (400 Hz) and 0.06–0.08 m on-ice along-track spacing of height measurements. (5) Cross-over differences of airborne laser altimeter data are −0.172 ± 2.564 m along straight paths, which implies a precision of approximately 2.6 m for ICESat-2 validation experiments in crevassed terrain. (6) In summary, the comparatively light-weight experiment setup of a suite of small survey equipment mounted on a Eurocopter (Helicopter AS-350) and kinematic GPS data analyzed in post-processing using CSRS-PPP leads to high accuracy repeats of the ICESat-2 tracks. The technical results (1)–(6) indicate that direct comparison of ice-surface heights and crevasse depths from the ICESat-2 and airborne laser altimeter data is warranted. Numerical evaluation of height comparisons utilizes spatial surface roughness measures. The final result of the validation is that ICESat-2 ATLAS data, analyzed with the DDA-ice, facilitate surface-height determination over crevassed terrain, in good agreement with airborne data, including spatial characteristics, such as surface roughness, crevasse spacing and depth, which are key informants on the deformation and dynamics of a glacier during surge.

**Keywords:** ICESat-2; laser altimetry; kinematic GPS experiments; glaciology; surge glaciers; svalbard; density dimension algorithm for ice surfaces; airborne validation of satellite data

#### **1. Introduction**

With the Advanced Topographic Laser Altimeter System (ATLAS), NASA's ICESat-2, launched on 15 September 2018, carries the first space-borne multi-beam micro-pulse photon-counting laser altimeter system [1,2]. The ATLAS system records returns from every single photon in the 532 nm range of the sensor, which facilitates determination of alongtrack surface heights at the 0.7 m resolution of the sensor (under clear-sky atmospheric conditions), using the Density Dimension Algorithm for ice surfaces (DDA-ice) [3]. This capability of high-density measurements constitutes a better than 200-fold improvement in resolution over the 173-m spacing of the Geoscience Laser Altimeter System (GLAS) of the ICESat Mission, which operated from 2003 to 2009 [4,5].

Such increased density of measurements requires a field validation with high geodetical accuracy to establish that the surface heights determined from ICESat-2 accurately represent the actual surface heights and morphologies of an ice surface, especially for the complex morphology of heavily crevassed ice surfaces. Capturing crevassed surfaces accurately is important, because crevassing occurs as a sign of fast-moving and accelerating glaciers. Glacial acceleration is one of the largest sources of uncertainty in sea-level-rise assessment, according to the Fifth Assessment Report (AR5) of the Intergovernmental Panel of Climate Change (IPCC) [6]. The different spatial characteristics of crevasse fields, including crevasse spacing, depth and surface roughness, yield information on the deformation characteristics and ice dynamics during glacial acceleration [7–11].

To assess the surface-height determination capabilities of ICESat-2 ATLAS, we carried out an airborne geophysical evaluation campaign over Negribreen, Svalbard, during surge in summer 2019. The objective of this paper is a geodetical and glaciological evaluation using laser altimeter, image and geodetical data collected during this campaign. Negribreen accelerated to more than 200 times its quiescent-phase ice speed in July 2017, with a maximal speed of 22 m/day at the height of the surge acceleration [3]. The glacier's rapid mass transfer to the Arctic Ocean during the surge provides a prime example of the relationship between dramatic mass loss through calving and the surge process [12]. A surge in a Svalbard glacier typically lasts 7–10 years [13]. Maximal velocities in summer 2019 were still around 7 = 6 m/day, based on our analysis of Sentinel-1 Synthetic Aperture Radar (SAR) data.

Negribreen during surge is an ideal validation region, because many types of crevasses and other complex surface types exist in close proximity (Figure 1). This heavily crevassed environment provides a set of challenges for surface-height measurement and surfaceheight determination, whose solution will indicate that surface heights can be determined accurately in other crevassed and smooth glaciated regions as well. Because crevasses are relatively small features, a high geodetical accuracy and a close repeat of satellite track by airborne tracks is important, both will be analyzed in this paper. As there is not much literature about kinematic GPS data collection from a helicopter for the purpose of satellite data evaluation, the mathematical geodetical results of this paper are expected to be of use for future airborne experiments.

**Figure 1.** Negribreen during surge in 2017–2019. The multitude of different surface types indicates that a glacier during surge provides an ideal test situation for evaluation of height determination over crevassed and other complex ice surfaces. Photographs by Ute Herzfeld. (**a**) Negribreen surge July 2017, overview. (**b**) Wavy crevasses and folded moraines in central Negribreen, July 2017. (**c**) Several different crevasse provinces in central Negribreen, July 2017. (**d**) Calving front, ice fingers advancing into the ocean (front right), with smooth surface of non-surging Ordonannsbreen in the back left, July 2017. (**e**) Crevasse province near calving front July 2018. (**f**) Blocky crevasses near RGT 594 region, July 2018. (**g**) Ordonannsbreen-Negribreen boundary near calving front, August 2019. (**h**) Folded moraines and surge crevasses, August 2019.

The precision of ICESat-2 ATL06 (version 003) [14,15] ice surface height data has been determined as better than 7.2 cm, with a bias of less than 3.3 cm for smooth and flat ice/ snow surfaces at 88◦ Southern latitude in Antarctica, which is the circle of convergence of the ICESat-2 ground tracks [16]. Ref [17] conducted a field evaluation along the Chinese Antarctic Expedition Route (CHINARE, a route driven with a Pisten Bully) near Amery Ice Shelf, using GNSS, in December 2019 to February 2020, finding that ICESat-2 ATL06 (version 003 data) are precise to 1.5 cm with 9.1 cm precision and surface heights of ATL03 photon events [18,19] are accurate to 4.3 cm with 8.5 cm precision. Since this is a route that was driven with a Pisten Bully, it is all smooth, uncrevassed ice. However, the ICESat-2 ATL06 data product is a 40 m along-track ice-surface height product with 20 m postings, which generally does not report crevassed and otherwise morphologically complex glaciated terrain, and the photon classification on the ATL03 product often fails over such terrain [3]. The study in this paper complements such determinations of measurement accuracy over smooth and flat ice surfaces by investigations of surface-height determination over crevassed and otherwise complex regions, derived using the DDA-ice. The surface-height-determination problem has been addressed in [3], here, we focus on the geodetical component of ICESat-2 validation campaign data over crevassed terrain.

In this paper, we describe the ICESat-2 ATLAS data, experimental and instrumental setup of the airborne geophysical field campaigns, the process for geodetical correction of the airborne data using GPS/GNSS data and IMU data, and the resultant accuracies. Results will include (1) the capability to repeat satellite tracks during a helicopter-based airborne validation campaign, (2) the type of geodetical corrections and associated accuracy, (3) spatial characteristics of airborne data (on-ice spacing of resultant altimeter data), (4) determination of cross-over heights from airborne laser altimeter data after correction as an indicator of vertical precision of the field data, (5) the resultant accuracy of surface-height determination and variability, derived from ICESat-2 ATLAS data using the DDA-ice, and (6) the resultant ability to derive spatial characterizations of the morphological properties of the ice surface, exemplified by a surface-roughness parameter.

#### **2. ICESat-2 ATLAS Data**

#### *2.1. ATLAS Instrumentation, Data Collection Geometry, Basic Corrections and ATL03 Data Product*

ICESat-2 ATLAS data are collected along-track, for 3 pairs of 2 beams, a strong beam and a weak beam per pair, where the weak beam has a quarter of the energy of the strong beam [1,2]. Across-track separation of the beams on the Earth's surface is 3.3 km between the centers of adjacent pairs and 90 m for the beams within each pair; ICESat-2 beam pattern and track geometry are illustrated in Figure 2. The sensor operates in the 532 nm wavelength (green light) with a pulse-repetition rate (PRF) of 10 kHz. The PRF results in a nominal 0.7 m spacing of the laser pulses on the Earth's surface, under clear-sky atmospheric conditions. ATLAS has a footprint diameter of less than 17.4 m at 85% encircled energy. The Field of View (FOV) of the receiver telescope is 83.8 μrad, equivalent to 45 m FOV on the surface of the Earth. Observatory performance has been assessed in [20], where the actual footprint is characterized as closer to 10 m in diameter. As the FOV moves along the satellite ground track and returns from every single photon (in the 532 nm wavelength domain of the sensor) are recorded, surface structures at much higher resolution than footprint diameter can be resolved (see Section 2.2 and [3]). The instrumentation and derivation of the recorded photon point cloud are described in [2], and any technical component of the instrument and data set necessary for understanding the work in this paper is found in [3].

**Figure 2.** ICESat-2 beam pattern and track geometry. For observatory orientations of ATLAS flying forward ("weak beams leading") and ATLAS flying backward ("strong beams leading"). Beam numbering is a feature of the transmitter array, while spot and ground track numbering is a feature of the ATLAS data products, see Table 1. (Figure from [3], Figure 3).

ATLAS records returns from every photon in the 532 nm wavelength domain of the sensor. This results in a so-called photon point cloud. The photon point cloud includes photon events from ambient light (background) and photon events that results from the active laser signal of the instrument; these types of photons are initially not distinguishable. Classification of photons into signal and background requires a sophisticated mathematical algorithm, since the classification problem is mathematically ill-posed. For the solution of the classification problem and the determination of surface heights, we apply the density-dimension algorithm for ice surfaces, the DDA-ice [3,21], briefly introduced in Section 2.2. Data in this paper stem from the ICESat-2 ATL03 data product (version 003) [19], described fully in the Algorithm Theoretical Base Document (ATBD) for Geolocated Photon Points [18]. There is a classification of photons provided on ATL03, however, the work in our paper utilizes the classification based on the DDA-ice (the ATL03 classification often fails over crevassed surfaces). Note that technically, the data sets (ICESat-2 ATLAS data products) are referred to by version number (here: version 3), while the algorithms used to create the data products are referred to by release number (here: release 3). Because the geodetical corrections have not changed significantly in a way that would affect the ATL03 data geolocation and especially the representation of crevasse characteristics, the field validation presented in this paper is expected to be valid for subsequent releases as well (to release 5, the last release at publication date of this paper) [22].

Of relevance for the airborne validation campaign of the ICESat-2 ATLAS data is the concept of Reference Ground Tracks (RGTs). ATLAS ground tracks are referenced to the RGTs, which are calculated ahead of time, prior to the mission.  Actual ATLAS ground tracks are derived after data collection and following detailed geolocation.

The ICESat-2 beam pattern and track geometry depend on the orientation of the observatory, which can be flown in two modes, forwards (+x) and backwards (−x) (Figure 2). All data in this paper were collected during an "ATLAS oriented backwards" phase in July and August 2019. The orientation of the observatory changes which ground track stems from a weak beam and which one from a strong beam (see, Table 1). Therefore, for each track used in this paper "strong beam" or "weak beam" are specified (tracks are labeled gt1r, gt1l, gt2r, gt2l, gt3r, gt3l where the number stands for the beam and the letter for "left" or "right". Beam 2 is the center beam and the RGT splits the distance between gt2r and gt2l by design. Reference pair tracks split the distance between (gt1r, gt1l) and (gt3r, gt3l). Both weak beam and strong beam data were evaluated during our validation campaign. For shortness, we use the terminology "predicted ATLAS tracks" for ground tracks calculated ahead of time from the RGTs and "actual ATLAS tracks" for ground tracks derived from ATLAS observations, corrected as summarized in the next paragraph.

After data collection, the geolocation process results in the association of latitude, longitude and height for each telemetered photon event in the point cloud. It uses Precision Pointing Determination (PPD) of the laser and Precision Orbit Determination (POD). ATLAS carries a Laser Reference System (LRS), two GPS antennas and Spacecraft Star Trackers. The process is described in short in [3], in more detail in [2] and in full in the geodetical part of the ATL03 ATBD for Receive Photon Geolocation ATL03g [22], see also [23]. The location of the "bounce point" of the laser altimeter on the ice surface is given in the International Terrestrial Reference Frame (ITRF, Petit and Luzum [24]) ITRF 2014 (ITRF14) as latitude, longitude and height w.r.t. the World Geodetical System WGS84 in realization G1150 with ellipsoid constants (*ae* = 6,378,137 m, <sup>1</sup> *<sup>f</sup>* =<PHONE_NUMBER>63) [2]. We will utilize the same referencing for the airborne laser altimeter data of our campaign. After geolocation, several geophysical corrections are applied to the photon points [2,3].

**Figure 3.** Negribreen ICESat-2 Airborne Validation Campaign 2019: Experiment setup, track maps, location map and satellite image. (**a**) Field team 2019 (f.l.t.r. Ute Herzfeld, Harald Sandal (pilot) and Gustav Svanström (technician), Matthew Lawson, in front of the Airlift/ NPI helicopter; photograph by Thomas Trantow), on the side of Negribreen, flight 1, 12 August 2019, (**b**) GPS base station on a ridge above Negribreen, 2018 (Connor Myers), photograph by Ute Herzfeld. (**c**) Planned and actual flight tracks and predicted ICESat-2 ATLAS ground tracks. Flight 1, 12 August 2019; Flight 2: 13 August 2019. Background Landsat-8 data (5 August 2019) [LC08\_L1TP\_215003\_20190805\_20190820\_01\_T1]. (**d**) Zoom into (**c**) for area of RGT 594 and RGT 450 analyzed in this paper. (**e**) Location map. (Figure panels (**c**–**e**) after [3]).

*2.2. The Density-Dimension Algorithm for Ice-Surface-Height Determination (DDA-Ice) in Crevassed and Other Complex Terrain—Description in a Nutshell*

The analysis of the ATLAS data utilizes the Density-Dimension Algorithm for ice surfaces (DDA-ice) [3,21]. The DDA-ice is an algorithm specifically developed for surfaceheight determination from micro-pulse photon-counting laser altimeter data, especially those of the ICESat-2 ATLAS instrument and similar sensors. The DDA-ice is introduced in [21] for data from the Slope Imaging Multi-polarization Photon-counting Lidar (SIMPL) [25,26], an airborne micro-pulse photon-counting laser altimeter system that was employed as an ICESat-2 predecessor instrument in test flights, and has been first applied to ICESat-2 ATLAS data in [3]. The data from the ATLAS instrument have the form of a point cloud of geolocated photons [2]. The DDA-ice utilizes a Gaussian radial basis function (*rbf*) for data aggregation. The *rbf* is evaluated for every single photon in the point cloud; this operation results in a density field. In the next step, an auto-adaptive threshold function is applied to separate signal and background photons; this function automatically adapts to the different background properties of daytime, twilight and nighttime, apparent surface reflectivity (ASR) and some instrument characteristics, such as the afterpulse of the lidar signal. The algorithm is termed density-dimension algorithm, because density is used as an additional dimension in the signal-background classification step. The resultant surface heights have the same resolution as the sensor, which is nominally 0.7 m along-track under clear-sky atmospheric conditions. The DDA-ice uses its own cloud flag algorithm (it does not rely on the cloud flag in the atmospheric data product ATL04 [27,28] ). The last step of the DDA-ice is the ground follower. The ground follower is a piece-wise linear interpolator, whose segmentation automatically adapts to surface roughness, resulting in 2.5 m spacing for crevassed (or generally, rough) surfaces and 5 m for uncrevassed (or generally, smooth) surfaces. The DDA-ice is controlled by a set of so-called algorithm-specific parameters, introduced to facilitate quick adjustment of the algorithm should the performance characteristics of the ATLAS lasers change. The values of the algorithm-specific parameters employed in the analysis in this study are given in Table 2. Input data for the DDA-ice are geolocated photons, as given in the publicly available ICESat-2 ATLAS data product ATL03 "Global Geolocated Photons", version 3 [19], described in Neumann et al. [2,29]. The

DDA-ice does not utilize the signal classification given in ATL03, but signal-background classification is an intrinsic step of the DDA-ice (for a comparison, see [3]).

**Table 2. DDA-ice parameters.** Strong beam, weak beam—ATLAS beams. Parameters are the same as used in [3], parameters used for strong beam are the same as parameters optimized for analysis of SIMPL\_green (532nm parallel polarized) channels. Variables Q and S are now algorithm-specific parameters for DDA-ice. Resolution of the ground follower is 2.5 m over rough/ crevassed surfaces and 5 m over smooth surfaces. See also [21].

#### **3. Experiment Setup for the Negribreen ICESat-2 Airborne Geophysical Validation Campaign**

#### *3.1. Geographic Region and Ice-Surface Types during Surge*

The Negribreen Glacier System (NGS) is located on Spitsbergen in Svalbard, it forms an outlet of the Filchnerfønna ice cap (see, Figure 3). The calving front of Negribreen was located at approximately (78.573◦ N, 19.083◦ E) in 2019. An overview of main areas of the Negribreen Glacier System, including Negribreen (in center), Ordonnansbreen (the glacier that joins Negribreen close to the calving front) and Akademikerbreen (the glacier that joins upper Negribreen from the North) is seen in the background Landsat-8 image in Figure 3c.

Ice surface types during the surge include many different crevasse classes (Figure 1). Figure 1a shows smooth ice of Ordonannsbreen contrasting with the highly fractured ice of Negribreen, which calves into the bay, part of the Arctic Ocean. Wavy crevasses and folded moraines, the indicators that a glacier is a surge-type glacier, are seen in Figure 1b. Figure 1c illustrates several types of crevasse provinces, delineated from neighboring provinces by relatively clear boundaries. Heavily crevassed ice that advances rapidly into the Arctic Ocean indicates the rapid mass loss during the surge in July 2017 (Figure 1d) and July 2018 (Figure 1e). Figure 1f shows a typical surge crevasse type, characterized by blocky crevasses with subordinate, thin multidirectional crevasses in the region of RGT594. The surface of slow-moving Ordonannsbreen (left) is only crevassed near the boundary to surging Negribreen and otherwise smooth (Figure 1g). Lastly, Figure 1h presents a picture of multi-generational crevasses resultant from several deformation events in an area where the folded moraines exist, taken in August 2019. While maximal ice speeds during the height of acceleration in July 2017 reached 22 m per day, the region of crevassed ice expanded during the evolution of the surge and the maximal speeds decreased. Ice-surface velocity magnitudes derived from analysis of Sentinel-1 SAR data in August of 2019 were about 6 m (5.96 m) near the front of the glacier and 1–2 m in the region of the surveys

analyzed in this paper (Figure 4a). High-resolution satellite imagery from Planet SkySat, collected on 2019, August 18, illustrates the type of complex crevassing that exists near the ice front, in the fastest-moving region of the glacier, over 3 years after the start of the surge in spring 2016 (Figure 4b).

**Figure 4.** Negribreen velocity magnitude and crevassing at time of validation campaign in August 2019. (**a**) Velocity magnitude, calculated from ESA Copernicus Mission Sentinel-1 SAR data for 9 August 2019 to 21 August 2019. Highest velocity magnitude was 5.96 m/day near the calving front. (**b**) Crevasse types in the fastest-moving region near the front, Planet SkySat image, collected 18 August 2019. [20190818\_150858\_ssc9\_u0002\_panchromatic\_dn.tif]. (Figure panel (**b**) after [3]).

#### *3.2. Experiment Setup*

We conducted airborne geophysical observation campaigns over Negribreen during surge in the summers of 2017, 2018 and 2019. The campaigns over Negribreen had two objectives, (1) study of the state of the surge in the Negribreen Glacier System and (2) validation of ICESat-2 data. The 2018 campaign was carried out in preparation of ICESat-2 validation, prior to launch of ICESat-2 (which was on 15 September 2018), while the 2019 campaign obtained synoptical airborne and satellite observations. Airborne data were collected from a helicopter (Eurocopter AS 350), operated by Airlift under the auspices of the Norwegian Polar Institute, chartered for our campaign by the University of Colorado Geomathematics Group (Figure 3a).

During the 2019 field validation campaign, airborne geophysical data were collected along predicted ICESat-2 ATLAS ground tracks (see, Figure 3c) on 12 August 2019 (Flight 1) and 13 August 2019 (Flight 2). The locations of ground tracks from our flight campaign are referred to as flight tracks, or flight ground tracks (with coordinates calculated as described in this paper). The analysis in this paper utilizes predicted ATLAS ground tracks, actual ATLAS ground tracks and flight ground tracks. A zoom into Figure 3c, seen in Figure 3d, shows the tracks along RGT 450 and RGT 594, which will be analyzed in this paper, their relative location and crevasse types surveyed. Planet SkySat satellite image data, specifically collected to support this project, provide additional information suitable for evaluation of ICESat-2 measurements. Figure 3f shows crevasse provinces near the calving front in a Planet SkySat image taken on 18 August 2019. Extensional wavy crevasses, indicative of the area of highest velocities at time of image capture, dominate the northern central part of Negribreen, north of the trace of a medial moraine in upstream Negribreen. The fast-moving region is flanked by shear crevasses in its north and south, then (in an outward direction from the center) more closely spaced crevasses in the north, followed by smooth, slow-moving ice provinces in the north and south.

#### *3.3. Instrumentation*

During the 2018 and 2019 airborne geophysical campaigns, we collected laser altimeter data and image data, which were geolocated using coincidently collected Global Positioning System (GPS) and Inertial Measurement Unit (IMU) data. The same instruments were used in both years, except for the GPS Rover, and included the following:

	- (2.1) Kinematic GPS
		- (2.1.1) GPS Base Station placed on the side of Negribreen (Trimble NetR9)
		- (2.1.2) GPS Rover, mounted on a skid of helicopter (Septentrio APS-3G in 2018, Trimble R10 in 2019) (10 Hz data)
	- (2.2) "Button GPS", a simple GPS system placed on the dashboard of the helicopter under the window (1 Hz data)

*(1) Laser Altimetry.* The LaserTech ULS is a class-1, eye-safe laser that operates at 905 nm (Near-Infrared) with a configurable pulse-repetition frequency (PRF), set to 4000 Hz with 10-point averaging to yield an effective PRF of 400 Hz during the Negribreen campaigns. Beam divergence is 3 mrad, which results in a 30 cm footprint at 100 m range (see, Table 3). The ULS was integrated by members of the University of Colorado Geomathematics, Remote Sensing and Cryospheric Sciences Laboratory (the first three authors) for altimetry use from small manned and unmanned aircraft, utilizing a so-called engineering kit provided by LaserTech. An earlier integration is described in [30], for use aboard the NASA Sierra UAV during the Characterization of Arctic Sea Ice Experiment (CASIE) in 2009. Findings from geostatistical analysis of the CASIE data are reported and related to sea-ice modeling in [31].

*(2) GPS data.* Two separate GPS units were used for redundancy. The first, called "Button GPS", is a low cost, plug-and-play type used for making quick, low-accuracy plots and serves as a backup for the second system. Button GPS was placed on the dashboard of the helicopter, where it had good sky view through the large front window. The second system is a high-cost high-accuracy two-part GPS system, which consists of a base station reference receiver and a helicopter mounted kinematic GPS unit (rover). The base station collects observations including atmospheric delay parameters and orbit information that allow for post-processing that increases the position accuracy of the kinematic GPS unit mounted on the helicopter. Sky view of the base station was unobstructed, because the base station was positioned on an elevated moraine edge or near a cliff edge on the southern side of Negribreen, several hundred meters above the glacier (Figure 3a,b). There were no large mountains obstructing the sky view of the base station. The rover was mounted on one of the skids of the helicopter, and thus the body of the helicopter and the helicopter rotors possibly interfered with sky view of the rover. The instrumentation for the differential GPS data collection was provided by UNAVCO in both years.

*(3) IMU data.* IMU data record the attitude of the helicopter during the flights. These data are collected to allow correction for attitude and vibration. The IMU is also owned by our Geomathematics Laboratory and integrated into the survey system. Our postprocessing includes coincidently collected IMU data to remove helicopter vibration with a low-pass filter, and to correct for roll, pitch, and yaw with a transformation matrix, taking into account the spatial offset between the mounting positions of the laser altimeter system and the GPS rover on the helicopter.

**Table 3. Technical specifications of the Universal Laser System (ULS) and settings during the Negribreen airborne geophysical observation and validation campaigns in 2017, 2018 and 2019.** FDA—United States of America, Federal Food and Drug Administration, CFR—Code of Federal Regulations Title 21, see https://www.accessdata.fda.gov/scripts/cdrh/cfdocs/cfcfr/cfrsearch.cfm ?FR=1040.10 (accessed on 3 October 2021).

*(4) GoPro time-lapse-imagery* was recorded at two images per second in linear imaging mode. The linear mode was selected, because it yields the least distortion of the imagery. Time-lapse imagery is better suited for comparison with discrete altimeter data than video data, results in much smaller data volumes and reduces battery energy usage.

*(5) Photographs* were taken using handheld cameras to provide overviews of the glacier system during surge, especially from high altitudes above ground level (a.g.l.), to collect documentation of crevasse provinces, calving events and other features of special interest to the study of the surge process. Several thousand images were taken during each flight.

#### **4. Comparison of Different Approaches to GPS-Data Collection and Processing: Real-Time Kinematic (RTK) Compared to Kinematic GPS Data Collection with Differential Correction Using Base and Rover Data**

#### *4.1. Scientific Motivation and Experiments*

In preparation of the ICESat-2 launch (15 September 2018), in summer 2018 we carried out experiments aimed at optimizing GPS data collection and resultant repeat accuracy of satellite ground tracks during airborne validation flights of ICESat-2 measurement capabilities. These include experiments with Real-Time Kinematic (RTK) data collection, which are motivated as described in the following paragraphs.

There is relatively little knowledge of the geodetical accuracy that can be obtained from kinematic data collection from a fast-moving platform, especially from a helicopter that is susceptible to wind and turbulence. In an approach that utilizes a base and rover with post-processed solutions, the challenge is to acquire sufficient satellite views to maintain a fixed solution in the GPS post-processing for at least 10 min, while the helicopter is flying. The faster the vehicle/rover moves, and the less linearly due to turbulence, the harder it is to find solutions. The mounting location of the GPS receiver and, for a helicopter, interference with rotor blade occultation, create additional unknowns [32].

The problem is compounded by the limited mounting options for a GPS rover on a helicopter. In our Svalbard experiments, the rover had to be mounted on one of the skids for safety reasons, where it had limited sky view due to the body of the helicopter and the rotors. Therefore, we included experiments of RTK data collection in our 2018 experiments, as an alternative to the post-processed base-rover GPS data analysis.

Combining the 2018 pre-launch experiments with 2019 experiments along near-time ICESat-2 tracks, we have results from three types of GPS experiments:

Type 1 GPS experiments involve setting up a GPS base station at the side of the glacier on an elevated moraine or near a cliff edge for the duration of the airborne survey, typically two hours over the Negribreen region, and collecting kinematic data with a rover GPS mounted to the helicopter. During these experiments, the pilot follows a predetermined flight path. Experiments of best flight path segmentation were carried out. After the field experiment, base station data, recorded at 10Hz, are used to differentially correct rover data, also collected at 10 Hz. Because of regulations of the governor of Svalbard, the GPS base cannot be left in the glacier region. While a 24 h recording time for the base would be preferred, the length of time of about 2 h has been found to be sufficient for the accuracy of the project objectives.

Type 2 experiments require a base station and a rover as well. But in contrast to type-1 experiments, during RTK data collection, the correction of the rover position is carried out during the flight, using information sent over by the base station, and displayed on a blue-tooth-enabled device. This enables a correction of the flight path during flight.

In 2019, the GPS base station malfunctioned or was set incorrectly, recording only 15-second data, thus the differential correction applied to the 2019 GPS rover data uses Spatial Reference System Precise Point Positioning (CSRS-PPP), a service provided by Canadian Geodetic Survey, Natural Resources Canada, for post processing [33,34]. This situation yields an unintentional third type of experiment, whose results may be of interest for flights where a base station cannot be placed.

#### *4.2. Results*

(1) Kinematic GPS data collection with differential correction of rover data with base station data in post-processing yields the highest accuracy of positioning, of our three types of experiments. Experiments on flight-path segmentation indicate that quarter-crossings of the glacier work best, aiming to fly a straight line across a 2–3 km distance. Full glacier crossings are similarly good with respect to resultant accuracy of track repeats and require fewer waypoints to be entered manually by the pilot into his on-board system. Notably, the entire survey setup requires good collaboration and communication between the pilot, the experiment lead (the PI) and the survey team who continuously monitors instrument functionality and height (altimeter range). As a result, full glacier crossings of about 10–20 km flight distance depending on location were used in 2019.

(2) Experiments with RTK indicated that the added accuracy from the real-timekinematic correction is more than offset by the inaccuracies introduced by the continued reactions of the pilot to adjust the flightpath and thus the aircraft attitude to the updated position information (from the RTK system). Flying straight towards points, predetermined by the expedition lead/ scientist with 2–4 points for a glacier crossing of 10–20 km track length and shared with the pilot prior to the flight, yields data from a more "quiet" helicopter flight and more consistent, linearly interpolatable data. Thus, the repeat accuracy of a pre-determined satellite track is not limited by the accuracy of the GPS solution, but instead by the maneuvering characteristics of a helicopter. The pilot cannot fly a straight path when the quick changes in the RTK-GPS directions call for unpractically frequent changes of the flight direction of the helicopter. A helicopter turn necessitates leaning the helicopter to the side, which then introduces additional distortions to the flight path and results in off-nadir laser pointing.

(3) Results from CSRS-PPP [33,34] are sufficiently accurate for validation of ICESat-2 surface height determination and morphological characteristics, as will be detailed in the following sections, while most likely significantly less accurate than positions differentially corrected with a base station. The CSRS post-processing uses precise satellite orbits, clock and bias corrections from a global network of receivers. Because of the much shorter base line resultant from a placing a base station in the field area, differential corrections using a locally placed base station are known to yield highest accuracy [35].

(4) Despite the problem of interference of the rotor blades and the helicopter body with view of GPS satellites, the amount of GPS rover data collected was found to be sufficient for accurate geolocation of the ULS data, despite the mounting location of the GPS rover on the skids.

(5) Outlook: Based on the 2018 and 2019 experiments (and previous experiments in similar types of field surveys), highest accuracy will be achieved by using differential GPS with a base station on the side of the glacier (as shorter baselines yield higher position accuracy) and a GPS rover on the aircraft. It may be easiest, especially for students new to GPS data collection, to utilize a rover and a base GPS instrument from the same brand, while not necessary to process the data (which can be converted to RINEX format independently of the receiver brand). Trimble has a proprietary format, but RINEX conversion remains a possibility.

The summer 2018 data are not used in the remainder of this paper, because they were collected prior to the launch of ICESat-2 on 15 September 2018.

#### **5. Surface-Height Determination from Synchronization of Laser and GPS and IMU Corrections**

*5.1. Processing of 2019 GPS Data*

Processing of 2019 GPS data was carried out using the following steps:

In essence, the CSRS software was able to provide a fixed solution for the rover data. To best match with the ICESat-2 data and their coordinate reference, we selected ITRF14 as the reference frame for the CSRS-PPP solutions. The coordinates used in this analysis are given as *(rover\_data\_set\_PPP\_v3)*, in csv format (human readable, including lat, lon, elevation, time, reference frame) with reference to ITRF14. The Canadian Geodetic Survey uses the GRS80 ellipsoid (*ae* = 6,378,137 m, <sup>1</sup> *<sup>f</sup>* =<PHONE_NUMBER>01) for the realization of ITRF14 in CSRS-PPP (Brian Donahue, Canadian Geodetic Survey, friendly pers. comm. 19 January

2022). For practical purposes, the ITRF14 corresponds to WGS84 (conversion from ITRF14 to WGS84 results in identical coordinates, for a more precise description, see Section 2). In summary, practically the same coordinate reference is utilized for the airborne data and the ICESat-2 ATLAS data.

#### *5.2. Comparison of Button GPS and Corrected Rover GPS Data*

In addition to the rover GPS, GPS data were also collected using a simple instrument, the Button GPS, located on board the helicopter. The Button GPS is a plug-and-play type inexpensive instrument, that is placed under the windshield of the helicopter and has proven to be very reliable during several years of airborne campaigns. The Button GPS records location at 1 Hz, while the rover GPS records location at 10 Hz. This Button GPS has lower accuracy, especially of height data, and lower frequency than the rover GPS, but is directly connected to the ULS laser system and records on the same computer. The rover GPS recorded only time-from-start-of-rover-time, it records on a GPS-internal computer and hence its times are independent of the ULS. Thus, we are using the time stamps from the Button GPS to align the corrected rover-GPS data with the ULS measurements.

This process is illustrated in Figure 5, which shows heights and times from the two GPS receivers during flight 2. Button GPS recordings appear to be "ahead of" the rover GPS with respect to distance, an effect that results from starting the Button GPS first. Since the Button GPS and ULS laser were being controlled from the same computer, we have exact time synchronization between them. In contrast, the rover GPS was controlled through its internal computer (and a handheld controller), which results in timing independent from that of the ULS laser. The first correction step involved time-shifting the rover GPS to align with the Button GPS/ULS laser (Figure 5).

#### *5.3. Derivation of Surface Height Using ULS and Differential GPS or Button GPS*

Ice-surface height is calculated by subtracting the laser-retrieved height from the GPSretrieved height. Two apparent surface height values result from pairing ULS data with rover GPS data and Button GPS data (Figure 3b). In addition to the time lags between the rover GPS and the ULS laser, the change in height of the aircraft and the surface roughness of the glacier contribute to the apparent differences in height. For example, when the GPS starts to move down (aircraft lowers), the laser-retrieved height would normally be lower as well. However, in our case, the laser is ahead of the GPS, so it has already started recording lower heights. The different start times contribute a small part of the large apparent height difference values between the rover GPS-ULS and the Button GPS-ULS data, the main source of the apparent height difference is the poor height determination of the Button GPS. ULS data with ranges greater than 500 m are excluded from the analysis, as the ULS does not provide enough returns at ranges larger than 500 m (see Section 7).

#### *5.4. Synchronization of Laser and Differential GPS and Button GPS*

Crevasse locations are used to further improve the determination of a time shift, performed visually to match. Time accuracy is 0.1 s. Results in Figure 5c show that crevasse morphology lines up in the recordings of the Button GPS and the differential GPS, thus using crevasse morphology to aid in determination of a time shift is a feasible approach. There are still height differences on the order of tens of meters, which are due to the inherent properties of the low-resolution Button GPS. Because absolute heights from Button GPS data are less accurate than heights from corrected rover GPS data, heights from corrected rover GPS data are used in the following analysis.

#### *5.5. IMU Correction*

An Inertial Measurement Unit (IMU) records the attitude of the survey platform, here the helicopter, during flight, at (in our case) a frequency of 1 Hz. As the laser operates from the helicopter, it actually measures the range from the helicopter to the Earth's surface, rather than the height above the surface, and along a vector to the ground that is not normal to the ground surface. The IMU correction is a change of reference system, from that centered at the moving helicopter to an Earth-centric system. Prior to the actual correction, a low pass filter (moving average filter) is applied to remove helicopter vibration from the IMU data. The attitude reported by the IMU sensor is utilized in the form of Euler angles for correcting lidar laser range returns, and an orientation matrix for correcting the location of the range measurement in (x,y) coordinates. Formally, the IMU correction utilizes a quaternion multiplication (introduced e.g., in [36]). The IMU correction is applied to the GPS-corrected ULS data (with GPS corrections as described in Sections 5.1–5.4). It should be noted that the IMU corrections do not only yield corrected locations on the ground, but also corrections of the range values (this is a matrix multiplication applied to all components of the ULS data).

Figure 6a,b shows the magnitudes of the components (roll, pitch and yaw) of the IMU corrections. Figure 6c,d illustrates the magnitude of range corrections. At the scale of absolute range, this is not resolved (in Figure 6c,d). The correction magnitude (Figure 6c,d) is less than 5 m on average, with some spikes up to 10 m. These usually occur when the aircraft is turning or during heavy winds. We experienced heavy winds near the end of flight 2, which explains the largest difference.

**Figure 5.** Illustration of the process of surface-height determination from synchronization of ULS and GPS data. Negribreen validation campaign. Flight 2, 13 August 2019. Corrected GPS data are used

(see Equation (1)). (**a**) Comparison of GPS heights and times from two GPS receivers. Button GPS (red), recorded at 1 Hz, started first, together with ULS and recorded on the same computer as the ULS (exact time synchronization of ULS and Button GPS). Rover GPS (blue), recorded at 10 Hz, started second, time recorded on the GPS-internal computer. The first correction step is time-shifting the rover GPS to align with the Button GPS/ULS laser. (**b**) ULS apparent surface heights, derived when pairing ULS data with Rover GPS data (blue) and Button GPS data (red). (**c**) ULS surface heights, after time-shift correction between ULS laser and the rover GPS. Now both elevation measurements and crevasse morphologies resemble each other.

**Figure 6.** IMU values and IMU-corrected ULS ranges. Negribreen validation campaign. Flight 2, 13 August 2019. (**a**,**b**) IMU values (roll, pitch and yaw) for flight2, part1 (**a**) and flight2, part2 (**b**). (**c**,**d**) Range correction (total corrected and uncorrected range and magnitude of correction) for flight2, part1 (**c**) and flight2, part2 (**d**).

#### **6. Analysis of Airborne Geophysical Data and Results: Ice-Surface Heights, Technical Information and Crossover Analysis**

In this section, we analyze the airborne geophysical data (without comparison to ICESat-2 data) to derive ULS range and ULS-based ice surface heights along flight tracks (Section 6.1). Section 6.2 is devoted to technical information that may be useful to design future airborne validation campaigns from small aircraft, including helicopter velocity, optimal ULS range and aircraft altitude above ground, resultant data frequency and on-ice spacing. These results generalize to other types of manned aircraft and to some extent, unmanned aircraft. Precision of the ULS height data in crevassed terrain is assessed using a crossover analysis (Section 6.3).

#### *6.1. Results Part 1: ULS Range and Ice-Surface Heights Based on GPS and IMU-Corrected Data*

Following the GPS and IMU corrections as described in Section 5, we now have the correctly geolocated coordinates for the ULS data. After application of the IMU and GPS geodetical corrections to the ULS data, we obtain ULS data with geographical coordinates

$$x\_{\\$^{\text{lcp}}\text{corr}} = (lat\_{\text{corr}}, lon\_{\text{corr}}, z\_{\text{corr}}), \tag{1}$$

where *zcorr* is the corrected surface height. As described above, the GPS coordinates are given in IRTF14 (as output from the CSRS-PPP software). The bounce point of the ATLAS signal is also given in the IRTF14 reference frame as latitude, longitude and height w.r.t. the WGS-84 (G1150) ellipsoid for each point in the photon cloud reported on the ATL03 data product [2,3,18,22] (Section 2.1). Therefore, corrected ULS data from the validation campaign are directly comparable to ATL03 data.

Lastly, ULS data *xgeogcorr* given with reference to the geographic coordinate system are converted to ULS data with reference to the Universal Transverse Mercator (UTM) coordinate system

$$x\_{utmorr} = (north\_{corr}, east\_{corr}, up\_{corr})\tag{2}$$

with reference to the WGS84 ellipsoid and UTM zone 33 (central meridian 15◦ E). We use *upcorr* = *zcorr*. We are now able to determine the range of the ULS to the ice surface at the laser measurement rate (Figure 7a), which combined with the aircraft position gives the surface height of the glacier (Figure 7b). The corrected surface heights can be used in several glaciological applications, including monitoring of surface elevation change during the surge, indicative of mass transfer during the surge, and derivation of the spatial characteristics of the crevassed ice surface (see, Section 7).

#### *6.2. Results Part 2: Technical Information: Helicopter Velocity, ULS Range, Data Frequency and On-Ice-Spacing*

The range of the ULS is an important experiment parameter to control, as the ULS functions for altitudes up to about 500 m and works best for shorter ranges up to about 300 m. Helicopter safety requires that the instrument be flown at sufficient distance above the surface, dependent on (a) visibility and surface reflectance, and (b) topographic relief, morphology and crevassity of the terrain. Following [11], we use the term "crevassity" as a measure of the void fraction in glacier ice as a material, akin to porosity. Crevassity is utilized here as a general term to describe the magnitude and complexity of the overall deformation at a given location, which for a surge glacier manifests itself in the crevasse characteristics.

All ICESat-2 validation flight lines were flown at this optimal height range, with tracks in the upper central Negribreen (inlcuding those for RGT 450 and RGT 594) collected with a range of 100-200 m and tracks traversing towards the upper inflows flown at 200–350 m range (Figure 7a). The ULS is originally designed for measurements at much shorter ranges and there is no practical lower limit of the measurement range (we test the instrument in the laboratory, for distances down to 1 m). The practical lower limit of the range is constrained by flight operations safety, as the helicopter pilot needs to maintain a safe height above the highly reflective, sloping and crevassed ice surface, at a flight speed that is efficient for data collection.

Flight velocity magnitudes, evaluated after GPS and IMU corrections, are given in Figure 7c and after along-track smoothing in Figure 7d. Speeds for the validation tracks (RGT 450 and RGT 594) are generally between 25 m/s (90 km/h) and 35 m/s (126 km/h), slower over steep ice falls. Averaged helicopter velocities were 17 m/s (averaged for the whole survey flight) and 30–35 m/s while flying over the ICESat-2 segments. Mean velocities were 34.426 m/s over the RGT 450 segment and 32.0451 m/s over the RGT 594 segment. The ULS functions well at higher velocities as well, as supported by the fact that frequency of recorded data is generally close to the 400 Hz nominal frequency of ULS

setting used during our experiments (Figure 7e, for ULS settings, see Table 3). Resultant on-ice along-track spacing was 0.06–0.08 m for the validation tracks (Figure 7f).

**Figure 7.** Results from ULS data collection. Negribreen validation campaign. Flight 2, 13 August 2019. The scale does not allow for identification of crevassed regions. (**a**) ULS range aircraft to ice surface. (**b**) ULS surface height (elevation). (**c**) Aircraft velocity magnitude. (**d**) Aircraft velocity magnitude, smoothed. (**e**) Data points per second. (**f**) On-ice spacing of ULS data. Background: Landsat-8 image from 5 August 2019 [LC08\_L1TP\_215003\_20190805\_20190820\_01\_T1].

Flight altitudes above ground are further constrained by the other instruments: Higher altitudes allow collection of better imagery, i.e., imagery that shows larger crevassed regions. Flight altitudes of 10 km or more are flown for collection of overview imagery. Altitude in transit from one observation target to another can be higher than 300 m and still yield altimeter data useful for science. For example, the blue stretch in Figure 7e is flown at 400 m above ground, resulting in lower data frequency and somewhat larger on-ice spacing of 0.08–0.14 m. ULS data with ranges greater than 500 m are excluded from the analysis, as the ULS does not work properly at this range, and the situation of large ranges only occurs when we flew high to collect overview imagery and the ULS was not turned

off. In conclusion, the ULS can be operated with full recording frequency for altitudes up to 300 m above ground level (a.g.l.). The advantage of flying at lower altitudes a.g.l. is that the on-ground spacing is smaller, as the spacing is a direct function of altitude above ground. Between 300 m and 500 m, the recording frequency can drop somewhat, but still yield useful science data. Lower altitudes of 100 m can be used for flights with unmanned aerial systems over sea ice, where risk is lower and the surface is flat, as during the Characterization of Arctic Sea Ice Experiment (CASIE) [30,31].

#### *6.3. Results Part 3: Crossover Analysis of ULS Heights*

The crossover analysis is performed on the corrected height estimates for the 2019 flight data in the upper glacier (see Table 4). Crossover height differences are calculated after averaging all points within a 10 m radius of the intersection location for each crossing flight segment, resultant in 14 crossover locations with a mean elevation difference of 0.111 ± 3.056 m. Similar values are achieved when the averaging radius is 1 m and 0.1 m. Largest differences (red and dark blue points in Figure 8) occur at locations where the aircraft is turning significantly, where laser ranges are large and IMU corrections are most relied upon, and at locations of large crevasses where height variations on the order of tens of meters occur. However, most of the crossover differences are near zero (green points) when the aircraft is flying straight over non-crevassed areas or areas with beginning crevassing. If we remove the 3 crossover locations that occur on significant turns, we obtain a mean elevation difference of −0.172 ± 2.564 m for the remaining 11 points. Since underflights of ICESat-2 tracks follow straight flight paths, the estimated the precision of our measurements is approximately 2.6 m in crevassed terrain. Accuracy estimates would require ground elevation truths independent of and coincident with our flight data, however collection of ground measurements on a glacier during surge is not feasible (dangerous). Crevasse depths ranging to 45 m were observed with ICESat-2 and determined with the DDA-ice [3].

**Figure 8.** Crossover elevations from ULS data collection. Negribreen validation campaign. Flight 2, 13 August 2019. There are 14 crossover locations. Removing the three crossover points on large turns, for the remaining 11 points the mean difference is −0.172 ± 2.564 m. Background: Landsat-8 image from 5 August 2019 [LC08\_L1TP\_215003\_20190805\_20190820\_01\_T1].

**Table 4. Crossover analysis for Flight 2, 13 August 2019, Negribreen ICESat-2 Validation Campaign 2019.** Estimates are given when using all crossover locations (seen in Figure 8) and for locations only along relatively straight flight segments, i.e., where the helicopter was not performing a significant turning maneuver. Crossover estimates are averaged within a radius specified in the last column of the table.

#### **7. Analysis of Geolocated ULS and ICESat-2 ATLAS Data**

#### *7.1. Overview*

Following the geodetical and technical results in the last sections, we now have established the technical basis needed to (1) quantitatively investigate the capability of a light-weight experiment using a small helicopter to repeat satellite tracks accurately and (2) compare surface heights and crevasse morphologies from ULS and ATLAS measurements. Results relating to the first topic will be summarized in Section 8 and results regarding the second topic will be reported in Section 9. Here, the analysis is carried out as a case study for each of five sets of tracks, to facilitate understanding the relationships between flightrepeat capability of our experiment setup dependent on a number of factors including wind, resultant track separation, crevasse representation in ULS data and in ATLAS data for weak and strong beams, crevasse morphology, crevasse characteristics and surface roughness.

ICESat-2 tracks for RGT 450 and RGT 594 were repeated, including both strong and weak beams, resultant tracks for ULS and ATLAS are shown in Figure 9a,b and listed in Table 5. A useful tool in this context is plotting the results in 1-km along-track segments, as shown in Figure 9. Figure 9a,b provides the 1-km segments for flight 2, parts 1 and 2, respectively, carried out on 13 August 2019 over Negribreen.

Data from flight 1 are not used in the analysis in this paper, because the GPS rover did not record data. Data from the Button GPS are still useful for analysis of surge evolution, reported elsewhere.

Flight 2 has two parts, each part includes an overpass of RGT 594 (gt1 l and gt1r) and RGT 450 (gt1l and gt1r) tracks. The track of the RGT 594 weak beam (gt1r) was overflown twice (both passes in part 2 of flight 2) in order to obtain information on the quantitative measures listed above from a repeat experiment. However, because the weather worsened significantly during flight 2, the first repeat occurred during approximately normal weather conditions and the second one during unusually strong winds of over 30 knots (60 km/h). The increasing winds led the pilot to end the data collection for this day soon after. The three tracks for RGT 594 were flown in the order (1) "RGT 594, gt1r, part 1", (2) "RGT 594, gt1l", (3) "RGT 594, gt1r, part 2", with other flight lines in between, hence wind speed increases in this order. Labels in the figures are as follows: "RGT 594, gt1r, part 1" refers to overpass 1 (less wind). "RGT594, gt1r, part 2" refers to overpass 2 (stronger wind). During strong winds, the pilot has to hold the helicopter at an angle relative to the direction of flight and tilted, which makes following the prescribed tracks a lot harder and increases the IMU values. The IMU still functions for this situation (see Figure 6). Resultant track separations for RGT 594, gt1r, part 2 are higher than for the other tracks, but crevasse comparisons are still useable for evaluation of ICESat-2 tracks, as the analysis will show.

The tracks of the RGT 594 gt1r and gt1l are the location of analysis of 2018 prelaunch Negribreen airborne campaign data and first post-launch ICESat-2 data from early winter 2018 reported elsewhere, which indicated that ICESat-2 measures summer's crevasses through winter's snow cover. Because the identification of crevasses under snow was based on non-synoptic airborne and satellite data with several months separation, an additional objective for the 2019 validation campaign was the collection of synoptic satellite and airborne data for the same crevasse field surveyed in summer 2018 with the ULS and in winter 2018 with ICESat-2 ATLAS.

In the following subsections, we analyze each of the tracks (1) RGT 450, gt1l (strong beam) (Figure 10), (2) RGT 450, gt1r (weak beam) (Figure 11), (3) RGT 594, gt1l (strong beam) (Figure 12), (4) RGT 594, gt1r (weak beam), overpass 1 (Figure 13), and (5) RGT 594, gt1r (weak beam), overpass 2 (Figure 14).

For each example, we show the location of the ground track sections, indicated by red and green tracks, superimposed on a Landsat image from 5 August 2019 (collected only 8 days before flight 2), in panel (a) of Figures 10–14. Zooming into the area, with the flight tracks plotted over the aforementioned SkySat image from 18 August 2019, gives and impression of the crevasse provinces under the flight tracks, and visually illustrates the across-track distances between the predicted and actual ATLAS tracks and the actual ULS tracks (panel (b) in Figures 10–14). The distances are then given in a graph that illustrates the variability of the distance between the ULS track and the actual ATLAS track, with the distance between the predicted and actual ATLAS tacks added in for comparison (panel (c) of Figures 10–14). Note that at time of flight, only the predicted ATLAS tracks were available, thus we aimed to underfly the predicted tracks for each beam. Track variability is exaggerated by the over-heightened plot style. Variability of the ATLAS-predicted to ATLAS-actual distance is shown in better resolution in panel (d) of Figures 10–14. In summary, capability for track repeat ranges from an average of 11.32 m to 24.87 m, with an outlier of 60m during strong storm conditions. Differences between ICESat-2 predicted and actual range from 6.51 m to 15.65 m for the average difference of a glacier crossing.

Figure panels (e) and (f) in Figures 10–14 aid in the comparative analysis of the surface heights of crevassed regions. Segment lengths in (a), (c), (d) and (e) are the same within each figure and of approximately 2–3 km length, while segment lengths in (b) and (f) show subsets for better representation of the crevasse morphologies. To further aid the comparison of panels (e) and (f), we have added a blue bar in panels (e) that shows the location and coverage of the DDA-ice results seen in figure panels (f) of the same figure. Surface heights of ICESat-2 ATLAS data are derived by application of the DDA-ice, which was briefly described in Section 2.2.

**Table 5. Average distances of ATLAS predicted vs. ATLAS actual and ATLAS actual vs. ULS actual. Values calculated for entire glacier crossings**. Note that Figures 9–13 highlight 2–3 km sections for each crossing.

**Figure 9.** Flight path segments, resultant from ULS data collection, and co-location with ATLAS ground tracks over crevasse fields. Negribreen validation campaign. Flight 2, 13 August 2019. (**a**,**b**) Flight paths, colored in 1-km segments to match the along-track distance vs. elevation plots in Figures 9–13. Background: Landsat-8 image from 5 August 2019 [LC08\_L1TP\_215003\_20190805\_20190820\_01\_T1]. (**a**) Flight 2, part 1. (**b**) Flight 2, part 2. (**c**) Example of co-located airborne ground track (ULS), left, and ATLAS ground track (RG 450 1l), right. Even with an 18 meter average separation, the airborne ground track (ULS) and the ATLAS ground track cross the same crevasse provinces. Background: Planet SkySat image from 18 August 2019 [20190818\_150858\_ssc9\_u0002\_panchromatic\_dn.tif].

**Figure 10.** ICESat-2 ATLAS vs. ULS data. Track repeat accuracy and results. RGT 450, gt1l (strong). Negribreen validation campaign. Flight 2, 13 August 2019. (**a**) Location of ATLAS segments over Negribreen and closest ULS data (note the tracks are so close that the green line disappears under the red line). For segment numbering, see Figure 9 (ATLAS DDA segments). Background: Landsat-8 image from 5 August 2019 [LC08\_L1TP\_215003\_20190805\_20190820\_01\_T1]. Segment lengths in (**a**,**c**–**e**) are the same, segments in (**b**,**f**) are subsets. (**b**) Track separation between our ULS flight path, the ATLAS predicted flight path and the ATLAS actual flight path for granule: ATL03\_20190727132129\_04500405\_003\_01.h5. Background: Planet SkySat image from 18 August 2019 [20190818\_150858\_ssc9\_u0002\_panchromatic\_dn.tif]. (**c**) Separations between ATLAS Track 450 (gt1l) actual and ULS flight path, and between ATLAS Track 450 (gt1l) actual and ATLAS Track 450 (gt1l) predicted. Calculated for actual ATLAS and ULS flight paths from data shown in (**b**). (**d**) Track separation between predicted and actual ATLAS 450 (gt1l) flight paths. (**e**) ULS surface heights and ICESat-2 ATLAS DDA-ice surface heights superimposed, after GPS and IMU correction of ULS paths (ground locations). (**f**) ICESat-2 ATLAS surface heights over crevassed region derived using DDA-ice, RGT 450, gt1l, segment 25. DDA-ice parameters as in Table 2. ATL03 data set given in caption for (**c**). Rel. 03 data.

**Figure 11.** ICESat-2 ATLAS vs. ULS data. Track repeat accuracy and results. RGT 450, gt1r (weak). Negribreen validation campaign. Flight 2, 13 August 2019. (**a**) Location of DDA segments over Negribreen and closest ULS data (note the tracks are so close that the green line disappears under the red line). For segment numbering, see Figure 9 (ATLAS DDA segments). Background: Landsat-8 image from 5 August 2019 [LC08\_L1TP\_215003\_20190805\_20190820\_01\_T1]. Segment lengths in (**a**,**c**–**e**) are the same, segments in (**b**,**f**) are subsets. (**b**) Track separation between our ULS flight path, the ATLAS predicted flight path and the ATLAS actual flight path for granule: ATL03\_20190727132129\_04500405\_003\_01.h5. Background: Planet SkySat image from 18 August 2019 [20190818\_150858\_ssc9\_u0002\_panchromatic\_dn.tif]. (**c**) Separations between ATLAS Track 450 (gt1r) actual and ULS flight path, and between ATLAS Track 450 (gt1r) actual and ATLAS Track 450 (gt1r) predicted. Calculated for actual ATLAS and ULS flight paths from data shown in (**b**). (**d**) Track separation between predicted and actual ATLAS 450 (gt1r) flight paths. (**e**) ULS surface heights and ICESat-2 ATLAS DDA-ice surface heights superimposed, after GPS and IMU correction of ULS paths (ground locations). Blue line: Section covered in f. (**f**) ICESat-2 ATLAS surface heights over crevassed region derived using DDA-ice, RGT 450, gt1l, segment 25. DDA-ice parameters as in Table 2. ATL03 data set given in caption for (**c**). Rel. 03 data.