Uncertainty assessment of a permanent long-range terrestrial laser scanning system for the quantification of snow dynamics on Hintereisferner (Austria)

A permanently installed terrestrial laser scanner (TLS) helps to investigate surface changes at high spatio-temporal resolution. Previous studies show that the annual and seasonal glacier volume, and subsequently the mass balance, can be measured by TLSs. This study systematically identifies and quantifies uncertainties and their sources of the permanent long-range TLS system at Hintereisferner glacier (Ötztal Alps, Austria) in order to assess its potential and limitations for detecting glaciologically relevant small-scale surface elevation changes, such as snowfall and redistribution events. Five uncertainty sources are analyzed: the registration method, the influence of the instrument and hardware limitations of the TLS, the effect of atmospheric conditions on the laser beam, the scanning geometry, and the uncertainty caused by rasterization. The instrument and hardware limitations cause the largest uncertainty to the TLS data, followed by the scanning geometry and influence of varying atmospheric conditions on the laser beam. The magnitude of each uncertainty source depends on the distance (range) between the TLS and the target surface, showing a strong decrease of the obtained spatial resolution and a concurrent increase in uncertainty with increasing distance. An automated registration method results in an uncertainty of ±0.50 m at grids of 100 by 100 m. After post-processing, a 0.1-m vertical accuracy can be obtained allowing the detection of surface changes of respective magnitudes and especially making it possible to quantify snow dynamics at Hintereisferner.

A permanently installed terrestrial laser scanner (TLS) helps to investigate surface changes at high spatio-temporal resolution. Previous studies show that the annual and seasonal glacier volume, and subsequently the mass balance, can be measured by TLSs. This study systematically identifies and quantifies uncertainties and their sources of the permanent long-range TLS system at Hintereisferner glacier (Ötztal Alps, Austria) in order to assess its potential and limitations for detecting glaciologically relevant small-scale surface elevation changes, such as snowfall and redistribution events. Five uncertainty sources are analyzed: the registration method, the influence of the instrument and hardware limitations of the TLS, the effect of atmospheric conditions on the laser beam, the scanning geometry, and the uncertainty caused by rasterization. The instrument and hardware limitations cause the largest uncertainty to the TLS data, followed by the scanning geometry and influence of varying atmospheric conditions on the laser beam. The magnitude of each uncertainty source depends on the distance (range) between the TLS and the target surface, showing a strong decrease of the obtained spatial resolution and a concurrent increase in uncertainty with increasing distance. An automated registration method results in an uncertainty of ±0.50 m at grids of 100 by 100 m. After post-processing, a 0.1-m vertical accuracy can be obtained allowing the detection of surface changes of respective magnitudes and especially making it possible to quantify snow dynamics at Hintereisferner. KEYWORDS topographic lidar, RIEGL VZ-6000, uncertainty assessment, terrestrial laser scanning, cryosphere, atmosphere hydrology and sea level rise. Furthermore, it is used to assess respective climate change impacts and is subjected to detection and attribution studies. In order to calibrate and validate highresolution distributed glacier mass balance models, detailed information of glacier surface changes is needed (Klok and Oerlemans, 2002;Hock and Holmgren, 2005;Machguth et al., 2006).
The climatic mass balance information on the surface is often acquired with the traditional "direct glaciological method" (Kaser et al., 2003;Cuffey and Paterson, 2010;Cogley et al., 2011) by measuring the relative surface elevation changes at 10 to 50 selected locations (points) on the glacier. The measurements are transferred to mass changes in combination with density measurements of snow and ice and are then extrapolated to the glacier-wide surface. This "glaciological method" is laborious and resolves the complexity of spatial mass balance patterns only coarsely and is only grossly able to mirror the driving processes. On the contrary, the "geodetic method, " carried out from terrestrial or airborne platforms (Geist and Stotter, 2007;Fischer et al., 2016;Klug et al., 2018), allows for detecting comparatively high-resolution surface elevation changes relative to the surrounding bedrock. Again, obtained volume changes need measured, modeled, or assumed density information for conversion into mass changes. However, the geodetically measured glacier mass balance contains, beyond the climatic signal, the internal and basal changes in the ice column. In addition, it also contains effects of ice flow divergence, which is the elevation change component originating from the ice flux to maintain mass continuity (Cuffey and Paterson, 2010, Sect. 8.5.5). The ice flow divergence is omitted by integrating the elevation changes over the entire glacier area (Kuhn et al., 1999;Zemp et al., 2010;Klug et al., 2018). The effect of ice divergence on the point glacier mass balance may be neglected on glaciers with moderate to low ice flow dynamics and over short time periods such as daily up to monthly time intervals (Cuffey and Paterson, 2010). If, as on many mountain glaciers, internal and basal mass changes can also be neglected, the geodetically obtained mass changes can provide spatially and temporally distributed information on the climatic mass balance.
For most temperate and land-terminating glaciers, the representation of the snow cover dynamics poses the central deficiency in distributed mass balance models (Greuell and Bohm, 1998;Machguth et al., 2006;Carturan et al., 2012;Molg et al., 2012;Gurgiser et al., 2013;Ayala et al., 2015;Prinz et al., 2016). It crucially impacts the spatially distributed glacier surface mass balance by modulating accumulation and, subsequently, also ablation patterns. Data quantity and quality both from ground measurements and remote sensor systems have not yet been sufficient to reflect the actual snow cover processes on a glacier scale so far, which also limits the calibration and evaluation of distributed mass balance models (Machguth et al., 2006;Carturan et al., 2012;Gardner et al., 2013;Zemp et al., 2013).
A possible way to overcome this deficiency and acquire more information on surface elevation changes over the glacier is the measurement by a terrestrial laser scanner (TLS). A TLS emits a laser beam as the active sensing carrier. The distance between the sensor and the surface target is derived from the travel time of the laser beam and point clouds are created with a high point density (>1 points m −2 , depending on distance and scan settings) (Wehr and Lohr, 1999). Several studies have already proven that the total annual and seasonal mass balance can be acquired by TLS (Fischer et al., 2016;Wang et al., 2018;Xu et al., 2019). Since only the processes larger than the uncertainty of the TLS can be measured, an uncertainty assessment of the TLS is mandatory.
Various studies (Soudarissanane, 2016;Friedli, 2020) have already stated the main sources for uncertainty in TLS point cloud data: 1) instrument and hardware limitations, 2) scanning geometry, 3) atmospheric conditions between TLS and target surface, 4) surface reflectance properties, and 5) post-processing including registration and georeferencing. Additionally, the rasterization of point clouds contributes to the uncertainty budget when working with a digital elevation model (DEM). Previous research assessed only selected uncertainty sources and at shorter distance ranges than the system at HEF (Soudarissanane, 2016;Friedli, 2020;Kuschnerus et al., 2021;Dong et al., 2020;Schaer et al., 2007). Here, we addressed the full spectrum of the main uncertainty sources of TLS point clouds and derived DEMs and at distances up to 4,500 m. This paper aims to: 1) identify and quantify the main sources of uncertainty for measuring surface changes with a permanent long-range TLS system and 2) assess the potential of the permanent long-range TLS system at HEF for detecting glaciologically relevant surface elevation changes.
The findings contribute to a better understanding of requirements and limitations of permanent TLS systems for observing glacier surface dynamics.

