**3. Results**

#### *3.1. LULC Change*

The analysis of almost 200 years of changes in LULC displayed distinct trends in the two municipalities (Figure 2).

Even though forest covered more than one third of the total land area in both case studies throughout the investigation period (Figure 2A), in Paldau, its area decreased by 113.2 ha (8%) from 1820 to 2015, while in Waidhofen, the forest area increased by more than 1000 ha (18%) from 1820 to 2015. In both areas, cropland extent declined while grassland and "settlement and other" expanded. In Waidhofen, in 2015 only 1% of the area was used as cropland (1820: 26%), while in Paldau, the share diminished from 43% to 37%. Areas classified as "settlement and other" increased in Waidhofen more than fivefold from 1820 to 2015 (246 ha to 1431 ha) and in Paldau more than sevenfold (104 ha to 750 ha).

Biomass extraction nearly doubled in both municipalities from one time cut to the next (Figure 2A; Paldau: 4969 to 10,630 to 18,674 kg/ha/a, Waidhofen: 4661 to 8410 to 16,773 kg/ha/a). Cereal yields increased most strongly in both municipalities (factors 9.8 and 7.9 in Waidhofen and Paldau, respectively). Grassland yields increased less strongly and reached their maximum in 1960. In both study areas, wood extraction declined slightly from 1820 to 1960 (Waidhofen: −6%, Paldau: −2%), but it more than doubled from 1960 to 2015 (Waidhofen: 5067 kg/ha/a, Paldau: 5667 kg/ha/a). In 1820, wood extraction had the largest share in biomass extraction in both municipalities.

The landslide distribution for the different LULC classes and time cuts showed contrasting patterns in Waidhofen and Paldau (Figure 2C). While in Waidhofen, the change pattern in landslide distribution and LULC distribution by LULC category was similar (Figure 2A,C), in Paldau the change pattern showed clear differences. There, 82.4% of the landslides were located in present-day (2015) forest areas, which also showed the highest landslide frequency ratio (FR) of 2.5 (Figure 2C); in Waidhofen, 42.2% of the landslides were found in present-day forested areas with a FR of nearly 1. In Waidhofen in 2015, grassland showed the highest FR of 1.14 (Paldau: FR of 0.78). Furthermore, in Paldau, 69% of the landslides were located in continuously forest-covered areas (28.7% for LULC change pattern) and in Waidhofen 20.7% (27.5% for LULC change pattern), respectively. In both study areas, landslides were least frequent on present-day cropland (Paldau: FR 0.05, Waidhofen: 0.31).

**Figure 2.** LULC change in the study areas for different time cuts. (**A**): Changes of LULC patterns (in %, left) and changes of biomass extraction (in kg/ha/a, right). (**B**): Cultural landscape persistence maps showing persistent historical structures in present-day LULC. (**C**): Landslide distribution for different LULC patterns (in %) and their changes (left), and landslide frequency ratios for the LULC pattern in 2015 (right).

*3.2. LULC Legacy Effects on Landslide Occurrence*

3.2.1. Model Performance and Transferability

We applied SpCV with the AUROC as performance measurement to gain information on the model's capabilities to discriminate landslides from non-landslides observations. The performance assessment showed distinct differences between the study areas (Figure 3, Table A4 in the Appendix B).

**Figure 3.** Model performance in SpCV and transferability assessment. (**A**): Study-area specific SpCV and transferability estimates. (**B**): SpCV estimates using the combined data from both municipalities.

For Waidhofen the median AUROC (mAUROC) in SpCV ranged between 0.78 for GAM-1820 and 0.80 for GAM-2015-Masked, i.e., acceptable to excellent discrimination capabilities. In contrast, in Paldau, GAM-Base had the lowest mAUROC of 0.88 and GAM-2015-Masked the highest mAUROC of 0.93, i.e., excellent to outstanding discrimination. AUROC estimates for Waidhofen were substantially more variable than for Paldau (interquartile ranges, IQR, 0.13 to 0.14 versus 0.03 to 0.06). Model transfer to the other study area resulted in a strong drop in model performance, with a decrease of 0.11–0.14 for model transfer to Waidhofen and 0.05–0.12 in Paldau. As in SpCV, GAM-2015-Masked performed best (Waidhofen: 0.69, Paldau: 0.85), followed by GAM-Base (Waidhofen: 0.65, Paldau: 0.83) in both study areas.

