Improving landslide inventories by combining satellite interferometry and landscape analysis: the case of Sierra Nevada (Southern Spain)

An updated and complete landslide inventory is the starting point for an appropriate hazard assessment. This paper presents an improvement for landslide mapping by integrating data from two well-consolidated techniques: Differential Synthetic Aperture Radar (DInSAR) and Landscape Analysis through the normalised channel steepness index (ksn). The southwestern sector of the Sierra Nevada mountain range (Southern Spain) was selected as the case study. We first propose the double normalised steepness (ksnn) index, derived from the ksn index, to remove the active tectonics signal. The obtained ksnn anomalies (or knickzones) along rivers and the unstable ground areas from the DInSAR analysis rapidly highlighted the slopes of interest. Thus, we provided a new inventory of 28 landslides that implies an increase in the area affected by landslides compared with the previous mapping: 33.5% in the present study vs. 14.5% in the Spanish Land Movements Database. The two main typologies of identified landslides are Deep-Seated Gravitational Slope Deformations (DGSDs) and rockslides, with the prevalence of large DGSDs in Sierra Nevada being first revealed in this work. We also demonstrate that the combination of DInSAR and Landscape Analysis could overcome the limitations of each method for landslide detection. They also supported us in dealing with difficulties in recognising this type of landslides due to their poorly defined boundaries, a homogeneous lithology and the imprint of glacial and periglacial processes. Finally, a preliminary hazard perspective of these landslides was outlined.


