remote Distinct Susceptibility Patterns of Active and Relict Landslides Reveal Distinct Triggers: A Case in Northwestern Turkey

: To understand the factors that make certain areas especially prone to landslides, statistical approaches are typically used. The interpretation of statistical results in areas characterised by com-plex geological and geomorphological patterns can be challenging, and this makes the understanding of the causes of landslides more difﬁcult. In some cases, landslide inventories report information on the state of activity of landslides, adding a temporal dimension that can be beneﬁcial in the analysis. Here, we used an inventory covering a portion of Northwestern Turkey to demonstrate that active and relict landslides (that is, landslides that occurred in the past and are now stabilised) could be related to different triggers. To do so, we built two landslide susceptibility models and observed that the spatial patterns of susceptibility were completely distinct. We found that these patterns were correlated with speciﬁc controlling factors, suggesting that active landslides are regulated by current rainfalls while relict landslides may represent a signature of past earthquakes on the landscape. The importance of this result resides in that we obtained it with a purely data-driven approach, and this was possible because the active/relict landslide classiﬁcation in the inventory was accurate.


Introduction
Data-driven models can be thought as empirical tools that extract functional relationships from past phenomena to estimate the expected behaviour of the same phenomena in a pre-defined (or ill-defined) future. This framework is commonly referred to as Hutton's uniformitarian principle, and is more commonly translated as the past is the key to the future [1][2][3]. Hutton first and subsequently Lyell helped to develop and spread the concept of uniformitarianism, replacing the then prevailing idea of catastrophism. Since then, this concept has formed the backbone of any landslide susceptibility study [3].
Landslide susceptibility models (LSM) can be used to predict the spatial occurrence of future landslides by assuming, consistently with the uniformitarian principle, that in any given area, slope failures will occur under the same circumstances and because of the same conditions that caused them in the past. However, this principle may not always hold true [4]. First-failure landslides and reactivations may have different controls, acting both on their triggers and kinematics: think of the peak and residual shear strengths, or the role of strong earthquakes as opposed to aftershocks or rainfall [5]. Changes in material properties also are reflected in morphological changes which, in turn, affect the process dynamics [6].
Uniformitarianism still represents a fundamental component of the literature [7][8][9][10][11][12], albeit the current climate change has led us to question its present-day validity [13]. In fact, data-driven models are generally built upon the effects of past events, which may not inform on the slope response to events of different nature in the future. This topic is hardly addressed in the literature as it requires an understanding not only of the evolution of the triggers, but also of the predisposing or preparatory conditions. These may behave differently in time [4], reflecting changes in mechanical [14], hydrological and hydraulic [15], or thermal conditions [13,16].
Moreover, the complexity of a seismically active area does not help perform a straightforward analysis and interpretation of both the predisposing and triggering factors. This is true, for instance, in Turkey (Figure 1a), where the huge variability and complexity of the territory makes the modelling difficult to tackle. The capability and feature of an existing landslide inventory [17], able to discern relict (termed "inactive" in the inventory) from active landslides, comes to our aid as it could enable the distinction of different processes in place and their drivers. Moreover, especially for active landslides, the definition of the main triggering factors is crucial for the evaluation of economic costs related with the frequency and magnitude of disaster events [18].
Seismically active areas may contain many landslide bodies. The seismic shaking likely triggers their first movement, while subsequent remobilisations become increasingly related to different triggers and predisposing factors as time goes by [16,19]. Different triggers can also produce distinct patterns in space [20]. In Northwestern Turkey, which we took as our study area, we can discern two sub-areas: the North Anatolian fault region in the southwest, characterised by a higher density of relict bodies, and the region close to the Black Sea in the northeast, richer in active ones (Figure 1b,c). These sub-areas display an attitude of surface processes to be related to distinct triggers: seismicity and rainfall (Figure 1d,e). However, these differences may be difficult to discern in an inventory in which landslide types or activity states or stages are not classified.
The effect of biases in susceptibility modelling has been explored in the literature [21][22][23]. The necessity to operate with an unbiased area [24,25] led us to focus on a specific sector in Northwester Turkey, rich in landslides but not too tectonically complex, and sufficiently geologically and geomorphologically homogeneous.
In terms of modelling approaches, the literature offers many options. We opted for the Generalised Additive Model (GAM; [15,26]), which can explain the spatial distribution of landslides via a family of Bernoulli exponential functions, in which the influence of the covariates can be captured via linear and nonlinear relationships. As such, the approach allows us to display the uncertainties in the estimations, which are intrinsically part of a Bayesian framework [27]. This statistical implementation is utilised here for the first time to investigate relationships between two distinct models covering the same area but differing by a categorical entity (active/inactive landslide). Furthermore, we decided to use Slope Units (SUs) as they are geomorphologically-consistent subdivisions that can be linked with landslide processes and are thus preferrable to grid-based subdivisions [28][29][30][31].

