Towards the Use of Land Use Legacies in Landslide Modeling: Current Challenges and Future Perspectives in an Austrian Case Study

Abstract

Land use/land cover (LULC) changes may alter the risk of landslide occurrence. While LULC has often been considered as a static factor representing present-day LULC, historical LULC dynamics have recently begun to attract more attention. The study objective was to assess the effect of LULC legacies of nearly 200 years on landslide susceptibility models in two Austrian municipalities (Waidhofen an der Ybbs and Paldau). We mapped three cuts of LULC patterns from historical cartographic documents in addition to remote-sensing products. Agricultural archival sources were explored to provide also a predictor on cumulative biomass extraction as an indicator of historical land use intensity. We use historical landslide inventories derived from high-resolution digital terrain models (HRDTM) generated using airborne light detection and ranging (LiDAR), which are reported to have a biased landslide distribution on present-day forested areas and agricultural land. We asked (i) if long-term LULC legacies are important and reliable predictors and (ii) if possible inventory biases may be mitigated by LULC legacies. For the assessment of the LULC legacy effect on landslide occurrences, we used generalized additive models (GAM) within a suitable modeling framework considering various settings of LULC as predictor, and evaluated the effect with well-established diagnostic tools. For both municipalities, we identified a high density of landslides on present-day forested areas, confirming the reported drawbacks. With the use of LULC legacy as an additional predictor, it was not only possible to account for this bias, but also to improve model performances.

1. Introduction

Landslides pose a threat to human lives and infrastructure. A changing climate and land use/land cover (LULC) alter the landslide risk and thus have societal consequences [1,2]. In Austria, landslides are relevant natural hazards preconditioned by factors such as lithology, geomorphology, tectonic structures and LULC, and are mainly triggered by long-lasting heavy rainfall and rapid snowmelt [3]. Therefore, understanding the factors that increase the chance of landslide occurrence is crucial for spatial planning in the face of ongoing and expected future climate and LULC changes.

LULC types, and their changes, are reported to have different hydrological and geomechanical effects controlling slope stability [4]. While forest is often considered to stabilize slopes [5], forest harvesting or road construction undercutting slopes may reduce slope stability [6,7]. However, in landslide studies, LULC has often been considered as a static factor representing solely the present-day LULC [8] (i.e., latest available LULC). Recently, more studies have begun to account for historical LULC in landslide analysis, assigning LULC dynamics an important role in explaining landslides [9,10,11,12]. Beguería [9] and Persichillo et al. [11] discovered a high landslide susceptibility on abandoned cultivated land, even after revegetation by shrubs or trees in the Spanish Pyrenees and in the Oltrepò Pavese (Italy), respectively. Gariano et al. [10] and Pisano et al. [12] found evidence that land management reduced landslide occurrences in Southern Italy (Calabria and Molise), supporting the importance of LULC changes in spatial planning practice. However, due to the availability of mainly remote-sensing products (aerial or satellite imagery), these landslide analyses were only able to consider historical LULC since the mid-20th century (e.g., since 1954 in [12] or since 1957 in [9]). To the authors’ knowledge, only one study used historical cartographic documents as additional sources (e.g., Napoleonic cadastral map) to assess the long-term legacy effects of LULC on mass movements [13].

In general, legacy effects describe the influence of past events or processes on later states, often spanning decades to centuries [14,15]. Long-term legacy effects of past LULC have been shown to exist in the context of socio-ecological dimensions such as contemporary forest structure, management and disturbance risk [16] or biodiversity [17,18]. For mass movements, Lopez-Saez et al. [13] revealed the potential of historical LULC changes in explaining the paradoxical observation of reduced rockfall hazards despite an increased urban exposure in the Grenoble conurbation since 1850. Especially the forest densification at the upper part of the slope was considered to contribute to the identified decrease in rockfall frequency and energy for volumes up to 5 m³ [13].