In the performance assessments for the combined data, the mAUROC and IQR values of the combined validation data fall in between the estimates of Paldau and Waidhofen (Figure 3B, Table A4 in the Appendix B). GAM-Base had the lowest mAUROC of 0.81 and GAM-2015-Masked the highest mAUROC of 0.84.

Regarding the effect of LULC legacy predictors on the performance, for Paldau the performance estimates were marginally, but ye<sup>t</sup> significantly higher with the inclusion of legacy information (order of significant "<" and non-significant "=" differences in AUROCs: GAM-Base < GAM-1960 < GAM-1820 = GAM-2015 < GAM-2015-Masked, Table A5 in the Appendix B). However, for Waidhofen, such a tendency was not identifiable (GAM-2015 < GAM-1820 = GAM-Base < GAM-1960 < GAM-2015-Masked, Table A5 in the Appendix B). Furthermore, for both study areas the landslide model excluding observations located in continuously forest-covered areas (i.e., GAM-2015-Masked) showed the highest estimates in SpCV (combined and non-combined data) and transferability assessment.

#### 3.2.2. Variable Importance

The assessment of variable importance showed differences in the variable ranking of the study areas, although slope angle was shared as a top-tier predictor (Figure 4, Table A6 in the Appendix B). In Waidhofen, the top three ranks were identical across all settings: 1. slope angle (mDD 7.45–8.63%), 2. lithology (mDD 3.72–6.38%), and 3. plan curvature (mDD 2.09–2.83%). In Paldau and across all settings, four different variables occurred in the top three ranks: slope angle (mDD 4.38–7.68%), profile curvature (mDD 3.05–6.37%), slope aspect (S-N, mDD 3.94%), and a LULC legacy variable (LULC 2015 or LULC legacy 1820/1960, mDD 3.64–3.92%). In all models, slope angle was the most important variable, except for Paldau's GAM-2015-Masked, where profile curvature was more important (mDD 6.37%).

**Figure 4.** The three most important variables in terms of decrease in deviance explained [%]. Refer to Table A6 in the Appendix B for an overview of variable importance for all predictor variables.

The top-ranked variable was always a land surface variable, specifically slope angle and in one setting profile curvature (Table A6 in the Appendix B). Soil variables were less important with a highest rank of 11 with 0.62% mDD in Waidhofen (hydraulic conductivity in GAM-2015-Masked), and of 7 with 1.3% mDD in Paldau (hydraulic conductivity in GAM-2015-Masked), respectively. LULC legacy variables were less important in Waidhofen (mDD 0.26–1.07%), while for Paldau, they showed a high mDD of 3.64–3.92%. Lithology, in contrast, was only important in Waidhofen (mDD 6.34–6.38%; Paldau: 0.25–0.39%).

#### 3.2.3. Predictor-Response Relationships

In both study areas, based on the transformation function plots in logit scale, we found a general agreemen<sup>t</sup> on the predictor-response relationships among all model settings, except for the GAM-2015-Masked (Waidhofen: Figure A2, Paldau: Figure A3). However, some differences in the predictor-response relationships for the three most important variables could be identified (Figure 5).

**Figure 5.** Comparison of predictor-response relationships of the three most important predictors (left column, (**<sup>a</sup>**–**<sup>c</sup>**) and (**g**–**i**)) and the LULC legacy variables (right column, (**d**–**f**) and (**j**–**l**)) for each study area. Grey: 95% Bayesian credible interval of GAM-Base. Lithological units in Waidhofen: Reference unit: 'Upper Austroalpine limestone', 0: talus and glacial deposits, 1: Inneralpine Neogene, 2: Klippen zone, 3: flysch zone, 4: Upper Austroalpine marls. Note: the y axes are plot-dependent, and the x axes of non-parametric transformation functions are limited to the 5th to 95th percentile range.

