Understanding Landslide Susceptibility in Northern Chilean Patagonia: A Basin-Scale Study Using Machine Learning and Field Data

: The interaction of geological processes and climate changes has resulted in growing landslide activity that has impacted communities and ecosystems in northern Chilean Patagonia. On 17 December 2017, a catastrophic ﬂood of 7 × 10 6 m 3 almost destroyed Villa Santa Lucía and approximately 3 km of the southern highway (Route 7), the only land route in Chilean Patagonia that connects this vast region from north to south, exposing the vulnerability of the population and critical infrastructure to these natural hazards. The 2017 ﬂood produced a paradigm shift on the analysis scale to understand the danger to which communities and their infrastructure are exposed. Thus, in this study, we sought to evaluate the susceptibility of landslides in the Yelcho and Rio Frio basins, whose intersection represents the origin of this great ﬂood. For this, we used two approaches, (1) geospatial data in combination with machine learning methods using different training conﬁgurations and (2) a qualitative analysis of the landscape considering the geological and geomorphological conditions through ﬁeldwork. For statistical modeling, we used an inventory of landslides that occurred between 2008 and 2017 and a total of 17 predictive variables, which are geoenvironmental, climatological and environmental triggers derived from volcanic and seismic activity. Our results indicate that soil moisture signiﬁcantly impacted spatial susceptibility, followed by lithology, drainage density and seismic activity. Additionally, we observed that the inclusion of climatic predictors and environmental triggers increased the average performance score of the models by up to 3–5%. Based on our results, we believe that the wide distribution of volcanic–sedimentary rocks hydrothermally altered with zeolites in the western mountains of the Yelcho and Rio Frio basin are highly susceptible to generating large-scale landslides. Therefore, the town of Villa Santa Lucia and the “Carretera Austral” (Route 7) are susceptible to new landslides coming mainly from the western slope. This requires the timely implementation of measures to mitigate the impact on the population and critical infrastructure