For landslide analysis and modeling, landslide inventories are a fundamental source to improve the understanding of the factors that precondition and trigger landslides. In the last decade, low-cost, airborne LiDAR-based high-resolution digital terrain models (HRDTM) became available area-wide for all federal provinces of Austria (≤10 m × 10 m) [19]. Many studies demonstrated the potential of LiDAR HRDTM and its derivatives to identify landslides, and thus to substantially improve conventionally created landslide inventories [20,21,22], especially underneath the forest cover, where passive remote-sensing sensors are of limited utility [23].

Generally, landslide inventories have an unknown level of incompleteness [24] and thus may be biased. Inventory biases have previously been studied for remotely-sensed, event-based or archival inventories, and can be caused by focusing solely on administrative boundaries, damage reports or on a single triggering event [25,26]. While most authors reported the usefulness of LiDAR-derived historical landslide inventories, only few analyzed drawbacks of the data source [26,27]. For example, it is very challenging or even impossible to determine a landslide’s exact extent, absolute age, trigger, and potential for reactivation when using only HRDTM derivatives [28,29]. Additionally, according to Petschko et al. [27], forest cover may have a conservation effect on the landslide morphology (i.e., “young” morphology of a “very old” landslide), while landslides may be underrepresented on agricultural land and near infrastructure due to land rehabilitation (i.e., “very young” morphology but no visibility on orthophoto) resulting in a landslide inventory that is substantially biased towards a high landslide density in recently forested areas. Analyzing the effect of systematically incomplete landslide inventories, Steger et al. [26] discovered that landslide susceptibility models emphasized bias-describing predictors (e.g., larger regression coefficients), and as a consequence the bias was directly propagated into the landslide predictions. While Petschko et al. [27] recommended to exclude inventory-biasing observations (i.e., “old” and “very old” landslides) or to drop the bias-describing predictor (i.e., present-day land cover) from modeling, Steger et al. [26] included the bias-describing predictor (e.g., forested area) as a random effect and used only the fixed effects to make model predictions.

In this study, we investigated the association between LULC legacies and landslide susceptibility using an airborne LiDAR-derived historical landslide inventory in two municipalities in Austria (Waidhofen an der Ybbs and Paldau). We addressed the following main questions: (i) are long-term LULC legacies important and reliable predictors of landslide susceptibility? And (ii) can LULC legacies help to understand and account for possible inventory biases in modeling present-day landslide susceptibility? Additionally, we analyzed the transferability of landslide models between study areas and the effect of dropping inventory-biasing observations.

For the analysis, we digitized and classified LULC patterns for three-time cuts comprising nearly 200 years using various spatial data sources. In addition, yields and livestock statistics were compiled from archival sources and statistical publications, and summarized as socio-ecological variables reflecting plot-level LULC legacies. For the assessment, we used generalized additive models (GAM) within a modeling framework considering different combinations of LULC legacy implementations while also accounting for land surface variables (e.g., slope angle, etc.) and lithological conditions as possible confounders. We evaluated the effect of LULC legacies using well-established diagnostic tools for model assessment and interpretation.

2. Materials and Methods

2.1. Study Area

The study was conducted in two municipalities in Austria: Waidhofen an der Ybbs (referred to as Waidhofen) in Lower Austria, and Paldau in Styria. The two municipalities represent different landscapes (Figure 1A,C).

Waidhofen is located in the Ybbstaler Alps from 14°39′ E 47°52′ N to 14°56′ E 48°01′ N, covers an area of 131 km² and has a population of about 13,000 inhabitants [30]. Its elevation rises towards the south with a relative relief of 54–651 m/km². In the limestone-dominated south, up to 1205 m above the Adriatic (m AA) are reached; towards the flysch zone in the north the mountains transition into gentle hills (302 m AA) [31]. In contrast, Paldau lies in the East Styrian Basin with a relative relief of 6–138 m/km² (282 to 465 m AA), characterized by unconsolidated sediments of the Neogene to Quaternary period [32]. The municipality is an agriculturally favorable region mainly with corn and pig farming. Paldau has a population of 3000 inhabitants [33], it extends from 15°43′ E 46°54′ N to 15°51′ E 46°59′ N and covers 39 km².