Study area
The research area is the HEF (Figure 1), a valley glacier located in the Rofental catchment in the Ötztal Alps (Austria). The glacier has a length of approximately 6,300 m, ranging between the Weißkugel mountain peak (3,739 m a.s.l., green triangle in Figure 1C) and the glacier tongue at 2,460 m a.s.l. (2018). The glacier has been the key research site of the University of Innsbruck for glaciological studies since the early days of glacier research (Blumcke and Hess, 1899) and is now a part of the wider Open Air Laboratory Rofental (Strasser et al., 2018). Continuous long-term mass balance observations dating back to the year 1952/53 (Kuhn et al., 1999) and velocity and ice thickness change measurements since 1895 (Span et al., 1997) are available. The glacier is classified as one of the key "reference glaciers" by the World Glacier Monitoring Service (WGMS) (Zemp et al., 2009 The catchment is instrumented with several permanently operated and automated weather stations and rain gauges. The most important station for this study is Im Hinteren Eis (IHE, 3,245 m a.s.l.), which is located on the orographic right side of the glacier on the ridge at the Austrian-Italian border (blue dot in Figure 1C) and was installed in September 2016. The main instrument at this location is the year-round permanently installed TLS, which was installed in September 2016 and is in operational use for monitoring surface elevation changes on the glacier since 2019. The TLS is mounted on a stable frame positioned in a small container. The scanning process was automated in 2020, along with the setup, as described in detail by Voordendag et al. (2021). The TLS is a RIEGL VZ-6000 (RIEGL, 2019a), which can scan up to 6,000 m. Due to its laser wavelength of 1,064 nm, the instrument is exceptionally capable of measuring snow and ice. To date, the RIEGL VZ-6000 at HEF is the only permanent installation of such a long-range system worldwide, but we are also aware of the permanently installed RIEGL VZ-2000 TLSs (ranging up to 2000 m) at coast lines (Vos et al., 2017;Anders et al., 2019;Deruyter et al., 2020) and in a boreal forest (Campos et al., 2021).
Furthermore, Im Hinteren Eis is equipped with an automatic weather station (AWS), which is located 50 m horizontally and 25 m vertically away from the TLS. The AWS provides all common meteorological data at a 1-min resolution and also includes a 3D sonic anemometer to measure turbulent fluxes. Last, two webcams are installed that take pictures of the glacier every 30 min 1 .
A second important measurement location is Station Hintereis (StHE, 3,026 m a.s.l.), built in 1966 and located at the orographic left side of the glacier (Strasser et al., 2018). This station presently hosts an AWS with a 10-min average of common meteorological variables.

Analyses of five uncertainty sources
We are mainly interested in the uncertainty of the gridded data, as they ease the comparison of scans to identify vertical changes at the glacier. The measurements are subjected to six potential uncertainty sources (Table 1). However, this study does not account for the surface reflectance properties, as RIEGL VZ-6000 delivers good reflectance from snow and ice surfaces. As we are mainly interested in surface elevation changes during the season where the glacier is covered in snow, we assume a good reflectance and a negligible respective influence on surface elevation changes (Kaasalainen et al., 2008;Fritzmann et al., 2011). Thus, five uncertainty sources remain, and they will be discussed individually in this study. The uncertainty sources σ in this study are calculated as uncertainties in vertical direction of a single laser pulse or gridded data and are given in meter.
The uncertainty that applies to a single-point measurement originates from the laser pulse source, the medium of the laser beam, and/or the target surface. The laser pulse source is the TLS, which has to deal with instrument and hardware limitations (Soudarissanane, 2016) (σ instrument , Sect. 3.3). Atmospheric perturbations due to 1 https://www.foto-webcam.eu/webcam/hintereisferner1/ fluctuations and changes in air pressure, air temperature, and relative humidity along the path of the laser beam influence the velocity and pathway of the laser pulse in the medium (Friedli, 2020) and thus the calculated distance of a laser pulse (σ atm , Sect. 3.4). The target surface is the footprint of the laser beam (Schaer et al., 2007;Sheng, 2008) and is elaborated as the scanning geometry in this study (σ geo , Sect. 3.5). These three uncertainty sources are assumed to be random, independent uncertainties. The assumption of the independence of these uncertainty sources is a simplification, as they are all affected by the distance from the TLS to the surface R (800-4,500 m), which is several orders of magnitudes larger than these uncertainties (see Sect. 3.3-3.5). So, the influence of a small uncertainty in R, say R ± δR, is negligibly small on σ instrument , σ atm , and σ geo , and R can be seen as an independent constant, which differs for every single laser pulse. After this simplification, the error propagation law for addition is applied in Eq. 1 to calculate the total error propagation of a single laser pulse.
Subsequently, the acquired point clouds are gridded to DEMs by taking the mean of all the points in the 1-m grid cell. The number of points per grid cell is called the point density and decreases with the distance from the TLS. Likewise, if the point density increases, the uncertainty decreases, as the mean is taken over a higher number of points. The point density can be adjusted by changing the scan settings, but in our study, the angular step width of the TLS is set to 0.01°and not changed. The point density ρ point can be seen as a repeated measurement and thus, due to the error propagation, the uncertainty of σ point must be divided by the square-root of the number of measurements (Joerg et al., 2012) to get the uncertainty of every grid cell σ gridcell (Eq. 2).
The calculated σ gridcell does not yet include the uncertainty caused by the registration (see Sect. 3.2). The registration of the point clouds introduces a systematic error and is a rigid transformation and translation applied to the point clouds. Opposed to the other random uncertainty sources, the registration uncertainty depends on the method the user applies to the data, independent of the previously calculated σ gridcell (Joerg et al., 2012;Fey and Wichmann, 2016). Therefore, the variable σ reg is added to Eq. 2 to acquire Eq. 3.
The registration uncertainty cannot be distinguished from the other uncertainty sources (Fey and Wichmann, 2016), and is therefore first assessed in Sect. 3.2.

Available data
For the uncertainty analyses, a period with stable surface conditions and as many TLS scans as possible are needed. Stable surface conditions are defined as a period where surface change due to snow drift or melt is unlikely to happen. Temperatures Frontiers in Earth Science 04 frontiersin.org have to be low and the snow grains should already be well bonded throughout the snowpack. These conditions were met on 5 and 6 November, 2020, when the average temperature at IHE was −2.9°C (at 1.5 m above ground), the average wind speed was 3.4 m s −1 (at 3 m above ground), and no precipitation had been recorded since 27 October. The webcams did not give any evidence for blowing snow. Overall, 28 scans, taken between 5 November 10:55 UTC and 6 November 2020 15:05 UTC, were available. The scan taken on 5 November at 10:55 UTC is defined as the reference scan, as this is the first scan of the hourly TLS acquisitions.

Registration of scans
Acquired scans need to be registered in a common coordinate system and preferably also georeferenced in a global coordinate system. Previous publications employing a RIEGL VZ-6000 (Gabbud et al., 2015;Fischer et al., 2016;Xu et al., 2019) use RIEGL RiSCAN PRO software (RIEGL, 2019b) to manually register and georeference the scans. Unfortunately, this processing software cannot be implemented in an automated time series analysis. In order to overcome this limitation, to make longer time series analyses feasible and to test the effect of the registration method on the data, we designed an automated georeferencing tool chain in the System for Automated Geoscientific Analyses (SAGA GIS) (Conrad et al., 2015;Laserdata, 2022) based on an iterative closest point (ICP) algorithm (Besl and McKay, 1992), which we name ICP HEF . In turn, the approach using RiSCAN PRO is called ICP RiSCAN . For running ICP HEF , a fix transformation matrix for the scanner's orientation and position (SOP) is derived for a single selected scan using the RiSCAN PRO module multi-station adjustment (MSA) (Prokop and Panholzer, 2009;RIEGL, 2019b). All further scans are then adjusted with this fix transformation matrix. ICP HEF assumes that the position of the TLS is stable, and thus, the same transformation matrix is applicable to all scans, whereas ICP RiSCAN slightly changes the transformation matrix after manually applying the MSA algorithm. The advantage of ICP HEF is that only one scan has to be registered manually to an existing project.
An alternative automated registration approach would be with the OPALS software (Pfeifer et al., 2014), where stable areas are selected at the surface of interest and a transformation matrix is based upon these stable areas. With this SOP matrix, the entire scan is referenced, but this approach is not possible at HEF, as most areas, except the walls of the research hut at StHE, get covered by snow. Another stable area is the mountain ridge in the direction of Weißkugel, but this is far away (>3,000 m) and, thus, barely reliable. Likewise, another possibility would be to install reflective targets, but this is not possible at HEF, as the footprints are large and thus large reflective targets are needed. These targets are also likely to get covered with snow or will move as they are not installed on a stable surface.
After registration and georeferencing, the point clouds are gridded to DEMs at 1 m resolution by taking the mean over all the points within a grid cell. A total of 20 areas surrounding and on the glacier are selected for the comparison between scans ( Figure 1C). The selected areas cover 100 by 100 m and contain 10,000 pixels each. The differences between the 28 scans and the reference scan are calculated, and the average difference of the areas relative to the reference scan are calculated and reported. This is done for both the registration methods with ICP RiSCAN and ICP HEF , respectively.
The mean differences relative to the reference scan of the 100 by 100 m areas with ICP RiSCAN are shown in Figure 2A. They range between −0.15 and +0.04 m in a vertical direction, even though zero change was expected. The largest deviations relative to the reference scan (red encircled in Figure 2) are found for the areas furthest away from the TLS, whereas the closer areas (yellow to orange colors in Figure 2A) have deviations closer to 0.0 m.
ICP HEF shows different results, as the deviations range between −0.62 m and +0.47 m ( Figure 2B). The deviations relative to the reference scan increase if the distance from the selected area to the TLS is increasing, so the deviations can strongly be assigned to the scanning distance. The changes observed in Figure 2B are compared to the measurements of the internal inclination sensors of the TLS. The mean roll and pitch of every scan has been calculated ( Figure 2C, D), and the pitch anomaly of the scans shows a similar pattern as the mean differences in Figure 2A. This shows that the assumption that the position of the TLS is stable during consecutive scans is not valid because the TLS moves slightly.
Taking a closer look at single scans with ICP RiSCAN (Figure 2A), it can be seen that the differences between the 20 areas range between 0.06 m (on 6 November at 4:30 UTC) and 0.19 m (on 5 November at 16:01 UTC). The smallest ranges of difference between the 20 areas are found at night, which might indicate differences caused by the atmospheric conditions or incoming solar radiation. The deviations have been compared to air temperature and wind data from the nearby AWS IHE. The temperature data showed an increase of approximately 4°C after 6 November, 6:00 UTC, but no correlation with the deviations in Figure 2A was found. Similarly, higher wind speed in the night between 5 and 6 November and a peak in the turbulence kinetic energy on 6 November at around 3:00 UTC did not show correlation with the deviations relative to the reference scan.
It can be concluded that the average vertical accuracy of the TLS data at grids of 100 m is ±0.10 m for scans registered and georeferenced with ICP RiSCAN and ±0.50 m for ICP HEF . The accuracy strongly depends on both the distance from the TLS to the surface and the applied registration method. However, the absolute magnitude of the registration uncertainty σ reg remains elusive as, in addition to the remaining registration uncertainty, the other uncertainty sources are also included in this calculation.

Instrument and hardware limitations
Preliminary results applying ICP HEF show stripes which extend radially from the position of the TLS to the surface in the DEMs, particularly visible in difference plots between two DEMs, i.e., at the glacier tongue, around areas 17-19, in Figure 1C. The magnitude of the stripes in the vertical direction is about ±0.10 m in a difference plot ( Figure 1C).
The radial stripes are possibly caused by wind, turbulence, and thermal changes of the mounting platform of the TLS during data acquisition (Kuschnerus et al., 2021;Voordendag et al., 2021). The exact source is not investigated in this study, as accurate measurements of wind, temperature, and turbulence were unavailable at the exact location of the TLS and, thus, the radial stripes are only investigated with the data of the internal inclination sensors. It is hypothesized that these stripes are caused by movements of the TLS at high frequencies. The internal inclination sensors measure at 1 Hz and have an accuracy of ±0.008°. The TLS measures the time-of-flight of a laser pulse and calculates the distance from the TLS to the surface with it. The position of this laser pulse is calculated with the beam direction and beam origin. The beam direction and origin of the laser pulse should be adjusted by using the measurements of the internal inclination sensors, if the TLS moves. The x-, y-, and z-coordinates in the Scanner's Own Coordinate System (SOCS) are called vertices in RIEGL's RiVLib (Eq. 4, where i denotes x, y, or z) (RIEGL, 2013).
Regardless of the adjustments made in the beam direction and origin, an uncertainty attributed to the radial stripes σ instrument remains for every single laser pulse. It is assumed that this depends on the distance R and the uncertainty of the internal inclination sensor as given by the manufacturer of ±0.008°(Eq. 5). We consider these radial stripes as the total uncertainty attributed to the instrument and hardware limitations. Figures 3A, B show the measurements by the inclination sensors in roll (x) and pitch (y) of the reference scan, respectively. The variation in the inclination data indicates high-frequency vibrations. The TLS at IHE is installed slightly tilted to adjust its fieldof-view downward to the glacier. A sinusoid is plotted through the inclination data with a least-squares method (Figures 3A, B).
The shapes of the curves of the beam direction (Figures 3E, G) and beam origin (Figures 3H, I) in the x-and y-directions depend on the rotation of the TLS around its tilted z-axis (Lichti and Shaloud, 2010). The beam direction and origin in the z-direction (Figures 3G, I) depend on the rotation of the mirror in the TLS and show the relief of the mountain ridge. Although the measured inclination data of the reference scan is noisy, Figures 3E, F show a rather smooth beam direction. The frequency of the inclination data is 1 Hz, and these data are not used to correct individual scan lines, but a mean value is calculated from the data [B. Groiss (Riegl), personal communication, 29 October 2020]. This value is used to level the instrument by means of an SOP matrix and shows the good correspondence between the inclination data and beam direction. This means that small movements in the scanner are not corrected in the individual scan lines and, thus, the individual scan lines show deviations in the form of radial stripes. These deviations become especially apparent if the difference between two DEMs with a grid cell size of 1 m is calculated.
Theoretically, it would be possible to correct the beam direction with the inclination measurement, but this is not feasible as the effective measurement rate of the TLS is 23,000 measurements per second (pulse repetition rate of 30 kHz) and the frequency of the inclination data is only 1 Hz. In addition, the starting time of the inclination sensors cannot directly be coupled with the start of the TLS measurement, as an approximate 5-second difference is found between the total scanning time and the number of inclination measurements.
The uncertainty related to the instrument and hardware limitations σ instrument , according to Eq. 5, is found between 0.09 and 0.64 m and only depends on the distance from the TLS to the surface (Figure 4).

Influence of atmospheric conditions on the laser beam
The distance R of a laser pulse is measured with the time-of-flight of a laser beam to the surface τ and the velocity of the laser beam in the propagating medium c g (Eq. 6) (RIEGL, 2019b).
The laser beam is delayed by the group index of refraction n g relative to the speed of light in vacuum c 0 . n g depends on the group refractivity N g at 0°C, 1,013.25 hPa, and 0% relative humidity (RH) and the wavelength λ (1.064 μm) (Eq. 7) (RIEGL, 2019b). The group refractivity N L at ambient moist air at temperature T (°C), pressure p (hPa), and vapor pressure e is calculated with Eq. 8 (RIEGL, 2019b), and this leads to the group index of refraction (Eq. 9) (RIEGL, 2019b The RIEGL VZ-6000 is able to correct for n g , but assumes constant temperature, pressure, and humidity along the path of the laser beam. This assumption is not valid in mountainous terrain at long ranges as air temperature and air pressure differ widely over this terrain. The atmospheric influence on the laser beam can only be assessed with spatio-temporal information at high-resolution of the relevant atmospheric variables (air temperature, RH, and air pressure) along the laser beam path. However, three-dimensional observations of the atmosphere over the glacier are not available. To overcome this data gap, numerical simulations with the Weather Research and Forecasting (WRF) model are performed over HEF and surroundings. A detailed model evaluation study (Goger et al., 2022) within a measurement campaign over the glacier (Mott et al., 2020) has shown that the WRF model is able to simulate the general atmospheric boundary layer structure over HEF successfully. The same high-resolution large-eddy simulation (LES) setup is used with a horizontal grid spacing of Δx = 48 m for 5 and 6 November, 2020. The model output is used as three-dimensional, realistic information of atmospheric variables above the glacier. A path along the laser beam through the atmosphere is interpolated at increments of approximately 13 m. It is assumed that the last 10 m of the laser beam travels through the glacier boundary layer. The temperature in these last 10 m is linearly interpolated between the simulated temperature at 2 m and the simulated surface temperature. The simulation of the air temperature, RH, and air pressure is compared against the observations, similar to Goger et al. (2022) and enables the calculation of the changing n g along the laser beam path.
The difference in path length ΔR relative to the entire distance R is consequently calculated with Eq. 10 (Friedli, 2020)

FIGURE 4
Distributed TLS uncertainty related to instrument and hardware limitations (σ instrument ).
The laser beam is also refracted due to the changing n g . The angle of refraction γ is calculated with the derivative of the refractive index perpendicular to the path dn(s) dq (Eq. 11) (Friedli, 2020). It is decided not to simulate the refraction of the laser beam in our study, as the derivative of the refractive index perpendicular to the path is needed, and this is not accurate in our high-resolution simulation with a horizontal grid spacing of Δx = 48 m.
The calculated change of the laser beam distance is not yet the vertical component of the uncertainty caused by the atmospheric conditions. The vertical component of the uncertainty as attributed to atmospheric effects σ atm is calculated with slope γ and incidence angle θ in Eq. 12 and clarified with Figure 5A. σ atm = ΔR ⋅ cos θ cos γ (12) Figure 6 shows a cross-section of the simulated air temperature for area 4 on 5 November at 12:00 UTC, and 6 November 00:00 and 07:00 UTC. This shows the high variability of the air temperature over space and time. The green star indicates the position of the TLS and the green line indicates the laser beam path in the simulation. The position of the TLS and the length of the laser beam path in the simulation slightly differ from the reality ( Table 2) because the model uses a coarser topography than that measured with the TLS (Figure 6). The simulated data show a decrease in the length of the laser beam between −0.0029 and 0.0194 m (Figure 7). The anomalies increase with increasing distances but are relatively small compared to the uncertainty caused by the instrument and hardware limitations and, as we will later see, the scanning geometry. However, temperatures having the largest influence on the velocity of the laser beam through the atmosphere (Friedli, 2020) were close to the default temperature of −2.0°C of the TLS at 5 and 6 November (−2.9°C on average). Therefore, the change in the laser beam distance with a temperature increase and decrease of 10°C was also calculated. This shows changes in distance as low as −0.0424 m (−10°C) and as high as +0.0210 m (+10°C), respectively. A similar sensitivity study was performed for RH (±25%) and air pressure (±10 mbar). The simulated air pressure only varies with the altitude and is relatively constant over time. The simulated relative humidity varies more, but uncertainties due to variations in the relative humidity are so small that they can be neglected (Friedli, 2020). Thus, changes in air pressure and RH have no significant influence on the change of distance measurements.

Scanning geometry
This section deals with the incidence angle and the consequent footprint of the laser beam on the target surface to quantify the contribution of the scanning geometry. The shorter the distance and the closer the incidence angle to the normal, the smaller the respective uncertainty and vice versa (Soudarissanane, 2016).
The footprint of the laser beam is an ellipse under the assumption that the illuminated surface is a plane, and the laser beam is not perpendicular to this plane. The incidence angle θ to the surface has been calculated with the terrain normal n and the laser beam direction l in Eq. 13 (Schaer et al., 2007).
A normal of the terrain was calculated with the slope (γ) and aspect (α) (Zevenbergen and Thorne, 1987) from the DEMs and Eq. 14.
The extent of the major axis M of this ellipse depends on the incidence angle with the surface from the TLS, the distance R, and  . 3.4). The extent of σ atm is calculated with slope γ and incidence angle θ. (B) AB is the major axis of the elliptical footprint, θ the incidence angle, and β the beam divergence.

FIGURE 6
Temperature plot for the cross-section from TLS to area 4 for 5 November 12:00 UTC, and 6 November 00:00 and 07:00 UTC. The black area is the topography as used in the simulation, with white areas indicated as the glacier area. The blue line is the topography as measured with the TLS. The green line is the laser beam in the simulation and the green star is the position of the TLS in reality.
Frontiers in Earth Science 10 frontiersin.org

FIGURE 7
Change in the distance traveled by the laser beam in mm for 20 different areas over time (Y-axis). The legend shows the distance of the TLS to the simulated coordinates. The vertical gray lines denote the time a new scan was started. The temporal resolution of LES simulation is 15 min.
the beam divergence β of the laser pulse. The extent of the minor axis m of the ellipse only depends on the distance and beam divergence. The beam divergence β of the RIEGL VZ-6000 is 0.12 mrad (RIEGL, 2019a). The major axis of the footprint is derived using Eq. 15 (Sheng, 2008), where H, OA, and OB comprise the geometry from Figure 5B; the extent of the minor axis is derived using Eq. 16 (Sheng, 2008).
Last, the footprint area A F , which is the illuminated spot on a planar surface from the laser beam, is calculated with Eq. 17.
To calculate the contribution of the scanning geometry to the entire uncertainty budget, the major axis is decomposed in x-, y-, and z-components. The z-component is the main interest, as we are interested in vertical changes over the glacier. The vertical component is approximately 1/3 of the vertical component of the major axis M Z (Schaer et al., 2007). We define this extent as the vertical uncertainty caused by the scanning geometry σ geo (Eq. 18).
The distance measured by the TLS and as used in the LES (Sect. 3.4), the footprint area, mean altitude above the sea level, and incidence angle for the boxed areas from Figure 1C are given in Table 2. The total footprint area (Eq. 17), which is the illuminated spot on the planar surface, increases up to an area of 1.33 m 2 . Area 20 has a very large σ for the footprint, which is caused by the incomplete coverage of the 100 by 100 m area and the high incidence angle, and thus, it is assumed to be an outlier in the dataset.
After the decomposition of the x-, y-, and z-components and Eq. 18, the contribution of the scanning geometry in the uncertainty budget is approximated between 0.02 and 0.25 m (Figure 8) in the vertical direction for a grid size of 1 m. σ geo is between 0.20 and 0.25 m close to Weißkugel and 0.02 and 0.05 m in the area between IHE and StHE.

Rasterization
In the previous three sections, all uncertainty sources that apply to a single pulse are described. These uncertainty sources are now

FIGURE 8
Distributed TLS uncertainty related to the scanning geometry (σ geo ).

FIGURE 9
Uncertainty of a single laser pulse distributed over the glacier (σ point ).
combined to form the total uncertainty σ gridcell , without accounting for the registration uncertainty. The uncertainty of a single laser pulse is found in Figure 9. As the uncertainty depending on the influence of the atmosphere on the laser beam over time is highly variable and depends on the distance from the TLS, we have taken a value of 0.01 m for σ atm in Eq. 1 as this is approximately the average of the change in distance from Figure 7. Furthermore, σ atm is an order of magnitude smaller than σ instrument and σ geo , and so, the influence of the atmosphere is relatively small. The uncertainty of a single laser pulse ranges between 0.09 m at the area closest to the TLS and 0.66 m in the accumulation zone close to Weißkugel. These values are close to the values found in Sect. 3.3 and show that σ instrument is the largest contributor to the uncertainty budget. Furthermore, all three uncertainty sources (σ instrument , σ geo , and σ atm ) strongly depend on distance R from the TLS.
The total uncertainty of the DEMs strongly depends on the point density. The point density found in Figure 10 is between 1 point m −2 in the accumulation zone, furthest away from the TLS and above 13 points m −2 directly at the glacier in the line between IHE and StHE.
The total uncertainty of the DEMs improves if the point density is taken into account (Eq. 2; Figure 11). The uncertainty at the glacier tongue closest to the TLS is only 0.012 m but remains at 0.66 m in the accumulation zone. This is caused by the low point density of one or two points m −2 in the accumulation zone, and thus, the uncertainty is close to the value found in Figure 9.
In the calculation of the total uncertainty of the DEMs, the registration uncertainty cannot be accounted for (Sect. 3.2).
Nevertheless, the registration uncertainty with ICP RiSCAN is in the order of a few centimeters, whereas it is in the order of decimeters for ICP HEF . As the uncertainty of the registration method is taken as an average over the 20 areas of 100 by 100 m (Sect. 3.2), the effect of the instrument and hardware limitations is evened out and is not visible anymore. Thus, in Figure 2B, the apparent uncertainty sources are the uncertainties caused by the atmospheric conditions, the scanning geometry, rasterization, and registration. The same accounts for Figure 2A but with a smaller contribution for registration. Nevertheless, the theoretically calculated σ cell in Figure 11 could be a few decimeters larger if the registration uncertainty from Figure 2B is taken into account, but this only applies for large anomalies in pitch measurements relative to a reference scan. Furthermore, the value of 0.008°in Eq. 5 is fairly generously chosen, and the contribution of σ instrument to the total uncertainty budget is probably smaller. The consequences of the different uncertainties obtained on the ability to detect surface changes at the glacier are discussed in Sect. 4.

Detecting phenomena at the glacier surface
The aim of the installation of the TLS is to measure glacier surface changes at HEF. This section discusses a series of different phenomena that involve glacier surface changes with their Frontiers in Earth Science 13 frontiersin.org appurtenant time scales, albeit having different magnitudes and can possibly be detected with a TLS. The ability to measure the annual geodetic mass balance with a TLS has already been proven in other studies (Fischer et al., 2016;Wang et al., 2018;Xu et al., 2019). The geodetic method is comparable to the glaciological method under the assumption that we integrate over all the processes that lead to elevation changes at any single point of the glacier over the total glacier area (Kuhn et al., 1999;Zemp et al., 2010;Klug et al., 2018). The geodetic mass balances derived from the automatically registered TLS data were compared to the glaciological mass balances for the hydrological years (between 1 October and 30 September) 2016/17 to 2019/20 for HEF. Differences of −117 kg·m −2 (mean) and 162 kg·m −2 (rootmean-square error) were found. This is within in the glaciological mass balance uncertainty, which is quantified at ±210 kg·m −2 for HEF (Klug et al., 2018). Thus, reanalysis (Zemp et al., 2013) or better registration is not needed and also shows the applicability of the automatically registered TLS data for yearly mass balance measurements at HEF.
The glacier mass balance can likewise be measured over seasons (Gabbud et al., 2015;Xu et al., 2019;Li et al., 2021), months, weeks, and potentially days with a TLS, given snow density information is available. This is especially interesting during the melt season, where melt evolves quickly, and snow is distributed as patches over the glacier. This can be useful validation data for atmospheric (Goger et al., 2022), snow cover (Mernild et al., 2006), or distributed mass balance models (Klok and Oerlemans, 2002;Hock and Holmgren, 2005;Machguth et al., 2006) and has magnitudes between centimeters and meters. Additionally, with daily measurements, the glacier loss day can be determined, which is the day in the hydrological year when the volume and, with similar snow distribution patterns, the mass is the same as that on 1 October. This means that after this date, the glacier will only lose mass in the present climate conditions.
On shorter time scales, the process with a magnitude of centimeters to meters at the glacier is snowfall. During snowfall events, snow is not evenly distributed over the glacier but is mostly driven by orographic precipitation in winter and by two precipitation processes that are mainly terrain and wind-driven, namely, 1) snowfall enhancement caused by the interaction with the local flow field and local cloud formation processes (Mott et al., 2014) and 2) pure particle flow interaction (preferential deposition of snowfall) . After snowfall, the snow is redistributed over the glacier. These snow drift phenomena have a vertical spatial scale of <1 m (Filhol and Sturm, 2015) and a horizontal spatial scale of 1-100 m (Marsh et al., 2020). Depending on the atmospheric conditions, snow drifts and settles after snowfall (Mott et al., 2018). During the day, snow fall and redistribution are still hard to model, but the TLS data can be a useful tool to validate models that take them into account. Snow fall, redistribution, and snow settling are difficult to distinguish from each other as they happen simultaneously, so these are partly combined in the investigation of snow distribution.