In both study areas, the chance of landslide occurrence was higher on steeper slopes. Among the model settings, the predictor-response relationship of GAM-Masked-2015 indicated the comparatively highest rOR compared to GAM-Base (e.g., for 15◦, Paldau: 1.4, Waidhofen: 1.6). Minor differences of the other models in reference to GAM-Base were mostly at lower and upper quartiles of the value ranges (Figure 5a,g), but fell within its 95% Bayesian credible interval.

In Waidhofen, for the variable plan curvature, the chance of landslide occurrence was higher on concave than on convex surfaces, but without substantial differences between the models (Figure 5b). Regarding the variable lithology (Figure 5c), each lithological unit showed a higher chance of landslide occurrences relative to Upper Austroalpine limestone (e.g., GAM-Base: OR of 13 for flysch zone). However, among the model settings, the rORs of GAM-Masked-2015 were generally lower by a factor of 0.6, and fell partially outside the Bayesian credible interval of GAM-Base.

In Paldau, for the variable profile curvature, the chance of landslide occurrence was lower on concave than on convex surfaces (Figure 5h). Substantial differences between the models were identifiable only for GAM-2015-Masked on convex surfaces with a curvature greater than 0.005 m<sup>−</sup><sup>1</sup> (e.g., landslides were 1.3 times more likely relative to GAM-Base). Among the model settings, south-exposed slopes were less susceptible to landslides than north-exposed slopes, with only marginal differences between models (Figure 5i).

Regarding the LULC variables, in Waidhofen LULC 2015 was significant only in the GAM-2015-Masked setting for "settlement and other" with an OR of 0.43 relative to forest areas (Figure 5d). Using additionally the biomass extraction as LULC legacy interaction term (i.e., GAM-1960 and GAM-1820), the predictor was significant and showed a higher chance of landslide occurrence with higher biomass extraction for each LULC class (Figure 5e,f). Compared to forest areas, landslides were more likely on cropland and less on grassland and "settlement and other". However, using the long-term legacy in GAM-1820, the contrast relative to forest areas showed higher ORs for cropland (e.g., for 10,000 kg FW/ha/a: OR of 2.95 for GAM-1820 and of 2.53 for GAM-1960), but lower OR for grassland (OR of 0.83 for GAM-1820 and of 0.46 for GAM-1960) and "settlement and other" (OR of 0.67 for GAM-1820 and of 0.43 for GAM-1960) compared to GAM-1960.

In Paldau, for the variable LULC 2015 in GAM-2015 and GAM-2015-Masked, all LULC classes were significant terms and also showed a similar pattern with minor differences with respect to forest area (reference level; Figure 5j): For GAM-2015, the chance of landslide occurrence on cropland was only 0.07 times as high as in forest areas (OR of 0.06 for GAM-2015-Masked), 0.32 times as high on grassland (OR of 0.33 for GAM-2015-Masked), and 0.26 times as high on "settlement and other" (OR of 0.25 for GAM-2015-Masked). Including the biomass extraction as LULC legacy interaction term (i.e., GAM-1960 and GAM-1820; Figure 5k,l), for both models, the coefficients of grassland and "settlement and other" were not significant anymore, while with a higher biomass extraction the chance of landslide occurrences on cropland was lower (e.g., for 10,000 kg FW/ha/a: OR of 0.16 for GAM-1960 and of 0.27 for GAM-1820) and higher on forest areas (OR of 3.4 for GAM-1960 and of 3.7 for GAM-1820), respectively. Using the long-term legacy in GAM-1820, the ORs of the LULC classes relative to forest areas were generally lower compared to GAM-1960 (e.g., for 10,000 kg FW/ha/a on cropland: OR of 0.05 for GAM-1960 and of 0.07 for GAM-1820).