The geological setting coupled with the characteristic very local, frequent and intense rain events in summertime create conditions that make both study areas particularly prone to landslides of different types and magnitudes [34,35,36]. In the last decades, landslide occurrences have caused substantial damage to settlements and infrastructure in both study areas [37,38].

2.2. Data

This study is based on land surface and landslide data from various sources and at different spatial resolutions. As our target resolution, we used 10 m × 10 m to account for the dependence of landslide susceptibility on local-scale topography.

2.2.1. Land Surface Data

For both study areas, airborne LiDAR-based HRDTMs of 1 m × 1 m resolution were available, which were provided by the GIS department of the Styrian government for Paldau (acquisition year: 2009), and by the provincial government of Lower Austria for Waidhofen (acquisition year: 2017), respectively (

Table A1. Sources of land surface data.

Study Area	Source Holder	Resolution
airborne LiDAR-based high-resolution digital terrain model
Waidhofen	provincial government of Lower Austria	1 m × 1 m, acquisition year: 2014
Paldau	GIS department of the Styrian government	1 m × 1 m, acquisition year: 2009
hydrologic and hydropedologic parameters
Waidhofen	Austrian Research Centre for Forests	50 m × 50 m, year: 2014
Paldau		100 m × 100 m, year 2017
geological basemaps
Waidhofen	Geological Survey of Austria	1:50,000
Paldau	GIS department of the Styrian government	1:50,000

For downscaling to the target resolution of 10 m × 10 m, bilinear interpolation was applied.

).

Datasets of hydrologic and hydropedologic characteristics were compiled by the Austrian Research Centre for Forests in 2014 for Waidhofen [39] (50 m × 50 m resolution) and in 2017 for Paldau [40] (100 m × 100 m resolution), respectively, and were made available by the respective federal states. The datasets include information on soil parameters such as soil type, total pore volume, and hydraulic conductivity of the topsoil (0–20 cm).

Geological basemaps at a 1:50,000 scale were provided by the Geological Survey of Austria for Waidhofen and by the Styrian GIS department for Paldau, respectively. In both geological basemaps, alluvial deposits were corrected in order to match valley floors visible in the HRDTM. Furthermore, we reclassified the geological units into a smaller number of relevant classes. In Waidhofen the reclassification was based on lithology and geomechanical properties which resulted in seven lithological units of (i) alluvial deposits, (ii) talus and glacial deposits, (iii) Inneralpine Neogene, (iv) Klippen zone, (v) flysch zone, (vi) Upper Austroalpine marls, and (vii, reference level) Upper Austroalpine limestone. In Paldau, we specified five geological units based on grain size distribution and age of origin: (i, reference level) Neogene formations with coarse-grained layers, (ii) Neogene formations dominated by fine-grained sediments, (iii) pre-Würmian Pleistocene formations, (iv) Würmian and Holocene sediments, and (v) other units.

2.2.2. LULC Legacy

By combining unique historical spatially explicit information on LULC in the two study areas and numerical information on land use intensity, we were able to create information on LULC legacies that could be used as an input in landslide susceptibility modeling. This is, to our knowledge, the first study that uses such a long, multi-temporal and spatially explicit LULC record to better understand regional-scale landslide susceptibility.