FIGURE 11
Distributed TLS uncertainty with the point density taken into account (σ gridcell ).
Last, special cases of snow redistribution are avalanches. Avalanches can occur within seconds but have magnitudes up to meters in the vertical direction and runout lengths of meters to kilometers. The avalanches are clearly visible in webcam pictures at HEF, yet the DEMs give quantitative information, which is of interest for the snow avalanche community as the input for model calibration (Prokop et al., 2008;Schaffhauser et al., 2008;Hancock et al., 2020).
Our study shows the possibility to investigate the mentioned phenomena. First, on longer time scales, processes with magnitudes in the order of meters can be detected with the automatically registered data.
Second, this automated registration is also helpful to select case studies and to qualify days when glacier surface changes have occurred. Subsequently, the scans taken at these days can be postprocessed with RiSCAN PRO, delivering an uncertainty of the TLS data of ±0.10 m. Moreover, the grid size in this study is set to 1 m, but to achieve the best combination between accuracy and spatial horizontal resolution, the grid size can also be set an order of magnitude larger. A coarser horizontal resolution increases the amount of measurements per grid cell and hence leads to a smaller contribution of the instrument and hardware limitations to the uncertainty budget as these are evened out in the data. Eventually, this leads to a more accurate automatically registered DEM, which is still usable in the generally coarser models, though having a lower horizontal spatial resolution.
Nevertheless, it is known that snow drift mainly occurs at mountain crests to leeward slopes (Mott et al., 2018), and these areas have an insufficient coverage and a rather low point density (1-2 points m −2 ) in the data of HEF. Therefore, future research could benefit from the coverage of the glacier accumulation zone with airborne laser scanning (ALS), which has already been done in the past at HEF at annual and seasonal temporal resolution (Helfricht et al., 2014;Klug et al., 2018). These studies can be improved with more frequent acquisitions, which might be feasible with uncrewed aerial vehicles (UAVs) (Zieher et al., 2019).

