3.1 Generation of the topographic data

3.1.1 Processing of airborne lidar and photogrammetric SfM–MVS point clouds

Several data sets were used for the reconstruction of the terrain surface for the entire catchment area. These include both current airborne lidar data and historical aerial image series (Fig. 2). Thus, the time periods are based on the availability and quality of the data.

Figure 2Type of topographic data set and the resulting time periods.

Download

To determine the recent morphodynamics at the respective sites, available airborne lidar data from 2004/2005 to 2019 were used. The 2004/2005 and 2006 data of the three valleys were provided by the Autonomous Province of Bolzano and the province of Tyrol (Table 3). The latest airborne laser scanning (ALS) data sets of each valley (2017 and 2019) were collected in the ALS flight campaigns of the Chair of Physical Geography at the Catholic University of Eichstätt-Ingolstadt (Table 3) (Stark et al., 2022). In this case, lidar data sets were collected using previously determined flight strips. Direct georeferencing (position and altitude) of the trajectories was determined by Global Navigation Satellite System (GNSS) rover antenna and an inertial measurement unit (IMU) (Applanix AP 20), both located in the laser scanner. In addition, GNSS correction data were acquired on the ground during the flight missions using a dGNSS antenna (Fig. 3).

Table 3ALS and DEM data and corresponding flight mission attributes. NA: not available.

Download Print Version | Download XLSX

Figure 3ALS data collection on 8 August 2019 in Horlachtal. Helicopter with nose-mounted VP1 laser scanner, as well as the ground station, which recorded the dGNSS raw data during the flight time (Stonex S9III).

Download

In order to extend the temporal scope of this study by several decades (until 1953), previously digitised (high-resolution) overlapping historical aerial images were processed into historical DEMs. Except for the 1959 Martelltal survey, camera distortion parameters and focal lengths were provided for all data with the respective camera calibration certificates (Table 4).

The digitised image series were processed with the Agisoft Metashape Professional Edition software package (Version 1.6.6; Agisoft), using structure-from-motion (SfM) photogrammetry with multi-view stereo (MVS) algorithms to generate high-resolution point clouds. The generation of point clouds from digitised historical (aerial) image series requires different preparation and processing steps. First, all images of each series were resized to a common image size (uniform number of pixels along the x and y axes) without changing the image content. This step was necessary so that the software can assign all images to the same camera (source) and was carried out using Adobe Photoshop (CS6). This is of enormous importance in order to be able to use the appropriate distortion parameters for the respective camera models for the calculation. Thereafter, the image sets were imported into single folders, and a common global coordinate system (ETRS89 UTM zone 32N; EPSG code 25832) was defined. Next, all images were masked to exclude the black edges or the frame (instrument stripes with the camera metadata) in order to avoid interference with the orientation of the cameras (Gomez et al., 2015). Before the initial processing of images, we defined the fiducial mark information and lens distortion parameters in order to set the metric dimension of images and lenses. This information was included and used for the alignment of single images (SfM).

Since a global exterior orientation requires a large number of precisely surveyed ground control points (GCPs) distributed throughout the area, we used highly precise ALS data sets with millimetre accuracy (2019 Horlachtal; 2017 Kaunertal) to extract these GCPs and to define the exterior orientation of all data. The selection and extraction of GCPs was based on clearly identifiable objects (e.g. rock formations) that were also considered to be stable (geomorphologically unchanged) over the entire observation period. If a calibration certificate was available, the film camera option was used, fiducial marks defined, and the focal length set and fixed. All other lens distortion parameters (C_x, C_y, k₁, k₁, k₁, p₁, and p₂) were estimated and adjusted fully automatic using the auto-calibration function. In case of missing camera calibration certificate, an auto-calibration (no film camera) was performed. Both options were proposed by Stark et al. (2022).

According to these pre-processing steps, the point clouds were generated by the (i) initial joint orientation of the images; (ii) selection of ground control points (GCPs); (iii) final camera orientation (bundle block adjustment) including scale definition; and (iv) calculation of dense point clouds.

The processing of the 1959 point cloud, which was used in this study, is already described in Betz et al. (2019).

Table 4Overview of acquired historical image series for point cloud generation and corresponding DEMs by photogrammetry SfM.

^aRefers to the exact sites. ^bRefers to the entire data set.

Download Print Version | Download XLSX

3.1.3 Data statistics

The presence of various uncertainties in differently generated DEMs (Hodgson and Bresnahan, 2004; Bakker and Lane, 2017) also leads to uncertainties in the resulting DoDs (Lane et al., 2003; Rolstad et al., 2009; Cavalli et al., 2017; Anderson, 2019). Therefore, an uncertainty assessment was carried out using the DoD values from stable areas near each site. The size of the stable areas varied between 25 % and 75 % of the size of the corresponding site. In addition to the estimation of the precision (SD) and accuracy (RMSE), the arithmetic mean, minimum, and maximum values were also determined (Fig. 4).

Figure 4Statistical analysis of the uncertainty assessment. (a) Size of the test sites, (b) size of the stable areas next to the test sites, (c) max values (stable areas), (d) min values (stable areas), (e) standard deviation (stable areas), (f) RMSE (stable areas), and (g) arithmetic mean.