To generate this dataset, we collected, digitized and harmonized data depicting different spatial patterns for the years 1820, 1960 and 2015 (

Table 1. Sources of land use legacy information of at least acceptable spatial and thematic quality identified for the case study areas.

Year	Study Area	Source	Source Holder	Q	Q-Explanation
Land use and land cover
1820	Wh	Maps of Franciscan Cadastre	Provincial Archive of Lower Austria	++	Sharp delimitation of utilization unit
	P		Federal Office for Calibration and Measurement	++
1962	Wh	Aerial images	Federal Office for Calibration and Measurement	~	Differentiation based on greyscale aerial photography
1953	P			~
2015	Wh, P	Orthophotos & IACS *	Open Data Austria	++	Parcel-sharp delimitation of arable land and grassland
Agricultural yields (cereal and grassland)
1820	Wh	Text records of Franciscan Cadastre	Provincial Archive of Lower Austria	+	based on two cadastral municipalities of Waidhofen
1820	P	Sandgruber [42]		~	average of Styria
1960	Wh, P	Agricultural statistics	Statistics Austria Library	++	data on municipality level
2015	Wh, P	IACS *	Open Data Austria	+	data of farms in municipality
Wood yields
1820	Wh	Text records of Franciscan Cadastre	Provincial Archive of Lower Austria	+	based on two cadastral municipalities of Waidhofen
1820	P	Gingrich et al. [43]		~	average of Styria
1965	Wh, P	Weiss et al. [44]		~	Austrian average
2015	Wh, P	Forest inventory	Federal Forest Office	~	state averages

Wh: Waidhofen, P: Paldau; * Data status of temporal extent of 2016; Quality (Q): ++: high spatial and thematic quality, +: high spatial or thematic quality, ~: acceptable spatial and thematic quality. Adapted from Knevels et al. [45].

). Our data sources included the Franciscan Cadastre of 1820, aerial photographs of 1960, and aerial orthophotos combined with Integrated Administration and Control System (IACS) data of 2015 (“present-day”). Numerical information on LULC intensity (e.g., agricultural yields and machinery numbers) was collected at communal, district or provincial levels from archival sources and both historical and recent statistics, and was up- or downscaled to match the municipal level. The data sources assigned to a time cut slightly differed in their year of origin between the municipalities (Table 1). The selection of the time cuts was based on characteristic socio-economic framework conditions of the Central European cultural landscape following Bender [41] (i.e., time cuts around 1850, 1914, 1960 and present age) as well as on the availability and temporal proximity of data.

We created a vector GIS database with LULC polygons from these sources using ESRI ArcGIS. We distinguished between forest area, cropland and grassland, and classified the remaining LULC as “settlement and other” (for details refer to

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

land Use/Land Cover Category	Time Cut 1820	Time Cut 2015
	LULC Types in the Franciscan Cadastre *	IACS, Orthophotos
forest (including forest pasture)	Hardwood forests, Coniferous forests, Mixed forests, Chestnut forests, Meadows with fruit trees	All forest types digitized from orthophotos
grassland	Dry meadows, wet meadows, pastures, community pastures, shrubs	IACS agricultural parcel: Grassland, alpine pastures, pasture
cropland	Orchards, vegetable gardens, vineyards, arable land (with fruit trees, trees and vines)	IACS agricultural parcel: arable land
settlement and other	Marshes, lakes, ponds, rivers and streams, wastelands and bare rocks, buildings (all types), trails (all types)	Remaining area, which includes, e.g., buildings, impervious surfaces, water bodies, excavation pits and quarries, urban green, near-natural areas

* 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].

in the Appendix A). Due to varying data quality (Q, Table 1), in the first step, the cadastral maps (1820) and orthophotos (2015) were georeferenced and digitized (EPSG: 32633) since these allowed us to map boundaries between distinct land uses. In a second step, we used these layers as base maps for digitizing the less distinct areas from greyscale aerial photographs (1960). Spatially persistent natural and built structures were kept unchanged for all time cuts. Furthermore, we consistently used a display scale of 1:1000 to keep the digitization error constant. The estimated overall positional error is up to 3–5 m for 1820, 5–10 m for 1960, and less than 3 m for 2015.

Based on the numerical information, we quantified biomass extraction as yields (kg fresh weight per ha and year, kg FW/ha/a) at the lowest possible administrative level and allocated it to cropland, grassland and forests. Biomass extraction represents the output from the agricultural production system and is thus an indicator of the LULC intensity during historical time periods [46]. For the biomass extraction from cropland, we employed the cereal yields as proxy; data on harvest from meadows was collected for grassland yields, and wood extraction was derived from forest yields.