We extracted reported relationships of comparative studies, in which the landslide dates were known (for Paldau: event-based landslide inventory of southeast Styria in [73]) or the relative landslide age could be approximated (Waidhofen: remotely-sensed landslide inventory based on [21] used e.g., by [74]). For the land surface variables, we generally found a good match in the shape of the predictor-response relationships in both study areas (Figure 6), with the exception of slope aspect (S-N) in Paldau. Contrary to our results, in Knevels et al. [73] the predictor slope aspect was less influential with south-exposed slopes being slightly more susceptible to landsliding than north-exposed slopes.

Regarding the LULC variables, for Waidhofen, Steger et al. [74] reported that landslides were 1.7 to 2.0 times more likely on arable land or pastures compared to forests, while in our study, LULC was only significant when biomass extraction was included (i.e., GAM-1960 and GAM-1820). For southeast Styria, Knevels et al. [73] found non-forest areas more than three times as susceptible to landsliding as forest areas (OR > 3.6). Yet, in Paldau's GAM-2015 and GAM-2015-Masked, we found contrasting results (OR < 0.33 outside forest).

**Figure 6.** Comparison of predictor-response relationships between event-based landslide and historical LiDAR-derived landslide inventories. (**A**): Event-based inventory of Waidhofen (left, Table 1 in [74]) and southeast Styria (right, [73]). (**B**): Historical LiDAR-derived inventory of Waidhofen (left) and Paldau (right). Note: the logit scale is variable-dependent, the intercept for slope angle in Waidhofen was artificially set, and the x axes are limited to the 5th to 95th percentile range.

## **4. Discussion**

#### *4.1. Initial Objective: Effect of LULC Legacy on Modeling and Biases*

We were able to successfully link the historical LULC characteristics covering a period of almost 200 years in two Austrian municipalities to landslide occurrences identified in an airborne LiDAR-derived HRDTM. The established landslide susceptibility models showed performances with acceptable to outstanding discrimination capabilities, confirming the proposed modeling approach. Landslide models including LULC legacy predictors had significantly higher (for Paldau) or at least equal (for Waidhofen) performances compared to a reference model without LULC (GAM-Base, Table A5 in the Appendix B), demonstrating the potential of LULC legacies for explaining landslide susceptibility today.

The use of LULC legacies may improve the understanding of LULC dynamics and landslide occurrences. The LULC legacy categories were generally more susceptible to landsliding with a higher biomass extraction (excluding cropland in Paldau; Figure 5e,f,k,l). In particular, present-day forested areas previously used for agricultural activities were more prone to landsliding than continuously forested areas, confirming findings for less extensive historical time periods by Beguería [9] and Persichillo et al. [11]. Thus, the effect of LULC legacies on landslide occurrences might be a good explanation of contradictory empirical observations such as landslide events in forests [23]. Additionally, for Waidhofen, we conclude that forest areas are more susceptible to landsliding than grassland, after accounting for biomass extraction. While Goetz et al. [7] reported that forest harvesting may temporally lead to open or semi-open forest types and thus reduced slope stability due to different soil hydrological and mechanical conditions (e.g., reduced rainfall interception and root cohesion), Tasser et al. [76] found managed grassland to be significantly less erodible than abandoned areas. However, we acknowledge that potentially landslide mitigating land managemen<sup>t</sup> strategies are not incorporated in the analysis, and that the difference between forest and grassland may also be related to inventory biases.