Download

To determine the uncertainty in the sediment volume change (total sediment output; Fig. 7), the error propagation method for uncorrelated, correlated, and systematic error according to Anderson (2019) was applied. No threshold has been set for the level of detection of the DoDs, as Anderson (2019) clearly recommends not using this for volumetric calculations as it leads to bias in the results. For the final determination of the total error, the following formula was applied Eq. (1):

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_v=\sqrt{{\mathit{\sigma }}_{v,\mathrm{re}}^{\mathrm{2}}+{\mathit{\sigma }}_{v,\mathrm{sc}}^{\mathrm{2}}+{\mathit{\sigma }}_{v,\mathrm{sys}}^{\mathrm{2}}},\end{array}$$

where σ_v, re is the uncorrelated error, σ_v, sc spatially correlated error, and σ_v, sys systematic error.

3.2 Derivation of the regression lines

In this study, we followed the sediment-contributing-area approach of Neugirg et al. (2015a, b, 2016) and Dusik et al. (2019), who applied this approach at the slope scale and replaced real sediment traps in the channels, as originally based on the work of Haas (2008) and Haas et al. (2011), using the sediment-contributing-area model, with so-called virtual sediment traps in modelled channels in a DEM (Fig. 5).

The sediment-contributing-area model uses an empirical relationship between the log. sediment-contributing area as the independent variable and the log. mean annual bedload sediment yield as the dependent variable. Thus, the sediment-contributing area can be used as a predictor of sediment delivery in alpine catchments.

Linear regression analysis was used to show this significant correlation, which is formulated as Eq. (2):

where y is the log. mean annual bedload sediment yield, and x is the log. sediment-contributing area.

This has already been confirmed in several studies in both small and large catchments (ranging from hectare to square kilometres) and in different regions such as the Northern Calcareous Alps (Haas, 2008; Haas et al., 2011; Sass et al., 2012; Huber et al., 2015) and the French Northern Alps/Prealps (Altmann et al., 2021). Finally, it can be stated that a linear dependency of two variables x and y on a log–log scale has a fundamentally different behaviour than a usual linear dependency. In our case, we have y=log (sediment yield) and x=log (sediment contributing area). Back-transformation of Eq. (2) using the exp function yields gives the following relation between sediment-contributing area and sediment yield, as seen in Eq. (3):