We finally intersected the three vector datasets for 1820, 1960 and 2015 successively to one final file, keeping all available attributes, and converted the result to raster format at the uniform target resolution. Biomass extraction was summed over the time cuts to obtain a cumulative land use intensity (Figure A1 in Appendix A). This is an innovative approach for this socio-ecological indicator in the context of landslide science as other authors calculated cumulative materials flows [47] or greenhouse gas emissions [48] in other application contexts. To avoid artifacts due to geometric inaccuracies inherent in the vector files and from sliver polygons, isolated pixels were identified, and the affected grid cells were excluded from the landslide sampling design described below (see Section 2.2.3). The estimated digitization errors (3–10 m) might partly be counteracted by using a target resolution of 10 m × 10 m (i.e., resolution corresponds to largest estimated error).

The created historical LULC legacy dataset was made publicly available [49].

2.2.3. Airborne LiDAR-Derived Landslide Inventories

For both study areas, we derived historical landslide inventories by mapping landslides visible in the HRDTM following the approach proposed by Schulz [22]. We included landslides in earth and debris materials, focusing on the slide type with possible transitions to complex slide flows according to the classification scheme by Cruden and Varnes [50]. Landslides were digitized as polygon features also distinguishing between landslide body and scarp (Figure 1B). Selected mapped landslides were inspected in the field for validation, and we corrected the inventory where necessary. In total, in Waidhofen, 621 landslides were mapped, covering 5.31% of the municipal area, and in Paldau, 418 landslides (4.14% of the area; Figure 1C;

Table A3. Summary of the landslide inventory data.

		Waidhofen			Paldau
Number	landslides		621			418
	scarps		829			469
	bodies		663			348
	samples		974			559
Total Area [m2] (%)	landslides	6,976,638 (5.31)			1,621,250 (4.14)
		min	mean	max	min	mean	max
Area [m2]	landslides	113	11,235	1,163,088	30	3879	206,842
	scarps	2	517	79,640	12	140	1518
	bodies	52	9876	1,083,448	24	3949	206,842
Perimeter [m]	landslides	44	451	7250	26	163	1821
	scarps	7	101	3013	18	71	373
	bodies	29	297	4548	28	168	1821

in the Appendix A).

Following Hussin et al. [51], we randomly sampled landslide presence points in the landslide scarp area using the recommended 50-m distance constraint, and attributed the point to the corresponding grid cell in our target grid (in total 974 and 559 landslide points in Waidhofen and Paldau, respectively). For sampling non-landslide points, we defined the landslide-free area by excluding the mapped landslides and a surrounding 50-m buffer to account for digitization errors. We furthermore excluded so-called trivial areas—areas considered as not susceptible to landsliding (e.g., floodplains, flat areas) [52]. Isolated grid cells (see Section 2.2.2.) and anthropogenic structures with similar geomorphometric characteristics as landslides (e.g., quarries) were also masked. For the landslide absence locations, we distributed random points in a 1:1 sampling ratio using a minimum nearest-neighbor distance of 50 m to reduce spatial autocorrelation.

2.3. Methods

The relationships between LULC legacies and landslide distribution were analyzed using GAMs [53,54] while also accounting for the local topography as an important preparatory factor for landslides [8]. GAMs have become popular in landslide susceptibility studies due to their ability to model nonlinear relationships while allowing for a separate interpretation of additive effects in terms of odds ratios and variable importance [55].

Our analysis was conducted in the free and open source computing environment R (R version 3.5.3) [56]. We used the GAM implementation of the mgcv package [54] and the mlr package [57] as the modeling framework. Furthermore, for terrain analysis we used System for Automated Geoscientific Analysis (SAGA) GIS 6.3.0 [58] through RSAGA [59] and Terrain Analysis Using Digital Elevation Models (TauDEM) 5.3 [60] via R system calls.

For downscaling the input data to the target resolution, we applied bilinear interpolation for resampling using SAGA GIS. However, we acknowledge that we are unable to capture local-scale patterns of geology or soil parameters.

2.3.1. Landslide Susceptibility Modeling Design