Introduction
Landslides represent one of the main natural hazards with a strong socioeconomic impact on a global scale (e.g. Kirschbaum et al. 2015;Froude and Petley 2018;Mateos et al. 2020). A good-quality landslide inventory map is necessary for assessing landslide hazard (van Westen et al. 2008). There are some global databases that are actively maintained, such as the NASA Global Landslide Catalogue (https:// gpm. nasa. gov/ lands lides/ index. html) and the Global Fatal Landslide Database (https:// www. un-spider. org/ links-and-resou rces/ data-sourc es/ global-fatal-lands lide-datab ase-gfld-unive rsity-sheff ield). There are also inventories over a more specific spatial scale within a region or country that resulted mainly from the compilation of landslides after catastrophic triggering events (e.g. Hervás and Bobrowsky 2009;Mateos et al. 2012). In the case of Spain, there is a national database of landslides (the Spanish Land Movements Database, BD-MOVES, http:// info. igme. es/ catal ogo/) published in 2016 by the Geological Survey of Spain (Instituto Geológico y Minero de España, IGME-CSIC). Traditionally, inventory maps have been produced through multitemporal aerial photo-interpretation and field surveys. However, it remains difficult and time-consuming to produce and update landslide inventories in most regions of the world, especially in mountainous areas with high extension and poor accessibility (Bekaert et al. 2020). Herrera et al. (2018) compared the European Landslide Susceptibility Map (ELSUS v1) (Günther et al. 2014) with the mapped landslides in each country to analyse where to expect more landslides than those already inventoried. For example, the completeness of the national landslide inventory in Spain (BD-MOVES) is less than 5% . The inventoried landslides are usually the most morphologically visible on the landscape while other typologies with more diffuse boundaries are often overlooked. Therefore, new technologies such as satellite remote sensing or advanced landscape analysis are gaining prominence to improve and optimise landslides' mapping at a regional scale.
Differential Synthetic Aperture Radar Interferometry (DInSAR) is a remote sensing technique that exploits radar satellite images to derive multitemporal displacement measurements of the ground surface. Among numerous applications, DInSAR is a powerful tool to map active landslides and produce inventory maps (e.g. Bekaert et al. 2020;Reyes-Carmona et al. 2020;Crippa et al. 2021). Thanks to the wide coverage (up to a 250 km swath width) and the high temporal resolution (up to 1 day) of the radar images, DInSAR makes analysing very large areas and constantly updating a landslide inventory possible. Some initiatives, such as the Geohazards Exploitation Platform (GEP), developed by the European Space Agency (ESA), aim to promote the use of DInSAR techniques in a user-friendly way. The GEP is a web-based platform (https:// geoha zards-tep. eu/ # !) that allows users to perform automated and independent DInSAR analysis, offering quick results in just 24-48 h. The GEP services have already been successfully applied to discover new landslides, between other natural processes (e.g. Manunta et al. 2016;Galve et al. 2017;Tapete and Cigna 2017;Foumelis et al. 2019;Reyes-Carmona et al. 2021;Gaidi et al. 2021).
Landscape Analysis techniques can also be used to identify (1) recent geological processes, such as active tectonics, fluvial captures or landslides; (2) local conditions, like lithological contrasts; and (3) the imprint of past geomorphic processes, such as glacial erosional features or high-elevation low-relief surfaces (e.g. Larue 2008;Pérez-Peña et al. 2010;Antón et al. 2014;Troiani et al. 2014;Subiela et al. 2019). These phenomena usually disturb the drainage Landslides 20 • (2023) Original Paper network and express themselves topographically on rivers by creating knickpoints or knickzones (Walsh et al. 2012). A knickpoint or knickzone is an abnormal increase of the gradient in a specific segment of a river. Punctual changes in the gradient are commonly known as 'knickpoints' while 'knickzones' are referred to gradient changes that affect a longer transect of a river. Both knickpoints and knickzones can be assessed by indexes that analyse river gradient, such as the normalised steepness index (k sn ): the higher the gradient change, the higher the k sn index. Knickzones are reflected as clearly higher values than the rest of values along the river profile, which are considered as anomalous values. The successful application of gradient-related indexes to detect landslides has already been proven in several studies (Pánek et al. 2007;El Hamdouni et al. 2010;Walsh et al. 2012;Troiani et al. 2014Troiani et al. , 2017Penna et al. 2015;De Palézieux et al. 2018;Ahmed et al. 2019;Subiela et al. 2019;Piacentini et al. 2020;Gu et al. 2021;Liu et al. 2021). In recent years, GIS platforms and high-resolution DEMs have facilitated the application of geomorphometric techniques in terms of time-consumption and cost-effectiveness (Troiani et al. 2014. These techniques also allow studying large areas accurately and efficiently to produce geomorphological maps (Weibel and Heller 1991;Pike 2000), including landslide inventories (e.g. Subiela et al. 2019).
In this study, we used a new combination of DInSAR and k sn index data to explore its effectiveness for landslide detection and mapping in a mountainous area. The southwestern sector of the Sierra Nevada mountain range (Southern Spain) was selected as the case study. We consider this sector a complex mountainous area mainly due to its accessibility, the high local relief, the homogenous lithology and the difficult recognition of landslides on its landscape. Through DInSAR techniques, we obtained the first ground displacement map of this sector of Sierra Nevada, where unstable areas were detected and related to active landslides. To perform the Landscape Analysis, we applied a new morphometric index: the double normalised steepness (k snn ) index. The k snn index was derived from the conventional normalised steepness (k sn ) index. It enabled us to identify landslide anomalies by filtering the general tectonic signal observed in the Sierra Nevada from the k sn values. Our results show that, despite their limitations, the combination of both techniques facilitated the identification and mapping of large landslides. The use of high-resolution DEM-derived products and field observations was also essential for the delimitation of landslides. With such a data combination, we provided an inventory with a higher degree of completeness than the previous one (the BD-MOVES). Our analysis also allowed us to identify, for the first time, the existence of Deep-seated Gravitational Slope Deformations (DGSDs) in the Sierra Nevada, as well as to contemplate their related hazard.

Study area
The study area contains the northeastern part of the Guadalfeo River Basin, located in the Province of Granada, Southern Spain ( Fig. 1). This area (378.5 km 2 ) includes the sub-basins of six tributaries of the Guadalfeo River that are, from west to east, the Lanjarón, Sucio, Chico, Seco, Poqueira and Trevélez rivers (Fig. 1). These rivers drain the southwestern side of Sierra Nevada mountain range, where they have excavated steep V-shaped valleys due to the high topographic gradients: 35 km from 3479 m.a.s.l. (Mulhacén Peak) to the coastline. The Sierra Nevada was declared a Biosphere Reserve by UNESCO in 1986, a Natural Park in 1989, and a National Park in 1999. It is a privileged and representative spot of the Mediterranean high mountain systems. This Alpine setting also has a rich cultural and historical heritage linked to several relict Berber villages known as 'La Alpujarra' (Fig. 1). This region comprises 25 small picturesque villages with a total population of around 25000 people, being the municipalities of Órgiva (5420 inhabitants) and Lanjarón (3720 inhabitants) the most populated. Moreover, the Sierra Nevada is plenty of uncoated ditches excavated in the ground (locally known as 'acequias de careo'), originally from the Middle Ages (Martín-Civantos 2010) with an important cultural and hydrological value. This irrigation system was designed to infiltrate the snow melt and runoff water in the wetter months to have spring water supply during the driest months (Martos-Rosillo et al. 2019). Nowadays, more than 700 km of acequias is still working in the Sierra Nevada as a sustainable groundwater recharge system.
From a geological perspective, the Sierra Nevada is located in the central Betic Cordillera, which is the western termination of the Mediterranean Alpine orogen linked to the broad-scale collision between Africa and Iberia (DeMets et al. 1994). The main outcropping geological units are the Nevado-Filábride Complex, the Alpujárride Complex, and the Neogene-Quaternary sedimentary rocks (Fig. 2). The internal structure of the Nevado-Filábride Complex is very heterogeneous, characterised by multiple transposed foliations and lineations as the result of a complex polyphase deformation story (e.g. Jabaloy et al. 1993;Martínez-Martínez et al. 2002;Puga et al. 2011;Aerden et al. 2013;Ruiz-Fuentes and Aerden 2018). The subdivision of the Nevado-Filábride Complex is still under scientific discussion (e.g. Puga et al. 2002;Martínez-Martínez et al. 2002;Gómez-Pugnaire et al. 2012;Sanz de Galdeano and Lopez-Garrido 2016;Santamaría-López et al. 2019), but in this study, we followed the classification proposed by Martínez-Martínez et al. (2002). According to these authors, there are two lithological units of metamorphic rocks (Fig. 2): the Ragua Unit and the Calar-Alto Unit. Black-colored graphitic schists entirely form the Ragua Unit (Paleozoic). The Calar-Alto Unit is subdivided into two formations: the Montenegro Formation (Paleozoic), formed by graphitic micaschists with alternation of quartzites, and the Tahal Formation (Permian-Triassic), formed by light-coloured schists with isolated amphibolites and marbles. These three lithologies outcrop in most of the study area (Fig. 2), and we consider the lithological sequence relatively homogeneous from the mechanical point of view. The Alpujárride Complex is formed by Permian-Triassic metamorphic rocks that include, from older to younger, graphitic schists, phyllites and quartzites, mica-schists and dolomitic marbles (Fig. 2). The Neogene sedimentary rocks are related to fan deposits: conglomerates with intercalated sandstones (Aldaya et al. 1979). Quaternary deposits include fluvial deposits, travertines and landslide bodies (Fig. 2). The latter are those included in the Spanish Land Movements Database (BD-MOVES, http:// info. igme. es/ catal ogo/). The contact within the Nevado-Filábride and Alpujárride complexes and the inferred limit between the two main units of the Nevado-Filábride Complex are inactive extensional detachments. These detachments and other high-angle normal faults have conducted an extension and consequent exhumation of the complexes since the Miocene (Galindo-Zaldívar et al. 1989;Martínez-Martínez 2006). One of these normal faults (the 'Lanjarón Fault' in Fig. 2) is considered to be probably active, although there are no clear signs of activity at present (Sanz de Galdeano et al. 2003). For this reason, this fault is catalogued as a debated fault in the Quaternary Active Faults Database of Iberia (QAFI, http:// info. igme. es/ qafi/). The overall structure of Sierra Nevada is a large-scale antiformal fold that coincides with the highest elevations of the mountain range. Despite the uplifting stated earlier, Pérez-Peña et al. (2010) inferred that the present-day drainage pattern started to develop in the Pleistocene. These authors analysed several geomorphic indexes to demonstrate that the Sierra Nevada is tectonically active nowadays and that the recent uplift is concentrated within the southwestern sector, where our study area is located.
From a geomorphological perspective, the current morphology of the study area is highly influenced by the uplifting of the western part of Sierra Nevada, which has conditioned a strong river incision (e.g. incision rates of 5-5.9 mm/year, according to Chacón et al. 2001;Reinhardt et al. 2007). Consequently, river incision produced over-steepened slopes prone to landslides (El Hamdouni et al. 2010). These landslides are mostly planar slides, earthflows and rotational slides that occur mainly in the Alpujárride phyllites, the Nevado-Filábride schists and the Neogene granular or slightly cohesive deposits (El Hamdouni et al. 2010). Chacón et al. (2007) carried out a landslide inventory of the whole Province of Granada in which they identified, just in our study area, a total of 67 landslides ( Fig. 2): 7 rockfalls, 36 slides and 24 surface processes such as creeping and solifluction. These landslides were later included in the BD-MOVES, in which the 24 surface processes were classified as flows. In the highest elevations of the range, glacial and periglacial morphologies are predominant, with magnificent examples within the upper part of the Poqueira and Trevélez valleys (Gómez-Ortiz et al. 2002). These are related to small valley and cirque glaciers that were the most southern in Europe during the Last Glacial Maximum (Gómez-Ortiz et al. 2012). Deglaciation took place around 14 ka ago, and rock glaciers were formed immediately after, affected by intense periglacial conditions until 7 ka ago (Palma et al. 2017). These authors also established an Equilibrium Line Altitude (ELA) for the last glaciation at 2650 m in the southern sector of Sierra Nevada.
The climate in the Sierra Nevada corresponds to a semiarid cold mountain climate (Dsc according to the Köppen climate classification). Mean annual temperatures are around 0 ºC on the summit areas, and the average annual precipitation is around 710-750 mm. Snow is usually present from early December to the end of May, being the snowline settled at 2100 m.a.s.l.

Methodology
The methodology of this work consists of the following main stages: (1) calculation of surface ground displacement, derived from DInSAR techniques through the Geohazard Exploitation Platform; (2) calculation of the double-normalised steepness (k snn ) index that we propose in this study for the first time; (3) examination of DIn-SAR results and identification of unstable areas; (4) interpretation of k snn anomalies; and (5) creation of an updated landslide inventory map by combining the DInSAR and k snn data with geomorphological observations. Essentially, our main interest was to evaluate the reliability of both techniques and their complementation to facilitate landslide mapping in a complex mountainous area such as the Sierra Nevada (Southern Spain).

Differential Synthetic Aperture Radar Interferometry (DInSAR)
To derive the DInSAR data, we used the Parallel Small Baseline Subset (P-SBAS) processing chain (Casu et al. 2014), which is based on measuring ground displacement in points that are Distributed Scatterers (DSs): small targets of similar radar amplitude that usually correspond to natural features (e.g. agriculture areas, open fields, bare soil, rock surfaces). For this reason, the SBAS methods are very suitable for analysing rural environments and arid areas with low vegetation and debris (Even and Schulz 2018), such as the Sierra Nevada Range. The P-SBAS processing chain has been implemented on the European Space Agency (ESA)'s Geohazards Exploitation Platform (GEP), as detailed in De Luca et al. (2015). It was possible to use the P-SBAS in a fully automated and unsupervised manner through the GEP web-portal thanks to a granted permission of access in the framework of the ESA Network of Resources (NoR) Initiative. For being an automated processing chain, we only had to select the desired input satellite images and decided on a few parameters that were latitude and longitude of the reference point, polarisation, processing mode, DEM type and coherence threshold.
As radar images are acquired in two different geometries (i.e. ascending and descending orbits), we carried out a processing per each acquisition geometry. We used 101 Sentinel-1B images for the ascending orbit processing with a temporal sampling of 12 days from the 30th of September 2016 to the 13th of March 2020. For the descending orbit, we used 241 Sentinel-1A and Sentinel-1B images with a temporal sampling of 6 days and covering a period from the 22nd of December 2014 to the 19th of March 2020. For both processing jobs, the following parameters were established: vv polarisation, SRTM-1 DEM, coherence threshold of 0.85 and reference point in Lat 36.848, Long −3.497 (WGS84 projection). In the ascending orbit, the satellite travels along the NNW-SSE direction and looks to the east, while in the descending orbit, the satellite travels along the SSW-NNE direction and looks to the west. The direction to which the satellite looks is named Line Of Sight (LOS) direction, and the DInSAR velocity is always calculated along this direction. For this reason, the detected movement of the ground is registered as approaching or distancing from the satellite: negative values indicate that points move away from the satellite while positive values refer to points moving toward the satellite. Therefore, the output of each processing was a set of points representing the LOS mean displacement or velocity in mm/year. The pixel size of each point was 90 m.
The DInSAR surface ground displacement rate (or mean annual velocity) maps are represented in equal intervals by establishing a threshold for discriminating stable from unstable points as two times the standard deviation of all the measured velocity points (Barra et al. 2017). Therefore, the stability range was set between 6 and −6 mm/year for the ascending orbit processing and between 5 and −5 mm/year for the descending orbit. It is important to remark that the stability range also represents the general noise of the results (i.e. the sensitivity of the technique). This fact implies that a point classified as 'stable' can be either truly stable or unstable, but the displacement cannot be detectable by the technique.

Landscape Analysis
The morphometric analysis of rivers was computed through the Python library 'landspy', freely available at https:// github. com/ geolo vic/ lands py. This computing tool is based on the stream-power model, which relates the local channel slope and the contributing drainage area upstream (Perron and Royden 2013) (Eq. 1): where S is the slope of the channel, A is the up-stream drainage area, k s is the steepness index, and θ is the concavity index.
The traditional way to analyse k s and θ is through linear regression in logarithmic area-slope river profiles. This procedure presents the problem of the high autocorrelation in both parameters, increased even by the logarithmic scale (Wobus et al. 2006;Kirby and Whipple 2012). As the θ index does not vary in high ranges, a solution to derive the steepness index is by using a fixed reference concavity (θ ref ) to obtain a normalised steepness index (k sn ). The most popular approach to derive the k sn index from a fixed reference concavity was proposed by Perron and Royden (2013), through the integration of Eq. (1) and the definition of the Chi index (χ). By applying this integration, the k sn index is calculated by linear regression between the Chi index and elevation (i.e. the slope of a Chi-elevation plot is the k sn index).
Even normalising the steepness index, the highest values still occur in high-relief areas. This fact makes it difficult to compare gradient changes in areas with prominent topographic differences. In these areas, compared to areas with low-to-moderate topography, Chi-elevation profiles are steeper and, thus, the k sn values are higher. This is the case of the Sierra Nevada Range, where active tectonics have generated high topographic gradients and high k sn (1) S = k s A values along the main rivers , complicating the identification of knickpoints unrelated to tectonics. To discard these topographic gradient trends resulting from active tectonics, we proposed a double-normalised steepness index (k snn index) by normalising the k sn index with the mean slope of the Chi-elevation plot (mean k sn ). This normalisation eliminates these trends and highlights knickzones that were not evidenced through the k sn index analysis.
To derive the k snn index for the study area, the only input required by the tool 'landspy' was a 10-m resolution digital elevation model (DEM), obtained from the Andalusian Environmental Information Network (REDIAM, https:// www. junta deand alucia. es/ medio ambie nte/ portal/ acceso-rediam). Following, channel gradients and the Chi, k sn and k snn indexes were derived for the six selected sub-basins of the Lanjarón, Sucio, Chico, Seco, Poqueira and Trevélez rivers. These indexes were calculated for rivers' segments of 700-m length, which was considered an appropriate value for the scale of our analysis. The Chi index was computed with a reference concavity of 0.45, which is a suitable value according to previous studies in the Betic Cordillera (Bellin et al. 2014;Azañón et al. 2015). Once obtained the K snn index, we defined five intervals by a natural break classification for a proper data visualisation. The k snn values higher than 7.4 were considered to be anomalous. This threshold is the cut value of the natural break intervals that is closest to one standard deviation of the data (7.6). Lastly, we focused on identifying anomalies (k snn values higher than 7.4) with a length equal to or longer than two channel segments (1400 m) along the trunk rivers and their tributary channels.

Landslide detection and mapping
Once we obtained the raw DInSAR and k snn results, we identified the unstable areas from the DInSAR ground displacement maps and the k snn anomalous values (i.e. knickzones) from the k snn map. Therefore, we inspected the spatial distribution of the unstable areas and knickzones in combination with the following existing data on a GIS environment: (1)  . Crossing all of this information, we associated unstable areas and knickzones with landslides.
Our main aim was to delimit the landslides' boundaries as accurately as possible. For it, the exhaustive examination of products derived from high-resolution DEMs was essential. We used 2-m and 5-m DEMs, freely available at the Spanish National Geographic Institute (https:// centr odede scarg as. cnig. es/ Centr oDesc argas/ index. jsp) to derive the hillshade, slope, aspect, rugosity and topographic openness maps. These products were also exported to Google Earth for a 3D, more accurate visualisation of landslide-related features. The hillshade model combined with the slope, aspect and topographic openness maps allowed recognising the slope breaks related to the head and lateral scarps. The hillshade model was also useful to Landslides 20 • (2023) Original Paper identify secondary scarps, benches and rock deposits within the landslides, as well to delimit the slide masses, which were expressed as an increase of rugosity and convexity of the ground surface. Moreover, we carried out field surveys for a further visual inspection of morphologies and rock deposits of the landslides, as well as for checking the observations made by the GIS analysis. All of this work conducted us to provide an updated landslide inventory of the SW sector of Sierra Nevada. We also made a classification of the mapped landslides into different typologies.

DInSAR displacement rate maps
The mean displacement rate or velocity maps in ascending and descending orbits are shown in Fig. 3a, b, respectively. A total number of 33 unstable areas were detected: 16 areas by the ascending orbit (polygons from 1 to 16 in Fig. 3a) and 17 areas by the descending orbit (polygons from 17 to 33 in Fig. 3b). Some of these areas are coincident in both geometries (1 and 19, 3 and 23, 4 and 24, 11 and 29, 15 and 31), what means that we detected 28 different unstable areas within the study area. The ascending orbit processing provides a better point coverage within the western slopes of the valleys, and the maximum LOS velocity recorded is −32 mm/year along the Trevélez River's valley (area 13 in Fig. 3a). On the contrary, the descending orbit provides the better point coverage within the eastern slopes, with a maximum LOS velocity recorded of −31 mm/ year along the Poqueira River's valley (area 26 in Fig. 3b). Figure 4a shows the k sn map, where it is hard to identify anomalies (or knickzones) as values are consistently high. This fact is due to the strong river incision related to the active uplift of Sierra Nevada . Such tectonic signal produces a steeper chielevation profile that can be described by its general slope (mean k sn ) (Fig. 4b). By normalising the k sn values (Fig. 4c), this tectonic signal was reduced, and we obtained the k snn index (Fig. 4d), which clearly evidence knickzones in the k snn map (Fig. 4e).

k sn and k snn maps
We detected 10 knickzones within the study area, named from number 1 to 10, which are distributed as follows (Fig. 4e): knickzones 1 and 2 along the Lanjarón River, knickzone 3 along the Sucio River, knickzones 4 and 5 along the Chico River, knickzone 6 along the Seco River, knickzones 7 and 8 along the Poqueira River and knickzones 9 and 10 along the Trevélez River. We also detected other 15 k snn anomalies, named with letters from 'a' to 'o' distributed along the tributary channels of the trunk rivers (Fig. 4e). There is just one anomaly along the tributary channels of the Seco and Chico rivers (anomaly 'a' and 'b' respectively). Other seven anomalies (from 'c' to 'i') along the Poqueira River and six anomalies (from 'j' to 'o') along the Trevélez River were also identified (Fig. 4e).

Landslide inventory map
The landslide inventory of the study area is shown in Fig. 5. The k snn values of the trunk rivers and tributaries (Fig. 4e), as well as the unstable points from both DInSAR geometries (Fig. 3), are also plotted to facilitate the correlation between these data and the mapped landslides. Through both DInSAR and k snn anomalies data together with geomorphological observations, we could delimit a total of 28 landslides. Such a mapping implies that 126.8 km 2 of the study area is affected by landslides, which means 33.5% of its total extension. Table 1 details the associations of the DInSAR unstable areas and K snn anomalies for each of these landslides. Their names and abbreviations were established according to the trunk river where they are located (Fig. 5, Table 1): 'L' for Lanjarón, 'Su' for Sucio, 'C' for Chico, 'Se' for Seco, 'P' for Poqueira, and 'T' for Trevélez (Fig. 5). They are also numbered from lowest to highest towards the headwater for each river.
From DInSAR results, unstable areas from 1 to 32 (Fig. 3) are associated with 25 different landslides ( Fig. 5; Table 1). Some of them are almost entirely active (e.g. the Lanjarón-1 or Sucio-1 landslides), while most show only active sectors within a larger landslide body (e.g. the Chico-1, the Poqueira-1 or the Trevélez-2 landslides). Regarding the k snn anomalies, we could confidently assume a dominant role of landslides and glacial morphologies on the knickzone generation after dismissing the influence of other phenomena. Anomalies linked to lithological contrasts are not expected as the valleys' slopes are formed mainly by schists from the Nevado-Filábride Complex (Fig. 2). Similarly, anomalies along the trunk rivers cannot be related to significant tectonic structures, such as clearly active faults. This is the case of the Lanjarón Fault (Fig. 2), which is not spatially correlated with knickzones 1 and 3 (Fig. 4e), what suggests that it is an inactive fault. Therefore, six of the trunk rivers' anomalies could be spatially associated with landslides (numbers 1 to 5 and 7), and from the tributary channels' anomalies, ten were related to landslides ( Fig. 4e; Table 1). Out of the remaining anomalies, numbers 8 and 10 and letters h, i, n and o, are linked to glacial and periglacial morphologies (Figs. 4e and 5), according to the mapping of Gómez-Ortíz et al. (2002). The origin of knickzone 6 is unknown, while knickzone 9 can be linked to either processes: the Trevélez-2 Landslide or active tectonics, according to the hypothesis from Azañón et al. (2015). Similarly, anomaly k results from a fluvial capture that cannot be certainly due to the Trevélez-2 Landslide (Figs. 4e and 5).
We also made a preliminary classification of the mapped landslides in four different typologies: (a) Deep-seated Gravitational Slope Deformation (DGSD), (b) rockslide, (c) earthflow, and (d) rock spreading. Figure 6 reclassifies the landslide inventory shown in Fig. 5, taking into account these typologies, and Table 2 summarises the main characteristics of DGSDs and rockslides (i.e. the two most common typologies).Most of the mapped landslides are of DGSD type (17 landslides), developed within the Nevado-Filábride schists (Fig. 2) and located along the lower and medium part of the valleys. They are large landslides of variable size, with areas from 1.4 to 31.6 km 2 , occupying 28.4% of the study area ( Fig. 6; Table 2). These DGSDs affect entire slopes of the trunk rivers' valleys, and most of their head or main scarps reach the valley ridges. However, the DGSDs do not show well-defined head scarps and/or lateral boundaries, which made their delimitation an intricate task (Fig. 7). Most of these DGSDs are compounded by smaller-size rotational slides or rockslides that generate multiple minor scarps and benches, which facilitated the delimitation of the slide masses, which were also well-evidenced by an increase of the slope rugosity in the hillshade (Fig. 7). Active movements, with LOS velocities up to −32 mm/year, are registered within punctual Landslides 20 • (2023)  Figs. 3 and 5). The longest knickzones (numbers 1 to 5 and 7 in Fig. 4e) along the trunk rivers are associated with these DGSDs as their downslope force generates stream stretches and deviations, which disrupt the river equilibrium profile. In some cases, other anomalies along tributary channels are linked to slope breaks of nested movements located at the lower part of the larger DGSDs (e.g. anomalies b, c, and d in Fig. 4e).
We also mapped eight other landslides that we classified as rockslides, according to the descriptions of Crosta et al. (2014) and Borrelli and Gullà (2017). They are mainly distributed in the upper part of the Trevélez valley, involving the Nevado-Filábride schists (Fig. 2). The Fig. 4 a k sn map of the study area. b Chi-elevation plot of the Poqueira River. The mean k sn is indicated with a red-coloured dashed line. c k sn profile of the Poqueira River. d k snn profile of the Poqueira River. The detected knickzones are also shown along the profile. e k snn map of the study area. Anomalies located along the trunk rivers are indicated with numbers from 1 to 10, and anomalies along tributary channels are indicated with letters from 'a' to 'o' Landslides 20 • (2023) rockslides occupy 4.3% of the study area, and they have a smaller size than the DGSDs (areas up to 4.9 km 2 ). These slope movements are recognisable by well-defined curved main scarps that are covered by debris (Fig. 8). They are characterised by multiple secondary movements as well as by a very high ruggedness in the hillshade model (Fig. 8a) and waviness of the ground surface (Fig. 8b, c). These slope movements imply a deep sliding of the bedrock together with a shallow sliding of debris generated from thermal alterations (gelifraction). We also found some nival deposits like protalus ramparts that are accumulated in benches of the larger landslide body (Fig. 8b). This fact proves the influence of gelifluction and gelifraction processes into the rockslides' kinematics, as well as the debris mobilisation by snow melting processes in a periglacial environment. These rockslides show LOS velocities up to −23 mm/year (Figs. 5 and 6). Some of them can be related to k snn knickzones of short length along the tributary channels (e.g. anomalies l and m in Fig. 4e), while none of them generate anomalies along the trunk rivers.
Lastly, two earthflows (the Seco-1 and the Trevélez-1 landslides) and one rock spreading (the Seco-2 Landslide) (Figs. 5 and 6) were also identified within the Alpujárride Complex and Neogene sedimentary rocks (Fig. 2). These three landslides entail the minority of the study area (0.9%). All of them are active landslides (maximum LOS velocities of −19 mm/year) (Fig. 3), and there are no k snn anomalies related to them (Fig. 4e).

Strength of the combination of DInSAR and k snn analysis for landslide detection
The DInSAR services of the Geohazards Exploitation Platform (GEP) have already been demonstrated to be effective tools for landslide detection (e.g. Galve et al. 2017;Reyes-Carmona et al. 2021;Gaidi et al. 2021;Cigna and Tapete 2021) as well as the k s and k sn index analysis (e.g. Pánek et al. 2007;Walsh et al. 2012;De Palézieux et al. 2018;Ahmed et al. 2019;Gu et al. 2021). Nevertheless, the advantages of combining both data sets have not been explored until the present study. From a quantitative point of view, 25 of the total 28 mapped landslides (89.3%) have been identified as active landslides through the DInSAR data, while 17 landslides (60.7%) have been attributed to k snn anomalies along the trunk rivers and/ or tributary channels. All the landslides were revealed at least by one of the methods, while 14 landslides (50%) were evidenced by both. These results are satisfactory enough as we provided a better landslide inventory than the previous one of the Spanish Database of Landslides (BD-MOVES). Our new inventory doubles the area affected by landslides (33.5%) in relation to the BD-MOVES (14.5%) (Fig. 6). It is also remarkable that most of the landslides inventoried in this work have larger dimensions than those previously

Original Paper
identified (e.g. landslides along the Lanjarón valley, the Poqueira-1 or the Trevélez-2 landslides in Figs. 5 and 6). Other large landslides were identified for the first time in this study, such as the Lanjarón-4 or the Poqueira-2 landslides (Figs. 5 and 6). We also dismissed several landslides at higher elevations due to the absence of ground motion, k snn anomalies and, most importantly, any recognisable landslide-related features. Instead, glacial and periglacial morphologies mask any other processes in these areas (Gómez-Ortiz et al. 2002;Palma et al. 2017).
The great potential of DInSAR is revealing the active landslides, which rapidly highlights the slopes that should be further investigated. However, in our study case, ground movements were generally detected only within some sectors and not in the whole landslides' bodies (Fig. 3). This fact implies that the large size of the landslides may be underestimated if the attention is focused on mapping just the active zones, which usually correspond to smaller nested movements. In this sense, the identification of k snn knickzones along trunk rivers was useful to reveal the complete Chico-2 C-2 20, 21 4, 5 -DGSD Chico-3 C-3 22 5 -DGSD  Fig. 5) or even where there is not registered DInSAR points (e.g. the Lanjarón-2 or Lanjarón-4 landslides in Fig. 5). The combination of these two datasets also provided two different temporal perspectives about the landslides: DInSAR shows a very short-term or recent activity, while k snn anomalies point out a long-term activity that shows that a landslide has being perturbing fluvial channels at century or millennial timescale. The scarcity of measurement points in some areas is due to some limitations that should be mentioned when applying DInSAR in mountainous areas. Some natural terrain properties usually scatter the radar signal of the satellite images, which introduces noise to the DInSAR processing and reduces the point coverage (Hanssen 2001). Some of these properties are steep slopes (e.g. the lower part of the valleys), dense vegetation (e.g. the Chico River valley's landslides), terrace fields (e.g. the Poqueira-1 or Trevélez-2 landslides) and snow cover (e.g. highest elevations of the Sierra Nevada) (Fig. 3). Other intrinsic limitation of DInSAR is the decrease of the radar sensitivity when true landslide displacement deviates from the satellite LOS (Schlögel et al. 2015), what makes slow movements to be not registered or underestimated. Typically, DGSDs have velocities of millimeter/year order, close to the usual stability range of DInSAR processing (around 5 mm/year). This fact implies that unstable points may have not been registered in the cases of the Lanjarón-2, Lanjarón-4 or the Trevélez-2 landslides (Fig. 6). Therefore, it is important to remark that these landslides could be either stable or unstable but not detected by our DInSAR processing. Combining ascending and descending data is usually helpful to deal with such limitations, as radar sensitivity varies in each geometry depending on the slope orientations. Some examples are the Poqueira-1 or the Trevélez-6 landslides in which unstable areas (5 to 9 and 13 in Fig. 3a) were detected just by the ascending processing. On the contrary, the descending processing detected the Lanjarón-1 or the Poqueira-2 landslides' unstable areas (18 and 26 in Fig. 3b, respectively). In this sense, the Geohazards Exploitation Platform (GEP) afforded us to obtain processing in both orbits in a very cost-effective way.
The great potential of the k snn analysis is revealing the true extension of large landslides. In the study area, the longest knickzones along the Lanjarón, Sucio, Chico and Poqueira rivers (numbers 1 to 5 and 7 in Fig. 4b) correspond to large DGSDs (Fig. 6). This type of landslides shows a relevant control on the evolution of drainage network along the valleys' bottom, where knickzones are commonly formed. The downslope force of DGSDs generates deviation and narrowing of the river channels, which may shift the focus of fluvial erosion (Korup 2006). It is also possible for

Original Paper
a channel bed to be uplifted by the thrust of the landslide mass, if the failure plane extends below the channel (Bartarya and Sah 1995). These actions originate anomalous changes in the gradient of the river profiles when landslides are active and their gravitational force is able to counteract fluvial erosion. For this reason, rivers cannot reestablish their equilibrium or steady-state profile, and knickzones are generated, which is reflected by anomalously high values of gradient-related geomorphic indexes. As the case of Sierra Nevada, other studies worldwide show the spatial coincidence of landslides, including DGSDs and other large rockslope instabilities, with anomalous values of gradient-related indexes (Korup 2006;El Hamdouni et al. 2010;Walsh et al. 2012;Troiani et al. 2014Troiani et al. , 2017Penna et al. 2015;Subiela et al. 2019).
The limitation of the k snn analysis is that an anomaly does not necessarily have to be formed. According to Troiani et al. (2014), the formation of knickzones by landslides is dependent on several factors, such as the landslide size, the amount of sediment delivered by the landslide (landslide activity) or the river capacity to incise the landslide deposit (erosion power). For example, the large knickzones of the Lanjarón, Sucio, Chico and Poqueira rivers were generated as these four rivers' erosive power may be lower than the landslide activity. These cases are contrary to the case of the Trevélez River, where there are no long knickzones associated with large DGSDs, such as the Trevélez-2 or the Trevélez-6 landslides (Figs. 5 and 6). Nevertheless, shorter knickzones along tributary channels and DInSAR data were essential to recognise some landslides such as the Trevélez-6 or the Trevélez-9 (Fig. 5). As many natural processes can contribute to a knickzone generation, another challenge of the k snn analysis is decoding which is the dominant process (Walsh et al. 2012). Frequently, a predominant process may mask other processes of our interest, such as landslides. Some examples are knickzones 8 and 10 (Fig. 4e), where glacial and periglacial morphologies generate strong anomalies that cannot be certainly attributed to nearby landslides (e.g. the Poqueira-6 or the Trevélez-11 landslides). Another example is knickzone 9 (Fig. 4e), whose origin can be controversial. According to Azañón et al. (2015), it can be related to the water gap after the Guadalfeo River's migration that resulted from the recent uplifting of Sierra Nevada. Nevertheless, we consider that the Trevélez-2 Landslide could also contribute to this knickzone generation (Fig. 5). To deal with these ambiguities, new tools need to be developed to unmask anomalous values related to a specific process. In this sense, the k snn index calculation made it possible to eliminate a considerable influence of tectonic uplift and the consequent topographic gradients in the most active sector of Sierra Nevada. Thus, the visualization of knickzones related to landslides was greatly facilitated by the k snn index (Fig. 4e), in contrast with the conventional k sn index (Fig. 4b).
Despite their limitations, we conclude that both DInSAR techniques (e.g. Notti et al. 2010;Bianchini et al. 2013;Bekaert et al. Walsh et al. 2012;Troiani et al. 2017;Liu et al. 2021) can optimise the landslide detection in mountainous areas. Our results demonstrate that both methods not only can be well complemented, but also limitations of each one can be compensated. In this sense, the visual inspection of k snn anomalies and DInSAR data rapidly spotlight the slopes of interest on which to focus to recognise geomorphological features for landslides' delimitation. It should also be remarked that both the Geohazards Exploitation Platform (GEP) and the Python library 'landspy' are very user-friendly tools for obtaining quick results of DInSAR and geomorphic indexes, respectively. This makes both initiatives promising to improve and update landslide databases not only for the scientific community but also for public administrations.

Prevalence of DGSDs among the landslides affecting the SW of Sierra Nevada
Although the Province of Granada, including the Sierra Nevada, has been analysed thoroughly by different research teams for 30 years (see Chacón et al. 2007 and compiled references therein), DGSDs and their prevalence have not been pointed out until the present study.
The main landslide research and inventories of the Sierra Nevada area were produced in the 1980s, 1990s and the early 2000s (e.g. Chacón and Sofia 1992;El Hamdouni 2001;Chacón et al. 2007) through photointerpretation, field work and basic GIS analysis, before the free availability of high-resolution DEMs such as those used in our research. Despite that more recent studies (e.g. Azañón et al. 2008;Jiménez-Perálvarez et al. 2011Jiménez-Perálvarez 2018) estimated a medium/high degree of landslide susceptibility in the SW sector of Sierra Nevada, the mapped landslides were scarce. Moreover, the attention and popularity of DGSDs in the landslide research community have progressively increased since the 1990s (Chigira 1992a, b;Dramis and Sorriso-Valvo 1994) and especially, during the last two decades (Agliardi et al. 2001;Korup 2005;Gutiérrez-Santolalla et al. 2005;Ambrosi and Crosta 2006;Agliardi et al. 2009;Crosta et al. 2013;Chigira et al. 2013a, b;Agliardi et al. 2013;Tsou et al. 2015; Della Seta Trevélez-11 (T-11). Red-coloured asterisks facilitate the comparison between both images of the slope valley. b Lateral view of a rockslide within the Trevélez-6 Landslide and its related features: main scarp (m.s) covered by debris, wavy surface (w.s), protalus rampart (p.r) and debris from gelifraction processes (d) covering a large part of the slope. c Panoramic view of the Poqueira-4 Landslide (rockslide) and its main features: curved main scarp (m.s) covered by debris and a wavy ground surface (w.s). A glacial cirque (c) and a rock glacier (r.g) located above the landslide are also indicated Landslides 20 • (2023Mariani and Zerboni 2020;Crippa et al. 2021). This type of landslide was already described before through different terms such as sackung (Zischinsky 1969), mass rock creep (Radbruch-Hall 1978) or rockflow (Varnes 1978). Therefore, it is not surprising that DGSDs in the Sierra Nevada were not mapped in previous inventories, as these slope movements can be difficult to identify if a surveyor is unfamiliar with them and/or due to the lack of high-quality topographic data. In our study area, the DGSDs recognition and delimitation were supported by (1) the knowledge acquired in other geological settings such as the Alps (e.g. Agliardi et al. 2001), Apennines (e.g. Di Luzio et al. 2004), Pyrenees (e.g. Gutiérrez-Santolalla et al. 2005 or the Carpathians (e.g. Pánek et al. 2011) that helped us in their identification by comparison with other known examples; (2) the high-resolution DEMs that offered us enough detail of the ground surface, even in forested areas, to identify morphological features of DGSDs; (3) a cutting-edge technology such as DInSAR that allowed to identify wide areas of ground motion along the slopes and (4) landscape analysis techniques that provided information about how large landslides perturb rivers and where we had to look up the hillside. In this way, we could focus our geomorphological research directly on the slopes where these techniques provided us data for then, recognising scarps, benches and slope convexities associated with DGSDs. Previously to this research, no one had ever had these resources to identify such large landslides. However, future research will for sure improve the mapping of DGSDs as they are always difficult to delimit accurately.
The challenges for DGSDs detection in the SW of Sierra Nevada should also be mentioned. As the Nevado-Filábride Complex is a homogeneous sequence of schists (Fig. 2), there are no clear key layers or lithological contacts, which usually facilitate the recognition of slope ruptures . In this context, only the surficial morphological features guided the recognition of these landslides and their boundaries, the latter being very diffuse and poorly defined. The general absence of well-developed DGSD-related morphologies (e.g. double ridges, open trenches, or counterscarps) (Fig. 7) also made mapping most of the DGSDs complex. Moreover, the presence of glacial and periglacial morphologies usually makes difficult the surficial mapping of landslides (Weidinger et al. 2014) due to the similarity of slope deposits with glacial moraines (Hewitt 1999) and between the scarps with glacial cirques (Turnbull and Davies 2006). Some examples are the Trevélez-4 and Trevélez-6 landslides ( Fig. 8a, b), where main scarps were mapped as glacial cirques by Gómez-Ortiz et al. (2002) and Palma et al. (2017). Similarly, we interpreted as protalus ramparts (Fig. 8b) some deposits that were mapped as moraine segments by Gómez-Ortiz et al. (2002).
It is worth noting that it was expectable to find DGSDs in the Sierra Nevada as they are widespread in other Alpine mountain ranges (Jarman et al. 2014;Del Rio et al. 2021a, 2021bCrippa et al. 2021) where tectonic exhumation controls topography , and the constant relief uplift produces a high fluvial incision of valleys (Korup et al. 2007;Tolomei et al. 2013;Tsou et al. 2015;Demurtas et al. 2021). DGSDs usually affect the entire length of high-relief valley flanks , and we consider that the local relief of these valleys is high enough (0.8-1 km) to cause the gravitational collapse of their slopes. Furthermore, the rocks that compose the Sierra Nevada are common materials (i.e. foliated metamorphic rocks) where DGSDs are prone to occur . For this reason, this article is the first, but it should not be the last to investigate and map DGSDs in other sectors of the Sierra Nevada. Detailed morpho-structural studies about the internal segmentation of specific DGSDs as well as the research of their predisposing or causal factors (e.g. Agliardi et al. 2001Agliardi et al. , 2013Ambrosi and Crosta 2006;Crosta et al. 2013;Crippa et al. 2021) could also be carried out for a comprehensive understanding of this phenomena and its integration into the relief evolution of the Betic Cordillera.

Human-slope interactions in the SW of Sierra Nevada: implications of the new landslide inventory
The fact that the Sierra Nevada is a Natural and National Park implies a special commitment to its management and protection, which includes a better knowledge of natural processes such as landslides and their related hazard and risk . Therefore, the newly inventoried landslides may have positive and negative implications on infrastructures and populations that should be taken into account.
Regarding the acequias de careo, landslides may positively affect their proper functioning. The fractured rocks of a slide mass may work better as a groundwater reservoir than an intact rock massif. However, the water infiltration from the acequias could be excessive and ineffective. Water infiltration could also trigger an acceleration of a landslide, and in turn, these accelerations could damage the acequias. As an example, we observed a transect of a waterproofed acequia ('Acequia de Bérchules') that runs across the Trevélez-10 Landslide (Fig. 9a, b), probably to avoid excessive water infiltration or to prevent damage to the infrastructure. Further research has to be carried out to determine the positive and negative influence of landslides on water infiltration along the acequias.
The high local relief could have represented an important limitation for the human settlement in the SW sector of Sierra Nevada. However, the convex profiles of the slopes and the abundance of benches probably facilitated the creation of villages such as Pampaneira, Bubión and Capileira (Figs. 1 and 7b), as well as terraced agriculture and livestock farming activities. These slope morphologies are related to DGDSs, and our inventory brings to light their importance in the historical human occupation of the area. Despite this, the well-known negative issues of landslides should also be considered in our study area. As an example, Pampaneira, Bubión and Capileira (the three most touristic and famous villages of La Alpujarra region, Fig. 1) are settled within a large DGSD: the Poqueira-2 Landslide (Fig. 7b). Bubión is located just on top of a secondary or nested movement within this larger DGSD (Fig. 9c). DInSAR data evidenced ground displacement there, with LOS velocities up to −31 mm/year (unstable area 26 in Fig. 3b). Several damages, such as collapses of dry walls or piping phenomena, were observed in some terraced fields (Fig. 9d), which also prove the ground activity in this area. The ground movement has been generating damage during several decades in the penstock of the hydroelectric plant of Pampaneira (Alonso et al. 2021), which runs through the Poqueira-2 Landslide (Fig. 7b). Other eight villages, such as Pitres, Pórtugos and Busquístar (Fig. 1), are settled within the largest DGSD of the study area: the Trevélez-2 Landslide (Figs. 5 and 6). Although DGSDs are slow slope movements, they are long-lived phenomena that can evolve into faster secondary movements, such as rotational slides or debris/rock avalanches, which are Landslides 20 • (2023)

Original Paper
potentially destructive and may generate major risks to infrastructures and human lives (Soldati 2013). According to these authors, if monitoring and mitigation measures are focused on a single secondary movement within a larger DSGD, they may result incomplete and noneffective. However, the presence of DGSDs helps to identify large slopes that may be susceptible to catastrophic landslides in the future (Tsou et al. 2015). That is why their recognition, research and monitoring should be a priority in the Sierra Nevada.
Overall, detailed geological studies of these large landslides should be performed to understand better their internal structure and kinematics  and to model possible evolution scenarios for a correct hazard assessment (Soldati 2013;Spreafico et al. 2021). Similarly, in situ monitoring such as inclinometers or extensometers (Corominas et al. 2000), Global Positioning Systems (GPS) (Brunner et al. 2003) or the exploration drilling and geophysical techniques (Rogers and Chung 2017) should also be carried out for a precise subsurface characterisation of the landslides. In this work, we produced an updated and more accurate landslide inventory that is the starting point to assess the landslide hazard over an area (van Westen et al. 2008). Our new research has important implications for such assessment in the SW Sierra Nevada because a larger area than that initially mapped is potentially unstable. This fact also evidences that it is necessary to review and update the existing inventories by combining classical methods with innovative techniques to elaborate a landslide inventory as completely as possible.

Conclusions
Our work emerges the potential of integrating data from DInSAR techniques and Landscape Analysis to detect large landslides in a mountain range. Both are well-implemented tools that, when combined, considerably facilitated the mapping and understanding landslides in the SW part of Sierra Nevada. We provided an updated inventory of 28 landslides affecting 33.5% of the total area, compared with the area previously mapped in the Spanish Landslide Database (14.5%). Regarding the Landscape Analysis, we first proposed the k snn index derived from the conventional k sn index for reducing the effect of active tectonics in the Sierra Nevada. The visual inspection of DInSAR ground motion data (28 active areas) and k snn anomalies along rivers (17 knickzones) rapidly spotlighted the slopes to focus for landslide research. We proved that the limitations of both techniques could be compensated (e.g. DInSAR data showed activity in punctual sectors of larger landslides' bodies, for what k snn anomalies were useful to reveal such large sizes). Highresolution DEM-derived products were also essential for accurately delimitating the landslides' boundaries.
We distinguished two main typologies of landslides that have not been described in the area until the present work: rockslides Fig. 9 a Panoramic photograph of the Trevélez-11 (T-11) Landslide. Its main scarp is drawn by a black-coloured dashed line. An acequia de careo, which runs across the landslide, and its waterproofed transect are also indicated. b Detail of the waterproofed transect (black-coloured plastic) of the acequia de careo. c Lateral view of a secondary (or nested) movement and its scarp within the Poqueira-2 (P-2) Landslide. The villages of Capileira and Bubión are also indicated. d Damages of terraced fields within the nested body of the Poqueira-2 Landslide. Piping phenomena, minor scarps and damaged dry walls are evidences of the active ground movement. The village of Pampaneira is also indicated and DGSDs, the latter being the prevailing ones. Overall, the recognition and delimitation of these landslides were challenging due to their large size and diffuse boundaries, what makes them usually difficult to envisage. The presence of glacial morphologies and the homogeneous lithology (schists) also hindered the recognition of landslides' features. Finally, we suggested the relevant role that DGSDs may have had in the landscape evolution of the Sierra Nevada, and we offered a preliminary vision of their potential hazard, as DGSDs are likely to evolve into faster secondary movements. Our new inventory has relevant implications as landslides are larger and more abundant than previously considered, but further geological research and monitoring are still necessary for a proper landslide hazard assessment.

Funding
Funding for open access publishing: Universidad de Granada/CBUA. This work was mainly funded by the projects B-RNM-305-UGR18, A-RNM-508-UGR20 and P18-RT-3632 ('RADANDALUS') from the European Regional Development Fund (ERDF)/Junta de Andalucía-Consejería de Transformación Económica, Industria, Conocimiento y Universidades. This work was also supported by the ERDF through the project 'RISKCOAST' (SOE3/P4/E0868) of the Interreg SUDOE Programme and the project 'MORPHOMED' (PID2019-107138RB-I00) from the Spanish Ministry of Science (MCIN)/State Research Agency (SRA). The work of Jorge Pedro Galve was also funded by the 'Ramón y Cajal' Programme (RYC-2017-23335) of the Spanish Ministry of Science. Access to the Geohazard Exploitation Platform (GEP) of the European Space Agency (ESA) was provided by the NoR Projects Sponsorship (Project ID: 63737). Open access charge was funded by the University of Granada/CBUA.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.