Thus, the relation between sediment yield and sediment-contributing area is a polynomial of the form $y=a\cdot x^b$. In particular, the slope in the log–log model represents the exponent of the polynomial in the standard model. The relation between sediment yield and sediment-contributing area is (nearly) linear if slope is (close to) one. In this case, the exponential of the intercept in the log–log model represents the slope of the linear relation in the standard model, meaning that independent of the actual size of the sediment-contributing area, 1 m² provides the same amount of sediment yield, given by exp (intercept). On the other hand, if the slope in the log–log model is considerably greater than one, then the standard model shows polynomial behaviour, meaning that in the same site, increasing the sediment-contributing area provides more sediment yield per square metre. The steps of the sediment-contributing-area approach of this study are composed as follows and were implemented in SAGA LIS and R. The elevation changes in DoDs (using no threshold) generated from multitemporal data were routed downslope and accumulated using the D8 algorithm (O'Callaghan and Mark, 1984). The resulting accumulated DoD values (accDoD) in every raster cell correspond to the net volume of the sediment balance within its contributing area. On steep slopes, accDoD will be negative and represents the sediment yield of this contributing area (Pelletier and Orem, 2014); if it is close to zero, it means that all eroded sediment has been re-deposited within the contributing area. As in the previous sediment-contributing-area studies by Neugirg et al. (2015a, b, 2016), the application of the parameters used in the original sediment-contributing-area model (Haas, 2008; Haas et al., 2011), which lead to the reduction in the hydrological catchment to the sediment-contributing area, is omitted because the sites and the modelled channels are consistently steep, uncovered, and have short slope lengths, making this reduction obsolete. Therefore, the sediment-contributing area is identical to the catchment area in this study.

In detail, channel initiation points were delineated using a threshold of 20 m² of the flow accumulation that was computed using the D8 algorithm (O'Callaghan and Mark, 1984). Channels that were shorter than 10 m were discarded. To ensure statistical independence through avoiding overlapping contributing areas, a stratified sampling scheme was adopted that included one randomly selected raster cell per channel. Pairs of values (sediment yield and the corresponding sediment-contributing-area size) were randomly extracted from the corresponding channels for each site, and a regression line was calculated accordingly. To quantify the uncertainty due to random selection, this sample was repeated 100 times, resulting in 100 regression models of sediment yield on the sediment-contributing area.

Furthermore, we added two conditions and further developed the sediment-contributing-area approach accordingly. In order to obtain more stable regression lines, the range of values of the sediment-contributing-area size was divided into quartiles (with an equal number of cells within the quartiles) to ensure a homogeneous distribution of the extracted values. Additionally, samples that contained points with a high leverage (greater than 0.5) in the regression model were discarded, and the sampling was repeated until a number of 100 samples was reached.

Figure 5Derivation of the sediment-contributing-area approach using the example of test site KG2. (a) Determination of the size of the catchment area, (b) determination of sediment yield, (c) modelling of the channels and exemplary placement of the virtual sediment traps, and (d) calculation of the regression lines.

3.3 Calculation of the sediment output

Additionally, the total sediment output volume, the mean annual sediment output (divided by the corresponding number of years) and the specific mean annual sediment output (additionally divided by the area of the site) were calculated for each site and time period.

The following equation was used for this (Eq. 4):

where Σ DoD is the sum of the corresponding DoD subset values, and L² is the cell size.

3.4 Generation of meteorological data

Using data generated with a regional climate model (RCM), the influence of the changes in climate forcing (air temperature and precipitation) on morphodynamics was investigated. For dynamical downscaling of climate data for the beginning of the study period until 2015, we used the Advanced Research WRF (ARW) model of the Weather Research and Forecasting (WRF) model (version 4.3), which is based on fully compressible and non-hydrostatic equations (Skamarock and Klemp, 2008). The Twentieth Century Reanalysis version 3 (20CRv3) data set (Compo et al., 2011; Giese et al., 2016; Slivinski et al., 2019), with a spatial and temporal resolution of 1°×1° and 3 h, respectively, was used as driving data (initial and boundary conditions). The simulation was performed in three nested domains, with grid spacing of 18 (Domain 1), 6 (Domain 2), and 2 km (Domain 3). For our simulations, we mainly used the physics and dynamics options proposed by Collier and Mölg (2020) and listed in Table 5. However, the Noah land surface model, prescribed eta levels by Collier et al. (2019), and the 24 United States Geological Survey (USGS) land use categories were used. The temporal resolution of simulated data in D3 is 1 h for temperature and 15 min for precipitation.

Table 5Overview of the Weather Research and Forecasting (WRF) Model configuration.

Download Print Version | Download XLSX

The simulated temperature and precipitation data were extracted at the location of each of the five glacier forelands (Fig. 1). These are the centres of the respective sites and represent the corresponding glacier foreland. In addition, a corresponding elevation correction of the climate data was applied for temperature.

For the analysis, we used the mean annual air temperature (2 m above ground), as well as the corresponding trends and the mean number of ice days (days with maximum temperature <0 °C). In addition, the number of warm air inflows from October to May was determined in order to identify corresponding snowmelt processes at the sites. A warm air inflow is defined as a period of at least 3 d in which more than 70 % is above 0 °C, following a previously colder period of 5 d (100 % below 0 °C). In addition, the precipitation patterns were analysed. For this purpose, the mean annual precipitation totals, the mean annual winter (October to May), and the mean annual summer precipitation totals (June to September), as well as the corresponding trends, were determined in order to identify seasonal changes. Furthermore, various continuing classes (4 mm for 1 h resolution and 10 mm classes for daily totals) were used to analyse corresponding changes in individual extreme events and daily precipitation totals. Individual precipitation events were defined as one event, regardless of length, if they were contiguous throughout and were separated if there was no precipitation for at least 1 h. To minimise the noise generated in the data, both data sets were also filtered for extremely small events by changing the values from <0.01 to 0 mm. The calculation of the mean annual winter precipitation was always carried out over the entire winter. For example, for the winter of 1953, data from October 1952 to May 1953 was included. The average summer precipitation was calculated accordingly from June 1953 to September 1953. Furthermore, precipitation was differentiated into snow and rain events. The determination of a threshold to distinguish rain from snowfall is very dynamic in mountainous regions and difficult to estimate. However, the difference between rain to snow depends mainly on surface air temperature, as well as air humidity, with snow occurring mainly between 0 and 3 °C (Froidurot et al., 2014), and the lower the humidity, the higher the probability of snowfall will be. In this study, the threshold from rain to snow was defined at ≤0 °C, as rain is almost excluded below this temperature (Froidurot et al., 2014; Fehlmann et al., 2018).

4.1 Sediment-contributing-area approach

All determined regression lines show a positive correlation between log mean annual sediment yield and log sediment-contributing area (Fig. 6; Appendix A), which means that sediment yield increases with the corresponding sediment-contributing area. In the following, only the median of the 100 regression lines (median slope) of all sites and time periods is used to qualitatively describe corresponding differences. Mostly, there is a decrease in sediment yield between the different time periods and a decrease in the slopes of the regression lines, which is due to a decrease in sediment yield per square metre of the sites. With regard to Sect. 3.3, a decreasing intercept together with an almost constant, although slightly decreasing, slope close to one can be seen over the different time periods in the log–log model, indicating that the relation between sediment-contributing area and sediment yield remains almost constant. The sites KG3, KG4, KM1, KM2, and KW2 show such a behaviour. On the other hand, the areas KG1, KG2, KG5, KW1, and MH2 show clearly larger differences in the time periods. In the earliest time period, the slopes considerably larger than one (in the log–log model) show polynomial behaviour, which means that in the same site an increasing sediment-contributing area provides clearly more sediment yield per square metre. In the later time periods, the slopes also tend towards one, so that the models of the different groups become similar. In addition to this general trend (10 sites), an increase in sediment yield and an increase in the slope of the regression line for site HG1 were observed, showing an increase in sediment dynamics over the time periods in this case, which is in contrast to the previous observations. Site MH1 shows a similar level of sediment yield (between the time periods) with higher slopes of the regression lines, also indicating an increase in sediment yield. Furthermore, slopes of the regression lines below 1 occur in all time periods but especially in the second and third.

Figure 6Relationships between log sediment-contributing area and log sediment yield for 100 samples of each site and the corresponding time periods. In addition, the median regression line (median slope) is represented by a slightly thicker and darker line.

Download

4.2 Volume calculations of the sediment output

The analyses of the sediment output of the sites confirm the results of the sediment-contributing-area approach (Fig. 7). In general, there is a clear and continuous decrease in mean annual sediment yield of 10 sites over the different time periods. In contrast to this trend, site HG1 shows a clear increase in mean annual sediment yield. Site MH1 also shows an increase but at a very low level, which can also be described as a geomorphic activity of a similar level. In total, the mean annual sediment output decreases across the different time periods. Nevertheless, there is also very high temporal and spatial variability in this change at the sites HG1 and MH1, which also shows a clear increase in geomorphic activity, as well as a slight increase (respectively activity at the same level).

Figure 7Bar plots of total sediment output (with error range according to Anderson, 2019; see Sect. 3.1.3 and 3.4), mean annual sediment output, and mean annual specific sediment output (erosion rate) of each site and time period.

Download

4.3 Meteorological regime

4.3.1 Air temperature

The mean annual air temperature (2 m above ground) of all selected positions of the glacier forelands shows a statistically significant warming trend over the entire study period of 60, 64, and 65 years (Fig. 8). Overall, there is a positive total change of + 1.75 °C (annual trend 0.03; p value <0.05; R² 0.39) for the Horlachtal/Grastalferner glacier forelands; + 1.68 °C (annual trend 0.03; p value <0.05; R² 0.38) for the Kaunertal/Gepatschferner glacier foreland; + 1.70 °C (annual trend 0.03; p value <0.05; R² 0.38) for the Kaunertal/Gepatschferner Münchner Abfahrt (MA) glacier foreland; for the Kaunertal Weißseeferner glacier foreland of + 1.64 °C (annual trend 0.03; p value <0.05; R² 0.36); and for the Martelltal Hohenferner glacier foreland of +2.23 °C (annual trend 0.04; p value <0.05; R² 0.45).

Figure 8Mean annual 2 m air temperature of the glacier forelands within the study periods (95 % confidence interval is included).

Download

The analysis of the mean annual ice days shows a decrease between the time periods, especially from the second to the third time period (in the Martelltal from the first to the second), with a decrease in ice days between 18.4 and 20.7 d, which corresponds to almost 3 weeks (Fig. 9). The analysis of the mean annual warm air inflows shows a decrease from the first to the second time period in the glacier forelands of Grastalferner, Gepatschferner, and Gepatschferner MA and a more pronounced increase from the second to the third time period. In the glacier foreland of Weißseeferner, there is a consistent increase, with the latter being equally more pronounced, whereas in Martelltal there is only a slight increase. Thus, the number of warm air inflows has generally increased.

Figure 9Mean annual ice days and mean annual warm weather inflows with the corresponding changes between time periods.

Download

4.3.2 Precipitation

Over the entire time period (60, 64, and 65 years), all study areas show a decreasing trend in mean annual, mean summer, and mean winter precipitation (with the exception of summer precipitation in the Martelltal, which shows a positive trend) (Fig. 10). However, only the changes in winter precipitation (entire study period), summer precipitation (second time period 2004/2005–2019) in the Martelltal (Hohenferner glacier foreland), and annual precipitation (third time period 2006–2019) in the Horlachtal (Grastalferner glacier foreland) are statistically significant, although the latter both cover only 13 and 14 years. In the Horlachtal, the first two time periods show a decreasing trend in precipitation, while the third time period shows an increase in precipitation, which is significantly more pronounced in summer than in winter. In the Kaunertal, winter precipitation shows a slight increase in the first time period and a stronger decrease in summer precipitation. The second time period shows a slight increase in summer and a slight decrease in winter. The third time period also shows a strong decrease in summer and winter precipitation. In the Martelltal, on the other hand, winter precipitation decreases significantly, and summer precipitation increases significantly, especially in the second time period, although time periods of different lengths are analysed.

Figure 10Mean annual, mean summer, and mean winter precipitation of the respective glacier forelands.

Download

In the following, the changes in different precipitation classes (as well as with a different temporal resolution) between the individual time periods are analysed. The calculated changes are based on Appendix B and C. In these tables, the number of all precipitation events of the corresponding time periods is shown and divided into corresponding precipitation classes. Using the number of years per time period, this results in an mean annual number of events per class. The calculated changes result from the comparison of the mean occurrence of the precipitation classes of the previous time period. Both precipitation events with a resolution of 1 h (Fig. 11 and Appendix B) and daily precipitation totals (Fig. 12 and Appendix C) were analysed. The highest temporal resolution (1 h) shows that the classes >0 to 4 and 4 to 8 are subject to the highest variations (Fig. 11); for example, precipitation events of the class >0 to 4 occur 7.63 fewer times in the glacier foreland of the Grastalferner in the second time period compared to the first. In general, it can be seen that the higher precipitation classes tend to decrease, albeit very slightly, but there are still changes with both an increase and a decrease in the different precipitation classes. The daily precipitation totals also show a high variation, with both a decrease and an increase over the different time periods (Fig. 12). In general, there are also very slight changes. Nevertheless, the decrease in the higher three classes predominates when comparing the third with the second time period.

Figure 11Change in precipitation classes between the different time periods with a 1 h resolution (extreme precipitation events).

Download

Figure 12Change in precipitation classes between the different time periods with a 24 h resolution (daily precipitation totals).

Download

The analysis of the mean annual rainfall and snowfall events shows that there is a consistent increase in rainfall and at the same time a consistent decrease in snowfall (except for the Weißseeferner glacier foreland) (Fig. 13).

Figure 13Mean annual rain and snowfall events with the corresponding changes between time periods.

Download

5.1 Assessment of the sediment-contributing-area approach

Using the relationship between accumulated sediment yield from DoDs and sediment-contributing area/catchment area (log–log model) of different sites and different time periods, we show a long-term monitoring of several geomorphologically active sections of the LIA lateral moraines. This is a clear difference compared to previous studies that used the space-for-time substitution (SFTS) approach in proglacial areas, in which studies used recent morphometric or morphodynamic differences between sites located along a gradient of deglaciation age to infer long-term changes in morphodynamics (Ballantyne and Benn, 1994; Curry, 1999; Curry et al., 2006). This means that long-term studies with quantitative data are rare (Schiefer and Gilbert, 2007; Betz et al., 2019; Altmann et al., 2020; Betz-Nutz, 2021; Betz-Nutz et al., 2023). The approach shown here provides reliable results and requires only a few input data (Neugirg et al., 2015a, b; 2016; Dusik, 2019; Dusik et al., 2019). The results mainly show a decrease in sediment yield, as well as a decrease in the slope of the regression lines (suggesting less sediment yield per square metre) over the different time periods, indicating a decrease in geomorphic activity at these sites. In some sites, we observe contrasting changes. There is an increase in sediment yield and an increase in the slope of the regression line at HG1 and almost no change in sediment yield on the y axis but an increase in the slope of the regression line at MH1, which can be described as an increase (HG1) and a constant geomorphic activity (MH1). Moreover, in the earlier time periods, a clearly higher variability in the sediment yield (slope of the regression line clearly higher than 1) was observed at the respective sites, which is no longer reached in the later ones. Thus, it is possible to describe two different types of change in sediment yield (size of sediment yield between time periods and variability in the sediment yield within a time period, which can also be compared with the other time periods). Slopes of the regression line below 1 could occur when spots appear within the area that are no longer active, which could be an indication of stabilisation, which occurs mainly in the second and third time period.

The p values of the coefficients are mostly below the alpha level of 0.05, so it is assumed that the relationships between sediment-contributing area and sediment yield are statistically significant in almost all cases (∼92 %) (Fig. 6). To determine the proportion of the variance in the dependent variable that can be explained by the independent variable, the R squares (R² or the coefficient of determination) of all regression lines were analysed (Fig. 6; Appendix A). The relationship between sediment-contributing area and sediment yield shows varying correlations within the site and the time periods. The median R² values range from 0.59 to 0.91 in the first time period, from 0.37 to 0.93 in the second time period, and from 0.3 to 0.94 in the third time period (Fig. 6; Appendix A). The number of channels modelled differed between the time periods at the same sites due to the different quality of the DEMs and the slightly different size of these. As in Heckmann and Vericat (2018), the accumulation of DoD values resulted in very small positive values at some sites. Such errors are due to the quality of the DoD, different bulk densities of eroded vs. deposited materials, and the inability of the flow-routing algorithm to fully reproduce sediment transfer in reality especially when flow directions changed within one time period. Where positive accDoD values occurred, they were small and manually corrected to the zero.

Nevertheless, the D8 algorithm simplifies complex sediment transport processes such as fluvial activity, landslides, and debris flows, which have different frequencies, magnitudes, and forms of erosion and accumulation. As the individual time periods cover several years, no reference can be made in this study to individual processes that can be attributed to extreme precipitation events or to seasonal differences. Therefore, we compare different time periods based on mean annual sediment yield, which includes all geomorphological processes. Accordingly, the aim was not to model individual erosion processes but to compute sediment yield of each cell. The individual sites have a slightly different area within the different time periods, which is mainly due to headcut retreat (Heckmann and Vericat, 2018; Betz-Nutz et al., 2023). The lateral boundaries also changed slightly due to the quality of the DEMs and geomorphological slope processes, while the lower boundary did not change.

By processing historical aerial photographs into DEMs (by SfM–MVS), the temporal aspect of sediment-contributing-area studies could be quickly and cost-effectively extended to several decades (up to the 1950s), which previously spanned only a few months or several years (Neugirg et al., 2015a, b, 2016; Dusik, 2019; Dusik et al., 2019). However, as Schiefer and Gilbert (2007) have already shown, the shorter the time intervals and the lower the quality of the aerial images, the more difficult it becomes to detect surface changes, so in the process of this study, several series of aerial images had to be sorted out that were actually available due to a poor data quality. Furthermore, it should be noted that the accuracy and precision of the historical DEMs strongly depend on the respective generation, for example, whether they were generated with or without a calibration certificate (as was the case, for example, with the 1959 aerial photo series in the Horlachtal/glacier foreland Hohenferner), which ultimately influences the sediment-contributing-area results and the calculated erosion rates (Stark et al., 2022).

5.2 Geomorphic activity

The geomorphic activity is directly related to the characteristics of the sites. The sites with the highest mean annual sediment output (> 100 m³ yr⁻¹) (such as KG1, KG2, KG4, KG5, MH2, KM2, and KW2) show strong gully formation and are overall characterised by larger areas, longer max slope lengths, and higher mean and max slope gradients (Table 1). In contrast, the sites with lower mean annual sediment output (<100 m³ yr⁻¹) (such as KW1, KG3, HG1, KM1, and MH1) show less gully incision. These sites tend to be characterised by smaller areas, smaller max slope lengths, and smaller mean slope and max slope gradients (Table 1). The strong influence of slope length and slope gradient on morphodynamics is also shown by previous studies (Ballantyne and Benn, 1994; Curry, 1999; Curry et al., 2006; Betz-Nutz et al., 2023). KG3 also appears to be somehow stabilised by bedrock in the lower part of the slope, which could mitigate the erosion of this site, as also shown by Jäger and Winkler (2012). Elevation and aspect, however, do not seem to have an influence on geomorphic activity, which is also shown in the study by Curry et al. (2006). Since only bare and sparsely vegetated areas were investigated, no findings on the influence of vegetation on morphodynamics can be made in this study. Solifluction processes could also not be observed, probably due to the composition of the moraine material. Presumably, the morphodynamics are still so high that the vegetation does not yet have the opportunity to develop accordingly. In general, we assume that debris flows are the most common process, as described, for example, by Ballantyne (2002a) and Curry et al. (2006). Thus, material stored in the gullies is transported downslope by debris flows mainly as rain or snow events in the spring or heavy rainfall events during rainstorms in the summer months (Ballantyne and Benn, 1994; Ballantyne, 2002b; Curry et al., 2006; Dusik et al., 2019).

However, the high mean annual sediment yield and corresponding erosion rate in the first (3471 m³ yr⁻¹, 465 mm yr⁻¹) and second (2922 m³ yr⁻¹, 245 mm yr⁻¹) time periods of site KG1 (Fig. 7) can probably also be attributed to individual landslides and deep-seated slope failures that are in some cases linked with melting dead ice bodies, as these processes are more likely to occur after deglaciation and are characterised by high magnitude and low frequency, which has also been shown by Blair (1994), Hugenholtz et al. (2008), and Cody et al. (2020). On the less incised slopes (e.g. MH1), small-scale processes such as fluvial erosion or snow drifts probably occur (Betz-Nutz, 2021), which ultimately show no clear trend in the increase or decrease in the morphodynamics but can be described as a constant geomorphic activity. In Betz-Nutz et al. (2023) and in this study, six similar lateral moraine sections (although other exactly defined sites) were investigated. The test sites KG1, KG2, KW1, KW2, and MH1 (in Betz-Nutz et al., 2023; GPF1, GPF2, WSF1, WSF2, and HF1) showed similar erosion rates and the same long-term trends. In the case of site MH2, different trends were determined (stagnation in Betz-Nutz et al., 2023, and a decrease in this study), which can be attributed to the differently defined site and the slightly different study period.

In the sense of a process–response system, it is noticeable that the first-mentioned group of sites with the higher erosion rates (except KG1) had considerable influence from melting dead ice in the lower slope area at least until 2006 (KG2, KM2, KW2, and MH2), or the glacier was still present at the bottom of the slope (KG4 and KG5), which could be identified by the interpretation of the DoDs. Melting of the dead ice can lead to destabilisation of the slope, which can enhance erosion processes of the upper slope areas, as the support is no longer present, the sediment becomes saturated, and there can be an increase in the slope gradient due to the subsidence of the lower part of the slope (Altmann et al., 2020; Betz-Nutz et al., 2023). However, the highest slope gradients are also present here, which also plays a major role. In addition, site HG1, where erosion is increasing, shows an undercutting of the slope by the adjacent stream, which leads to a destabilisation or lowering of the erosion base and a typical formation of a debris cone and alluvial fan, with a successive reduction in the slope gradient that is missing (Fig. 14). It can be assumed that individual strong rainfall events in the second and third time periods in combination with changing flow pathes due to the retreat of the Grastalferner acted here as an impulse and affected both the site itself and the adjacent stream.

Figure 14Overview of the DoDs of the corresponding time periods of site HG1. (a) DoD 1954–1973, (b) DoD 1973–2006, (c) DoD 2006–2019, (d) orthophoto 2020 (provided by the province of Tyrol), and (e) photo of the site from 2019 by Anton Brandl.

5.3 Paraglacial landscape adjustment

5.3.1 The sediment activity concept

The finding of mainly decreasing geomorphic activity of the LIA lateral moraines in this study is largely consistent with previous model-based studies describing the paraglacial landscape adjustment with a decrease in geomorphic activity in proglacial areas over time, such as the theoretical model “paraglacial concept” of Church and Ryder (1972) or the “sediment exhaustion model” of Ballantyne (2002a, b). The geomorphic activity of gully systems is given as a few decades to centuries (Ballantyne, 2002a, b). Furthermore, it is stated that there is a high temporal and spatial variability in this development. The model provides an appropriate approximation (Ballantyne, 2002a, b). Within the paraglacial adjustment process, different geomorphological processes result in different durations of occurrence. Furthermore, different land systems react at different rates and on different spatial scales. Thus, external perturbations can occur, leading to secondary peaks and time delays (Ballantyne, 2002a, b). While 10 out of 12 sites fit the model descriptions, two test plots show opposite morphodynamics, which can be described as a delay of the paraglacial adjustment process or that the response systems can run counter to such an adaption.

To estimate the changing morphodynamics, we therefore propose the following simplified description of the landscape evolution using the sediment activity concept, based on the results of this study. Due to the study design, this concept is only valid for geomorphologically active areas (in this case, the upper lateral moraine section) on the LIA lateral moraines and until about 170 years after the end of the LIA (Table 6; Fig. 15). The concept distinguishes between an earlier and a later phase. The earlier phase (mainly 1950s to 1970s) is characterised by a wide range between areas with high and low variability in sediment yield within the area, as well as high and low mean annual sediment yield (erosion rate/volume). In contrast, the later stage (mainly 1970s to 2000s and 2000s to the end of the survey 2017/2019) shows a decrease in this range. Although the decrease in the morphodynamics predominates, there are also increases in morphodynamics. The time of ice release was not integrated here, so that the time periods refer to the actual time. In addition, we give two examples in Fig. 15. While example A shows high variability (polynomial behaviour) within the area and high sediment yield (erosion rates/sediment output), example B shows low variability (constant/linear behaviour) and low sediment yield (erosion rates/sediment output). Ultimately, the relationship between sediment yield and the size of the catchment has changed so that erosion within the area is more constant today.

Figure 15The sediment activity concept. Description and illustration of the change in sediment activity over time (a) and two corresponding examples (b).

Download

Table 6Tabular summary of the simplified conceptual model.

Download Print Version | Download XLSX

Through introducing the sediment activity concept of this study, we present a different description of the paraglacial adjustment process which is based on the actual sediment yield. However, the concept is only valid until about 170 years after the end of the LIA and at the sites of this study. The sediment activity concept presented here is also compatible with the results of other studies; for example, Betz-Nutz et al. (2023) mostly show a decrease in erosion rates over a similar time period but also a remaining at similar levels and an increase. In addition, the studies by Church and Ryder (1972), Ballantyne and Benn (1996), Curry (1999), Ballantyne (2002a), Curry et al. (2006), and Schiefer and Gilbert (2007) show a decrease in geomorphic activity over time, which is also consistent with the model presented here. Also, we assume that the concept can also prove its validity in further proglacial active areas and partly over an even longer period of time.

5.3.2 Erosion rates

The methodologies for determining long-term average erosion rates in proglacial areas are based on gully volume estimates (Ballantyne and Benn, 1994; Curry, 1999; Curry et al., 2006) and sediment volume calculations due to surface changes using DoDs, as shown by Betz-Nutz et al. (2023) and this study (Table 7). In glacier forelands in western Norway, this amounts, for example, to minimum estimates of 50–100 mm yr⁻¹ (max estimates of min 200 mm yr⁻¹) (Ballantyne and Benn, 1994) and a minimum of 5.5–169 mm yr⁻¹ (in different glacier forelands) (Curry, 1999). A further study in the Swiss Alps shows erosion rates of min 49–151 mm yr⁻¹ (in different glacier forelands) (Curry et al., 2006). The work of Betz-Nutz et al. (2023) and this study show erosion rates over several decades and distinguish between different time periods, making it possible to show differences between them. Both studies show that the mean erosion rates in the individual time periods decrease (Table 7), although in individual cases there is also a constant and an increase in erosion rates over time. Although there has been a clear decrease in geomorphic activity, stabilisation of the sites is not yet apparent, which means that the paraglacial adjustment is still ongoing. Within this study, we observed that the sites still show a high geomorphic activity even after they have been deglaciated for 76–159 years. A stabilisation of the gully systems, as shown by Curry (2006), cannot be observed. Other studies such as Lane et al. (2017), Dusik (2019), Altmann et al. (2020), and Betz-Nutz et al. (2023) also show the still ongoing paraglacial adjustment processes. Comparing the long-term erosion rates of gully systems from the different studies ultimately shows high variability in the adjustment, as these studies were also conducted using different methods, on different timescales, and in different regions. Differences are probably mainly due to the different local conditions, such as the geomorphological settings; these include, for example, the different characteristics of the lateral moraine sections, such as slope gradient, slope length, time of ice exposure, dead ice influence, and the development of vegetation. Furthermore, the lateral moraines have different sedimentological characteristics related to their genetic origin. In addition, different meteorological conditions prevail in the different regions.

Table 7Studies on long-term erosion rates (several decades) of gully systems on the LIA lateral moraines in different glacier forelands.

Note: HT is for Horlachtal, KT is for Kaunertal, and MT is for Martelltal.

Download Print Version | Download XLSX

5.4 Meteorological drivers

The decrease in the mean annual number of ice days and the increase in the number of warm spells over the different time periods and the associated potential increase in snowmelt on the slopes could also lead to an increase in morphodynamics, as these processes represent important preparatory steps for erosion processes in spring (Haas, 2008), such as increased saturation of the slope du to snowmelt, loosening of the upper sediment layers or the delivery of material by snow slides or small wet avalanches that is then available for debris flows in the summer months (Dusik et al., 2019). Klein et al. (2016), for example, also show an increase in the frequency and intensity of snowmelt in the Swiss Alps. Mean annual precipitation decreases slightly across time periods but is not statistically significant (except for winter precipitation for the entire study period (1959 to 2019) and summer precipitation in the second time period from 2005 to 2019 in the Martelltal). Other studies also show that the decrease in precipitation in the European Alps is low (Brugnara et al., 2012), that there is no clear trend in precipitation (Hock et al., 2019), or that it is mainly subject to regional influences and decadal variations (Mankin and Diffenbaugh, 2015). Extreme precipitation events (1 h resolution) and daily precipitation totals also show only minor changes. Differentiation of precipitation, on the other hand, shows a clear increase in rainfall and a decrease in snowfall, which is also shown by Serquet et al. (2011), Beniston et al. (2018), and Hock et al. (2019), who found that the rainfall on snow events in spring as a preparatory factor for the erosion processes in the summer months increase. The simulated meteorological data generally show lower temperatures and larger precipitation amounts when compared to three automatic weather stations operated by TIWAG (Tyrolean Hydropower AG, Innsbruck, Austria). These stations are located in the vicinity of our sites. The simulated mean annual temperatures extracted at the location of the weather stations at Horlachalm (1987–2015) (approx. 6.5 km linear distance to the site in the Grastalferner glacier forefield) and Weißseeferner (2007–2015) (approx. 500 m linear distance to the sites in the Weißseeferner glacier forefield), covering the same time period indicate a difference of −1.05 and of −0.87 °C, respectively, after accounting for differences in elevation. However, at Gepatschalm (2010–2015) weather station (approx. 2.5 km linear distance to the sites in the Gepatschferner glacier forefield), the difference between the simulated and observed mean annual temperatures is 0.13 °C, indicating that the magnitude of the discrepancies depends on the station data used for the comparison. The simulated precipitation, however, is generally larger with mean annual precipitation sums of 1531, 1655, and 1820 mm at the location of Horlachalm (1990–2015), Gepatschalm (2010–2015), and Weißseeferner (2007–2015), respectively, while the weather stations recorded values of 803, 1086, and 924 mm, indicating large discrepancies, especially when compared to Horlachalm and Weißseeferner weather stations. The data sets from which the temperature and precipitation were extracted are both based on coarsely resolved data, which makes a comparison with measurement data in the field difficult, although the corresponding trends are able to be used well for this study. The large difference between simulated and recorded precipitation is mainly due to winter precipitation (Fig. 11) when the weather stations are not always able to record total snowfall accurately; additionally, fog precipitation or precipitation in combination with stronger winds are not recorded correctly.

The weather and climate study periods are based on the predefined study time periods given by the availability and quality of orthophotos and not on the usual climate periods. This results in large differences in the length of the different time periods, which must be taken into account.

There are several sources of uncertainty in the simulated data, amongst them the dynamic initial and boundary conditions, as the forcing data have their own sources of uncertainties. Furthermore, the choice of the reanalysis data used for forcing the model has an influence on the final results. Additionally, for such a long simulations, an updated sea surface temperature (SST) is recommended. Since there are no SSTs available for the 20CRv3, we have generated SST fields from the skin temperature. Other sources of uncertainty are the static boundary conditions like the fixed land use categories and topography, as well as model simplifications and choices in the parameterisation of the physics and dynamics. In our simulations, we have used spectral nudging in order to keep the model from large deviations from the forcing data. Short test runs indicate that the use of spectral nudging improves the simulated data, especially with respect to precipitation. However, the strength of nudging also has an influence on the final results. Since the purpose of this study is not to test how strongly to nudge, we have used the default values in WRF.