For landslide susceptibility modeling, we related land surface variables, soil parameters, lithological units and LULC legacies as predictors to landslide occurrences (

Table 2. Predictor variables for landslide susceptibility modeling.

Variable(s)	Software	Setting	Method
land surface variable
convergence index (100 m, 500 m)	SAGA GIS	r = 100 m, 500 m	[61]
curvature (plan, profile)	SAGA GIS		[62]
flow accumulation, D-Infinity	TauDEM	log-transformed	[63]
normalized height	SAGA GIS	w = 5; t = 2; e = 2	[64]
slope angle	SAGA GIS		[62]
slope angle, catchment area	SAGA GIS		[65]
slope aspect (S-N, W-E)	SAGA GIS	cosine, sine transformed	[62,66]
topographic position index (TPI)	SAGA GIS	r = 500 m	[67]
topographic wetness index (SWI)	SAGA GIS		[65]
soil
total pore volume		up to 20 cm depth, median
hydraulic conductivity		up to 20 cm depth, median
lithology
geology		* ref: Waidhofen (vii), Paldau (i)
land use/land cover legacy
LULC 2015		ref: ‘Forest’
biomass extraction (1820, 1960)		sum since 1820 and 1960

Setting, scale-dependent parameters: r: radius; w,t,e: Parameters in SAGA GIS module Relative Heights and Slope Positions; * ref = reference level (see Section 2.2.1).

for overview).

Our model design enabled us to explore relationships between LULC legacy and landslide distribution, and to improve the understanding of the potential biases in airborne LiDAR-derived landslide inventories. We created landslide susceptibility models with different sets of input variables: (i) The baseline model ‘GAM-Base’ excludes LULC legacy variables; for the assessment of the LULC legacy effect, we built (ii) a GAM using the present-day LULC as an additional predictor (GAM-2015), (iii) a GAM based on the LULC legacy information from 1960 to 2015 (GAM-1960), and (iv) a GAM based on the LULC legacy information since 1820 (GAM-1820). Moreover, we tested (v) a GAM using the setting of GAM-2015, but excluding potentially inventory-biasing observations following the recommendations of Petschko et al. [27] (GAM-2015-Masked; i.e., all observations located in continuously forest-covered areas). Furthermore, we allowed modeled landslide occurrences to be dependent on the combined effect of the historical biomass extraction and present-day LULC class rather than modeling LULC legacy variables as additive terms. Thus, the predictor of LULC legacy information in GAM-1820 and GAM-1960 was implemented as a parametric, linear interaction term between the LULC of 2015 and the historical biomass extraction (sum since 1820 or 1960). Moreover, we tested the model’s transferability between study areas, but excluded predictors that are specific to each area (i.e., lithology).

For both study areas, the airborne LiDAR-derived historical landslide inventory and the created landslide susceptibility models are available as Supplementary Materials.

2.3.2. Assessment of the Effect of Land Use Legacy

The empirical effect of LULC legacy on landslide occurrence was assessed in terms of model performance, variable importance, odds ratios (OR) of the LULC classes, and transformation function plots of the three most important predictor-response relationships.

For the model assessment, we applied a k-fold spatial cross-validation (SpCV) to achieve independent test areas and thus a bias-reduced predictive performance as a measure of model generalization [55]. For SpCV, we partitioned the data into five disjoint folds using k-means clustering of the coordinates (k = 5), and repeated this procedure 100 times. In each repetition, subsequently, four folds were used as training data while the remaining fold was used for validation until each fold was used once for validation (i.e., k models per repetition, 500 models in total). Furthermore, we ensured comparability between the models’ performance estimates using identical training and validation data for each study area and repetition.

The area under the receiver operating characteristic curve (AUROC) was computed as the performance measure. The AUROC is a quality measure suitable for binary response data and is a common evaluation tool for landslide susceptibility models [55]. AUROC values lie between 0.5 (no discrimination) and 1.0 (perfect discrimination), and were interpreted following the recommendations of Hosmer et al. [68]. Additionally, differences in model performances were tested using Wilcoxon signed-rank tests ($\alpha$ = 0.05; R coin package) [69,70], and p-values were adjusted for multiple comparisons according to Benjamini and Hochberg [71].

