The potentials of high-resolution photogrammetry for analyzing glacier retreat in the Ötztal Alps, Austria

Glaciers all over the world experience an increasing mass loss during recent decades due to change in the global climate. This leads to considerable environmental consequences in the densely populated Alps and many other mountain ranges in the world. We used high-resolution aerial photogrammetry within the AlpSenseBench project to investigate glacier retreat in great spatial and temporal detail in the Ötztaler Alps, a significant glacier area in Austria. Long-term in situ glaciological observations are available for this region, and a multitemporal time 15 series of digital aerial images with a spatial resolution of 20 cm acquired over a period of 10 years exists. Glacier retreat of all 25 glaciers in the region, including the Vernagtferner, was analyzed by investigating glacier extent and surface elevation changes, derived from the aerial images by digital surface model (DSM) generation. Due to different acquisition dates of the large scale photogrammetric surveys and the glaciological data, a correction was established using a dedicated unmanned aerial vehicle (UAV) survey across the major part of the Vernagtferner. 20 This allowed us to compare the mass balances from geodetic and glaciological techniques, which reveals the potentials of the combination of these two techniques for gaining a better insight into glacier changes and its spatial distribution. The results show a clear increase of glacier mass loss for all glaciers in the region, including the Vernagtferner over the last decade. Additionally, the influence of debris-cover on mass balance, as well as the magnitude of dynamic processes, could be quantified. The comparison of geodetic elevation differences and the 25 interpolated glaciological data reveals that there exists a high potential in detecting local peculiarities of mass balance distribution and for correcting small scale deviations, not revealed in the interpolated glaciological information. The availability of high resolution multi-temporal digital aerial imagery for most of the glaciers in the Alps will provide a more comprehensive and detailed analysis of climate change-induced glacier retreat. 30 https://doi.org/10.5194/tc-2020-263 Preprint. Discussion started: 27 October 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
The impacts of climate change are widespread and clearly visible in the Alps (Rogora et al., 2018) but particularly evident in the dwindling glacier resources (Sommer et al., 2020;Zekollari et al., 2019). Over the past 100 years, the temperature in the Alps has increased almost twice as fast compared to the global average, resulting in nearly 2 °C higher mean air temperatures (Auer et al., 2007;Marty and Meister, 2012). By the end of this century, mean 35 air temperatures are expected to further rise by several degrees Celsius (Gobiet et al., 2014;Hanzer et al., 2018).
Due to this ongoing climate evolution, alpine glaciers may lose half of their volume by 2050 compared to 2017 (Zekollari et al., 2019). The response of glaciers to climatic variations is directly related to the mass and energy balance at the glacier surface. Internal deformation and basal sliding redistribute ice from regions with mass gain (accumulation) to regions with mass loss (ablation), compensating the local imbalance with some lag in the 40 response. Changes in glacier geometry, however, can be directly related to the mass balance if integrated across the entire glacier. Thus, the geodetic, glacier mean mass change (considering a volume to mass conversion) can be directly compared to glaciological mass balance measurements, based on local ablation and accumulation data.
We focused on a study site within the Ötztaler Alps, Austria, including the Vernagtferner, one of the reference glaciers in the World Glacier Monitoring Service (WGMS) system. Mass-balances have been determined here 45 using the glaciological method since 1965, while a series of historical maps back to 1889 demonstrates the longterm glacier evolution over more than a century.
Digital photogrammetry has been established as the standard method for cadastral surveys during the first decade of the new millennium. It is also used by BEV (Bundesamt für Eich-und Vermessungswesen, Austria) for their periodical digital elevation updates. Consequently, an immense imagery database has been generated, which 50 allows three-dimensional reconstruction and mapping of vast areas in the Alps with a resolution in the order of decimeters. Thereby horizontal and vertical changes can be detected and analyzed on a large scale and with unmatched precision. Combining high-resolution photogrammetry with the extensive knowledge data base for the Vernagtferner acquired over more than a century, this study presents a unique combination of methods that allows the extraction of relevant glaciological information (glacier retreat and mass balances) with greater accuracy from 55 aerial imagery. Data suitable for photogrammetric processing were available in three epochs over a period of 9 years (2009,2015,2018), and appropriate in situ observations were used to evaluate and improve the remote sensing information.
The workflow was applied for all 25 glaciers within the study area, including the Vernagtferner. Past studies already demonstrated the potential of spatiotemporal change analysis of alpine glaciers using photogrammetric 60 data (Fugazza et al., 2018;Legat et al., 2016;Rossini et al., 2018). Though, geodetic mass balances from photogrammetry are often restricted by acquisition times and, therefore, not comparable with fixed date glaciological measurements (Fischer, 2011). In this study, a correction due to differences in acquisition dates between photogrammetric and glaciological data is achieved by employing an additional unmanned aerial vehicle (UAV) survey. This calibration allows a parametrization between photogrammetry and glaciology, revealing 65 further potentials for these two techniques for accessing a more accurate insight into glacial retreat. https://doi.org/10.5194/tc-2020-263 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.