Comparison to other studies
The uncertainty sources discussed in this study were separately assessed in previous studies and in different contexts. The scanning geometry has been evaluated for a TLS on short ranges (≤50 m) inside a building (Soudarissanane et al., 2011). The focus in this study was on point cloud quality, when incidence angles are unfavorable and scanning ranges are long. Similar to our study, Soudarissanane et al. (2011) found that larger incidence angles lead to larger noise levels and a higher uncertainty in the position of the laser beam. Likewise, the extent of the footprint for an ALS acquisition was investigated (Schaer et al., 2007). The estimation of the uncertainty of the footprint in our study is based on the study by Schaer et al. (2007) (Eq. 18), even though the scanning geometry in an ALS acquisition is more favorable due its more perpendicular and closer range (<1,000 m) to the surface than the permanent TLS setup. On the other hand, TLS has no positioning errors related to the GPS/INS navigation of the airborne carrier.
The instrument and hardware limitations of TLS systems are also assessed in other studies (Lichti and Licht, 2006;Griebel et al., 2015), but they use different systems. The radial stripes found in our study have only been mentioned by Kuschnerus et al. (2021), but the cause of the stripes was not investigated.
The uncertainty caused by the atmosphere has previously been modeled in an idealized situation by Friedli (2020). This study investigated the velocity difference and refraction of the laser beam due to atmospheric conditions. The differences in the laser beam found by Friedli (2020) range between −0.11 and 0.27 m at a longest simulated distance of 1,841 m and, thus, an order of magnitude smaller than that found in our study. Friedli (2020) also stated that the change strongly depends on the temperature profile in the boundary layer. We assumed a linear temperature profile in the last 10 m of the laser beam path based on the simulated temperature at 2 m and the surface temperature. Different temperature profiles for these last meters of the laser beam path have also been tried, but no significant influences on the laser beam properties were found. Furthermore, our setting is not an idealized situation, and the temperature changes at one point in time are smaller than that in an idealized situation. In our study, the refraction of the laser beam has not been investigated, as the interpolation of the derivative of the refractive index perpendicular to the path dn(s) dq (Eq. 11) is not accurate in a simulation with a horizontal grid spacing of 48 m. The maximum change in the refraction angle found by Friedli (2020) is 6.3 mgon, but this strongly depends on the distance and temperature profile in the boundary layer. The refraction angle change of 6.3 mgon equals 0.0057°, which is approximately the same order of magnitude as the uncertainty added by the instrument and hardware limitations. Under the assumption that the boundary layer temperature profile is weak on a snow-covered glacier, the contribution of the angle of refraction is likely less strong than that found by Friedli (2020) and makes the laser beam hit the surface within the calculated footprint. This also validates neglecting the influence of refraction of the laser beam in our study. Eventually, our study is able to enhance the results of Friedli (2020), as the data of the atmospheric variables are based on a highresolution simulation and variable over time. Therewith, we show the novelty of our study, where the uncertainties as assessed by different studies are combined and applied to an actual permanent long-range TLS system.