As a measure of variable importance, we extracted for each variable the mean decrease in deviance explained (mDD, %) under the consideration of all SpCV models. The mDD indicates the explanatory contribution of a variable to the overall explained deviance of the corresponding model [72,73]. The higher the mDD value, the greater is the contribution of a variable, and thus its importance. To compute the mDD, we left the variable of interest out while fixing the remaining smoothing parameters (mgcv::gam sp argument) during model (re-)training, and subsequently measured the percentage differences of the deviance explained [73].

Transformation function and OR plots were used to explore the relationships between landslide distribution and the three most important predictors as well as LULC legacy of each model setting. Additionally, we extracted comparative predictor-response relationships reported in studies of the same area [73,74]. A transformation function plot shows the predictor-response relationship as a parametric (linear) or non-parametric smoothing function by using the additive structure of a GAM. We visualized predictor-response relationships on the logit scale (i.e., linear predictor scale). An OR indicates the chance that an outcome occurs given a specific exposure, relative to a reference exposure [75]. An OR less than one means an exposure with lower odds of the outcome while an OR higher than one shows an exposure with higher odds of the outcome, while accounting for the other variables in the model; an OR of one is associated with no influence of the exposure [75]. We calculated ORs for the LULC classes with ‘forest’ as reference level. Additionally, we derived ratios of ORs (rOR) by dividing a model’s predictor-response relationships by the corresponding relationship in GAM-Base. rOR enables a more sensitive comparison between models.

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).

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. Results of the model assessments within the study area with SpCV, and assessment of model transferability between both areas.

Model	$\tilde{\mathit{x}}$	Min	Max	IQR	Transfer
A: Waidhofen
GAM-Base	0.79	0.7	0.91	0.14	0.65
GAM-2015	0.79	0.69	0.91	0.14	0.64
GAM-2015-Masked	0.8	0.68	0.91	0.13	0.69
GAM-1960	0.79	0.7	0.91	0.14	0.65
GAM-1820	0.78	0.7	0.91	0.14	0.65
B: Paldau
GAM-Base	0.88	0.83	0.93	0.03	0.83
GAM-2015	0.89	0.85	0.94	0.06	0.8
GAM-2015-Masked	0.93	0.88	0.98	0.04	0.85
GAM-1960	0.89	0.85	0.94	0.06	0.78
GAM-1820	0.89	0.85	0.94	0.06	0.79
C: Combined
GAM-Base	0.81	0.74	0.91	0.07
GAM-2015	0.81	0.74	0.91	0.07
GAM-2015-Masked	0.84	0.75	0.93	0.07
GAM-1960	0.82	0.74	0.91	0.07
GAM-1820	0.81	0.74	0.91	0.07

AUROC Statistic: median ($\tilde{x}$), minimum (Min), maximum (Max), interquartile range (IQR), transferability estimate (Transfer).

in the Appendix B).

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 yet 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. Wilcoxon signed-rank tests for AUROC performance differences within each study area estimated with SpCV.

Model	mAUROC	N *	Z	p-Values	$\mathit{r}$
A: Waidhofen
GAM-2015	0.79
<GAM-1820	0.78	36	2.7	0.01	0.45
=GAM-Base	0.79	35	0.93	0.18	0.16
<GAM-1960	0.79	36	4.24	<0.001	0.72
<GAM-2015-Masked	0.80	36	2.24	0.02	0.37
B: Paldau
GAM-Base	0.88
<GAM-1960	0.89	66	5.01	<0.001	0.62
<GAM-1820	0.89	66	4.35	<0.001	0.54
=GAM-2015	0.89	66	1.54	0.06	0.19
< GAM-2015-Masked	0.93	66	7.03	<0.001	0.87

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

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. Variable importance measured as mean decrease in deviance explained (%), rank of variable in parentheses.

Variable	Study Area	GAM-Base	GAM-2015	GAM-2015- Masked	GAM-1960	GAM-1820
land surface variable
convergence index, 100 m	Wh	1.23 (5)	1.27 (5)	1.3 (6)	1.15 (5)	1.17 (5)
	P	0.91 (7)	0.97 (8)	0.71 (11)	0.96 (8)	0.98 (7)