Study Area
The geology and geomorphology of Turkey is unique and extremely complex, owing to both past and ongoing processes in place. Various studies exist [32][33][34][35][36], in which the national settings are dissected per geological history and geomorphological processes. Figure 1 displays the large-scale geomorphological and geological features.
The diversity of morphologies derives from a geodynamic environment that is still quite active and determines a variety in outcropping lithologies [35]. Three main landslidedominated landscapes are recognised, corresponding to the tectonostratigraphically-distinct Western, Central, and Eastern Pontides [37]. The lithological units forming the Pontides vary along the belt, featuring west-east-oriented sub-parallel bands of sedimentary, metamorphic, and igneous rocks. While the western portion is richer in Triassic to Paleogene sedimentary and medium-grade metamorphic outcrops, the central zone comprises Eocene volcaniclastic and sedimentary rocks and Palaeozoic metamorphic rocks, and the eastern zone features Paleogene and Cretaceous plutonic and igneous formations underlying Eocene and Neogene sedimentary and volcanic formations [38].
The portion chosen as our study area ( Figure 1) corresponds to the Zonguldak quadrangle in the Western Pontides. Here, the North Anatolian Fault System (NAFS) produced a landform dominated by mountain belts and plateaus. The NAFS is an over 1600 km long, right-lateral strike-slip, active transform fault running along Northern Anatolia in the E-W direction, that also separates the study area in a southern and a northern sector (Figure 1a).
Climatically, the area belongs to the Black Sea climatic region in the north and the continental inner Anatolian climatic region in the southeast [39]. The former receives rainfall throughout the year (>1000 mm mean annual precipitation, up to 2300 mm in its eastern portion) [39]. The north-facing slopes of the coastal mountain belt are comparatively wetter as they intercept the weather fronts, and this is reflected by a relative abundance of landslides ( Figure 1e). Precipitation decreases southward, where the Palaeocene-Eocene flysch and Palaeocene-Middle Miocene volcanics are the most landslide-prone units [39].

Mapping Units
We used SUs as terrain partitions characterised by similar hydrological and geomorphological conditions [42]. Each SU has a distinct shape given by the interplay between lithotypes and morphometries, and thus offers morphological and lithological characteristics that can be analysed statistically. An SU-based subdivision is not the only possible choice. In fact, most contributions in the literature opt for a regular lattice or pixel-based subdivisions [3]. These, however, even though they can be expressed at a fine to very-fine resolution, do not reflect any natural characteristics. Conversely, SUs can, better than pixels, represent geomorphological processes (e.g., [43,44]) and, at the same time, reduce the computational burden, especially in models covering large areas (in our case, SU-based calculations are~100 faster than pixel-based ones).