Conclusion
The potential of the permanent long-range TLS system at HEF for detecting glaciologically relevant surface elevation changes has been assessed. Five uncertainty sources are analyzed: the registration method, the influence of the instrument and hardware limitations of the TLS, the effect of atmospheric conditions on the laser beam, the scanning geometry, and the uncertainty caused by rasterization. The uncertainty sources are investigated separately, showing a strong dependence on the distance from the TLS to the surface. The registration method is tested with an automated ICP approach (ICP HEF ) and with commonly used RiSCAN PRO software (ICP RiSCAN ). This results in an average vertical accuracy of the TLS data of ±0.50 m with ICP HEF and ±0.10 m with ICP RiSCAN at grids of 100 by 100 m. The precision of the inclination sensors of the scanner cause an uncertainty between 0.09 and 0.64 m. The influence of the atmosphere on the velocity of the laser beam strongly depends on air temperature and to a lesser extent on RH and air pressure, but only has an order of magnitude of millimeters. The contribution of the scanning geometry to the uncertainty budget is approximated between 0.02 and 0.25 m. The instrument and hardware limitations cause the largest uncertainty to the TLS data. The total uncertainty with the point density of the TLS taken into account is between 0.01 m at the area between the TLS and StHE and 0.66 m in the areas close to Weißkugel.
For the annual and seasonal mass balance, the automated approach has proven to be sufficient, as the decimeter-range uncertainty is comparable to the uncertainty of the glaciological mass balance method. Phenomena smaller than a few decimeters, such as snowfall and snow distribution, can be studied after increasing the registration accuracy with post-processing. The uncertainty of the TLS system is an order of magnitude smaller than the processes of interest. Thus, snow drift phenomena with a vertical spatial scale >0.10 m and a horizontal spatial scale of 1-100 m can still be investigated with the data at hand. This also shows the suitability of the system setup to validate high-resolution atmospheric models that compute snow (re)distribution by wind.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.