Decadal Surface Changes and Displacements in Switzerland

Multi-temporal, high-resolution, and homogeneous geospatial datasets acquired by space- and/or airborne sensors provide unprecedented opportunities for the characterization and monitoring of surface changes on very large spatial scales. Here, we demonstrate how an off-the-shelf, open-source image correlation algorithm can be combined with SwissALTI3D LiDAR-derived elevation data from different tracking periods to create country-scale surface displacement and vertical change maps of Switzerland, including Liechtenstein, with minimal computational effort. The results show that glacier displacement and ablation make up the most significant fraction of the detected surface changes in the last two decades. In addition, we identify numerous landslides and other geomorphic features, as well as manmade changes such as construction sites and landfills. All produced maps and data products are available online, free of charge.


Introduction
The increasing availability of high-resolution and homogeneous geospatial remote sensing datasets offers unforeseen possibilities for long-term monitoring of the Earth's surface. These datasets include products derived by orbital and airborne optical, multispectral, thermal, laser scanners, and radar sensors (e.g., Emery & Camps 2017;Olsen 2016). In particular, LiDAR (Light Detection And Ranging) plays an important role in environmental and surface observation, characterization, and monitoring due to the accuracy (on the order of centimeters) and ground sampling distance (on the order of meters or lower) of its products (Tarolli 2014). In Switzerland, the Federal Office of Topography (swisstopo) has initiated a long-term program for the acquisition, processing, and release of airborne LiDAR-based, high-resolution, high-fidelity DTM (Digital Terrain Model) products at the country scale, known as SwissALTI3D (SwissALTI3D 2013 & 2019). SwissALTI3D represents a unique opportunity to study and characterize a variety of small-(meter) and large-scale (kilometer) surface changes that occurred over a period of nearly 20 years (as of today), covering two countries: Switzerland and Liechtenstein.
In this work, we showcase how the SwissALTI3D dataset and derivatives of it can be used to characterize and monitor surface changes on the country-scale using modern, off-the-shelf (opensource) digital image correlation algorithms. We systematically process DTM data from two different tracking periods to identify potential changes in vertical elevation and horizontal displacement. We provide a set of Switzerland-and Liechtenstein-wide maps that can be used as base for further investigations at regional and local scales and briefly discuss how our approach could be adopted to study other regions of the world. In the following sections, we describe the applied data and methods, their characteristics, and limitations, and present results for a few selected areas.