convergence index, 500 m	Wh	0.87 (8)	0.86 (9)	0.99 (9)	0.83 (10)	0.86 (9)
	P	0.93 (6)	1.16 (6)	3.25 (4)	1.2 (6)	1.18 (6)
curvature, plan	Wh	2.09 (3)	2.1 (3)	2.83 (3)	2.17 (3)	2.17 (3)
	P	2.38 (4)	1.43 (5)	0.84 (10)	1.46 (5)	1.49 (5)
curvature, profile	Wh	1.45 (4)	1.58 (4)	1.78 (4)	1.65 (4)	1.57 (4)
	P	3.32 (3)	3.2 (3)	6.37 (1)	3.08 (3)	3.05 (3)
flow accumulation	Wh	0.42 (13)	0.43 (12)	0.48 (14)	0.5 (13)	0.46 (12)
	P	0.07 (15)	0.2 (16)	0.46 (13)	0.22 (16)	0.21 (16)
normalized height	Wh	0.02 (15)	0.02 (16)	0.05 (16)	0.04 (16)	0.03 (16)
	P	0.95 (5)	0.95 (9)	3.15 (5)	0.93 (9)	0.93 (9)
slope angle	Wh	8.01 (1)	7.7 (1)	8.63 (1)	7.45 (1)	7.61 (1)
	P	7.68 (1)	4.38 (1)	4.62 (2)	4.77 (1)	4.64 (1)
slope angle, catchment area	Wh	1.11 (7)	1.15 (6)	1.1 (8)	1.07 (7)	1.15 (6)
	P	0.55 (11)	0.47 (12)	0.2 (15)	0.44 (12)	0.43 (12)
slope aspect, S-N	Wh	1.12 (6)	1.07 (7)	1.11 (7)	1.06 (8)	1.12 (7)
	P	3.94 (2)	2.3 (4)	1.25 (8)	2.36 (4)	2.33 (4)
slope aspect, W-E	Wh	0.61 (11)	0.58 (11)	0.84 (10)	0.52 (12)	0.56 (11)
	P	0.32 (13)	0.29 (14)	1.68 (6)	0.27 (15)	0.27 (15)
TPI	Wh	0.86 (9)	0.89 (8)	1.64 (5)	0.89 (9)	0.91 (8)
	P	0.7 (9)	0.84 (11)	0.14 (16)	0.84 (11)	0.84 (11)
SWI	Wh	0.66 (10)	0.66 (10)	0.39 (15)	0.68 (11)	0.67 (10)
	P	0.66 (10)	0.39 (13)	0.95 (9)	0.35 (13)	0.35 (13)
soil
total pore volume	Wh	0.35 (14)	0.31 (14)	0.51 (13)	0.29 (15)	0.31 (15)
	P	0.47 (12)	0.86 (10)	0.68 (12)	0.87 (10)	0.87 (10)
hydraulic conductivity	Wh	0.44 (12)	0.41 (13)	0.62 (11)	0.35 (14)	0.38 (14)
	P	0.9 (8)	0.98 (7)	1.3 (7)	0.98 (7)	0.97 (8)
lithology
lithology/geology	Wh	6.38 (2)	6.37 (2)	3.72 (2)	6.37 (2)	6.34 (2)
	P	0.25 (14)	0.27 (15)	0.39 (14)	0.27 (14)	0.28 (14)
land use/land cover legacy
LULC 2015	Wh		0.26 (15)	0.6 (12)
	P		3.64 (2)	3.92 (3)
biomass extraction, 1960 *	Wh				1.07 (6)
	P				3.81 (2)
biomass extraction, 1820 *	Wh					0.43 (13)
	P					3.85 (2)

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

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%).

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 agreement 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).

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⁻¹ (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).

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 management 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 suggest 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 forest-covered 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), yet the OR of the LULC categories in reference to forest were similar to GAM-2015. However, since the relationships are in general agreement 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; Figure A2 and Figure 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 agreement 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 yet 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 suggest 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.