Landslides
The landslide inventory was compiled by the General Directorate of Mineral Research and Explorations of Turkey (MTA) and published in 1:25,000 scale for the entire national territory (http://yerbilimleri.mta.gov.tr/, accessed on 1 February 2022). It is a general polygon-based-inventory that carries some complexities because the mapped phenomena are related to a plurality of causes. Thus, the identification of unstable areas should benefit from the use of as many variables as possible to discern the main factors affecting slope stability.
The catalogue [39], within our study area, comprises 4084 active landslides and 1140 inactive landslides. The former are those that were moving at the time of mapping (1997)(1998)(1999)(2000)(2001)(2002)(2003), while the latter could be classified as relict according to the UNESCO Working Party on World Landslide Inventory (1993) [48]. Qualitatively, we notice an abundance of active phenomena in the northeastern area where the climatic influence of the Black Sea is stronger. A less dense but significant presence of landslides can be seen in the southwestern sector ( Figure 2).
We extracted the highest point in the landslide polygon to better represent the source material and/or lithology [49,50]. Subsequently, separately for each model, we reported the count of active or inactive landslides within each SU, attributing the presence of landslides to those SUs containing at least one point. Finally, we identified 2822 SUs (X% of SUs) with active phenomena, and 983 (Y%) with inactive phenomena.

Modelling Strategy
In a binomial GAM, the data (y) are assumed to be conditionally independent given the linear predictor η: where p i is the binomial probability. Here, we assume N i = 1 for all i because we have binary data. The η i as a function of p i is called the link function, and we describe it using a logit, but we note that other link functions are possible. The linear predictor η is where we put the additive model: where β j are the fixed (or linear) effects, with weak priors, describing the linear relationship of the covariates x j . Each f represents a random (or nonlinear) effect with and τ is a constant. For the two f, we use a spline model, also referred to as Random Walk of the first order (RW1; [51]) A RW1 induces adjacent class dependence among mean slope and precipitation bins, respectively [52]. The whole implementation makes use of the INLA framework [46,47].
To quantitatively compare the inventories of active and inactive landslides, we relied on the effects of the selected covariates. Active and inactive phenomena should present distinct correlations with physical variables owing to changes in predisposing and triggering conditions as well as their values. For instance, active and inactive phenomena could be characterised by distinct distributions of slope angles or vegetation coverage. However, while morphological and climatic factors can change rather rapidly, geological, lithological, and structural factors are not expected to vary over human time scales (in absence of catastrophic events).
Possible multicollinearity issues among covariates ( Figure S1) were eliminated by discarding those showing more than 0.75 collinearity with another covariate [53][54][55]. The final list of covariates is reported in Table 1. We preferred covariates that are well known in the literature [3,56], and analysed their linear effects in most cases. We investigated the nonlinear effect of slope and precipitation to better capture the role of gravitativehydrological processes in landsliding.
For model fitting, we used the whole sets of active and inactive landslides (separately). For validation, we used a tenfold Cross Validation (CV) with mutually exclusive subsets, implying that no SUs are repeated across CV replicates, and thus, there is no autocorrelation [57]. We used the Area Under the Receiving Operating Characteristic Curve (AUC) and the confusion matrix to evaluate the model performance. This is not the only possibility. In fact, new articles suggest a spatial CV with connected packages [58][59][60]. However, this solution is still under discussion [61], which is why we preferred to use a pre-consolidated methodology.

Results
We report the outcomes of the so-called fitting (within sample) and cross-validation (out of sample) procedures. The former is used to interpret the patterns of the explanatory variables, while the latter is a tool for model validation. We then produce two susceptibility maps (that is, two spatial probability maps), one for the active landslides, and one for the relict landslides. These maps are compared with those of a number of explanatory variables, to seek common patterns. It should be stressed that the susceptibility maps are intended as descriptions of past/present phenomena and not as a temporal prediction tool.

Distinct Patterns of Explanatory Variables
The main tool that we can use to evaluate the extent to which the classification into inactive and active landslides relates to distinct controlling factors is the analysis of the effects of these factors, represented by sets of covariates in the statistical model. We believe, in fact, that inactive and active landslide phenomena should be spatially correlated with physical variables in a distinct way. Predisposing factors that cause landslides must be congruent in inactive as well as in active landslides. Hence, what should stand out the most is that the predisposing factors can change over time.
The posterior marginal distributions of the linear effects of each covariate in the models constructed with inactive and, separately, active landslides are displayed in Figure 3. Notably, about half of the covariates exhibit distinct effects in the two models, suggesting that the phenomena featured in the two classes may be controlled by different processes.
Geomorphologically, reasonable patterns are described. The role of Northness is consistent with the distribution of precipitation, which comes from the Black Sea, north of the study area (Figure 1). Negative values of RSPµ are observed in active landslides. Seemingly counterintuitively, PGAµ is positively correlated with inactive landslides. However, this is consistent with the observation [39] that these landslides may be related to historical earthquakes in the NAFS. The elongation of the SUs has a negative effect on both active and inactive landslides, as elongated slopes offer less room for large, deep-seated landslides to form. The average slope negatively correlates with the presence of active landslides, while no correlation is found with inactive ones. For active landslides, this suggests that phenomena on very steep slopes are unlikely; in fact, unstable bodies are quickly removed from steep slopes, while landslide inventories tend to better capture movements over more gentle slopes, which can remain active for a longer time. Hydrological covariates (TWIσ and SPIµ) seem to exert a negative effect in both active and inactive landslides, suggesting their preferential occurrence in the upper portions of catchments. The effect of lithology ( Figure S2, Table S1) also shows some differences between the models. For instance, granitoid areas are negatively correlated with active landslides, while the opposite holds true for carbonate rocks.

Distinct Landslide Triggers
The behaviours of the nonlinear effects of covariates for inactive and active phenomena also are different, as shown in Figure 4. These nonlinear effects clearly point to distinct triggers for inactive and active landslides. In fact, for inactive landslides, neither slope nor precipitation exert significant effects (95% confidence level). For active landslides, positive effects are seen within a certain range of slope angles (10-20 • ) that could be related to specific materials capable of sustaining prolonged landslides, and for large amounts of annual precipitation (>800 mm/year), capable of frequently triggering or sustaining a variety of movements. On the contrary, in areas with less precipitation, the effect is slightly negative. In order to validate the result, the out-of-sample performance of the model is investigated. This is done in two steps: the first one involves the use of AUC for each of the considered landslide classes, whereas the second one maintains the same structure but focuses on the summary metrics of confusion matrices. Figure 5a shows the ROC curves and their AUC values for ten cross-validations for the active and inactive landslide models. The AUC values (~0.8) can be deemed satisfactory and are consistent across the replicates, indicating robustness of the model [66]. The figure (bottom row) also shows the confusion plots of the two models, which are rather similar in both the high ability to detect true positive cases (~90%) and the lower ability to identify the true negative ones (47-50%). Consequently, the error rate (bottom right) is also of similar magnitude (44-50%), and it seems to fail on the stable conditions, together with the ratio between Predicted True Negatives and Observed Negatives.
However, we should remind that the SUs had been calculated only for areas with slope topography, as we excluded the flat areas that are obviously not susceptible to gravitative movements. This could be interpreted as a weakness of the model, but actually facilitates its ability in recognising the instable areas.

Distinct Susceptibility Maps
In Figure 6, we show the resulting susceptibility maps for the inactive and active landslides. The maps display markedly different spatial patterns that are consistent with the qualitative (Figure 1) and quantitative (Figures 3 and 4) observations that the distribution of active landslides better correlates with the annual precipitation and that the distribution of inactive landslides better correlates with the peak ground acceleration.
What is more, the two maps are not "one the negative of the other". The patterns that emerge are, in fact, distinct. Indeed, we do not observe a specular negative effect of precipitation in the distribution of inactive landslides (in fact, we do not observe a significant effect at all). However, we do see a negative effect of the peak ground acceleration in the distribution of active landslides, but the magnitude of this effect is smaller than that seen for inactive landslides.
The absence of correlation between the two maps is demonstrated quantitatively in Figure 6d, where the Pearson correlation shown is 0.5, indicating a random dependence. Similarities in the two maps mainly exist in areas with low density of landslides, independently of their state of activity, such as the northwestern and southernmost portions.

Controls and Fate of Active Landslides
From our analysis, it emerged that, consistently with the definition of relict landslides (used to classify the inactive landslides), the conditions that caused their occurrence in the past are distinct from those that are responsible for the active movements in the present. Moreover, the locations of the inactive landslides point to areas with high seismicity, suggesting that they may be earthquake-induced phenomena, now stabilised and insensitive to hydrometeorological forcing. Conversely, the distribution of active landslides reflects the pattern of present-time annual precipitation, suggesting the rainfall-induced nature of these phenomena. If we interpret the active landslides as slow-moving mass movements (that is, processes that remain active for a comparatively long time, and thus more easily captured by inventories), the correlation between the rainfall pattern and the spatial distribution of landslides makes very much sense. There is a large literature showing that the mobilisation and acceleration of slow-moving landslides along slopes is mostly governed by rainfall events causing an increase in pore water pressure and thus a reduction in the available shear strength [67]. On the other hand, in the absence of significant shifts in hydro-meteorological patterns, the stabilisation of these landslides should mostly be related to their transition from steeper to gentler slopes and/or to plains or valleys without significant fluvial erosion at the toe). In other words, active landslides could be described as meta-stable hillslopes materials experiencing creep and consolidation processes while the ratio between driving and resisting forces fluctuates over time mostly under the control of hydro-meteorological factors.
Here, we should stress that slow-moving landslides could rapidly turn into catastrophic landslides (and thus rapidly stabilise) if at some point the driving forces dramatically exceed their resisting counterparts. Various factors including seismicity or precipitation itself could trigger catastrophic landslides. Yet, this may seem a more likely scenario for a region exposed to intense seismic external forces rather than precipitation because, overall, even relatively low-intensity ground shaking may be more destructive than intense precipitation at triggering landslides [68]. In this context, it is not surprising that the relict landslides are mostly distributed closer to the North Anatolian Fault zone, whereas the active landslides concentrate far from it.