Study Area
The Ötztaler Alpen are located in the central-eastern Alps and represent one of the most extensive glaciated regions of Austria. The study area ranges from an altitude of 1.700 m to 3.768 meter above sea level (m.a.s.l.) and covers an area of more than 250 km² (Fig. 1). It combines the upper regions of the drainage basins of Rofental, Pitztal, 70 and Kaunertal. The location in the inner part of the Alps leads to relatively low precipitation amounts (Fliri, 1975), which for example, reach mean values of 660 mm yr -1 at the valley station Vent in 1969-2006(Abermann et al., 2009. For some of the 25 glaciers within the study site exist periods of mass balance measurements. One of the longest series of measurements is at the Vernagtferner, where regular monitoring by the Bavarian Academy of Sciences and Humanities (BAdW) began in 1965. This glacier is characterized by several sub-basins, which were 75 connected to a single glacier tongue in former times. Today the glacier covers an area of almost 7 km² within an altitude range between 2860 m and 3570 m. The mass balance is determined by the glaciological method, using measurements at ablation stakes, manual and geophysical depth soundings, and retrieving information from snow and firn pits. The annual and winter mass balances are determined independently of each other by measurements on the fixed dates of May 1 st and September 30 th , the dates of the glaciological balance year (Cogley, 2010;Cogley 80 et al., 2011). Besides, there is a long history of geodetic mapping at the Vernagtferner dating back to 1889 (Mayer et al., 2013a).

Photogrammetric data acquisition
After the introduction of the first large format digital aerial photogrammetric camera in 2000, this type of sensor quickly became the standard equipment for aerial surveys. The replacement of analog to digital equipment not only dramatically reduced the surveying costs, but also enabled the development of advanced three-dimensional 90 high-resolution reconstruction algorithms based on photogrammetry (Heipke, 2017;Hirschmüller, 2019). This kind of reconstruction works best with horizontal overlaps of more than 80%. Digital sensors allow surveys with this overlap because the high costs for photographic negatives no longer exist. Since the late 2000s, European state surveying agencies have been carrying out cadastral surveys every two to four years with digital sensors and horizontal overlap of 80%, thus creating an immense multitemporal database of aerial imagery to be explored both 95 commercially and scientifically.
In Austria, cadastral aerial surveys are conducted by the BEV with a nominal resolution of 20 cm. We use a BEV survey from 2015 as a basis for our investigations. Besides, surveys performed by 3D RealityMaps in 2009 and 2018 with the same or higher ground resolution are investigated. Table 1 shows the most relevant information on the conducted air surveys. 100 In addition to the airplane based surveys, a smaller test site was covered by UAV flights to retrieve high-resolution data at another acquisition date closer to the maximum of ablation in 2018 (Table 1). The processed area covers 6 km² containing most the glacier area of the Vernagtferner, including all glacier tongues. The glacier was almost snow-free at the time of the acquisition, which provides optimal processing conditions. 105

Glaciological data acquisition
Glaciological mass balance data are acquired at the Vernagtferner since 1965. Information is gathered as stake readings from ablation stakes for estimating the ice melt across the glacier. At the same time, snow depth and last 110 season firn deposits are determined from mechanical depth soundings with metal probes, which are then combined with density information from snow and firn pits to calculate the ice per water equivalent of the remaining snow and firn cover. While stake readings only require two length measurements per stake with an uncertainty of typically about 1 cm, more significant errors are included in the direct accumulation measurements due to uncertainties in the sample volume (about 5%) and the determination of density by using spring scales (about 4%). 115 https://doi.org/10.5194/tc-2020-263 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License. Therefore water equivalent accuracy is about 6% (Zemp et al., 2013). The typical number of stakes used for the annual mass balance measurements at the Vernagtferner is about 35, while some 4-5 accumulation measurements are collected at the end of the glaciological year, September 30th. The equilibrium line, which represents a mass balance of zero, is derived by comparing oblique terrestrial photographs of the transient snow line and firn extent with optical remote sensing information close to the date of the field measurements. The derived equilibrium line 120 has a horizontal location accuracy of about 10 m. Glacier boundaries for delineating the spatial mass balance distribution are derived from aerial surveys repeated roughly each decade, which are updated in the ablation region by annual GNSS (Global Navigation System Services) measurements of the glacier tongue geometry. The spatial error of these measurements is usually better than 1 m. The information from the stake readings, the depth soundings, the snow pits, and the location of the equilibrium line are combined to interpolate the spatial distribution 125 of the glacier mass balance into a raster file. Due to the sparse information in the accumulation region, it is necessary to manually correct the interpolation results in this region by the knowledge of the long term accumulation patterns, which are rather persistent. Errors introduced by uncertainties in the accumulation region, however, are rather small, especially during the recent decade where the accumulation area ratio is usually well below 30%. While ablation varies between 0 and up to 4500 mm in the ablation area, accumulation only varies 130 between 0 and about 300-400 mm in the accumulation area. Within this study we assumed the error of the interpolated raster to be 100 mm, which is in accordance to Zemp et al., 2013. It must be noted, that this relatively large error within the accumulation area will only affect the final mass balance by less than 2%.

Photogrammetric workflow 135
In order to determine geodetic glacier mass balances from aerial and UAV images, we use a workflow consisting of two main modules: the data processing and the vertical change analysis (Fig. 3). The main goal of the data processing module is the reconstruction from raw imagery data into 3D point clouds. Digital surface models (DSM) and true orthophoto mosaics (TOM) are then derived from these point clouds.   The first step in data processing is the aerotriangulation, which consists of the orientation of the aerial imagery to the real terrain. Modern photogrammetric survey systems deliver highly accurate position (GNSS) and orientation (IMU, inertial measurement unit) information for each image, which are also included in the aerotriangulation for 150 the georeferencing. Finally, at least 20 tie points were manually identified for each image block to enhance the orientation accuracy. For the large format imagery, the aerotriangulation was performed using the state of art software Match-AT by Trimble. For the UAV imagery, the software Metashape Professional from Agisoft was used.
The second stage of data processing is the generation of point clouds through three-dimensional reconstruction. 155 For this purpose, the state of the art Semi-global-matching algorithm (Heipke, 2017;Hirschmüller, 2019) was used. This algorithm is implemented within the software SURE (nFrames), which generates DSM and TOM from the point clouds. Positioning accuracies of TOM and DSM were computed based on ground control points from the BEV. The horizontal shift lies between 10 and 20 cm depending on the acquisition year and thus within the ground resolution of the images. 160 In order to calculate the DSM differences, existing systematic height shifts were corrected by using stable areas (mostly large rock formations) outside the glaciers. The 2015 DSM was chosen as the reference because it was derived from the official Austrian cadastral survey and is referenced to the Austrian national survey system. The mean shift in elevation compared to the reference DSM was computed based on 50 identified stable points, e.g., on solid rocks, for each DSM and 21 stable points for the UAV-survey. Based on this mean vertical shift over 165 stable ground, all DSMs except for the reference DSM were adjusted in height relative to the reference DSM of 2015. Subsequently, the DSM differences 9/ 2018-8/2015, 8/2015-9/2009, 9/2018-9/2009, and 9/2018-8/2018 were computed and named after their acquisition dates. https://doi.org/10.5194/tc-2020-263 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License.

Vertical change analysis
For the vertical change analysis (Fig. 2, orange part), glacier outlines were digitized visually using the TOMs at a 170 scale of 1:2.000. The geodetic glacier mass balance (B) for a single glacier is determined by integrating the DSM difference over the glacier outline at the beginning of the respective period. The result is divided by the mean glacier area of the period. More details can be found, for instance, in Fischer et al., 2015b. For the conversion in mm water-equivalent (w.e), the density assumption proposed by Huss, 2013, 850 kg m -3 ± 60 kg m -3 , was used. To determine altitudinal dependencies of the glacier mass balances, the glacier polygons were divided into 10 m 175 altitudinal intervals additionally. For each altitudinal interval, the mean absolute height change in meters was determined as well. For the conversion in mm w.e., an altitude-related density functions was used. It represents the gradual change from ice density (900 kg m -3 ) in the ablation region to firn density (550 kg m -3 (Cogley et al., 2011)) in the accumulation region, by using a linear transition zone of ±50 m across the equilibrium line altitude (ELA). 180

Comparison with glaciological data
To be able to compare the mass balances determined glaciologically and photogrammetrically, they must have the  (Table 2). To account for the temporal differences, which are limited 185 to the period of maximum ablation in August and September, the photogrammetric DSM difference 08/2018-09/2018 was used. From this, a correction function (sigmoid, see Fig. 3) was derived. For the regression, only those altitude levels were considered, which were at least 40 % covered by the UAV survey. The SD of the regression is 0.07 m ice a -1 (Fig. 3). Above 3180 m.a.s.l., thus, for the accumulation areas, the correction function is not based on any data and is therefore error-prone. Additionally, for the correction of the time periods, it had to 190 be assumed that the ablation process was constant during the period between August 21st and September 21st, 2018, and that the magnitude was comparable to the years of 2009 and 2015. Subsequently, the DSM differences were adjusted according to the number of days the investigation periods differ using a correction factor, and the correction function ( Table 2). The correction factor converts the number of 200 deviating days into months, thus represents the factor of which the monthly ablation (correction function) is added to the respective DSM difference.

205
Subsequently, for the density conversion of the adjusted DSM differences 2009-2018, 2009-2015, and 2015-2018, the altitude related density assumption (4.1.2) was used. Finally, the adjusted DSM differences were subtracted from the corresponding accumulated glaciologically derived rasters (sect. 3.2). The resulting Variation Rasters show the spatial deviations between the two methods. Using the Variation Raster, local methodical 210 deviations (e.g. debris cover or crevassed areas) were quantified. Therefore, the respective areas were manually digitized, and the mean value of the Variation Raster was calculated. This value can be related to the total mass balance to estimate the magnitude of methodological errors.

Vertical accuracy assessment
To assess the error distribution within the study area after the coregistration, 1.5 km² of ice-free, stable terrain were 215 analyzed representing a wide range of slope, aspect and altitude. Therefrom, the mean shift as well as the standard deviation within those areas was calculated and their relation to topography was investigated, following Nuth and Kääb, 2011. To estimate the error when averaging over extended areas, we followed Rolstad et al., 2009, by assessing the spatial covariance of the elevation differences with the use of semivariograms. For the application of the method, 220 we assumed that elevation difference are constant in space, no large-scale trends and that there is no significant variation of the variance in space.
Basic error propagation was implemented (Nuth and Kääb, 2011) in order to determine the compound error of DSM differences, density conversion (7 %, sect. 4.1.2, Huss, 2013), the correction function of the acquisition dates (SD = 0.07 m ice a -1 , sect. 4.2) as well as the error associated to the glaciological interpolation raster (sect. 225 3.2) for all presented results. The error within this study is indicated by the 95% confidence interval. Using the photogrammetrically derived glacier outlines, TOMs, as well as the DSM differences, the first results 230 are obtained by a visual interpretation for the entire study area. In general, glaciers were losing mass during the    The accuracy assessment conducted (sect. 4.3) allows an estimation of potential error related to the DSM differences and derived products. In general, the mean vertical error of the DSM differences ranges from -0.16 m to 0.1 m with SD not exceeding 0.42 m (Fig. 6a). As expected, SD increases with the slope (Fig. 6b). A clear 250 increase of the SD can also be found in lower altitudes (˂ 2600 m.a.s.l., (Fig. 6d), most likely due to the increasing influence of vegetation. A relation between aspect and the vertical error was found (Fig. 6c), that can be attributed to the horizontal shift of the DSMs (see sect. 4.1.1) (Nuth and Kääb, 2011).
Based on the error assessment, the confidence interval for all results within this paper was determined. Therefore, semivariograms were computed for the 1.5 km² bedrock areas for all observation periods, from which the range-255 value (100 m) was conducted and converted to the correlated area value (31416 m²). For more detailed information of this method see Rolstad et al., 2009.  (Fig. 7, bottom). However, it needs to be considered that the investigated periods do not cover full mass balance years, which influences the annual mean values considerably. Therefore height changes for the Vernagtferner were corrected (as described in sect.

Comparison with glaciological measurements
The photogrammetric, as well as the glaciological mass balances of the Vernagtferner, were further compared to detect discrepancies in spatial distribution or magnitude between both methods presented. Therefore, the DSM differences, as well as the overall geodetic mass balances, were compared to the equivalent accumulated 290 glaciologically derived mass balances. At the same time, the correction of the acquisition dates is evaluated. As already noted, the overall mass balance for the Vernagtferner is more negative in the period 2015-2018 than in the period 2009-2015. This can be seen in both, the glaciological and the geodetic mass balances (Fig. 8). The  (Table 2), the corrected photogrammetric data lies within the error bars of the glaciologic data. 300

Figure 8: Comparison between photogrammetric mass balances (blue) and glaciological mass balances (orange) for the investigation periods. Red bars visualize corrected photogrammetric data.
The spatial differences between both methods were further investigated using the Variation Raster described in 305 sect. 4.2. The Variation Raster represents spatially distributed deviations of the interpolated rasters based on the glaciological field measurements from the corrected DSM differences based on photogrammetric data (sect. 4.1.1, Fig. 10).
High spatial variations occur in debris-covered areas (red in Fig. 10), where photogrammetrically derived mass balances are higher than glaciological mass balances. Since no ablation stakes are located in supra-glacial debris 310 at the Vernagtferner, the glaciological mass balances for these areas are interpolated from surrounding information on clean ice. Debris cover above a certain thickness protects the glacier against incoming heat fluxes (Östrem, 1959) and therefore reduces the ablation. For the entire Vernagtferner, the neglect of debris-covered areas within a -1 (0.8 ± 0.6 %). 315 Further analysis of the Variation Raster (Fig. 10) emphasizes that photogrammetric and glaciologically derived mass balances do not only differ spatially due to debris cover or other external effects. It is also evident that the deviations depend on the altitude. Thus, in the lower reaches of the glacier, photogrammetric mass balances are higher than the glaciological mass balances. For higher altitudes, photogrammetric mass balances are lower than the glaciologically derived ones. This pattern is usually attributed to the ice dynamic component of elevation 320 change contained in the geodetic differences (submergence and emergence of ice and firn). Additionally, Fig. 9 clearly shows that comparing the different observation periods, the bias between both methods becomes larger within the accumulation areas.

345
Black outlines represent the different accumulation areas of the Vernagtferner; The mean submergence is given for each accumulation area in m w.e. a -1 , having an associated error of ± 0.25 m w.e.

Other glaciers within the study site
The major advantage of the photogrammetric mass balance assessment is that large areas can be investigated in detail. Therefore, we investigated the geodetic mass balance of all glaciers in the study area. The mean height 350 difference over the glacier surface was determined and converted into m w.e. as described above. The mean mass balance of all glaciers was -0.38 ± 0.05 m w.e. a -1 for 09/2009-08/2015, while it quadruples to -1.61 ± 0.22 m w.e. a -1 between 08/2015 and 09/2018. It must be noted that there was no correction applied for the acquisition times. Accordingly, the mass balances do not refer to identical annual periods. However, if a time correction is assumed to have a similar influence on the annual mass balances of other glaciers as it had for the 355 Vernagtferner (-27% for 09/2009-08/2015 and +21% for 08/2015-09/2018, Fig. 8), the mean mass balances of all glaciers can be estimated. Accordingly, the mean annual glacier mass balance within the study area would be -0.48 ± 0.06 m w.e. a -1 for the period 2009-2015 and -1.26 ± 0.06 m w.e. a -1 for the period 2015-2018. The comparative values for the Vernagtferner exceed those values by 51% and 21%, respectively (Fig. 8). Figure 11 compares the altitude-related height losses of the Vernagtferner to the mean of all glaciers of the study area and 360 the Hintereis-and Hochjochferner (Fig. 4 and Fig. 5).
We found that, for the Vernagtferner, glacier mass loss accelerated within the last decade. At the same time, the absolute surface height changes showed a maximum at low altitudes for all glaciers (Fig. 7). As it was shown for 370 the Vernagtferner, peak volume loss occurs at higher altitudes, compared to the altitude where maximum height losses occur, since the volume is directly linked to the area-height-distribution of a glacier. These findings correspond to general observations (i) for the Vernagtferner (BADW, 2019; Escher-Vetter, 2015) and (ii) for most of the glaciers in the Alps (Fischer et al., 2015b;Fischer et al., 2015a;Huss, 2012). This study showed the great potential of using photogrammetry for glacier studies. Thus, height changes can easily be investigated for large 375 areas with high spatial resolution and accuracy. As a result, for instance, we were able to estimate the decrease in the study site. Furthermore, dead ice bodies can be mapped very well using the high-resolution photogrammetric 380 data. For instance, on the north-facing side of the valley, below the glacier tongue of Hintereisferner, a large dead ice body is visible and ice loss could be determined for this feature, which is not included in the glaciological mass balance measurements (Fig. 4). Another result derived from the DSM differences is the vertical surface height changes as a function of elevation (Fig. 7). Even though they cannot be compared to glaciologically derived vertical balance profiles directly, they illustrate (i) the relation of height losses to altitude (Kaser et al., 2003;Pellikka and 385 Rees, 2010) and (ii) a shift between the two observation periods.
Moreover, the comparison of the uncorrected photogrammetric total mass balances reveals the great potential of photogrammetry for glacier analysis. A large number of glaciers may be analyzed with a high spatial resolution and for the exact same time. For the Ötztaler Alps, based only on the photogrammetric DSM differences, we found that the elevation change (i) is on average lower compared to the Vernagtferner, (ii) is correlated to the altitude of 390 the glacier tongue and (iii) increased between 09/2009-08/2015 and 08/2015-09/2018 (Fig. 11).
The error assessment conducted reveals the high precision of the photogrammetric data with SD not exceeding 0.42 m for all DSM differences (Fig. 6). Moreover, the variation of the vertical error is normally distributed. A relation was found for the vertical error regarding the aspect, which is presumably based on the vertical shift of the DSM (see sect.4.1.1). The error caused by the relation to the aspect can cause a larger error than the value we 395 suggest, especially for glaciers with homogeneous orientations. Nevertheless, the 95% confidence interval used in this study covers most of these variations. (Fig. 6). Future studies, that want to deduct this error, may follow the methodology suggested by Nuth and Kääb, 2011. Semivariograms determined in order to assess the error for averaged values prove that this assumption was possible, since semivariance approximates a clear sill (further information on the method can be found in Rolstad 400 et al., 2009).
We used the density assumption of Huss, 2013 to convert height changes to water equivalent. It must be noted that the DSM difference 08/2015 -09/2018 does not fulfill the prerequisite for this assumption, being a period shorter than three years. However, since the accumulation area ratio at the Vernagtferner is relatively small, the density of the lost volume will be close to the density of ice (900 kg/m³), this seems to be a good assumption. 405 To adjust photogrammetric mass balances to the acquisition time of the glaciologic measurements (September 30th), a correction was applied. An additional UAV survey was conducted, which made it possible to assess vertical height changes during the last month of the ablation period in 2018 (Fig. 3). Therefrom, a regression function (sigmoid) was determined to correct the offsets of the acquisition dates in 2009, 2015, and 2018. The correction includes an uncertainty for the accumulation areas since the UAV survey did not cover the entire glacier 410 area. Additionally, the intra-annual variation of melting and snowfall introduces another unknown, but important, source of error. More precisely, the transfer of the correction function from September 2018 to another year (here September 2009 and 2015) is strongly dependent on the meteorological conditions of the specific period. For instance, a snowfall event shortly before the survey would introduce an unknown error. The retrospective correction of the time periods therefore is not more than an approximation to the glaciological values. This 415 uncertainty must be taken into account, in particular for the shorter observation periods 2009-2015 and 2015-2018, where the correction has a higher impact (Table 2). This might explain that the deviations to the glaciological data increase for shorter periods (Fig. 8). Thus, we focused our analysis of the submergence on the period 2009 -2018, where the impact of correction is lowest (Fig. 10).
For transferability, individual correction functions from UAV surveys must be determined for each glacier, as 420 these are directly related to glacier-specific slope gradients, orientation, the height of the glacier tongue, and areaheight distributions (Fig. 11). Additional UAV surveys, which are low-cost and more flexible, but cover smaller areas, may be conducted to determine the correction function of any glacier within the study site. This allows the adjustment of the respective DSM difference to glaciological periods and therefore a comparison with the magnitudes of the glaciological data. Best results can be expected if the correction function represents the ablation 425 processes in the period that needs corrections. If the correction function is applied to photogrammetric mass balances retrospectively, intra-annual variation must be considered. In that case, if a more detailed analysis of the annual variations is needed, a mass balance model would produce more reliable results.
Despite the above-mentioned potential causes of errors, we showed that the presented workflow to correct the photogrammetric elevation changes to full multi-annual mass balances provides good results. These results agree 430 well with the glaciological mass balances (Fig. 8), but show characteristic deviations in the elevation distribution ( Fig. 9). Thus, they are (i) suitable for detecting systematic errors in the glaciological method (Fischer, 2011) and (ii) the method is transferable to other glaciers. This result is of great importance for future studies. Large-area, cadastral aerial surveys are usually recorded at any time in autumn and are therefore useless for verification of glaciological measurements. 435 For our case, the Vernagtferner, comparing the corrected photogrammetrically derived DSM differences with interpolated rasters from glaciological measurements, revealed small-scale differences as well as general offsets due to ice dynamics (Fig. 10). Such deviations usually occur in highly crevassed areas, in debris-covered areas, due to temporal changes in the snowpack within the accumulation area, internal and basal melting and interpolation errors of the glaciological method (Fischer, 2011). For crevassed areas, the deviations found at the Vernagtferner 440 are negligible. With regard to the influence of debris cover (about 3.8 % of the total glacier surface) on the glacier mass balance, the effect results to 0.8 ± 0.6 % of the annual total mass balance, which is not represented by the spatial interpolation of the glaciological measurements at the moment. Debris-cover might play an increasing role when debris-covered glacier areas increase due to further glacier decline (Scherler et al., 2018). This could increase the need to conduct high resolution geodetic observations in addition to glaciological measurements because point 445 measurements on debris-covered parts are not always representative. Surprisingly, the differences between geodetic and glaciological volume change in the accumulation areas do not reveal large spatial deviations, even though the glaciological results rely on very sparse in situ data. This demonstrates that the long-term accumulation conditions at Vernagtferner, which are known from former detailed investigations (Mayer et al., 2013a;Mayer et al., 2013b), are spatially relatively stable and can be used for scaling the point measurements. Even though this 450 comparison of geodetic and glaciological results in the accumulation region has the potential to improve the representation of spatial variability across these areas.
Assuming that internal and basal melting are negligible, and the influence of temporal snowpack differences is small due to the long period of investigation and the rather small ratio of accumulation area, the remaining differences result from the large scale dynamic processes. They become quantifiable when looking at the mean 455 differences per altitude of both methods (Fig. 9), whereby the high spatial resolution and accuracy is absolutely crucial. It is noticeable that the switch between emergence and submergence is about 100 to 150 m below the ELA of the Vernagtferner (BADW, 2019). While the ELA increased between the two periods (2009 -2015 vs. 2015 -2018) by about 60 m, the change from submergence to emergence decreased by about 90 m. The reasons for this https://doi.org/10.5194/tc-2020-263 Preprint. Discussion started: 27 October 2020 c Author(s) 2020. CC BY 4.0 License. development are unknown but might indicate a further deviation from balanced conditions during recent years. 460 Even though the derived submergence is depending on the correction function, the values for the Vernagtferner and its accumulation areas correspond to literature values and show differences between the individual accumulation areas, which support glaciological observations. (Lambrecht et al., 2011;Mayer et al., 2013b) Future aerial image acquisitions aiming at calculating glacier mass balances must be as close as possible to each other (same date within the year) and preferably at the end of the ablation season. Ideally, these acquisitions are 465 taken close to the standard glaciological mass balance data of September 30 th , to allow a direct comparison with available field information. Since aerial image surveys require cloud-free weather, it will not always be possible to acquire the images on this exact date. For this reason, temporal corrections, as successfully implemented within this study, will play an essential role in future studies. Additionally, future photogrammetric data acquisitions should be as large as possible, but at least cover the entire glacier area. 470

Conclusion
This study demonstrates the high potential of aerial images and the resulting DSM's for analyzing glacier retreat in great spatial and temporal detail. For all 25 glaciers within the study area, height changes were analyzed with a 20 cm spatial resolution and compared to each other. It was shown that the glacier retreat does not only take place in low altitudes, but that even high accumulation basins are meanwhile affected due to non-compensated ice flow. 475 The altitude of the most significant glacier height changes, as well as the amount of maximum volume loss, has increased within our observation period. If all aerial surveys that are available for large parts of the Alps were used, the majority of all alpine glaciers could be analyzed by the presented method in the future. The spatial resolution of the respective analysis would be significantly higher than in other recent satellite-based studies. The photogrammetric results were compared with glaciological in-situ measurements. The required correction of 480 acquisition dates was successfully applied by using an UAV survey. This shows that the main disadvantage of using photogrammetry, not being able to survey at the right time can be reduced by using additional UAV surveys, which can be better targeted in time. This study also shows that results from the glaciological method can be greatly complemented with photogrammetric analyses, which increases the accuracy of the glaciological mass balance series, reveals regions of anomalous mass balance conditions, and allows estimates of the imbalance 485 between mass balance and ice dynamics.
Data availability. Except for the aerial imagery and derived products all datasets are available on request.