SwissALTI3D
This work is based on the SwissALTI3D releases 2013 and 2019 (R13 and R19, respectively) provided by swisstopo via GeoVITE and the ETH library (SwissALTI3D 2013 & 2019). These releases contain DTM data from two different tracking periods. The R13 contains data acquired over the years 2000-2012, while R19 data are over the years 2012-2018. The minimum time difference of tile pairs between the two releases is 3 years and the maximum time difference is 17 years (SwissALTI3D 2013(SwissALTI3D & 2019. R13 contains DTM data that are either based on DTM-AV (airborne) LiDAR, a combination of LiDAR and stereophotogrammetry, and on stereophotogrammetry alone (for elevations > 2000 m). R19 contains DTM data that are based on "new" or "old" (DTM-AV) generation LiDAR in combination with 3D stereo measurements, and on a combination of stereophotogrammetry and 3D stereo measurements (for elevations > 2000 m) (SwissALTI3D 2013 & 2019). The used 2019 elevation data is showcased in Fig. 1 and an overview of the temporal difference distribution is displayed in Fig. 2.
The spatial resolution for all products used in this study is 2 m/pixel, with vertical accuracies of ± 0.5 m for the DTM-AV, ± 0.3 m for the new generation LiDAR, and between ± 0.1 and 3 m for the stereo-photogrammetric data (SwissALTI3D 2013(SwissALTI3D & 2019. For R13, parts of the stereo-photogrammetric data were acquired over an extended period from 2008 to 2011. Differences in snow and ice cover during this period can cause local artefacts in SwissALTI3D, such as linear breaks (edges), vertical offsets, or seams along the edges of snowfields. A similar problem exists when merging stereo-photogrammetric data and LiDAR data, as the respective acquisition flights are usually performed at different points in time. These artefacts can affect the image correlation routines and could result in noise in the final displacement velocity and change maps, especially at higher topographic elevations (> 2000 m).

DTM Patching and Pre-processing
We divided the SwissALTI3D data into 248 individual subpatches per release, i.e., to a total of 496 patches (248 pairs),  Fig. 1). Every patch has been stored as a 32 bit (floating point).geotiff with dimensions of 8750 × 6000 pixel.
The DTM patch pairs from the 2013 and 2019 releases were imported into Matlab R2017a for further processing. We customized and implemented the "hillshade.m" function by Hebeler (2022) (https:// de. mathw orks. com/ matla bcent ral/ filee xchan ge/ 14863-hills hade) to assess whether hillshades or DTMs are more suitable for characterizing and monitoring surface displacements (the input parameters for the hillshade function were: solar azimuth of 315° and solar elevation of 75°) (see, e.g., Fey et al. 2015).

Digital Image Correlation (DIC)
Digital image correlation (DIC) is a technique that correlates multi-temporal image/sensor acquisitions of a given region, identifying displacements that occurred over time, i.e., in between individual acquisitions. DIC is a powerful remote sensing tool that can be applied to a large variety of input data (e.g., optical, thermal, radar images; Dematteis & Giordan 2021;Cai et al. 2017;Walter et al. 2013) and for a variety of applications (landslide, volcano, material testing monitoring; e.g., Ayoub et al. 2009;Bickel et al. 2018).
We use an established, validated, and open-source DIC code that uses a Fourier transform-based algorithm (DIC_FFT, Bickel et al. 2018). The accuracy of that particular DIC algorithm has been tested using GNSS displacement data and ranges from 8 to 12% of the observed maximum surface displacement, outperforming all other tested DIC algorithms (Bickel et al. 2018). The performance of the algorithm has been further confirmed by independent studies (Dematteis & Giordan 2021).
DIC is generally not able to provide results at locations where substantial or extremely rapid changes in surface texture occur between individual acquisitions (i.e., a catastrophic loss of signal). However, while not providing quantitative information about surface displacement at those locations, the decorrelation indicates that a substantial change occurred, i.e., represents a qualitative measure for surface change. Similarly, DIC struggles to correlate extremely smooth surfaces with an overall lack of texture, such as lakes, plain bedrock, and glaciers/snow/ice; such regions might suffer from decorrelation and noise as well.
For this study, we separately ran DIC using (A) the original SwissALTI3D DTM data and (B) hillshades generated from those, using the input lighting parameters specified above. In the process, we observed that using hillshade pairs provided significantly better results, i.e., better correlation and spatial coverage of the displacement field than the application of DIC_ FFT to DTMs. Thus, we focus on hillshade-derived results in this study.
We performed a DTM and hillshade pair co-registration without image split and oversampling prior to DIC using the FFT-based image alignment function from Bickel et al. (2018) to correct for potential residual rigid shifts between pairs. Previous studies have shown that median and Wallis filtering (Bickel et al. 2018) prior to DIC can help to increase the signal and spatial coverage of the derived displacement field. However, in this work, we discarded the use of such pre-processing steps as the impact on the results was negligible while the computational cost increased significantly (by a factor ~ 2 to 8). The used DIC window size (wi) was 128 × 128 pixels, as it offered the best compromise of spatial coverage, spatial resolution, processing speed, and correlation quality (the used metric is the RMSE, root mean square error). We applied no correlator oversampling. We used a search window x-and y-direction stride (or skip)-that characterizes the shift of the DIC correlation window (or template/kernel) within the image during the correlation-of wi/2 (Bickel et al. 2018). Wi/2 represents the smallest stride that still returns a true, physical representation of the displacement field (wi/2 is the Nyquist frequency). No post-processing of the results has been performed in order to maintain the full interpretability of the DIC results, with one exemption: we applied a 3 × 3 pixel median filter to one of the georeferenced preview products (hillshadebased, wi 128 × 128 pixel) to improve its visual representation. Filtering can be applied by the end-user according to any individual, specific applications, and requirements.
After the correlation, the original.geotiff DTM data has been used to calculate the vertical difference between R13 and R19 for each correlation cell to derive the absolute vertical change that occurred between the different tracking periods. In order to derive horizontal displacement and vertical change velocities, we normalized the derived 2D and elevation changes (in meters) using the absolute acquisition time difference (in years) of the respective tile pairs, as reported by the SwissALTI3D release reports (2013 & 2019). The resulting maps report the mean 2D displacement and elevation change velocities in meters per year referred to the overall period of observation, i.e., assuming continuous displacement.
All products are stored and provided as .geotiff (preview) and .csv (science) files, enabling a detailed analysis in GIS software. Figure 3 visualizes the complete DIC processing pipeline and the outputs; all products are described in more detail in the Data chapter below. All used algorithms and produced datasets are available online, free of charge, in ETH Zurich's Research Collection (https:// doi. org/ 10. 3929/ ethz-b-00055 7355) and on swisstopo's map.geo.admin.ch webpage (https:// map. geo. admin. ch/) via 'Advanced tools/ Import' using the following link: https:// ogc. glamos. ch/ dsc-d.

Results and Discussion
We produce a total of 7 Switzerland-and Liechtensteinwide maps: (1) E-W horizontal surface displacement, (2) N-S horizontal surface displacement, (3) local vertical difference (elevation difference), (4) horizontal surface displacement magnitude (2D resultant), (5) average horizontal surface displacement velocity, (6) average vertical elevation change velocity, and (7) correlation RMSE. Figures 4 and 5 present some of the map products for the Aletsch region (swisstopo tile number 1269, as indicated in Fig. 1), while other map products are showcased on the maps associated with this paper. The technical details of all maps and related files are outlined further below. The processing time for the entire dataset amounts to ~ 2.5 h using a regular desktop PC (4 CPUs and 32 GB RAM). This highlights that fully automated, DIC-based change detection and monitoring on country scales can be time-and cost-effective.
We observe maximum horizontal and vertical displacement/change magnitudes up to ~ 180 and ± ~ 100 m in the Aletsch, Fiescher, and Rhone glacier regions, respectively. Generally, glaciers make up the dominant fraction of displaced features in Switzerland and Liechtenstein. The maximum mean displacement velocities are observed at the Aletsch, Gorner, and Fiescher glaciers with values of up to ~ 40 m/year. These values are generally lower compared with the velocities retrieved by recent studies in the Aletsch region when considering satellite radar products in the period 2011-2019 (e.g., Leinss et al. 2021). Disparities between the results can be due to differences in the time span considered, noise, as well as in the data sources and the methods applied. We note that a number of slopes in regions such as Ticino, Valais, and Grisons appear to feature a faint yet continuous and wide-spread horizontal displacement component, which is likely related to artifacts in the underlying DTMs and the sensitivity of the algorithm.
We observe the largest negative local elevation difference (up to ~ − 100 m) in mountainous glaciated regions, which appears to be predominantly related to ice/glacier ablation and retreat. In turn, we also observe regions that appear Fig. 3 Visualization of the DIC input, processing workflow, and output. All outputs (hillshade-based outputs only) are further described in the "Results and Discussion" section and Data chapters. DTM credits to swisstopo to experience positive local elevation differences, such as mountainous regions that feature seasonal ice and snow accumulation. Locally, we observe small-scale elevation differences that are not directly linked to physical features or processes-those might be related to (random) noise remaining in the SwissALTI3D products, e.g., related to the presence and/or growth of vegetation. We additionally observe significant changes associated with artificial reservoirs and natural lakes, such as the Sihlsee, Lago del Sambuco, and Lago di Luxxone. Those signatures might represent water level differences or might be caused by systematic vertical offsets in the used DTMs, a potential by-product of the swisstopo measurement campaigns.
Our products directly or indirectly reveal a series of known past or active surface displacements and mass wasting events. For example, we observe displacement and decorrelation artefacts that are caused by the Piz Cengalo rock avalanche events that occurred between 2009 and 2015 in Val Bondasca (Walter et al., 2020), as well as an apparent downslope displacement of the deposited scree in the direct vicinity of the Plan Lo and Laret sections of Val Bondasca (Canton Graubünden). Furthermore, we identify a signature of the large Moosfluh instability in the direct vicinity of the Aletsch glacier (Figs. 4 and 5). The mean surface velocities derived for the Moosfluh landslide align well with the values derived via in situ as well as remote sensing observations . We also identify a displacement of the toe of the currently highly active Brienz/Brinzauls landslide (Canton Graubünden, Kenner et al. 2022) that shows a displacement of up to ~ 14 m between 2009 and 2015. We note that there might be an unknown number of currently active, yet unknown slope displacements in the Switzerland-and Liechtenstein-wide displacement dataset we produced. We further detect the DIC decorrelation signatures of a number of large construction and landfill sites, such as the General Guisan Quai construction site in Zurich, the ICTR Baragge di Sopra construction site in Bellinzona, and the In Divan landfill site in Sigirino (Fig. 6).
The correlation RMSE values range from 0 to ~ 0.4 m for the hillshade-based products, indicating a reasonably low overall error. We identify larger RMSE values exclusively in regions that experienced substantial, rapid change (loss of correlation, e.g., due to construction or catastrophic slope failure), and in mountainous regions, usually above 2000 m, usually related to snow and ice-related DTM and hillshade artefacts (surface texture changes related to ice and snow deposition or ablation, as described earlier).

Fig. 4
Detailed information for swisstopo tile number 1269 (Aletsch glacier region). A RGB aerial image, B DIC_FFT correlation RMSE, C horizontal surface displacement magnitude (pixels with values below 2 m are omitted as they are either stable or their displacement is below the sensitivity of the algorithm -DIC_FFT can resolve displacements of at least one pixel, i.e., 2 m; RGB aerial image in the background, and D vertical elevation difference (negative is downwards-loss, positive is upwards-gain). Note that in panel C white regions (e.g., in the center of the Aletsch glacier) are not decorrelated but were mismatched due to the smoothness of the glacier's surface, erroneously indicating "stable" ground. Panel C represents a raster product; panels B and D represent vector products; all were derived from a hillshade pair using a wi of 128 × 128 pixels and a wi/2 skip, i.e., feature a spatial resolution of 128 m/pixel. The white box in panel A denoted with an * denotes the Moosfluh landslide area that is displayed in greater detail in Fig. 5. Grid facilitates cross-comparison of features among the panels. Image and DTM credits to swisstopo We note that our DIC-driven approach (DIC_FFT; Bickel et al. 2018) can be deployed in any location -world-wide -that features a consistent and multi-temporal geospatial remote sensing dataset derived by sensors such as LiDARs, RaDARs, and optical/multispectral cameras. Here, the scale (space and time) of discernable surface change depends on some of the key characteristics of the available dataset, such as the acquisition mode, spatial resolution, and temporal baseline.

Produced Data
All generated outputs are made accessible as two product types (raster .geotiff and vector .csv). There are three products in total: (1). Hillshade-based vector (128 m/pixel, no resampling), floating point ( Product (1) consists of a single file that contains all different displacement and change maps in uncompressed, tabular form, where each column represents one displacement and change map-this product is supposed to be used for scientific analyses. (2) and (3) consist of multiple georeferenced rasters with compressed (and binned values)-these products are supposed to be only used as a thumbnail/preview product. Currently, we only provide one type of raster (the horizontal surface displacement magnitude in meters, rounded), but might add other displacement and change maps in the future. It is important to note that all displacement values contained in those rasters have been rounded to full integers before they were written into 8bit files, i.e., they do not represent fully accurate displacement valuesthese are preview/quick-view files only.
Product (1) contains columns that represent the following parameters (column name/header is shown in brackets):

Fig. 5
The Moosfluh instability (outline approximated by dashed line), located right next to the current Aletsch glacier terminus (RGB aerial image) (A). The Moosfluh instability is characterized by distinct vertical elevation differences (B) and horizontal displacements (C). Panels B and C highlight how the surroundings of the glacier are modified due to the retreat of the ice (see, e.g., Storni et al. 2020). The Moosfluh displacement field shows compartments with different displacement and vertical change magnitudes. Panels B and C are derived using hillshades, a wi of 128 × 128 pixels, and a wi/2 skip, i.e., feature a spatial resolution of 128 m/pixel. Grid facilitates crosscomparison of features among the panels. Image and DTM credits to swisstopo <12> "vertical elevation change velocity (m/a)" (VELOZ).
Product (1) features a precision of two decimals, i.e., values are rounded to two decimal places; columns < 10 > and < 12 > are not rounded to preserve very small displacement velocities. Products (2) and (3) feature a precision of zero, i.e., are rounded to zero decimal places (full integers, preview files).

Conclusions
We used a validated and open-source digital image correlation algorithm (DIC_FFT) in combination with a countryscale, high-resolution & high-fidelity, multi-temporal DTM dataset (SwissALTI3D) to create-for the first time-vertical change and horizontal surface displacement maps of Switzerland and Liechtenstein. The produced maps characterize and quantify meter-scale surface change that occurred between the years 2000 and 2018, containing glaciers, landslides, construction sites, landfills, mining sites, and other geomorphic and manmade changes. This work demonstrates the cost-and time-effective potential of DIC for the identification, quantification, and monitoring of surface change on country-type spatial scales. The analysis of the produced displacement maps is ongoing-we note that an unknown number of small-and large-scale changes and surface displacements may yet to be discovered. All used algorithms and produced datasets are available online, free of charge, in ETH Zurich's Research Collection (https:// doi. org/ 10. 3929/ ethz-b-00055 7355) and on swisstopo's map.geo.admin.ch webpage (https:// map. geo. admin. ch/). On map.geo.admin.ch, please add the data via 'Advanced tools/Import' using the following link: https:// ogc. glamos. ch/ dsc-d.