Accuracy of the Active/Inactive Landslide Classification
The binary classification into active and inactive landslides in the inventory was performed well. This is demonstrated by the fact that it resulted in the production of two distinct and uncorrelated susceptibility maps. In other words, in addition to suggesting differences between the conditions responsible for landslides in the past and in the present in the study area (supporting a non-uniformitarian view in this highly dynamic context), this result also suggests that we are dealing with a well-done classification. Logically, active and relict landslides should not be difficult to discern, but the point here is that, if a bias existed in this classifier, it would have resulted in less distinct (and thus more spatially correlated) susceptibility maps. Seeing this the other way round, if a classifier is expected to define distinct regions of space and this does not occur, the severity of the classification bias could be quantified from the degree of correlation between the maps generated, separately, for the distinct values of the classifier.

Conclusions
The analyses presented in this work aimed at investigating differences in the spatial patterns of relict and active landslides in a landslide-rich geomorphological context. These differences, expected in the light of qualitative observations on the possible landslide triggers and predisposing factors, were expressed quantitatively using a purely data-driven approach, confirming the validity of such methodology, suggested in the literature [69], and the accuracy of the classification operated in the inventory. The result that the susceptibility patterns of relict and active landslides in the study area are spatially distinct and correlate with distinct explanatory variables suggests that, while current rainfall patterns may explain the distribution of active landslides, seismicity may have had an impact on the relict landslides.
Overall, we believe our work can represent a summary of good practices in the definition of landslide susceptibility mapping and hopefully serve as a reference standardised assessments in both common and specific applications. It also brings novelty as it presents a general slope unit-based susceptibility model through a Bayesian approach in a study area, namely the Turkish Northwesternmost Sector, so far not investigated with this technique.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14061321/s1. Figure S1: Collinearity between the variables used in the model; Figure S2: Fixed effects expressed as marginal distributions for each lithological unit for landslide type; Table S1: List of lithologies used in the model.  Data Availability Statement: Elaborated data are presented in the manuscript. Raw experimental data can be provided by the authors upon reasonable request.

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.