Inventory biases are a known issue in landslide modeling [25,26]. We sugges<sup>t</sup> that the identified substantial difference in the predictor-response relationship for LULC compared to other studies, is due to the underlying airborne LiDAR-derived historical landslide inventory. In particular, we assume that in Paldau, the high number of landslides in continuously forest-covered areas (69%, FR of 2.5 in 2015) can be explained by an inventory bias, which ultimately produces biased predictor-response relationships (OR of 0.33 for outside forest relative to forest areas vs. OR of >3.6 in Knevels et al. [73]). A bias was also detectable in Waidhofen's inventory, albeit weaker (about 20% of the landslides in continuously forested-covered areas). This is consistent with the reported drawbacks of such landslide inventories [26,27]. Unreflecting model interpretation might thus lead to the conclusion that forest areas are more susceptible for landslides than other LULC categories (i.e., contrasting the findings of [4,5,13,73]). We assume that the inventory bias results from the fact that landslide traces on agricultural land are quickly "tilled away" (Figure 1D), and that landslides affecting the built environment are removed during reconstruction, while in forests, they are preserved for an extended period of time (see Figure 12 in [28]).

For bias mitigation, we explored three approaches. Following the recommendations of Petschko et al. [27], (i) dropping the bias-describing predictor from modeling (i.e., LULC; GAM-Base) had only minor effects on the predictor-response relationships compared to the other models (excluding GAM-2015-Masked). We assume that confounding effects may explain these minor differences (e.g., slope angle: forest area more likely located on steeper slopes [26]). Besides, (ii) excluding inventory-biasing observations (i.e., continuously forestcovered areas; GAM-2015-Masked) led to substantial differences of predictor-response relationships for some variables (e.g., slope angle in both study areas, profile curvature in Paldau, lithology in Waidhofen, Figure 5), ye<sup>t</sup> the OR of the LULC categories in reference to forest were similar to GAM-2015. However, since the relationships are in general agreemen<sup>t</sup> with other research results [73,74] (Figure 6), further research might help to evaluate and identify the "true" relationship. Moreover, reducing the models' training data may not always be suitable (e.g., sparse inventories), and might even decrease model performances [77,78]. (iii) Using LULC legacy information (i.e., biomass extraction, GAM-1820 and GAM-1960), the predictor-response relationships matched GAM-Base except for the minor differences mentioned above, and we could reduce the bias of a higher landslide chance on forest area relative to other LULC categories. Additionally, we identified a tendency towards relatively lower OR in forests using long-term legacies. Albeit the landslide inventory bias still had an effect on the landslide models for all approaches, for studies based on airborne LiDAR-derived landslide inventories we recommend to exclude bias-describing predictors (e.g., LULC) in the first place, or to enrich the feature space with a bias-reducing information predictor (e.g., LULC legacies).

Moreover, when considering LULC as a predictor, it should be kept in mind that each study area has its unique legacy. The transferability assessment showed a poorer transferability of landslide models that were fit only in one study area (AUROC loss up to 0.14, Figure 3). While the transformation functions of slope angle, plan curvature and the LULC legacy category forest matched between study areas, for other important variables they were partly (e.g., concave profile curvature) or even completely opposed (e.g., slope aspect, S-N and other LULC legacy categories; Figures A2 and A3). Moreover, landslide models using the combined data did not improve performances compared to both area-specific models. However, for event-based landslide modeling with especially sparse inventories, the opposite may be true (see Figure 4 in [73]).

Regarding the variable importance assessment, the results are in agreemen<sup>t</sup> with findings in other studies despite different landslide inventories: For Waidhofen, Steger et al. [74] also identified the slope angle as the most influential predictor, while present-day land cover ranked last. In the landslide model of Knevels et al. [73] fitted in southeast Styria (where Paldau is located), the present-day LULC (categorized into forest classes) and slope angle ranked second and third, respectively. We conjecture that the lower influence and different shape of the slope aspect relationship (Figure 5) could be due to effects such as pre- or post-failure HRDTM or the specific rainfall event in the study of [73].

#### *4.2. Study Data: Challenges and Requirements*