Introduction
Earthquakes and volcanic eruptions have intensely disturbed the landscape of the Patagonian Andes. In fact, in the last 60 years, the three following significant earthquakes have occurred: Valdivia, 9.5 Mw [1]; Chiloé, 7.6 Mw [2]; and Aysén, 6.2 Mw [3]. The Patagonian Andes were also the epicenter of one of the most explosive volcanic eruptions of modern times. In 2008, the Chaitén volcano spread pyroclastic material over a large area on a continental level [4]. However, in recent years, the glacial and periglacial relief of the Patagonian Andes has also been changing at an accelerated pace due to climateinduced changes [5]. The interaction of these processes has generated a growing landslide activity that has impacted communities and ecosystems in northern Patagonia, Chile [3,[6][7][8].
In the Province of Palena, the National Geology and Mining Service (SERNAGEOMIN) identified a total of 2533 landslides [9]. Morales et al. [7] identified 722 landslides in the last 19 years. These geomorphological processes have demonstrated their destructive capacity. For example, on 16 December 2017, a landslide of volcanoclastic rocks generated debris and mudflows, destroying 3 km of the southern highway (Route 7) and 50% of the town of Villa Santa Lucia, causing the deaths of 22 people [6,9].
Given the impacts of landslides on mountain communities, it is essential to understand them spatially. Currently, predictive landslide models have incorporated a wide variety of statistical techniques, including machine learning algorithms that have been widely used in different geomorphological studies around the world [10][11][12] . Among these algorithms, classification methods based on logistic regression and support vector machine (SVM) stand out [11,12]. On the other hand, field observation data have played a fundamental role in providing geospatial information for this type of analysis. Field observations have been used both to create landslides databases after extreme precipitation events and/or earthquakes and to incorporate the different predictors needed for the statistical modeling of these natural hazards.
The following theoretical bases support geomorphological evolution linked to landslide activity. (a) In areas where landslides exist, future events may occur; (b) these areas share similar characteristics, such as geology, geomorphology, hydrology, hydrogeology and geotechnical characteristics and (c) a similar geo-environment generates similar patterns in the evolution of landslides [13]. In this sense, the susceptibility analysis is a widely used tool to understand the spatial probability of landslides based on certain conditioning factors [11,14]. Generally, susceptibility models are analyzed from a stationary perspective since landslide databases do not often have temporality due to information limitations. This makes it difficult to understand how different endogenous and exogenous elements of the evolution of the landscape interact. For this reason, in our study, we try to minimize this bias by considering a temporal database developed for the north of Patagonia in Chile during the period 2008-2017.
According to the area's background, landslides are recurring processes in the region with great destructive power; yet landslide susceptibility models have not been implemented to understand these dangers and mitigate the risk. The objective of our study was to evaluate the susceptibility of landslides in the Yelcho and Río Frío basins, the point of intersection where the alluvium that destroyed the town of VSL originated. We used two approaches, namely, (1) geospatial data in combination with machine learning methods, using different training configurations (see Sections 3.1-3.3); and (2) a qualitative analysis of the landscape considering the geological and geomorphological conditions (see Section 3.4). With these approaches, we attempt to uncover the variables which govern the different levels of susceptibility in this area of Patagonia and how they are distributed spatially. According to our field observations, how do geological conditions work together to generate these levels of susceptibility? Is our statistical approach sufficiently robust to capture the importance of the geological and geomorphological conditions of the landscape? Finally, we wish to provide information to technical institutions to better manage the danger of landslide processes in this portion of Chilean Patagonia.

Location
The study area includes the Yelcho and Rio Frio basin (43°05 S-43°48 S, 72°37 W-72°4 W) in a segment of the Andes mountain range in Patagonia, Chile ( Figure 1). Here, during the summer of 2017, a flood partially destroyed the town of Villa Santa Lucía (VSL) and part of Route 7 or the southern highway, the main overland connection route in Chilean Patagonia. The catastrophic result of this alluvium has made it one of the most important landslides in the historical record of southern Chile. The study area is a critical connection point (Route 7 and CH-235) for important cities in northern Patagonia. The history of the area and its characteristics reflect the importance of understanding the danger of landslides in these mountain basins, which, to date, has not been studied.

Geological and Geomorphological Setting
In the study area, the main geomorphological features are associated with the Liquiñe-Ofqui Fault System (LOFS), an intra-arc fault system parallel to the trench that extends for more than 1200 km along southern Chile, between 38°and 47°S in a NE direction [15][16][17]. This structural system aligns with an active volcanic arc where the Chaitén-Michinmahuida volcanic complex is located, whose last eruption strongly impacted southern Chile in 2008. The study area is characterized as having an irregular and rugged relief. The morphological features below are recognized.
The development of a valley of north-south tectonic origin [18] eroded by glaciofluvial action and filled by sedimentary deposits (Figure 1).
The Yelcho Lake basin, which is located in a NW direction covering more than 116 km 2 . The mountain ranges that make up these valleys constitute the second morphological feature, formed by intrusive rocks from the Norpatagonic Batholith, metamorphic rocks from the Principal Mountain Range metamorphic complex and volcanic rocks from the Pleistocene grouped in the Cordón-Yelcho volcanic complex [19][20][21] (Figure 1).
The third feature is glacial morphologies, such as cirques, troughs, hanging valleys and glacial lagoons from the Yelcho glacier, located to the west of the area [22,23].

Methodology
To analyze the susceptibility to landslides in the Yelcho and Rio Frio basins, we used geospatial data in combination with machine learning methods in an experimental setting. We also carried out a qualitative analysis of the landscape considering the geological and geomorphological conditions through fieldwork to understand the landscape elements that control landslides. Figure 2 shows the different methodological processes addressed in our work and which are explained in the following sections.

Landslide Database
For the training of the statistical models, we generated two landslide databases. Database with the presence of landslides: We used the landslide database developed by Morales et al. [7], which includes a total of 28 landslides in the study area from 2008 to 2017 (Figure 1). This database was created using temporal changes in the Normalized Difference Vegetation Index (NDVI) from Landsat images. To implement our analysis, we polygonized all the landslides located in the area of interest from the escarpment to the deposit area using QGIS 3.10. We relied on the mapping services Google Earth (January 2018) and Bing Maps (September 2013) using a mapping scale of 1:1500. The locations of the landslides were validated during two field surveys. Subsequently, to apply a grid-based classification method [11,24], all the pixels within the body of the landslide were sampled as slide pixels using a spatial resolution of 30 m, reaching a total of 6695 pixels. We used this method because of its widespread use in models of landslide susceptibility in various geographical areas of the world [24].
Database with an absence of landslides: We randomly chose 6695 pixels throughout the study area without landslides during the analysis period. We used this sampling mechanism for the negative class because it is one of the most used in the literature for the spatial prediction of landslides [11].

Climate
We used hydrometeorological variables that show a relationship with landslides activity according to previous work in the study area [7]. There is consensus that the erosive power depends on the regimes of precipitation and extreme events such as high 0 isotherm conditions with heavy rainfall, which have increased as a consequence of climate change, inducing important variations in the magnitude and frequency of landslides [25,26]. To represent the effects of climate and, mainly, precipitation, we used the total precipitation from the CHIRPS Pentad product, Climate Hazards Group InfraRed Precipitation with Station Data, with spatial resolution of 0.05°and temporal resolution of five days [27]. We also used total climate water deficit (CWD) and soil moisture from TerraClimate, Monthly Climate and Climatic Water Balance for Global Terrestrial Surfaces (University of Idaho), which has a spatial resolution of 0.04°and monthly temporal resolution [28]. We used the cloud processing tools provided by the climate engine [29] to calculate these climate variables from 2008 to 2017.

Geo-Environmental Information
We prepared a set of variables that represent the morphometric and hydrological characteristics of the study area (see Appendix A). We derived them from a 30 m spatial resolution Shuttle Radar Topography Mission (SRTM) DEM using the geoprocessing tools from QGIS 3.10. SRTM DEM has been widely recommended by recent research study [11,30,31]. The morphometric variables correspond to elevation, slope, aspect, curvature and TPI (topographic position index). We also calculated the excess of relief present in the study area using the algorithm provided by TopoToolbox [32]. The hydrological variables correspond to distance to drainage and drainage density. A set of geological variables was also considered, including distance to fault lines, fault line density and litology. The geological data (lithology and lineaments) were vectorized from the geological maps developed for the study area ( Figure 1) [9,20,21] using QGIS 3.10. The lithological classes were coded as numerical data by adding a field, to generate a single categorical variable subdivided into classes. Then, we transformed the variable into raster data using the information from the field contained in the initial vector file. A similar procedure was developed for the guidelines, for the subsequent estimation of density and distance using the processing tools of QGIS 3.10 .

Environmental Triggers
According to the historical record, major earthquakes have affected the area, such as Valdivia Mw 9.5 (1960), the largest recorded in history and Chiloe Mw 7.6 (2016). Earthquakes generate significant disturbances in the landscape, causing stress on rocks that can last for decades [33]. Our work represented this stress in the landscape considering the historical accumulation of PGA (peak ground acceleration) and PGV (peak ground velocity) between 1900 and 2017. For this, we used the seismic events that impacted the study area from the USGS earthquake catalog (https://earthquake.usgs.gov/earthquakes, accessed on 15 October 2021), reaching a total of twelve events. The variables PGA and PGV correspond to isolines that spatially represent a seismic event that occurred in a certain geographic area. To incorporate this information, we transformed the numeric field of the original vector file containing this information to a raster-type file to later estimate the PGA and PGV accumulated between the years 1900 and 2017 (see Appendix A). Volcanic eruptions can also generate critical environmental triggers, mainly affecting sedimentation patterns and influencing landslide activity. The study area was less than 30 km from the Chaitén volcano, which generated one of the most significant recorded eruptions in the Patagonian Andes. In this work, we used the isopacs that represent the thickness (in centimeters) of the pyroclastic material associated with the eruption of the Chaitén volcano from [4]. To incorporate this information, we vectorized the thicknesses associated with the isopaches and then transformed the numerical information field of these isolines to a raster file format (see Appendix A).

Multicollinearity Analysis
Before applying the statistical method, we ruled out the duplication of highly correlated variables (Table 1). We used a variance inflation factor (VIF) multicollinearity analysis. Based on what was suggested [34,35], the variables that met a VIF of less than five were used as predictors.
The following equation corresponds to the formulation of the VIF analysis: where R 2 i corresponds to the coefficient of determination for the ith independent variable, corresponding to a correlation measure of the independent variables [34].

Machine Learning Models and Parameter Estimation
In this work, we implemented seven different configurations (M) to analyze susceptibility sensitivity considering the presence of different estimators. M1-all estimators; M2-all except earthquakes; M3-all except weather; M4-all except eruptions; M5-all except weather-eruptions; M6-all except weather-earthquakes; M7-all except earthquakeseruptions; M8-all except weather-earthquakes-eruptions.
For statistical modeling, we used two machine learning algorithms implemented in the scikit-learn library that runs on the Python programming language (https://scikitlearn.org, accessed on 1 October 2021). The first algorithm used was logistic regression (LR); this statistical method has become one of the most used techniques in the evaluation of landslide susceptibility [36], since it allows the regression values to be expressed in terms of probability [37]. In the case of the model LRy i , the probability that the dependent variable y i (presence or absence of landslide) is equal to 1 is given by where n corresponds to the number of predictors used (X k ), β k are the coefficients of the function and β 0 is the intercept. For the implementation of the model, we used the hyperparameter adjuster GridSearchCV provided by the scikit-learn library with 10 crossvalidations. GridSearchCV generates a deep search of specified parameter values for a model. This search is optimized by cross-validation. The parameters adjusted by Grid-SearchCV correspond to Solver, which consists of different types of solvers tasked with finding the weights of the parameters that minimize the cost function. Penalty is a regularization method implemented to avoid overfitting, for which a penalty term is added, which, in this case, corresponds to the Ridge or L2 type. Tolerance is the convergence threshold value of the optimizer; when the improvement between iterations is less than the specified threshold, the algorithm stops returning the current model. The last adjusted parameter is Iterations, which is the maximum number of iterations carried out so that the solvers converge. The parameters of the model that fit the best are shown in Table 2. The second model used corresponds to the support vector machine (SVM). This machine learning algorithm maps the data set in a high-dimensional space, where it creates an optimal hyperplane to separate the characteristics or classes [38]. The algorithm finds the hyperplane when the separation margins between the classes are maximum; these margins are called support vectors [39]. Due to its high levels of performance, the SVM method has been widely used in landslide susceptibility analyses [40]. SVM classifications are mainly affected by the type of kernel function selected. There are different types, such as linear, polynomial, sigmoid and radial basis function (RBF) [39]. In our study, we applied the linear kernel since this allowed us to obtain the importance of the estimators for the prediction of the susceptibility probabilities. Similar to the previous model, we used the GridSearchCV hyperparameter adjuster with 10 cross-validations. The best performing regularization parameter C is shown in Table 2. This parameter (C) controls the balance between a decision limit and the correct classification of the training points, determining the model's fit to the training data.

Model Validation
To validate the models, we included a set of recommended performance metrics [41,42] that considered pixels classified as positive that exist (TP ), positive cases wrongly classified as negative (FN), cases wrongly classified as positive (FP) and negative cases correctly predicted as negative (TN).
Accuracy is the proportion in which the model can correctly classify all the positive and negative samples.
Precision is the proportion of identified landslides correctly classified by the model.
Recall is the proportion of landslides correctly detected by the model.
The F 1 -Score is the harmonic average weighted by Precision and Recall.
The ROC AUC represents the relationship between the rate of false positives and true positives.
True positive rate = TP TP + FN False positive rate = FP FP + TN Finally, we defined a Score that corresponded to the average of all the previously described metrics to determine each model's general performance.

Fieldwork
One of the largest landslides in the historical record of southern Chile occurred in the study area, the alluvium that partially destroyed the town of Villa Santa Lucia (VSL) in the summer of 2017. Considering this background, we carried out a field survey to evaluate the different susceptibility levels predicted by machine learning models in the area surrounding VSL. We analyzed the model results from a critical perspective, mainly considering the geological and geomorphological aspects. VSL is located in a strategic area for the connectivity of Chile. The valley covers the only land connection route through the national territory from north to south in Chilean Patagonia via Route 7.
Fieldwork was divided into a prospective survey, and a geological data collection and landslide characterization survey. We carried out a general description of the areas affected by landslides, considering the geomorphological context, the structural geology and a deep macroscopic analysis of the lithology involved in the scarp areas.

Importance of Predictors
We used 17 predictors to model the spatial susceptibility of the landslides that occurred between 2008 and 2017. According to the logistic regression model, when we used all the variables (LR1), the most predictive factors were soil moisture, geology, PGV, drainage density and curvature (ordered in decreasing order of importance) (Figure 3). In the case of the SVM model (Figure 4), the variables soil moisture, geology, drainage density, PGA and curvature were positively related to the susceptibility to landslide and showed greater importance in the prediction. By contrast, the remaining variables in the SVM model negatively influenced the probability of landslides. We observed that, in general, the variables thickness of pyroclastic deposits, TP, distance to drainage, distance to lineaments and elevation showed the lowest prediction coefficients for both the LR and SVM models. In this sense, it can be stated that the order of importance of the variables (without considering the value of the coefficient) is quite similar in the two models. We observed that, when climatic variables were not included, PGA became the most critical predictor in both models (LR3, LR5, SVM3 and SVM5). In general, the weight of the topographic, geological and hydrological variables did not present significant changes among the different configurations, which may be an indicator of the relatively constant importance of these variables under the experimental conditions developed in the present work.

Models' Performance and Landslide Susceptibility
We implemented eight configurations with both machine learning models. In general, we noted that, when the climatic variables were absent (M3, M5, M6 and M8), the average prediction capacity (score) in both models decreased by up to 3-5%. When we only used the terrain variables, the score decreased by 3% in both models; even so, they reached acceptable prediction conditions. In all the configurations implemented, the SVM model performed the best, except for the M5 model (Table 3).
Although good yields characterized the prediction results, the configurations used in M1 and M2 showed the best score in both models. In the case of the LR model, the M1 configuration was superior to the rest, with a score of 0.83 (Table 3). While the SVM model showed no differences in the M1 and M2 configurations, in the latter, we did not include seismic data. Table 3. Performance of the machine learning models for each configuration. Considering the M1 statistical model, we observed that SVM and RL showed similar patterns in the spatial distribution of susceptibility ( Figure 5). The results suggest that the probability of landslides is mainly restricted to valley areas with a high drainage density. As elevation increased, the probability of landslides began to be distributed mainly in morphologies of glacial origin, such as hanging valleys or glacial cirques. The north part of the study area (Yelcho basin) stood out with a high probability of landslides, concentrated mainly in the central valley. The same was observed in the headwaters of the Río Frio basin to the southwest of the area, where areas of high susceptibility also stood out.

Metrics
In the western mountains that delimit the Rio Frio and Yelcho basin, the danger was mainly restricted to the areas of the valleys connected to the Yelcho glacier; this includes areas of discharge from lagoons and in front of the glacier, whose drainage networks are connected to the central valley. Consequently, for example, both the town of Villa Santa Lucia and the Carretera Austral were found to be in an area of great danger. On the other hand, considering the temporal range analyzed, the eastern mountains of the cold river basin did not present greater danger. High susceptibility was restricted to the NE of the study area to the valleys that discharge into Lake Yelcho and the steep slopes surrounding the lake.

Field Observations
During the field survey, we analyzed the area around Villa Santa Lucia that corresponds to the intersection of the Rio Frio basin that drains to the south and the Yelcho basin that drains in a northerly direction. The area visited included a section of the central valley separated by two mountain ranges in an NS direction. The Yelcho Glacier covers the western range; it presents hanging glacial valleys that provide areas of sediment accumulation from landslides that occurred on the walls of the range (Figures 6b and 7). In addition, there is an extensive foothill made up of reworked moraine, alluvial, glaciolacustrine and volcanic sedimentary deposits (Figure 7c). In the northern zone, metamorphic and intrusive rocks emerge, constituting the western range's basement rocks. In this zone, the mountainous range is crowned by volcanoclastic rocks of the Cordón-Yelcho volcanic complex, a unit that is both covered by and juxtaposed with the Yelcho Glacier (Figure 7a). In the lower part of this mountain range, at the height of Cuesta Moraga and for approximately 4 km, a sedimentary sequence composed of quartz sandstones, mudstones and shales is partially exposed, in addition to a volcanic deposit (Figure 5a). In these units, multiple landslides and debris flows occur with areas of influence of a few meters, directly affecting the southern highway (Route 7) (Figure 6a). In addition, in the northwestern mountain range, the tuffs of the Cordón-Yelcho volcanic complex involve a deep landslide associated with the Villa Santa Lucía alluvium (Figure 6c). This part of the tuff formation comprises 500 m thick levels deposited on intrusive rocks (Figure 6c). The base of the landslide escarpment is located in front of one of the tongues of the Yelcho Glacier, where part of the sliding material and moraine deposits of the same lithology remain. In the upper part of the escarpment, there is a marked subvertical fracture that currently remains unstable (Figure 6b). The volcanoclastic rocks (tuffs) from the escarpment show an intense hydrothermal alteration; in this sense, the presence of zeolite fillings of 1-2 cm that are readily identifiable in hand samples drew our attention. Along the path of the flow through the bed of the Burritos River, we observed the intercalation of sedimentary deposits and wetland areas, which are also visible in satellite images (Figure 7e).
The eastern mountain range is composed of intrusive rocks (Figure 6f,g). It is characterized by a sharp relief that prevents the accumulation of material. At the foot of this mountainous range, there are narrow valleys covered by alluvial deposits of a few meters (Figure 7d). On the relief of steep slopes, we observed the development of a layer of soil of a few centimeters, which favors surface landslides particularly when vegetation is present. Northeast of Villa Santa Lucia, along route 235 (Figure 6g), landslides and flows occur from both slopes that enclose the valley directly, involving route 235 in its affected area. The occurrence of these landslides has, as a common factor, high angle slopes and heavily fractured rocks (Figure 7b). In general, we observed that the pebbles are angular and of variable size; they present flat ridged faces with epidote skates, characteristics that are consistent with the strongly jointed rocky walls. Similar patterns of landslides and intrusive rock flows were noted towards the south of Villa Santa Lucia. We identified that, in both the eastern mountain range and in the southwestern mountain range (Figure 6e), rock landslides are shallow and primarily associated with soil and vegetation cover movements, being restricted to secondary drainage networks.

Data Limitation and Approach
The methodology implemented in the present study allowed us to develop a susceptibility model for the period 2008-2017 using a set of predictor variables which considered the relief, climatology and environmental triggers such as earthquakes and volcanic eruptions. The climatic products used had resolution of 3.4-4 km. Better resolution should be able to improve the predictive capacity of the model's results, particularly the differences in the spatial distribution of susceptibility, as was shown recently in an experimental investigation [43].
The landslide inventory considered 28 events in the 2008-2017 period, which were limited by the spatial resolution of the satellite images from which they were derived. The existence of events not identified during this period would generate changes in the prediction levels of the configurations used and, therefore, in the results presented here. We stress that there is a much higher number of landslides in the study area without associated temporality that are not included in our study. The models implemented have good prediction capacity; however, better results could be obtained using non-linear machine learning models such as decision trees, assembly models and neural networks [12]. Susceptibility models are sensitive to class sampling techniques. In order to apply the grid-based classification method, in our work, we used all the pixels within the body of the slide; this sampling method is widely used for its simplicity [11,24]. In this sense, we believe it is significant to mention that applying other sampling configurations could potentially induce changes in yield and spatial susceptibility [11,44].

Prediction Capacity of the Variables
Some of the most important variables in the implemented models were the geological variables considered. As observed during the field survey, the type of rock is very important for generating landslides. For example, the volcanoclastic rocks of the western mountains of the Yelcho and Burritos basin resulted in an alluvium of 7 million m 3 in 2017. We believe that the zeolites identified in the tuffs of the Cordón-Yelcho volcanic complex in the western mountain range increase susceptibility to deep landslides. In this respect, various investigations carried out on rocks with zeolite have indicated that hydrothermal alteration with this mineral affects the stability of the rocks under humid conditions [45,46]. At present, the escarpment where this mega-slide originated still shows signs of gravitational instability, which suggests that an event of large magnitude cannot be ruled out. Although we believe that the weight assigned by the models to this predictor is consistent, our results do not agree with those previously obtained by [7] in northern Patagonia using different machine learning models. Differences in the analysis approach explain these contrasting results. While [7] considered the point-scale temporality for each year of analysis to study the spatial-temporal behavior of the landslides, we took all the variables as static elements with the idea of obtaining the spatial probability of landslides in a much more limited area. Therefore, there are significant differences produced by different learning conditions of the models, size of the study area and data set. In this sense, although the database used has time associated with it, we did not develop a time-dependent analysis, so our approach may underestimate the variability of susceptibility through time and space, as suggested by recent research [31,47].
According to our field observations, the structural geology of the area also impacts the development of landslides; however, in the implemented models, it did not stand out as a significant predictor. During fieldwork, we observed a wide development of landslides in areas weakened by geological faults, consistent with previous studies in the study area, which makes us think that the models we implemented here did not capture the importance of these elements. The importance of the drainage density does not surprise us; during the field survey, we observed that the landslides originated in areas of high slopes close to perennial and seasonal river networks that can transport soil-rock eroded by landslides to valley areas.
The set of PGA/PGV predictors, associated with the seismic activity that affects the study area, also positively influenced the prediction of susceptibility. The role of earthquakes in the evolution of the landscape is well known. Chile, located in one of the areas with the highest seismic activity on the planet, has suffered significant earthquakes in the past that affected this segment of Patagonia. Despite the critical role of seismic activity in the geomorphological evolution of the Andes, its effects on mountains have not been studied in depth considering different time scales. However, from studies carried out in other mountain ranges, we know that there is a dual role of earthquakes in the activity of landslides, as a trigger mechanism and as an agent of stress with long-term effects on the landscape [33,48], due to damage to rocks and the expansion and propagation of cracks that accelerate the weathering processes and subsequent erosion [33].
The pyroclastic deposits derived from the recent eruption of the Chaitén volcano (2008) dispersed in the study area show a low prediction capacity for the models implemented in our work. This does not agree with recent research that has shown the importance of this factor in the activity of landslides [7,8]; however, we believe that our results are a consequence of the distance of the basins studied from the eruptive center. Therefore, in areas closer to the volcano, the prediction capacity of this variable should increase.
Concerning the ability of the climate variables to predict, soil moisture was the most relevant. This result was observed in all the model configurations where this variable was used. This is consistent with research that has analyzed the role of soil moisture in landslide activity, which has even suggested the implementation of this parameter in early warning systems [49,50]. Although the remaining climatic variables did not stand out for their prediction capacity, our results show that, when they were not used, it reduced the predictive capacity (configuration M3, Table R3), exceeding the models that did not consider environmental triggers. We believe that the generation of landslide inventories with the date of the event could help better understand the role of dynamic predictors (environmental triggers and climatology) in the susceptibility distribution in complex terrains such as Patagonia.

Interpretation of the Results
According to the susceptibility model, we can highlight the following interpretation from our results.
There is a high susceptibility in the western mountain range, mainly associated with chain erosion processes due to the topographic gradient generated by hanging valleys and well-developed drainage networks. There is also a high susceptibility on the Villa Santa Lucía plain and the surrounding western area; the map shows that high susceptibility is associated with eventual flows from the sources of the Burritos and Orlando rivers. In this area, glaciolacustrine deposits and wetlands are essential sediments that can be transported by flows coming from this direction. In the western mountains of the Rio Frio basin, there are four glacial lagoons confined to slopes with high susceptibility. This suggests that there is a danger associated with flooding due to the overflow of glacial lagoons. Such processes are a growing danger in the Patagonian Andes associated with the retreat of mountain glaciers increased by climate change [51]. For this reason, we suggest evaluating the danger associated with these processes in this segment of the Andes through competent methodologies, mainly because they have caused havoc to the south of the study area [23,52,53]. Based on our results, we recommend that the authorities carry out engineering works to mitigate disasters associated with floods in the study area and promote the development of an early warning system for landslides in the Patagonian Andes. We also propose this study as a territorial planning tool, in order to avoid population settlement in high-risk areas.
The results for the eastern mountain range in the Río Frío basin correlate with what was observed in the field survey. The map shows that the landslide susceptibility of the eastern mountain range does not represent a danger for Villa Santa Lucía ( Figure 5). We believe that this occurs because the elevation of the land (approximately 500 m a.s.l.) between the eastern mountain range and the town (approximately 230 m a.s.l.) acts as a natural barrier that can help mitigate the effects of landslides from that direction (BB' profile in Figure 6). On the other hand, on the slopes of Lake Yelcho, there are areas with high-moderate susceptibility ( Figure 5); therefore, it is crucial to evaluate the susceptibility in all the walls that enclose the lake basin and to assess the danger of tsunamis induced by landslides that could impact infrastructures and nearby locations. In southern Chile, earthquake-induced landslides have produced damming of lakes (Valdivia earthquake, 1960) [54] and even tsunamis (Aysén earthquake, 2007) that cost human lives [3].
Field surveys are of utmost relevance to validate the results of the models. They can help identify possible prediction errors inconsistent with the geological and geomorphological context of the studied area and complement the strengthening of the analysis.

Conclusions
In this work, we analyze the susceptibility to landslides in the Río Frío and Rio Yelcho basins in northern Patagonia in Chile during the period 2008-2017, integrating machine learning methods and field observations. For statistical modeling, we used 17 predictive variables associated with topographic, hydrological, geological, climatic parameters and environmental triggers such as earthquakes and volcanic eruptions. Our results indicate that soil moisture significantly impacts susceptibility, followed by geology, drainage density and seismic activity. Additionally, we observed that the climatic and environmental disturbance predictors increased the predictive capacity of the models by up to 3-5% in the average performance score.
Based on our results, we believe that the wide distribution of volcanic sedimentary rocks, with zeolites in the western mountains of the Yelcho and Rio Frio basins, is highly susceptible to generating large-scale landslides. The eroded material can travel through drainage networks that open over an extensive foothill to the central valley. Our field results show an essential control of the structural geology that contributes to the weakening of the rocks due to the predisposition of landslide planes; however, this condition was not captured by the statistical models implemented.
As a result of the high level of susceptibility in the western cord, we recommend monitoring the slope displacement of volcanic-sedimentary rocks partially covered by the Yelcho glacier. Additionally, due to the instability conditions presented by the Villa Santa Lucia alluvial escarpment, we do not rule out that an event similar to that of 2017 could occur again at any time. For this reason, we recommend that the authorities carry out engineering works to mitigate disasters associated with floods in the study area. In addition, we believe that this study can be used by decision makers as a territorial planning tool, in order to avoid population settlement in high-risk areas.
Finally, we believe that the scientific community should advance in the development of an early warning system in the Patagonian Andes. In this sense, the generation of comprehensive landslide inventories would be essential to assessing the variability of the threat caused by landslides at different temporal and spatial scales in a constantly evolving relief.