We encourage the use of long-term LULC legacy data in landslide studies. The presented LULC legacy predictor is not only a simple and comprehensible way to highlight historical LULC dynamics but may also support spatial planning in preventive disaster reduction regarding future LULC changes. With historical cartographic documents of several countries (e.g., Habsburg, Franciscan or Napoleonic cadastral map) [79] and information on long-term biomass extraction [80,81] becoming increasingly available in a pre-processed form, there is a real opportunity to further explore LULC legacies in future studies.

For the creation of LULC legacy data, we used data sources of different spatial, temporal and thematic quality and resolution (Table 1). For the time cut of 1960, the available greyscale aerial photographs for the municipalities were nine years apart (Paldau 1953 and Waidhofen 1962). Thus, LULC changes in the context of the "economic miracle" (e.g., mechanized agriculture, urban growth and sprawl) beginning in 1955 may not ye<sup>t</sup> be present in Paldau [82]. Also, the interpretation of the greyscale aerial photographs of low image quality was particularly challenging (especially grassland vs. cropland). Additionally, the thematic resolution of the data sources used for the biomass extraction (i.e., wood and agricultural yields) was very heterogeneous in terms of the aggregated spatial statistical reference unit (e.g., the data availability for wood yields was for time cut 1960 on a national scale and for 2015 on for federal states). Nevertheless, we sugges<sup>t</sup> that the data still allows conclusions to be drawn about long-term trends [43,45,81]. We acknowledge that the digitization of LULC legacy data was a time-consuming process and historical imagery at an appropriate scale may not be available everywhere.

The landslide inventories were separately created for each study area by local experts. There is evidence that landslide inventories are not only incomplete to an unknown degree [24], but also that landslide inventories compiled by different experts may have positional mismatches up to 70% [23]. Therefore, in this study, consistent rules were agreed upon to ensure or at least strive for a homogeneous inventory quality (e.g., uniform digitization scale, order of digitization: first landslide scarp, then body, etc.). To overcome the drawbacks of the inventory (e.g., bias in forest areas) and to enhance the explanatory power of the LULC legacy variables in modeling landslides, we encourage to collect and analyze additional historical aerial photographs, orthophotos or satellite images in order to estimate landslide ages [4,11,27].

#### **5. Conclusions and Outlook**

In this study, LULC and legacies of its change could successfully improve the explanation of landslide distribution as we identified the LULC variable as a meaningful predictor in landslide susceptibility modeling. Higher biomass extraction resulted in higher landslide susceptibility (excluding cropland in Paldau), explaining different risk levels in areas with the same present-day LULC. Furthermore, we could confirm that airborne LiDAR-derived inventories may be biased towards currently forested areas (e.g., high landslide density in Paldau with a FR of 2.5). Using long-term LULC legacy variables accounting for changes for almost 200 years, we could successfully reduce the effects of LULC-related inventory bias. However, without any information on the failure date or at least an approximated time slice, the additional landslides preserved in forest areas may still lead to an unknown bias in the model, and thus lead to potentially contradictory relationships between landslide occurrences and LULC legacies (e.g., OR < 0.33 outside forest relative to forest areas in Paldau). Other bias-avoiding strategies such as removing inventory-biasing observation led to improved model performances, but also to different predictor-response relationships. Thus, the implementation and interpretation of LULC as a bias-describing predictor must be carefully considered. We highly encourage further research using event-based landslide inventories with known landslide ages or age ranges to avoid such biases.

The construction of historical LULC datasets is relevant for future use both from methodological and empirical perspectives. Methodologically, the approach chosen here can be used in future analyses in cases where both geographic and statistical information on LULC are available and can be combined. Empirically, the dataset established may be used for future analyses of different issues related to long-term LULC change, in particular those related to sustainable LULC intensification (e.g., cultural landscape change [45]). With the publication of the historical LULC legacy dataset [49], we highly encourage future research with this data and a replication of the analysis, especially once a reliable event-based landslide inventory for the study areas is available.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2073-445 X/10/9/954/s1: Landslide inventories and landslide susceptibility models of both study areas.

**Author Contributions:** Conceptualization, R.K., A.B. and H.P. (Helene Petschko); data curation, T.L., H.P. (Herwig Proske), G.H., R.K., C.P. and S.G. formal analysis, R.K.; funding acquisition, P.L., S.G. and H.P. (Helene Petschko); methodology, R.K., A.B. and H.P. (Helene Petschko); supervision, H.P. (Helene Petschko) and A.B.; writing—original draft, R.K.; writing—review & editing, R.K., A.B., S.G., G.H., T.L., P.L., C.P., H.P. (Herwig Proske) and H.P. (Helene Petschko). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was conducted within the Integrating Land use Legacies in Landslide Risk Assessment to support Spatial Planning (ILLAS) project funded by the Austrian Climate Research Program (ACRP), gran<sup>t</sup> number KR16AC0K13226. The contribution of S. Gingrich was also financially supported by the European Research Council (ERC-2017-StG 757995 HEFT). We acknowledge support by the German Research Foundation and the Open Access Publication Fund of the Thueringer Universitaets- und Landesbibliothek Jena Projekt-Nr. 433052568.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The long-term spatially explicit information on land use 1830–1960– 2015 are available at https://doi.org/10.5281/zenodo.4896571, accessed on 30 June 2021. The landslide inventories and susceptibility models of both study areas are accessible through the Supplementary Materials.

**Acknowledgments:** We are grateful to the GIS department of the Styrian governmen<sup>t</sup> and the Provincial Government of Lower Austria for providing the high-resolution digital elevation models, orthophotos, and the datasets of hydrologic and hydropedologic characteristics. Furthermore, we thank the Geological Survey of Austria and the Styrian GIS department for providing the geological base maps.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Descriptive Summary of Input Data**


**Table A1.** Sources of land surface data.

**Table A2.** Overview of the classification and recording of land use classes according to time cut.


\* LULC types not listed here are not present in the study area; LULC classes of the 1960 time cut were digitized in aerial photographs, modified from Knevels et al. [45].



**Figure A1.** Biomass extraction [kg FW/ha/a] as an indicator of historical land use intensity of study areas.

#### **Appendix B. Summary of Model Assessment Results**

**Table A4.** Results of the model assessments within the study area with SpCV, and assessment of model transferability between both areas.


AUROC Statistic: median (*x* ), minimum (**Min**), maximum (**Max**), interquartile range (**IQR**), transferability estimate (**Transfer**).

**Table A5.** Wilcoxon signed-rank tests for AUROC performance differences within each study area estimated with SpCV.


**mAUROC:** median AUROC of SpCV; **N:** Number of observations per group; \* tied observations were removed; **Z**: Z score;alternative hypothesis: *greater*, *α* = 0.05; <: Differences significantly greater from previous model, =: No significant difference to the previous model.


**Table A6.** Variable importance measured as mean decrease in deviance explained (%), rank of variable in parentheses.

> \* Biomass extraction in interaction with LULC 2015; Wh: Waidhofen, P: Paldau.

**Figure A2.** Predictor-response relationships of landslide models in Waidhofen. Grey: 95% Bayesian credible interval of GAM-Base. Reference level of LULC 2015: 'Forest'; Lithological units: Reference: 'Upper Austroalpine limestone', 0: talus and glacial deposits, 1: Inneralpine Neogene, 2: Klippen zone, 3: flysch zone, 4: Upper Austroalpine marls. Note: the y axes are plot-dependent, and the x axes of non-parametric transformation functions are limited to the 5th to 95th percentile range.

**Figure A3.** Predictor-response relationships of landslide models in Paldau. Grey: 95% Bayesian credible interval of GAM-Base. Reference level of LULC 2015: 'Forest'; Lithological units: Reference: 'Neogene formations with coarse-grained layers', 0: 'Neogene formations dominated by fine-grained sediments', 1: 'pre-Würmian Pleistocene formations'. Note: the y axes are plot-dependent, and the x axes of non-parametric transformation functions are limited to the 5th to 95th percentile range.
