**Geohazards: Risk Assessment, Mitigation and Prevention**

Editors

**Ricardo Castedo Miguel Llorente Isidro David Moncoulon**

Basel ' Beijing ' Wuhan ' Barcelona ' Belgrade ' Novi Sad ' Cluj ' Manchester

*Editors* Ricardo Castedo Department of Geological and Mining Engineering Universidad Politecnica ´ de Madrid Madrid Spain

Miguel Llorente Isidro Department of Climate Change and Geological Risks Spanish Geological Survey Madrid Spain

David Moncoulon R&D and Modeling Caisse Centrale de Reassurance (CCR Group) ´ Paris France

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: www.mdpi.com/journal/applsci/special issues/ geohazards risk).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

Lastname, A.A.; Lastname, B.B. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-9039-4 (Hbk) ISBN 978-3-0365-9038-7 (PDF) doi.org/10.3390/books978-3-0365-9038-7**

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license. The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons Attribution-NonCommercial-NoDerivs (CC BY-NC-ND) license.

## **Contents**



Reprinted from: *Appl. Sci.* **2023**, *13*, 5295, doi:10.3390/app13095295 . . . . . . . . . . . . . . . . . **400**


Reprinted from: *Appl. Sci.* **2023**, *13*, 7551, doi:10.3390/app13137551 . . . . . . . . . . . . . . . . . **572**

## **About the Editors**

## **Ricardo Castedo**

Ricardo Castedo is an Associate Professor at the Geological and Mining Engineering Department (School of Mines and Energy; Universidad Politecnica de Madrid). He obtained his B.Sc. (2008) in ´ Geological Engineering, M.Sc. (2014) in Teaching and Education, and Ph.D. (2012) in Engineering (Universidad Politecnica de Madrid). He had a one-year research stay at New Mexico Tech (USA) ´ and three research stays at the University of Leeds (UK), the University of Windsor (Canada), and Politecnico di Torino (Italia). He is the author of more than 35 publications in scientific journals and more than 40 contributions to national and international conferences. Reviewer and editor of scientifically indexed (JCR) journals and international conferences. His research interests include rock mechanics, geomorphology, modeling, the dynamic behavior of structures, safety and security engineering, testing in structures with explosives, explosives behavior and characterization, and the safety and modeling of explosives.

## **Miguel Llorente Isidro**

Miguel Llorente is Head of Scientific and Technical Projects in the Department of Climate Change and Geological Risks (Spanish Geological Survey, IGME-CSIC). He obtained his B.Sc. (2004) in Geology (University of Salamanca), M.Sc. (2005) in General and Applied Hydrology, and Ph.D. (2015) in Engineering (Universidad Politecnica de Madrid). He has two decades of work experience at the ´ Head Research Unit in Galicia, as a Head of the Group in Applied Geosciences Research Unit, Director of Geology with the Camara Minera del Peru, EU Expert for Geohazards, and as a Tsunami Expert ´ for the National Plan on Volcanic, Seismic, and Other Geophysical Process Monitoring. He is the author of more than 70 publications and contributions to national and international conferences and journals. He is also an Ad Honorem Profesor (Universidad Politecnica de Madrid). He is a reviewer ´ and editor of scientifically indexed (JCR) journals and international conferences. His research interests include geohazards, geological risks, numerical modeling, floods, tsunamis, earthquakes, volcanoes, engineering, mapping, and GIS.

## **David Moncoulon**

David Moncoulon was head of R&D and modeling at Caisse Centrale de Reassurance (CCR ´ Group) until his death in July 2023, when he was 50 years old. He obtained his B.Sc. (1998) in Agronomical Engineering at L' Institut Agro Montpellier and his Ph.D. (2014) in hydrology and risk modeling at the Universite Paul Sabatier Toulouse III. He worked three years at SEM Agriculture, ´ four years at the French National Centre for Scientific Research (CNRS), and more than fifteen years at CCR. He is the author of several publications in scientific journals and at national and international conferences. His research interests include hydrological modeling, agricultural yield loss modeling, damage modeling, climate modeling, flood modeling, prevention, risk management, and insurance.

## **Preface**

May this preface be a tribute to our colleague Dr. David Moncoulon, our fellow co-editor for this Special Issue (SI). He passed away too soon, too young, and still had much to contribute to research and development. We join his family, friends, and colleagues in mourning his loss.

This SI aims to provide a common forum to share knowledge about the growing threat of geohazards such as landslides, earthquakes, coastal changes, and drought, among many others, due to the expanding size of cities and urban areas, the critical use of agricultural lands, and the effects of climate change. A total of thirty-one papers (twenty-eight research papers, one review paper, and one technical note) on different topics like coastal changes, landslides, earthquakes, and others are presented in this reprint. The published works have been carried out in different parts of the world, such as Europe (i.e., Spain, France, Italy, etc.), Asia (i.e., China, South Korea, India, etc.), or America (i.e., Canada and Mexico).

Among the diverse research that the reader will find in this SI related to climate change and human intervention over coastal areas are the works of Fernandez-Hern ´ andez et al. [1], Bio et al. [2], ´ or Oh et al. [3]. These works use different data, such as aerial photographs (in some cases, more than 90 years ago), digital elevation models (DEMs), and video monitoring systems.

The contributions to this SI are mainly related to landslides (Saha et al. [4], Affandi et al. [5], Leonardi et al. [6], Gao and Zhang [7], Wu et al. [8], Liu et al. [9], Franco et al. [10], Cobos et al. [11], Niu et al. [12], and Gutierrez-Mart ´ ´ın et al. [13]), especially in populated areas, and focus on the development of prediction models, susceptibility maps, analysis of critical factors (e.g., rainfall, soil types, slope, etc.) or remediation techniques. These works include the use of geographic information systems, geotechnical flow modeling software, or electrical resistivity tomography.

Additionally, the contributions on earthquakes or seismic risks, including quantification of effects, stochastic generators, or analysis for insurance purposes (Cao et al. [14], Szabo et al. [15], Fidani [16], Gouache et al. [17], Hui et al. [18], Issadi et al. [19], Issadi et al. [20], and Perez-Moreno et ´ al. [21]), are noteworthy.

The reader will also find different laboratory studies on the influence of rainfall on slope erosion, the effect of dry-wet cycles on loess, or debris flow mobility (Tian et al. [22], Liu et al. [23], and Chang et al. [24]). To finish, readers will find articles on the effects of extreme droughts, forest fires, floods, or marine gas hydrate (Kapsambelis et al. [25], Gualdi et al. [26], Garrote et al. [27], Lee and Kim [28], Leon et al. [29], Jami et al. [30], and Li et al. [31]). ´

Finally, the editors would like to thank all the authors and reviewers, as well as the MDPI staff (especially Amy An), for their valuable contributions to this reprint.

## **Ricardo Castedo, Miguel Llorente Isidro, and David Moncoulon**

*Editors*

## *Article* **Anthropic Action on Historical Shoreline Changes and Future Estimates Using GIS: Guadarmar Del Segura (Spain)**

**Marta Fernández-Hernández 1,\* , Almudena Calvo <sup>2</sup> , Luis Iglesias <sup>1</sup> , Ricardo Castedo <sup>1</sup> , Jose J. Ortega <sup>1</sup> , Antonio J. Diaz-Honrubia <sup>3</sup> , Pedro Mora <sup>1</sup> and Elisa Costamagna <sup>4</sup>**


**Abstract:** A good understanding of historical change rates is a key requirement for effective coastal zone management and reliable predictions of shoreline evolution. Historical shoreline erosion for the coast of Guardamar del Segura (Alicante, Spain) is analyzed based on aerial photographs dating from 1930 to 2022 using the Digital Shoreline Analysis System (DSAS). This area is of special interest because the construction of a breakwater in the 1990s, which channels the mouth of the Segura River, has caused a change in coastal behavior. The prediction of future shorelines is conducted up to the year 2040 using two models based on data analysis techniques: the extrapolation of historical data (including the uncertainty of the historical measurements) and the Bruun-type model (considering the effect of sea level rises). The extrapolation of the natural erosion of the area up to 1989 is also compared with the reality, already affected by anthropic actions, in the years 2005 and 2022. The construction of the breakwater has accelerated the erosion along the coast downstream of this infrastructure by about 260%, endangering several houses that are located on the beach itself. The estimation models predict transects with erosions ranging from centimeters (±70 cm) to tens of meters (±30 m). However, both models are often overlapping, which gives a band where the shoreline may be thought to be in the future. The extrapolation of erosion up to 1989, and its subsequent comparison, shows that in most of the study areas, anthropic actions have increased erosion, reaching values of more than 35 m of shoreline loss. The effect of anthropic actions on the coast is also analyzed on the housing on the beach of Babilonia, which has lost around 17% of its built-up area in 40 years. This work demonstrates the importance of historical analysis and predictions before making any significant changes in coastal areas to develop sustainable plans for coastal area management.

**Keywords:** climate change; beach erosion; dune cliff; land management; data analysis

## **1. Introduction**

Coasts are valuable and vulnerable environments that provide numerous ecosystem services and support human populations in different ways [1]. It is estimated that sandy beaches make up approximately one-half of the world's ice-free coastline [2]. According to the United Nations, approximately 10% of the world's population currently lives in coastal areas that are less than 10 m above sea level, and about 40% lives within 100 km of the coast [3]. These areas are residentially, touristically and economically attractive; thus, they are becoming densely populated, but the concentration of infrastructures due to intense anthropic activity (i.e., ports, cities and resorts) puts additional stress on the coastal environment, sometimes pressing it to its limits [4].

**Citation:** Fernández-Hernández, M.; Calvo, A.; Iglesias, L.; Castedo, R.; Ortega, J.J.; Diaz-Honrubia, A.J.; Mora, P.; Costamagna, E. Anthropic Action on Historical Shoreline Changes and Future Estimates Using GIS: Guadarmar Del Segura (Spain). *Appl. Sci.* **2023**, *13*, 9792. https:// doi.org/10.3390/app13179792

Academic Editor: Athanasios Sfetsos

Received: 2 August 2023 Revised: 18 August 2023 Accepted: 28 August 2023 Published: 30 August 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Coastal zones are particularly vulnerable to natural erosive processes due to factors like wave action or storm events [5]. These processes can lead to the loss of land, degradation of coastal ecosystems, and increased risks for human settlements and infrastructure. The problem is emphasized when anthropic constructions modify coastal dynamics in some way, even though their primary objective is to protect the coast from erosion due to natural phenomena [6–8]. While these structures (i.e., harbors, channels or breakwaters) may provide immediate protection against erosion, they can disrupt the natural movement of sediment along the coastline, leading to unintended consequences [9,10]. Many of these hard structures were built from the 1950s to the 2000s, without accurate studies of the effects on coastal sediment dynamics and flora and fauna [1]. In contrast, in the last two or three decades, management strategies have focused on soft methods such as artificial beach nourishment, revegetation of dunes or bioengineering [11,12]. For example, some researchers have found a relationship between increased shoreline erosion in areas around the mouths of dammed rivers and in areas with high population densities in the Gulf of California [13]. Research conducted on the East Coast of South Korea found that the artificial structures constructed along the coast have not completely solved or stopped erosion but shifted it from one location to another [14]. On the contrary, work by Skilodimou et al. [15] shows how coastlines have been enhanced by 40% over the last 76 years due to human interventions, in this case, earth filling.

The effects of climate change on coastal regions are critical due to their exposure to rising sea levels, increased storm intensity and changes in oceanic and atmospheric conditions [16]. In general, rising sea levels can be a significant factor in shoreline changes [17]. Higher water levels mean that waves and storm surges can reach farther inland, eroding coastlines and flooding low-lying areas. This can result in the loss of beaches, dunes and other coastal landforms, which can lead to feedback loops, e.g., as coastal areas erode, sediment may be lost, reducing natural coastal barriers and exacerbating erosion. This, in turn, can increase vulnerability to coastal hazards. Changes in ocean wave characteristics, including alterations in wave height, period and direction, can further impact coastal erosion, sediment transport and stability. In addition to this, storms, including hurricanes and tropical cyclones, can cause significant changes due to intense erosion and the transportation of large amounts of sediment along the coastline [18]. An example is the work of Simeone et al. [19], who demonstrate that the beaches studied in Western Sardinia (Italy) show greater changes when faced with consecutive storms and not with an isolated one. These climate change impacts can lead to the loss of coastal land, damage to infrastructure and the displacement of communities. For these reasons, the analysis of the long-term effects of anthropic actions in past cases, in combination with information on the adverse effects of climate change (both at the global and regional levels), is fundamental to ensure the long-term sustainability and resilience of the new actions to be taken.

The aims of this work carried out in Guardamar del Segura, Alicante (Spain) are as follows: (1) analyze the historical evolution of the shoreline from 1930 to 2022, taking into account different anthropic constructions; (2) study the effects of these constructions on the erosion rates of the shoreline and a cliff dune; (3) evaluate the loss of built surface (houses) due to changes in coastal dynamics; (4) estimate the shoreline situation for the year 2040 with the information from available data and the use of two prediction models based on historical recession rates and changes in sea level; and (5) compare, with the use of these computational models, the positions that the coastline would have had in the years 2005 and 2022 without the constructions carried out in the area. The results and methodology provided by this work can be used as the basis for coastal planning and management, especially in areas with key infrastructures.

## **2. Study Area**

This work focuses on the coastline of Guardamar del Segura located in the southeast of Alicante, Spain, on the Mediterranean coast (see Figure 1). The area has a population of more than 16,000 inhabitants, and the main economic activity in the region is summer

tourism. The studied coast includes the beaches of Los Viveros (corresponding to zone B in Figure 1c) and Babilonia (zone C in Figure 1c). These two beaches have an actual length of 2.7 km and an average width of 12 m, reaching a maximum of 35 m. They are open beaches of golden sand with a grain size of around 0.27 mm [20], which is between fine and coarse sand. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 4 of 20

**Figure 1.** (**a**) Location of the study area. (**b**) Orthomosaic of the coast in 1930 and details of the transects used in the study. (**c**) Orthophoto of the coast in 2022. Note that the coordinates used are UTM ETRS89 H30N. **Figure 1.** (**a**) Location of the study area. (**b**) Orthomosaic of the coast in 1930 and details of the transects used in the study. (**c**) Orthophoto of the coast in 2022. Note that the coordinates used are UTM ETRS89 H30N.

Guardamar del Segura coast is a fetch-limited environment with a predominance of sea waves. This coast is an environment with small astronomical tides (as is common in the Mediterranean) that oscillate about 0.3 m (SIMAR data from Puertos del Estado [21]). However, the so-called meteorological tides can reach up to 0.45 m. The dominant wave direction is northeast, creating a net N-S-oriented coastal current. Significant wave heights in this area, from 1958 to 2023, are 50.80% of the time between 0 and 0.5 m, 37.52% between 0.5 and 1 m, 8.7% between 1.0 and 1.5 m, 2.06% between 1.5 and 2.0 m and the remaining 0.92% between 2.0 and 3.5 m [21]. It should be noted that according to the available data [21], the maximum monthly values up to 1970 did not exceed 3 m of significant height. In the 1980s, there were already some peaks that exceeded 4 m, a trend that remained constant with peaks between 2.5 m and 4 m until 2010, but with an average value of around 2 m. However, in the last 13 years, the average has risen around 3 m with some peaks in 2019 of almost 7 m. The latter peak is due to the large Isolated Depression at High Levels (DANA in Spanish) that occurred in the area in 2019 [22]. As for the peak period, the histogram is centered on the 5–6 s interval (23.94% of the time), distributed on the 2–4 s interval (29.12%) and then on the 4–5 s interval (19.59%), 6–7 s interval (17.89%) and 7–10 s interval (17.43%).

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 4 of 20

According to the latest data published by the Intergovernmental Panel on Climate Change (IPCC), the trend of accelerating sea level rise in the Mediterranean is consistent [23]. However, there are important differences depending on the methods and time horizons used in the analyses. In general, the accepted sea level rise referred to in the 20th century was 1.4 <sup>±</sup> 0.2 mm yr−<sup>1</sup> [24]. However, by 2050, it is estimated that the sea level may rise between 0.52 and 1.22 m above the mean relative sea level between 1996 and 2014. The relative sea level rise from 1927 to 2012 in Alicante was 1.28 <sup>±</sup> 0.5 mm yr−<sup>1</sup> [25].

This coastline has a long history of anthropogenic actions since the 18th century. The dune system in this area was left unfixed due to the massive felling of trees for the fishing industry during the 18th century. This meant that for a long time, the town of Guardamar was "threatened" by the movements of the dunes. From 1900 to 1930, the engineer Mira y Botella directed the reforestation of the area, which reduced the mobility of the dunes and, therefore, prevented the burial of the town. The houses of Babilonia (zone C in Figure 1) did not respond to the predatory model of the urbanism of the second half of the last century around the Spanish Levante. On the contrary, they had an environmental function and fit perfectly with the coastal physiognomy of the municipality in the 1930s. The number of houses increased until the end of the 1970s and remained constant until the 2000s, but since then, houses started to naturally and progressively collapse due to marine action. See Figure 2 for some photographs of the area since 2010. **Figure 1.** (**a**) Location of the study area. (**b**) Orthomosaic of the coast in 1930 and details of the transects used in the study. (**c**) Orthophoto of the coast in 2022. Note that the coordinates used are UTM ETRS89 H30N.

**Figure 2.** (**a**) View from north of Babilonia Beach on 2 April 2010. (**b**) View of Los Viveros artificial dune on 2 April 2010. (**c**) View from north of Babilonia Beach on 5 May 2021. (**d**) View from north of Los Viveros Beach and artificial dune on 13 August 2021.

The main period of anthropic actions began in the mid-1970s and lasted about 20/25 years. In the 1970s, a small groin was built in the southern part of the mouth of the Segura River, which was completed in the mid-1980s with the construction of a larger groin in the northern part. In the period between 1987 and 1988, a large dredging of the mouth of the Segura River was carried out. Between 1990 and 1994, to reduce the risk of flooding, its mouth was channeled by means of a 525 m breakwater with an east and northeast orientation (zone A in Figure 1c), the opposite of all those built in the Spanish Levante [26]. This orientation stops the longitudinal transport of sediments from the north and blocks the outflow of sediments transported by the Segura River [7]. This forces the "Confederación Hidrográfica del Segura (CHS)" to schedule periodic dredging due to sediment accumulation that, if not executed, could cause breakage of the breakwater [27]. A volume of 14,000 m<sup>3</sup> of material extracted in the works of the early 1990s was used to raise Viveros beach (zone B in Figure 1c), and the rest of the material was used to raise other beaches to the south of the area of this study. Between 1996 and 1999, the Guardamar marina was built. In this project, a volume of sludge of approximately 200,000 m<sup>3</sup> was obtained, basically sand, with a low proportion of fines, very similar to that which makes up the current dune system. This sediment was dumped on a foredune originating an artificial dune (with the appearance of a small cliff) about 500 m long and of variable height (zone D in Figure 1) with a trapezoidal section, ranging from 8 m in the sector closest to the riverbed to 4.5 m in the southern part [6,28]. From 2002 to 2011, an attempt was made to eliminate the dumped anthropogenic materials, leaving the existing dune sands uncovered, and native species were planted to fix them [29]. Currently, this small cliff dune has areas up to 10 m high and slopes of up to >45◦ [30].

## **3. Methodology**

The explanation of the methodology is divided into three blocks: shoreline evolution data, historical data analysis using the Digital Shoreline Analysis System (DSAS) [31] and prediction of future shorelines.

## *3.1. Dataset Preparation*

Twenty aerial images covering a period of 93 years were used, with images in the following years: 1930, 1946, 1956, 1977, 1985, 1989, 1997, 1999, 2003, 2005, 2007, 2009, 2012, 2014, 2017, 2018, 2019, 2020, 2021 and 2022 (see Table 1 for details). They were obtained from different sources: National Geographic Institute (Instituto Geográfico Nacional (IGN) [32]) and Valencian Cartographic Institute (Institut Cartogràfic Valencià (ICV) [33]) photo libraries. The used reference system is ETRS89, UTM projection zone 30. In all images, the sea is calm, and therefore, all available data can be used for a comparative study.

**Table 1.** Main information from the aerial imagery used in this research. Note that GSD is the ground sampling distance: the distance between two consecutive pixel centers measured on the ground. This concept is equivalent to spatial resolution in remote sensing.


The photographs of the first flight used (Ruiz de Alda 1929–1930) were treated with AgiSoft Metashape to obtain an orthomosaic of the study area since the photos were not taken systematically [34]. Control points identified on the 2022 orthomosaic were used to make the 1929–1930 orthomosaic. It seems clear that, although the 1930 orthomosaic has an error (Table 1), it can be considered within acceptable limits when working with data 50 years old or more [35], especially when you consider that it provides an important piece of information that is worth considering. The 1946, 1956, 1977, 1985, 1989, 1999 and 2017 aerial photographs were rectified using easily identified points on the 2022 orthomosaic, which was taken as a reference for the entire work. The remaining used images were orthomosaics produced by the official cartographic agencies and downloaded from their data servers.

Once the 20 mosaics with the orthophotos of all the years of study were uploaded to GIS, the next step was the vectorization of the shorelines, the cliff top and the houses of Babilonia Beach. The shoreline was drawn manually by the same operator instead of using automatic techniques more commonly used in remote sensing [36], following the last wet tide mark on the beach profile (see Figure 3 as an example). This line defines the boundary with the backshore and can be visually identified in orthophotos [37,38]. This methodology can be used in sedimentary littoral formations exposed to the open sea (i.e., beaches). To profile the cliff dune, the top (crest) of the slope is used, which is a visually perceptible feature [39]. The calculation of the number of dwellings and occupied area is based on the most current 2011 cadastral data [40]. On these cadastral data, the remaining 19 layers were manually modified to adapt them to each orthophoto from 1930 to 2022. All this work was carried out by the same operator to ensure a uniform criterion in the choice of geomorphological features. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 7 of 20

**Figure 3.** (**a**) Orthomosaic of the coast in 1930. (**b**) Orthomosaic of the coast in 1930 with the drawing of the shoreline. (**c**) Detail of the shoreline in the area of transects 4 and 8 (north of zone B in Figure 1). (**d**) Detail of the shoreline in the area of transects 25 and 29 (south of zone C in Figure 1). **Figure 3.** (**a**) Orthomosaic of the coast in 1930. (**b**) Orthomosaic of the coast in 1930 with the drawing of the shoreline. (**c**) Detail of the shoreline in the area of transects 4 and 8 (north of zone B in Figure 1). (**d**) Detail of the shoreline in the area of transects 25 and 29 (south of zone C in Figure 1).

Digital Shoreline Analysis System (DSAS) version 5.1 [31] was used to analyze the

The first analysis of the historical data was performed considering the complete time series from 1930 to 2022, except for zones A and D. In zone A, the analysis was conducted up to 1989, as this area was removed during the construction of the north breakwater and the channeling of the Segura River. Zone D was analyzed from 1999 to 2022 after human intervention. In addition, to improve the analysis and interpretation of the results and consider the important anthropic actions in zones B and C, an analysis of two-time intervals was also carried out: 1930–1989 (before the channeling of the Segura River) and 1997–

The classical linear regression model was selected to estimate the rates of shoreline change. This technique is based on accepted statistical concepts, includes all the data of the series and provides the necessary data for the prediction models used in this research. In detail, the data provided by the DSAS software with this calculation method are as

three parts, in addition to the cliff area (see Figure 1c). In total, there are 40 transects, 3 of which are in zone A (numbered 1 to 3), 14 are in zone B (numbered 4 to 17), 12 are in zone C (numbered 18 to 29) and 11 are in the zone D or cliff zone (numbered 1d to 11d). The transects are spaced 50 m apart in zones A and D, while in zones B and C, they were spaced 100 m apart. Each transect is a line perpendicular to the coastline, which joins the intersection points of the coastline of each of the available images (dates) with said per-

2022 (after the construction of the marina and artificial dune).

*3.2. DSAS—Historical Shoreline Analysis* 

pendicular line.

follows:

## *3.2. DSAS—Historical Shoreline Analysis*

Digital Shoreline Analysis System (DSAS) version 5.1 [31] was used to analyze the historical evolution of the coastline. For its analysis, the study region was divided into three parts, in addition to the cliff area (see Figure 1c). In total, there are 40 transects, 3 of which are in zone A (numbered 1 to 3), 14 are in zone B (numbered 4 to 17), 12 are in zone C (numbered 18 to 29) and 11 are in the zone D or cliff zone (numbered 1d to 11d). The transects are spaced 50 m apart in zones A and D, while in zones B and C, they were spaced 100 m apart. Each transect is a line perpendicular to the coastline, which joins the intersection points of the coastline of each of the available images (dates) with said perpendicular line.

The first analysis of the historical data was performed considering the complete time series from 1930 to 2022, except for zones A and D. In zone A, the analysis was conducted up to 1989, as this area was removed during the construction of the north breakwater and the channeling of the Segura River. Zone D was analyzed from 1999 to 2022 after human intervention. In addition, to improve the analysis and interpretation of the results and consider the important anthropic actions in zones B and C, an analysis of two-time intervals was also carried out: 1930–1989 (before the channeling of the Segura River) and 1997–2022 (after the construction of the marina and artificial dune).

The classical linear regression model was selected to estimate the rates of shoreline change. This technique is based on accepted statistical concepts, includes all the data of the series and provides the necessary data for the prediction models used in this research. In detail, the data provided by the DSAS software with this calculation method are as follows:


## *3.3. Prediction of Future Shoreline*

While it is challenging to make precise projections on future shorelines, researchers use available data and models to understand potential scenarios [41,42]. Calculated data on historical shoreline changes generally serve as the first, and perhaps most important, input for prediction models. However, these data alone are valid for extrapolation, if in the future it is assumed that the system remains in the same dynamic equilibrium as in the historical period analyzed. Otherwise, some variability must be added to the historical data (i.e., standard deviation) or, in more complex models, the elements in each zone that are most likely to be affected by future changes must be included.

In this paper, two simple but well-known models were used to estimate the position of the shoreline in the future: Lee–Clark [43,44] and Leatherman [45–47] models. The Lee–Clark model is based on the extrapolation of historical data (LRR). In this model, the standard deviation (LCI) of the measured erosions is incorporated to consider the variability in the rates over time. Erosion calculation is simple, and the result in year *T* would be: LCe = (*LRR* ± LCI) × *T* years. The Leatherman model (also known as historical trend analysis) is basically based on the "Bruun rule" to estimate the potential shoreline retreat (*R*<sup>2</sup> measured in m/yr) resulting from a rise in sea level as: *R*<sup>2</sup> = *S*2(*LRR*/*S*1), where *S*<sup>1</sup> is

the sea level rise rate (m/yr) during the analyzed period and *S*<sup>2</sup> is the projected sea level rise rate (m/yr) in the study area.

## **4. Results and Discussion**

The results and their discussion were structured by zones, starting with the evolution of the shoreline, then the cliff dune and ending with a section dedicated to the comparison of the coastline in 2005 and 2022 with and without the effects of anthropogenic actions.

## *4.1. Shoreline Evolution (Past and Future)*

Table 2 shows the rates of coastline change for the three periods analyzed (1930–1989; 1997–2022; and 1930–2022) in the three areas (A, B and C) that have or have had shorelines. In zone A, between 1930 and 1989, the erosion is important, with an average rate in the zone of 1.13 m/yr (Figure 4). As can be seen in Figure 5, most of this erosion is due to the period 1930–1956. Even in transect 1, the 1977 shoreline was ahead of the 1956 shoreline, reflecting a small accretion in the area, probably due to the construction of the small groin on the south side of the river mouth. It seems evident that in this area, there was some erosion near the natural mouth of the Segura River before it was channeled.

**Table 2.** Erosion rates for profiles 1 to 29, including zone A, B and C for three different study periods: 1930–1989, 1997–2023, 1930–2022. Negative sign means erosion, while positive sign means accretion for NSM and LRR.


29 −7.77 −0.08 2.00 0.08 −11.39 −0.27 0.71 0.43 −28.25 −0.17 0.17 0.34

29 −7.77 −0.08 2.00 0.08 −11.39 −0.27 0.71 0.43 −28.25 −0.17 0.17 0.34

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 10 of 20

**Figure 4.** The upper graph shows the NSM of all transects in zones A, B and C in the periods 1930– 1989 and 1997–2022. The lower graph shows the LRR values of each transect for the same zones and periods. **Figure 4.** The upper graph shows the NSM of all transects in zones A, B and C in the periods 1930–1989 and 1997–2022. The lower graph shows the LRR values of each transect for the same zones and periods. **Figure 4.** The upper graph shows the NSM of all transects in zones A, B and C in the periods 1930– 1989 and 1997–2022. The lower graph shows the LRR values of each transect for the same zones and periods.

**Figure 5.** (**a**) Shorelines and transects 1 to 3 from 1930 to 1989 (background image from 1989). (**b**) Image from 2022 of the same area. **Figure 5.** (**a**) Shorelines and transects 1 to 3 from 1930 to 1989 (background image from 1989). (**b**) Image from 2022 of the same area. **Figure 5.** (**a**) Shorelines and transects 1 to 3 from 1930 to 1989 (background image from 1989). (**b**) Image from 2022 of the same area.

In the northern area of zone B, near the breakwater (transects 4–8, Table 2 and Figure 6), the erosion during the period 1930–1989 is significant, with an average of 0.89 m/yr. It is worth remembering that between 1990 and 1992, the beach was regrown by about 50 m due to the creation of the marina [4]. As shown in Figure 4, this material added to the beach eroded rapidly, as the 1997 shoreline was behind the 1989 one. However, for the period 1997–2022, the values shot up to an average of 1.29 m/yr. It seems more than evident that the construction of the breakwater and the channelization of the river has accelerated a natural process (Figure 4). From transects 9 to 13, erosion rates only reach 0.36 m/yr for the period 1930–1989. This indicates that this zone, already farther from the mouth of the river, was less affected by the natural shoreline erosion. In contrast, after construction works in the area, the rate skyrockets to 1.22 m/yr, almost one meter more in a much shorter period of only 25 years. The erosion data for the area between transects 14 and 17 present a low LR2 for the first epoch analyzed. This is because the shorelines in this area are very

close together, having accretion and erosion epochs, which makes it difficult to fit the linear regression [48]. In any case, the net movement of shoreline erosion in these transects was 5 m for the 60-year period (1930–1989). However, for the period 1997–2022, the mean NSM for these transects is 31.71 m, which seems to confirm the increase in erosion south of the breakwater. Looking at the entire period of 92 years, the mean erosion of all transects in zone B is 0.92 ± 0.24 m/yr. It happens that throughout zone B, erosion rates have increased after anthropic actions, reflecting their negative effects, especially in the southern part (Figure 4). *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 11 of 20

**Figure 6.** Zone B. (**a**,**d**) Shorelines and transects from 1930 to 1989 (background image from 1989). (**b**,**e**) Shorelines and transects from 1997 to 2022 (background image from 2022). (**c**,**f**) Shorelines and transects from 1930 to 2022 (background image from 2022). **Figure 6.** Zone B. (**a**,**d**) Shorelines and transects from 1930 to 1989 (background image from 1989). (**b**,**e**) Shorelines and transects from 1997 to 2022 (background image from 2022). (**c**,**f**) Shorelines and transects from 1930 to 2022 (background image from 2022).

In zone C, for the period 1930–1989, something similar to zone B occurs with some of the linear regression, where the LR2 value is low. In this zone, some of the transects show a net accretion for this period (see Figure 4 and Figure 7a). Erosion at Babilonia Beach was very small until 1989, with only 8.72 m in 60 years, whereas in the last 25 years, 12.09 m, about 40% more, eroded in less than half the time. This has caused the loss of beaches in front of houses and, therefore, has left them directly exposed to sea waves (see Figure 2). However, as can be seen in Figure 7, the houses of Babilonia Beach act as an artificial barrier (or artificial shoreline), and therefore, there can be no further erosion except for the destruction of the houses themselves. If zone C is analyzed for the entire available time set, it can be observed that the fits are better and that the NSM value amounts to 33.86 m (0.32 m/yr). As in zone B, Figure 4 shows how erosion rates have increased after anthro-In zone C, for the period 1930–1989, something similar to zone B occurs with some of the linear regression, where the LR2 value is low. In this zone, some of the transects show a net accretion for this period (see Figures 4 and 7a). Erosion at Babilonia Beach was very small until 1989, with only 8.72 m in 60 years, whereas in the last 25 years, 12.09 m, about 40% more, eroded in less than half the time. This has caused the loss of beaches in front of houses and, therefore, has left them directly exposed to sea waves (see Figure 2). However, as can be seen in Figure 7, the houses of Babilonia Beach act as an artificial barrier (or artificial shoreline), and therefore, there can be no further erosion except for the destruction of the houses themselves. If zone C is analyzed for the entire available time set, it can be observed that the fits are better and that the NSM value amounts to 33.86 m (0.32 m/yr). As in zone B, Figure 4 shows how erosion rates have increased after anthropogenic actions.

pogenic actions. The effect of erosion on the houses of Babilonia Beach was also analyzed. During the 1930s and up to the 1970s, the built-up area increased, with 8533 m<sup>2</sup> (55 houses) in 1930; 12,277 m<sup>2</sup> (80 houses) in 1946; 17,976 m<sup>2</sup> (117 houses) in 1956; and 19,463 m<sup>2</sup> (126 houses) in 1977. This last value remained stable until 1999. From then on, and as a natural response to the values already discussed, the built-up area has been in continuous decline. Until 2005 the loss of floor area was a mere 131 m<sup>2</sup> (one house). But from 2005 to 2012, the loss increased by 2218 m<sup>2</sup> (14 houses) and has not stopped until 2022 (last data available), with a loss of another 1003 m<sup>2</sup> (6 houses). One should consider the numerous attempts to protect these houses with the use of rock dikes, Hesco barriers filled with rock and sand, etc., that have managed to somewhat slow down the destruction of these houses [49].

**Figure 7.** Zone C. (**a**,**d**) Shorelines and transects from 1930 to 1989 (background image from 1989). (**b**,**e**) Shorelines and transects from 1997 to 2022 (background image from 2022). (**c**,**f**) Shorelines and transects from 1930 to 2022 (background image from 2022). **Figure 7.** Zone C. (**a**,**d**) Shorelines and transects from 1930 to 1989 (background image from 1989). (**b**,**e**) Shorelines and transects from 1997 to 2022 (background image from 2022). (**c**,**f**) Shorelines and transects from 1930 to 2022 (background image from 2022).

The effect of erosion on the houses of Babilonia Beach was also analyzed. During the 1930s and up to the 1970s, the built-up area increased, with 8533 m2 (55 houses) in 1930; 12,277 m2 (80 houses) in 1946; 17,976 m2 (117 houses) in 1956; and 19,463 m2 (126 houses) in 1977. This last value remained stable until 1999. From then on, and as a natural response to the values already discussed, the built-up area has been in continuous decline. Until 2005 the loss of floor area was a mere 131 m2 (one house). But from 2005 to 2012, the loss increased by 2218 m2 (14 houses) and has not stopped until 2022 (last data available), with a loss of another 1003 m2 (6 houses). One should consider the numerous attempts to protect these houses with the use of rock dikes, Hesco barriers filled with rock and sand, etc., that have managed to somewhat slow down the destruction of these houses [49]. Having historical erosion data in different transects, it is possible to apply the prediction models explained in Section 3.3. An adequate time horizon, preferably shorter than the available observation period (92 years), should be chosen to consider the data variability and, at the same time, to reduce the impact of data uncertainty on the quality of the Having historical erosion data in different transects, it is possible to apply the prediction models explained in Section 3.3. An adequate time horizon, preferably shorter than the available observation period (92 years), should be chosen to consider the data variability and, at the same time, to reduce the impact of data uncertainty on the quality of the results [39]. Therefore, the simulation is extended for 18 years, up to the year 2040 (see Table 3). Uncertainties have been incorporated for each model, including lower/upper bounds, having a surface area that delimits the possible solution (Figure 8). In the Lee– Clark model, uncertainty in the calculation of future erosion (LCe) is incorporated by adding (upper bound) and/or subtracting (lower bound) the LCI to the LRR value for each transect (see Table 2). For the Leatherman model, the limits are introduced using a minimum and maximum estimate of the rate of sea level rise in the future (*S*2). According to local estimates, the future rates will be a maximum (upper bound) of 0.5 mm/yr, with a reasonable minimum (lower bound) value of 0.25 mm/yr [25]. As mentioned in Section 2, the historical relative sea level rise (*S*1) has been 1.28 mm/yr.

results [39]. Therefore, the simulation is extended for 18 years, up to the year 2040 (see Table 3). Uncertainties have been incorporated for each model, including lower/upper bounds, having a surface area that delimits the possible solution (Figure 8). In the Lee– Clark model, uncertainty in the calculation of future erosion (LCe) is incorporated by adding (upper bound) and/or subtracting (lower bound) the LCI to the LRR value for each transect (see Table 2). For the Leatherman model, the limits are introduced using a minimum and maximum estimate of the rate of sea level rise in the future (*S*2). According to local estimates, the future rates will be a maximum (upper bound) of 0.5 mm/yr, with a reasonable minimum (lower bound) value of 0.25 mm/yr [25]. As mentioned in Section 2, As can be seen in Figure 8 and Table 3, the prediction models overlap in most cases. In this situation, it can be said that the shoreline is most likely to be located at some intersection surface of the two models. The Lee and Clark model has a larger difference between the lower and upper limits than the Leatherman. This is probably because the variability in historical erosion rates is greater than the effect of sea level changes for this relatively shortsimulation period. In the northern part of zone B, the shoreline is predicted to reach the area currently occupied by the cliff dune that, therefore, would also be eroded. In zone C, the houses continue to act as a barrier, although in the northern zone, erosion will increase and the loss of built-up area is likely to increase following the actual trend.

### the historical relative sea level rise (*S*1) has been 1.28 mm/yr. As can be seen in Figure 8 and Table 3, the prediction models overlap in most cases. *4.2. Cliff Dune Evolution (Past and Future)*

In this situation, it can be said that the shoreline is most likely to be located at some intersection surface of the two models. The Lee and Clark model has a larger difference between the lower and upper limits than the Leatherman. This is probably because the variability in historical erosion rates is greater than the effect of sea level changes for this relatively short simulation period. In the northern part of zone B, the shoreline is predicted The historical evolution of the cliff dune shows episodic erosion with major events between 2005 and 2007 and, more recently, between 2021 and 2022 (see Table 4 and Figure 9). In general, this is due to periods of major storms and an increase in wave height in the region [4,50]. In any case, the mean erosion value (including the 11 transects) for this area is 1.12 ± 0.66 m/yr. In this part of zone B, erosion data are 1.11 ± 0.20 m/yr for the period

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 13 of 20

1930–2022 and a little bit higher, 1.31 ± 0.77 m/yr, for the period 1999–2022. This shows, as is logical, that both cliff dunes and shorelines are closely related. In zone C, the houses continue to act as a barrier, although in the northern zone, erosion will increase and the loss of built-up area is likely to increase following the actual trend.

to reach the area currently occupied by the cliff dune that, therefore, would also be eroded.

**Figure 8.** Shoreline model predictions (Lth—Leatherman; LC—Lee and Clark) by year 2040: (**a**) general view of zone B; (**b**) detailed view of zone B; (**c**) detailed view of zone C; (**d**) general view for zone C. **Figure 8.** Shoreline model predictions (Lth—Leatherman; LC—Lee and Clark) by year 2040: (**a**) general view of zone B; (**b**) detailed view of zone B; (**c**) detailed view of zone C; (**d**) general view for zone C. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 15 of 20

16 −6.36 −16.81 −13.85 −16.11 17 −4.20 −15.12 −11.55 −13.44 **Figure 9.** (**a**) Cliff dune erosion since 1999, zone D. (**b**) Cliff dune erosion predictions (Lth—Leatherman; LC—Lee and Clark) by the year 2040. **Figure 9.** (**a**) Cliff dune erosion since 1999, zone D. (**b**) Cliff dune erosion predictions (Lth—Leatherman; LC—Lee and Clark) by the year 2040.

15 −6.53 −16.63 −13.84 −16.10

C 18 −1.34 −12.09 −8.03 −9.34 19 −5.58 −15.94 −12.86 −14.96 20 −2.59 −9.05 −6.95 −8.09 21 −2.96 −8.46 −6.83 −7.94 22 −3.30 −8.04 −6.77 −7.88 *4.3. What Would Have Happened in the Area without Large-Scale Anthropogenic Actions?*  With the data obtained in the previous sections, it is possible to simulate where the shoreline would be without the anthropic actions and compare it with the real situation. For this purpose, the Lee–Clark model is used in zones B and C for the pre-works period, i.e., from 1930 to 1989, as shown in Table 2. The behavior of the Lee–Clark and Leatherman model prediction results are very similar to those obtained for the shoreline, with greater variability in the LC model and the Leatherman model mostly bordering the upper bound of the LC model. With both models, the prediction for the year 2040 is that erosion will be such that the paved road that parallels part of the marina will be eroded.

In zone B (Figure 10), it can be seen how erosion has been advancing faster than natural erosion would have. This effect is especially noticeable in 2022 (Figure 10c,d), about

the shoreline), while in 2022, that difference increased to 20.85 m. However, this effect is felt more in the central sector of zone B; at transect 12, there was a difference of 23 m in 2005 and of 36 m in 2022. It appears evident that the effect of upstream constructions has accelerated the erosive process on the shoreline, especially in the central section of zone B

(Figure 10b,d).


**Table 3.** Total shoreline erosion by 2040 in meters.

**Table 4.** Erosion rates for profiles in zone D (cliff area) for the period: 1999-2022. Negative sign means erosion, while positive sign means accretion for NSM and LRR. Also, the total shoreline erosion with both models by 2040 in meters is shown.


## *4.3. What Would Have Happened in the Area without Large-Scale Anthropogenic Actions?*

With the data obtained in the previous sections, it is possible to simulate where the shoreline would be without the anthropic actions and compare it with the real situation. For this purpose, the Lee–Clark model is used in zones B and C for the pre-works period, i.e., from 1930 to 1989, as shown in Table 2.

In zone B (Figure 10), it can be seen how erosion has been advancing faster than natural erosion would have. This effect is especially noticeable in 2022 (Figure 10c,d), about 32 years after the start of the main works. At some points, such as transect 6, there was a difference of up to 10.85 m in 2005 (measured from the upper bound of the prediction to the shoreline), while in 2022, that difference increased to 20.85 m. However, this effect is felt more in the central sector of zone B; at transect 12, there was a difference of 23 m in 2005 and of 36 m in 2022. It appears evident that the effect of upstream constructions has accelerated the erosive process on the shoreline, especially in the central section of zone B (Figure 10b,d). *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 16 of 20

**Figure 10.** Shoreline and area of possible erosion predicted with the LC model: (**a**) north of zone B by 2005; (**b**) middle sector of zone B by 2005; (**c**) north of zone B by 2022; (**d**) middle sector of zone B by 2022. **Figure 10.** Shoreline and area of possible erosion predicted with the LC model: (**a**) north of zone B by 2005; (**b**) middle sector of zone B by 2005; (**c**) north of zone B by 2022; (**d**) middle sector of zone B by 2022.

In Babilonia Beach, the erosion behavior is strongly conditioned by houses. An important characteristic in zone C is that in 2005 (Figure 11a,b), the houses had an important beach area that protected them from the direct action of the sea. In the 2022 images (Figure 11c,d), this beach does not exist anymore. Figure 11a,c show how a significant number of buildings have been destroyed from 2005 to 2022 and how these houses have been replaced by sand. This has left the way open for erosion, which has been increasing in contrast to what the model predicts (transect 19). According to the projections of the model used, erosion in some areas would have reached the housing area even without the works carried out, but these have undoubtedly had an important effect on the acceleration of this In Babilonia Beach, the erosion behavior is strongly conditioned by houses. An important characteristic in zone C is that in 2005 (Figure 11a,b), the houses had an important beach area that protected them from the direct action of the sea. In the 2022 images (Figure 11c,d), this beach does not exist anymore. Figure 11a,c show how a significant number of buildings have been destroyed from 2005 to 2022 and how these houses have been replaced by sand. This has left the way open for erosion, which has been increasing in contrast to what the model predicts (transect 19). According to the projections of the model used, erosion in some areas would have reached the housing area even without the works carried out, but these have undoubtedly had an important effect on the acceleration of this phenomenon.

phenomenon.

**Figure 11.** Shoreline and area of possible erosion predicted with the LC model: (**a**) north of zone C by 2005; (**b**) middle sector of zone C by 2005; (**c**) north of zone C by 2022; (**d**) middle sector of zone C by 2022. **Figure 11.** Shoreline and area of possible erosion predicted with the LC model: (**a**) north of zone C by 2005; (**b**) middle sector of zone C by 2005; (**c**) north of zone C by 2022; (**d**) middle sector of zone C by 2022.

## **5. Conclusions**

**5. Conclusions**  In general, the results presented and the methodology shown can be used to improve short- and medium-term decisions in areas where major construction is to be carried out, as well as where and when to intervene in coastal erosion processes. The conclusions in In general, the results presented and the methodology shown can be used to improve short- and medium-term decisions in areas where major construction is to be carried out, as well as where and when to intervene in coastal erosion processes. The conclusions in response to the objectives set out in the study are as follows:


mum erosion limit of the former. In this area, the effect of the rate of sea level change

not seem to be very important in shoreline erosion since this value was already quite high during the last century.

5. Comparing the actual shorelines in the years 2005 and 2022 with those simulated using the Lee–Clark model, it can be said that anthropogenic actions exhibit their effect by increasing erosion rates. This effect is more noticeable in the central zone of this study (southern part of zone B and northern part of zone C) and more limited in the marina zone (northern part of zone B) and in the houses on Babilonia Beach (southern part of zone C).

**Author Contributions:** Conceptualization, M.F.-H. and R.C.; methodology, A.C. and R.C.; software and data-based model learning, L.I., A.C. and A.J.D.-H.; validation, L.I., A.C. and M.F.-H.; formal analysis, R.C., J.J.O. and E.C.; investigation, M.F.-H., A.C. and P.M.; data curation, A.C. and L.I.; writing—original draft preparation, E.C., M.F.-H., R.C., J.J.O. and A.J.D.-H.; writing—review and editing, M.F.-H., A.C., L.I. and R.C.; visualization, A.C. and M.F.-H.; supervision, M.F.-H.; funding acquisition, R.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been partially funded by the Universidad Politécnica de Madrid with grants to PDI and PhD researchers for international research stays of one month or more (resolution of 30 May 2022, of the Rector of the Universidad Politécnica de Madrid).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not available to the public, as some of the images used must be requested from the Cartographic and Photographic Center of the Spanish Air Force and cannot be disseminated without permission.

**Acknowledgments:** The authors of this work would like to thank all the institutions that provide access to their material through official websites. The authors would like to thank the two anonymous reviewers for their contributions that have improved the final manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Linking Short- to Medium-Term Beach Dune Dynamics to Local Features under Wave and Wind Actions: A Northern Portuguese Case Study**

**Ana Bio 1,\* , José Alberto Gonçalves 1,2 , Isabel Iglesias <sup>1</sup> , Helena Granja <sup>1</sup> , José Pinho <sup>3</sup> and Luísa Bastos 1,2**


**Abstract:** Many coasts suffer from prevailing erosion, with them being particularly vulnerable to predicted climate change impacts, threatening coastal ecosystems, their services, infrastructures and populations. Understanding coastal morpho-sedimentary dynamics is thus essential for coastal management. However, coastal vulnerability may differ locally, depending on exposure/protection and local geological and morpho-hydrodynamical features, suggesting that a local approach to erosion risk assessment is needed to identify and understand local patterns. Digital elevation models of a 14 km long coastal stretch in northern Portugal that were extracted from aerial surveys obtained between November 2008 and February 2019 were analysed to quantify changes in shoreline position and sediment budgets, both for the whole study area and for distinct beach segments. The observed dynamics were subsequently analysed by considering prevailing wave and wind intensities and directions. Overall and during the decade analysed, the beach–dune system of the studied stretch slightly increased in volume (0.6%), although the shoreline retreated (by 1.6 m on average). Temporal variability in coastal dynamics was observed at all of the temporal scales considered—from seasons to 5-year periods—with them being related to variability in ocean and wind patterns. There was a trend from accretional to erosional conditions, with the first 5-year period showing a mean increase in the beach–dune system's volume of 0.6% and a mean shoreline progradation of 1.5 m, followed by 5-years with 0.0% volume change and 3.1 m shoreline retreat. Locally, the dynamics were very variable, with shoreline dynamics ranging from 24.0 m regression to 51.5 m progradation, and sediment budgets from 213.8 m<sup>3</sup> loss to 417.0 m<sup>3</sup> gain, per segment and for the decade. Stretches with relatively stable morphologies and others with erosional or accretional trends were found, depending on the beach type, shoreline orientation and the presence of defence structures. Rocky beaches were the least dynamic and sandy beaches the most dynamic, with mean shoreline position changes of 0.0 m and <sup>−</sup>3.4 m, respectively, and mean sediment budgets of <sup>−</sup>1.1 m<sup>3</sup> and <sup>−</sup>2.9 m<sup>3</sup> per linear meter of coastline, respectively, for the studied decade. The observed dynamics showed how local conditions interacted with meteo-ocean conditions in shaping local morpho-sedimentary dynamics, stressing the importance of a local approach to coastal erosion monitoring and risk assessment.

**Keywords:** coastal morpho-sedimentary dynamics; erosion/accretion patterns; remote sensing; beach types; meteo-ocean effects

## **1. Introduction**

Coasts are land–ocean interfaces of high environmental and economic value; they provide important ecosystem services, from coastal buffering and inland protection to nutrient

**Citation:** Bio, A.; Gonçalves, J.A.; Iglesias, I.; Granja, H.; Pinho, J.; Bastos, L. Linking Short- to Medium-Term Beach Dune Dynamics to Local Features under Wave and Wind Actions: A Northern Portuguese Case Study. *Appl. Sci.* **2022**, *12*, 4365. https://doi.org/ 10.3390/app12094365

Academic Editors: Miguel Llorente Isidro, Ricardo Castedo and David Moncoulon

Received: 31 March 2022 Accepted: 22 April 2022 Published: 26 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

cycling, biodiversity, and recreational and cultural environments. As a consequence, coasts are often densely populated and modified, suffering increasing anthropogenic pressure [1]. This raises concerns about their increasing vulnerability [2], particularly in the light of potential climate change impacts, such as sea-level rise and changes in ocean wave characteristics [3,4]. Furthermore, many coasts suffer from erosion due to a sediment deficit, mostly caused by human actions that disrupt sediment fluxes and retain sediments in reservoirs and constructions [5]. More than 20% of the European and 30% of the Portuguese coastline are estimated to suffer from coastal erosion [6]. Erosion causes land loss, threatens ecosystems and infrastructures through increased exposure to wave impacts, and increases flood risk [7], and it is expected to be aggravated by climate change [8]. For Portugal, for instance, the projected sea-level rise is 1.14 m (ranging between 0.39 m and 1.89 m, with a 95% probability) by 2100 [9]. Erosion and shoreline retreat are therefore likely to become an even more relevant problem in the coming decades. There is hence an urgent need for coastal zone management and maritime and coastal planning.

The retreat of coastlines can be prevented by hard structures, such as groynes, seawalls and breakwaters, or through soft methods, such as beach nourishment, dune stabilization and bioengineering [10–12]. Given the high implementation and maintenance costs of hard defence structures and their negative impacts on ecosystems, as well as on neighbouring hydro-morphodynamics, soft solutions have gained in popularity [12–14]. This trend is also visible in Portugal, where hard defence structures were implemented from the late 1950s onward, often producing unexpected or unwanted effects and exacerbating downdrift erosion problems [15,16]. Over the past two decades, management strategies have therefore favoured soft defence approaches through artificial beach nourishment, placement of fences, construction of footbridges and revegetation of dunes [17,18].

Furthermore, also after 1950, important dams were built in the main river basins of the Iberian Peninsula, and large harbours were constructed at the river mouths. These two anthropic interventions, next to other significant soil-use changes in the river basins, greatly contributed to diminishing the volumes of sediments that are transported to the coastal platform or feed the alongshore dominant drift. The consequent sediment starvation at the coast implies that the design principles of groynes and breakwaters are not satisfied, making these defence solutions inefficient in most coastal environments.

Sustainable coastal management requires an integrated evaluation of coastal protection through hard defence structures or soft interventions, such as beach nourishment. Longterm cost-benefits and feasibility depending on local conditions, vulnerabilities and values, as well as sediment availability, need to be assessed to decide on mitigation measures or the alternative of a managed retreat, sacrificing areas that are considered less valuable [11].

Coastal morphology reflects the local natural and man-made conditions. Beach morphological changes occur due to the influence of natural phenomena, such as ocean waves and coastal currents, wind and river-flow effects, as well as due to human activities, such as urbanisation, infrastructures and coastal defence intervention [19–21]. Coastal vulnerability may therefore differ locally, depending on exposure/protection and local geological features, as well as local hydrodynamic and meteorological conditions, suggesting that a local approach to monitoring and management is needed [22]. Spatial and temporal fine-scale coastal sedimentary morphodynamics studies showed that the trends observed in largescale studies may be the opposite of those observed in certain beach segments [23]. The consideration of local specificities, such as coastline orientation, exposure and natural or man-made structures, is thus crucial for the assessment of coastal risk and decision-making regarding protective measures for coastal populations and infrastructures.

Analogously, the characteristics of forcing variables, such as the prevailing wave and wind directions and intensities, have to be considered, as these are likely to influence the local impacts [24,25]. The importance of wave direction for the variability of coastal morphodynamics has been widely studied and acknowledged, with wave direction affecting sediment transport, beach rotation and morphology [26]. While, on open-coast, non-embayed beaches, morphology is considered to be primarily controlled by wave direction and energy, sediment transport will vary locally along a non-straight coastline due to uneven wave attenuation [27]. Lemke and Miller [22], who studied the impact of storm erosion on beach morphology, furthermore concluded that due to the propensity for beach conditions to change over short spatial scales, it is important to assess impacts on a local scale.

Another important parameter that affects beach morphology is the wind. Winds affect wave formation, as well as sediment transport on land, with this being particularly important for dune formation, as well as erosion. The aeolian sediment transport, which can be considered as the primary driver of the sediment flux not influenced by oceanic transport, is dependent on the wind velocity [28]. Wind speed needs to exceed a certain threshold value to start the sediment particle movement [29]. However, the wind direction and the angle of wind approach are also important for the transport of sediments, particularly regarding dune formation or erosion.

Other parameters that control sediment transport are the moisture content because moisture increases the inter-particle cohesion and, consequently, reduces the overall rate of transport [30]; vegetation, which helps to fix the sediments; and the availability of sediments for dune formation, which is also dependent on the sediment supply from waves, coastal drifts and rivers [29].

For sound coastal erosion management and mitigation, an analysis of the relationship between local conditions and forcing variables, as well as local morphodynamics, is thus needed to quantify and, ultimately, model these effects, allowing for much-needed prediction of future erosion trends and likely changes, for instance, following anthropogenic interventions (e.g., implementation of coastal defence structures such as breakwaters and groynes, or structures, such as promenades and buildings) or as a result of the predicted climate change impacts.

This can only be achieved through monitoring at temporal and spatial scales that are adequate to capture relevant variabilities and trends. However, high-temporal- and high-spatial-resolution studies of coastal morphology are rare, particularly in countries with a lower budget available for regular monitoring, although technological developments in remote sensing systems and methodologies nowadays allow for regular surveys at adequate spatial resolutions and affordable costs [31]. The present work aimed to show how a combination of digital elevation models, extracted from aerial photography, and in situ data of beach characteristics could provide the necessary information to characterize beach– dune morphodynamics at a local scale. This constituted a novel, local and comprehensive approach to the study of coastal sedimentary morphodynamics.

A case study is presented, where the morpho-sedimentary dynamics of a coastal stretch were analysed both temporally and locally and related to local features, as well as to the regional meteo-ocean conditions. Digital elevation and terrain models of the Vila Nova de Gaia coast in Northern Portugal, extracted from aerial photographic surveys carried out between November 2008 and February 2019, were analysed to quantify changes in shoreline position and sediment budgets for the study area as a whole and distinct beach segments. Changes were analysed at seasonal, yearly, five-yearly and decadal time scales. The observed dynamics were subsequently related to relevant parameters, including beach type/geology, the presence of a groyne and a breakwater, and dominant wave and wind directions and intensities.

## **2. Materials and Methods**

## *2.1. Study Area*

The study area covers a coastal stretch of about 14 km length in northern Portugal, between the Douro river mouth in the north and the city of Espinho in the south, corresponding to the coastline of the Vila Nova de Gaia municipality (Figure 1). This municipality is the third most densely populated in Portugal, with a population density in the coastal zone (up to 5 km from the coastline) of about 2500 inhabitants per square kilometre (data from the National Institute of Statistics, referring to 2011; www.ine.pt, accessed on

15 March 2022). Consequently, its coastline is densely occupied and urbanised. The coast is characterised by rocky beaches in the north (Figure 1a) and mostly sandy beaches with rocky outcrops in the centre and south. Being exposed to the high-energy ocean climate of the North Atlantic and suffering from sediment depletion, the Portuguese coast has several sectors that are prone to erosion [18,32]. Therefore, about 14% of the Portuguese coast is currently defended by artificial structures [33], although these structures, which mitigate sediment loss locally, often increase downdrift erosion [18]. The studied coastal stretch comprises a groyne (Figure 1b), which was built to retain sediments and to fix a submarine outfall, as well as a breakwater (Figure 1c) that was meant to be detached but developed a tombolo during construction, connecting it to the coast [15]. The beaches near Espinho have been managed with repeated artificial beach nourishment [34]. *Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 5 of 27

**Figure 1.** The Vila Nova de Gaia coast (left panel) and its location in Northern Portugal (red rectangle), and closeups of the rocky sector in the North (**a**), the groyne at Praia de Canide (**b**) and the breakwater at Aguda beach (**c**) (image: GoogleEarth, Data SIO, NOAA, U.S. Navy, NGA, GEBCO). **Figure 1.** The Vila Nova de Gaia coast (left panel) and its location in Northern Portugal (red rectangle), and closeups of the rocky sector in the North (**a**), the groyne at Praia de Canide (**b**) and the breakwater at Aguda beach (**c**) (image: GoogleEarth, Data SIO, NOAA, U.S. Navy, NGA, GEBCO).

The regional atmospheric climate varies seasonally. Mediated by the Azores Anticyclone [39], predominant winds from northerly and north-westerly directions occur during the whole year, but with the highest amplitude during the summer months [40]. During autumn and winter, southerly, south-westerly and westerly winds become dominant, though northerly winter winds continue to occur. The north-western Portuguese Atlantic coast is highly energetic, presenting mean significant wave heights of 2–3 m offshore, and mean wave periods of 8–12 s [35,36]. Waves usually come from the NW, which induces a longshore drift current from north to south. This current is in some areas inverted due to the presence of obstacles (such as breakwaters, jetties, groynes, ebb tidal deltas and bars) that promote wave diffraction. The continental shelf is 30–40 km wide [37] and the local tides are dominated by a semidiurnal regime [38], with amplitudes ranging between 2.5 m and 3.8 m.

Series of overlapping photos were taken from a small manned airplane carrying a digital photogrammetric camera (ZI-DMC, 7680 × 13,824 pixels, for the 2008–2010 surveys, and Vexcel UltraCam Falcon, 9420 × 14,430 pixels, for the 2018–2019 surveys) and flying at about 1000 m and 1600 m heights for the 2008–2010 and the 2018 and 2019 surveys, respectively, to provide high-resolution images with 10 cm and 12–13 cm groundsampling distance, respectively. Images were directly georeferenced using an on-board GNSS/INS system. For postprocessing, the GNSS relative positioning mode was used. For each image, the position of the camera projection centre and the attitudinal angles were obtained, and a boresight alignment with ground control points was carried out to correct slight systematic effects in the attitude angles of about 0.02 degrees [31,41,42]. In situ

2.2.1. Surveys and Elevation Models

*2.2. Data* 

The regional atmospheric climate varies seasonally. Mediated by the Azores Anticyclone [39], predominant winds from northerly and north-westerly directions occur during the whole year, but with the highest amplitude during the summer months [40]. During autumn and winter, southerly, south-westerly and westerly winds become dominant, though northerly winter winds continue to occur.

## *2.2. Data*

## 2.2.1. Surveys and Elevation Models

Aerial photos from surveys commissioned in the scope of several research projects were used to obtain orthomosaics and digital elevation models (DEMs) for the study area. Series of overlapping photos were taken from a small manned airplane carrying a digital photogrammetric camera (ZI-DMC, 7680 × 13,824 pixels, for the 2008–2010 surveys, and Vexcel UltraCam Falcon, 9420 × 14,430 pixels, for the 2018–2019 surveys) and flying at about 1000 m and 1600 m heights for the 2008–2010 and the 2018 and 2019 surveys, respectively, to provide high-resolution images with 10 cm and 12–13 cm ground-sampling distance, respectively. Images were directly georeferenced using an on-board GNSS/INS system. For postprocessing, the GNSS relative positioning mode was used. For each image, the position of the camera projection centre and the attitudinal angles were obtained, and a boresight alignment with ground control points was carried out to correct slight systematic effects in the attitude angles of about 0.02 degrees [31,41,42]. In situ measurements in built-up areas were used for calibration, as these contain features that can be used as accurate ground control points.

Surveys were taken on 14 November 2008, 23 April and 11 November 2009, 5 May 2010, 17 May 2018 and 20 February 2019 during spring low tides. DEMs with 1 m resolution were computed from the pairs of stereoscopic aerial images after extracting correlated points through stereo-matching using the Agisoft software [43]. Final DEM accuracies were in the order of 10 cm and, therefore, close to the image resolution [42].

Additionally, a Digital Terrain Model (DTM) from 2014 (kindly provided by the Direção Geral do Território, available at http://mapas.dgterritorio.pt/inspire/atom/downloadservice. xml (accessed on 15 March 2022)) with a 2 m resolution was used to complement the elevation data.

## 2.2.2. Wave Conditions

Wave data were obtained from a DATAWELL directional wave buoy located about 32 km NW of the study area, off the coast of Leixões harbour (41◦19.000 N, 008◦59.000 W), at 83 m water depth (processed data supplied by the Portuguese Hydrographical Institute—IH; https://www.hidrografico.pt (accessed on 15 March 2022)). Average significant heights (Hs) and wave peak directions for 3 h intervals were used (Figure 2). Wave roses were produced to show the wave direction distribution depending on the wave height (Figures 3–5), using 16 cardinal directions and frequencies for 3 significant-wave-height classes: above 4.5 m, between 3 m and 4.5 m, and above 4.5 m. These classes were based on the 85th and 95th percentiles of the mean significant wave heights (Hs) recorded during the study period, which were 3.1 m and 4.4 m, respectively. Three-hourly wave data were available for 89.5% of the total study period, with data missing for January 2010, January 2016, May and June 2018, and February 2019 due to equipment failure (mostly during winter) or maintenance (in spring/summer). Notice, that failures due to the damage or breakdown of the equipment during storm events meant that some data of extreme conditions, which are relevant for sediment transport, were missed.

conditions, which are relevant for sediment transport, were missed.

**Figure 2.** Significant wave heights (Hs) measured at the Leixões wave buoy and estimated wind velocities for the study area during the study period; survey dates are marked with red vertical lines. **Figure 2.** Significant wave heights (Hs) measured at the Leixões wave buoy and estimated wind velocities for the study area during the study period; survey dates are marked with red vertical lines. *Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 7 of 27

measurements in built-up areas were used for calibration, as these contain features that

2022)) with a 2 m resolution was used to complement the elevation data.

Surveys were taken on 14 November 2008, 23 April and 11 November 2009, 5 May 2010, 17 May 2018 and 20 February 2019 during spring low tides. DEMs with 1 m resolution were computed from the pairs of stereoscopic aerial images after extracting correlated points through stereo-matching using the Agisoft software [43]. Final DEM accuracies were in the order of 10 cm and, therefore, close to the image resolution [42]. Additionally, a Digital Terrain Model (DTM) from 2014 (kindly provided by the Direção Geral do Território, available at http://mapas.dgterritorio.pt/inspire/atom/downloadservice.xml (accessed 15 March

Wave data were obtained from a DATAWELL directional wave buoy located about 32 km NW of the study area, off the coast of Leixões harbour (41°19.00' N, 008°59.00' W), at 83 m water depth (processed data supplied by the Portuguese Hydrographical Institute— IH; https://www.hidrografico.pt (accessed 15 March 2022)). Average significant heights (Hs) and wave peak directions for 3 h intervals were used (Figure 2). Wave roses were produced to show the wave direction distribution depending on the wave height (Figures 3–5), using 16 cardinal directions and frequencies for 3 significant-wave-height classes: above 4.5 m, between 3 m and 4.5 m, and above 4.5 m. These classes were based on the 85th and 95th percentiles of the mean significant wave heights (Hs) recorded during the study period, which were 3.1 m and 4.4 m, respectively. Three-hourly wave data were available for 89.5% of the total study period, with data missing for January 2010, January 2016, May and June 2018, and February 2019 due to equipment failure (mostly during winter) or maintenance (in spring/summer). Notice, that failures due to the damage or breakdown of the equipment during storm events meant that some data of extreme

can be used as accurate ground control points.

2.2.2. Wave Conditions

**Figure 3.** Distribution of 3-hourly wave directions for three significant-wave-height classes (**left plots**) and of hourly wind directions for three velocity classes (**right plots**) for the last decade and the two 5-year periods. **Figure 3.** Distribution of 3-hourly wave directions for three significant-wave-height classes (**left plots**) and of hourly wind directions for three velocity classes (**right plots**) for the last decade and the two 5-year periods.

**Figure 4.** Distribution of 3-hourly wave directions for three significant-wave-height classes (**left plots**) and of hourly wind directions for three velocity classes (**right plots**) for the three yearly **Figure 4.** Distribution of 3-hourly wave directions for three significant-wave-height classes (**left plots**) and of hourly wind directions for three velocity classes (**right plots**) for the three yearly periods.

### periods. 2.2.3. Wind Conditions

Wind velocity and direction were obtained from the ERA5-Land reanalysis dataset [44]. ERA5-Land provides hourly high-resolution gridded climate variable estimations from 1950 to present with a resolution of 0.1◦ × 0.1◦ (i.e., about 11.1 km in the N–S direction and 8.3 km in the E–W direction at 41.1◦ N latitude), with the closest available location for the study region being at 41.1◦ N, 8.6◦ W. Hourly wind velocities (in m/s) were extracted and wind directions were computed from the respective estimates of the horizontal speed of air moving towards the east (u-component) and toward the north (v-component) at a height of ten metres above the surface of the Earth (Figure 2). Analogously to the processing of the wave data, wind roses were produced to show the wind direction distribution depending on the wind velocity, using 16 cardinal directions and frequencies for 3 wind-speed classes: more than 6.5 m/s, between 5 m/s and 6.5 m/s, and less than 5 m/s (Figures 3–5). The classes were based on the 85th and 95th percentiles of the hourly wind velocities recorded during the total study period, which were 5.2 m/s and 6.6 m/s, respectively. Notice that these data were model simulations and, although ERA5-Land wind data may be generally acceptable for coastal locations with moderate variability in topography, they should be used with caution for coastal zones because the ocean–land discontinuity affects model stability conditions [45].

**Figure 5.** Distribution of 3-hourly wave directions for three significant-wave-height classes (**left plots**) and hourly wind directions for three velocity classes (**right plots**) for the 3 seasons monitored (from top to bottom: winter, summer, winter). **Figure 5.** Distribution of 3-hourly wave directions for three significant-wave-height classes (**left plots**) and hourly wind directions for three velocity classes (**right plots**) for the 3 seasons monitored (from top to bottom: winter, summer, winter).

### 2.2.3. Wind Conditions *2.3. Analyses*

Wind velocity and direction were obtained from the ERA5-Land reanalysis dataset [44]. ERA5-Land provides hourly high-resolution gridded climate variable estimations from 1950 to present with a resolution of 0.1° × 0.1° (i.e., about 11.1 km in the N–S direction and 8.3 km in the E–W direction at 41.1° N latitude), with the closest available location for the study region being at 41.1° N, 8.6° W. Hourly wind velocities (in m/s) were extracted and wind directions were computed from the respective estimates of the horizontal speed of air moving towards the east (u-component) and toward the north (v-component) at a height of ten metres above the surface of the Earth (Figure 2). Analogously to the processing of the wave data, wind roses were produced to show the wind direction distribution depending on the wind velocity, using 16 cardinal directions and frequencies for 3 wind-speed classes: more than 6.5 m/s, between 5 m/s and 6.5 m/s, and less than 5 m/s (Figures 3–5). The classes were based on the 85th and 95th percentiles of the hourly wind velocities recorded during the total study period, which were 5.2 m/s and 6.6 m/s, The adopted procedure is schematised in Figure 6. The DEMs and DTM, on a regular 1 × 1 m grid and 2 × 2 m grid, respectively, were mapped and analysed in a GIS tool (using ArcGIS 10.6 and its Spatial Analyst and 3D Analyst modules). Differences between elevation models were computed to quantify sediment budgets and changes in shoreline position for different periods. The data series allowed for assessing changes that occurred during the following periods: an approximate decade, between April 2009 and February 2019; two approximate 5-year periods, between April 2009 and 2014 and between 2014 and February 2019 (we assumed that the 2014 surveys took place on the 30th of June, but precise dates are unknown); three periods of a year, between November 2008 and November 2009, between April 2009 and May 2010 and between May 2018 and February 2019 (not a complete year); and three seasons, between November 2008 and April 2009 (winter), between April 2009 and November 2009 (summer) and between November 2009 and May 2010 (winter) (Table 1).

respectively. Notice that these data were model simulations and, although ERA5-Land wind data may be generally acceptable for coastal locations with moderate variability in

**Figure 6.** Scheme of the adopted procedure, where digital elevation models (DEMs) and the digital terrain model (DTM) were analysed in terms of the spatial and temporal morphodynamics, which were subsequently related to local features and wave and wind patterns. **Figure 6.** Scheme of the adopted procedure, where digital elevation models (DEMs) and the digital terrain model (DTM) were analysed in terms of the spatial and temporal morphodynamics, which were subsequently related to local features and wave and wind patterns.

**Table 1.** Changes in volume of the beach–dune system and in the shoreline position (positive and

topography, they should be used with caution for coastal zones because the ocean–land

The adopted procedure is schematised in Figure 6. The DEMs and DTM, on a regular 1 × 1 m grid and 2 × 2 m grid, respectively, were mapped and analysed in a GIS tool (using ArcGIS 10.6 and its Spatial Analyst and 3D Analyst modules). Differences between elevation models were computed to quantify sediment budgets and changes in shoreline position for different periods. The data series allowed for assessing changes that occurred during the following periods: an approximate decade, between April 2009 and February 2019; two approximate 5-year periods, between April 2009 and 2014 and between 2014 and February 2019 (we assumed that the 2014 surveys took place on the 30th of June, but precise dates are unknown); three periods of a year, between November 2008 and November 2009, between April 2009 and May 2010 and between May 2018 and February 2019 (not a complete year); and three seasons, between November 2008 and April 2009 (winter), between April 2009 and November 2009 (summer) and between November 2009

discontinuity affects model stability conditions [45].

*2.3. Analyses* 

and May 2010 (winter) (Table 1).

negative values indicate seaward and landward movements, respectively) that were observed between subsequent surveys (upper rows), and for the yearly, 5-yearly and the decadal periods of the time series analysed. **Dates Period ∆ Volume (m3) ∆ Shoreline (m) Table 1.** Changes in volume of the beach–dune system and in the shoreline position (positive and negative values indicate seaward and landward movements, respectively) that were observed between subsequent surveys (upper rows), and for the yearly, 5-yearly and the decadal periods of the time series analysed.


To assess the local patterns, as well as the effect of beach exposure to wind and wave impacts, the studied coastal stretch was divided into segments. Therefore, the contour line of 1 m above the mean sea level (MSL) of the first survey (November 2008) was drawn and generalised (simplified) using the Douglas–Peucker simplification algorithm [46] with a specified maximum offset tolerance of 30 m. Each segment of this generalized line represented a beach segment, and for each segment, its length and orientation (facing direction) were extracted.

Considering the isoline 1 m above MSL as representative of the shoreline, shoreline dynamics were obtained by calculating the difference between surveys in 2D area per linear meter, using the simplified isoline to determine the coastal stretch or segment length. Given that the inside limit of the study area (and of the segments) is fixed, the change in area per linear meter corresponded to the mean change in shoreline position.

Volumes, differences in volume and differences in shoreline position between surveys were computed for the beach–dune system of the entire study area and per segment. For the segments, volumes were obtained using a buffer for each segment, with straight limits perpendicular to the beach line (Figure 7). To obtain comparable results for the different beach segments, which vary in length, changes in volume were presented per linear meter of segment length. Furthermore, the slope of the more dynamic part of the beach was computed per segment by considering the beach zone with elevations between 1 and 4 m only, approximately corresponding to the upper beach face. Mean slopes per segment were obtained based on the segment area for this elevation range and the segment length. Finally, the results were analysed by considering beach types, beach orientation and local wind and wave climates. *Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 12 of 27

**Figure 7.** Segmentation of the coastal stretch for local analysis (only part of the stretch is shown): November 2008 DEM (**a**); contours (1 m above MSL contour in red) (**b**); and the simplified contour line (in black) providing the segments used for analysis, with their length and facing direction (**c**). **Figure 7.** Segmentation of the coastal stretch for local analysis (only part of the stretch is shown): November 2008 DEM (**a**); contours (1 m above MSL contour in red) (**b**); and the simplified contour line (in black) providing the segments used for analysis, with their length and facing direction (**c**).

#### **3. Results 3. Results**

#### *3.1. Overall Morphodynamics 3.1. Overall Morphodynamics*

retreated 1.6 m on average.

*3.2. Local Morphodynamics* 

The beach and dune system of the studied coastal stretch had an approximate area of 1,400,000 m2 and a volume of 790,000 m3 above the MSL. Overall, the analysed time series started with a year showing marked accretion and shoreline progradation (Table 1), followed by less accentuated erosional and accretional periods. Shoreline dynamics did not always reflect sedimentary budgets, with some periods showing landward migration when the volume increased and seaward movements when the volume decreased. Looking at the morphodynamics for approximate seasons, years, 5-year periods and The beach and dune system of the studied coastal stretch had an approximate area of 1,400,000 m<sup>2</sup> and a volume of 790,000 m<sup>3</sup> above the MSL. Overall, the analysed time series started with a year showing marked accretion and shoreline progradation (Table 1), followed by less accentuated erosional and accretional periods. Shoreline dynamics did not always reflect sedimentary budgets, with some periods showing landward migration when the volume increased and seaward movements when the volume decreased.

decade of the data series (Table 1), different behaviours were observed for comparable periods. Seasonal changes, which could only be studied for the first one and a half years, showed accretion during the first winter and the summer period, followed by a nearly stable second winter. There was also inter-annual variability, with the early years (2008/2009, 2009/2010) showing increases in volume and shoreline progradation, whereas, during the last (approximate) year (2018/2019), volume was lost, though the shoreline position remained on average stable. Considering the changes observed during 5-year periods, the first (2009/2014) showed a slight increase in volume and shoreline progradation, while the second (2014/2019) showed stable volume and shoreline Looking at the morphodynamics for approximate seasons, years, 5-year periods and decade of the data series (Table 1), different behaviours were observed for comparable periods. Seasonal changes, which could only be studied for the first one and a half years, showed accretion during the first winter and the summer period, followed by a nearly stable second winter. There was also inter-annual variability, with the early years (2008/2009, 2009/2010) showing increases in volume and shoreline progradation, whereas, during the last (approximate) year (2018/2019), volume was lost, though the shoreline position

regression. For the approximate decade (2009/2019), there was a slightly accretional trend in terms of volume but an erosional trend in terms of the shoreline position, which remained on average stable. Considering the changes observed during 5-year periods, the first (2009/2014) showed a slight increase in volume and shoreline progradation, while the second (2014/2019) showed stable volume and shoreline regression. For the approximate decade (2009/2019), there was a slightly accretional trend in terms of volume but an erosional trend in terms of the shoreline position, which retreated 1.6 m on average.

## *3.2. Local Morphodynamics*

Of the 52 beach segments obtained through beach line simplification (Figure 8), 3 were not analysed because they represented hard defence structures without a beach (i.e., the southern face of the Praia de Canide groyne (Figure 1b), segment 14, and the two faces of the Aguda breakwater (Figure 1c, segments 41 and 42). In the northern part, segments 1 to 7 were predominantly rocky (Figure 1a), with slopes between 7.6◦ and 20.9◦ , with a mean slope of 10.6◦ . To the south, from segment 8 to segment 52, the beaches were sandy with rocky outcrops and slopes between 2.1◦ and 7.9◦ , with a mean slope of 5.1◦ .

Considering a decade (Table 1), sediment budgets and shoreline dynamics varied spatially (Figure 9). Segments showed sedimentary budgets ranging from losses of 213.8 m<sup>3</sup> to gains of 417.0 m<sup>3</sup> per linear meter of coast and shoreline dynamics ranging from 24.0 m regression to 51.5 m progradation. Stretches with relatively stable morphology and others with erosional or accretional trends were found, depending on beach type, shoreline orientation and the presence of defence structures. Rocky beaches (segments 1–7) were the least dynamic. Sandy beaches were the most dynamic, with mean shoreline position changes of 0.0 m and <sup>−</sup>3.4 m, respectively, and mean sediment budgets of <sup>−</sup>1.1 m<sup>3</sup> and <sup>−</sup>2.9 m<sup>3</sup> , respectively, per linear meter of coastline and for the studied decade. However, two zones of sandy segments could be distinguished as segments with varying erosional trends in the north and centre of the study area (8–39) and segments with accretional trends in the south (40–52). Patterns differed, however, for the two 5-year periods comprising the decade, except for the northern rocky segments, which remained relatively stable throughout. The accretion in the most southern segments (40–52) and the erosion in the northern sandy segments (12–20) took place during the first 5 years. The erosion of the central-southern segments (24–40) occurred during the second 5 years.

The three yearly periods analysed showed a spatial variability similar to the decadal period, with sediment budgets between <sup>−</sup>243.7 m<sup>3</sup> and 377.9 m<sup>3</sup> per linear meter and shoreline changes between −29.3 m and 30.7 m. There was a marked interannual variability (Figure 10). The first year (2008/2009) showed accretion for most segments, except for the rocky northern segments and the two most southern segments, which lost volume and retreated. The behaviour in terms of volume was not always analogous to the shoreline change. The second (2009/2010) and third (2018/2019) years analysed showed patterns that were more similar to the decadal pattern. The northern rocky segments were followed by predominantly eroding segments and a few accretionary southern segments, although there was some variability in the central part during the second year.

The seasonal analysis (Figure 11) showed a first overall accretional winter with a pattern similar to that of the first year, followed by a rather stable summer season, marked by accretion in some southern segments (40–44, 50–52). The second winter, however, displayed much more variability, with some central-zone segments losing a lot of volume. Seasons showed segments with sediment budgets between <sup>−</sup>264.3 m<sup>3</sup> and 301.6 m<sup>3</sup> per linear meter and shoreline changes between −26.1 m and 27.1 m.

Analysis of the segments per beach type and segment orientation (Figures 12–14), confirmed that rocky beaches presented less sedimentary dynamics in general than sandy beaches, as expected, independently of their orientation. Over the 10-year period and on average, rocky segments lost 1.1 m<sup>3</sup> of volume per linear meter (SD = 1.7) and did not change their shoreline position (change = 0.0 m, SD = 1.7). Sandy segments lost 2.9 m<sup>3</sup> of volume per linear meter (SD = 117.6) and retreated 3.4 m (SD = 10.7). Segments neighbouring defence structures, particularly the segment to the north of the Aguda breakwater, showed variable behaviour. For the decadal period (Figure 12), erosion and

shoreline retreat tended to be higher for westward oriented segments compared to northwestward oriented segments, although some segments (particularly the most southern segments of the study area) displayed marked accretion despite facing approximately west. The behaviours of the two 5-year periods were distinct, with a second 5-year period characterized by more intense erosion and shoreline retreat, which seems to become more severe as the segments tended towards a western orientation. Yearly patterns (Figure 13) showed accretion in most segments in the first year, independent of their orientation, and higher erosion in the second year, particularly in segments oriented towards the W. For the third analysed yearly period (corresponding to the last year of the data series), the shoreline dynamics were apparently related to the segment's orientation, similar to the behaviour observed for the last 5-year period. Seasonal patterns (Figure 14) were also distinct. Segment orientation seemed of little importance for the first winter period, which was marked by accretion, and for the following summer period, which showed less intense dynamics. The second winter, on the other hand, showed again a tendency towards higher erosion in the segments oriented towards the west. As seen earlier, volumetric changes did not always correspond to shoreline changes. *Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 14 of 27

**Figure 8.** Orthomosaic of the aerial photographs, with the analysed coastal segments (simplified 1 m isoline in red, segments labelled with their ID numbers) and the delimitation of the beach–dune system (yellow line); segment 14 marks the Canide groyne, segments 41 and 42 mark the Aguda breakwater. **Figure 8.** Orthomosaic of the aerial photographs, with the analysed coastal segments (simplified 1 m isoline in red, segments labelled with their ID numbers) and the delimitation of the beach– dune system (yellow line); segment 14 marks the Canide groyne, segments 41 and 42 mark the Aguda breakwater.

**Figure 9.** Decadal and 5-yearly changes in volume and shoreline position per segment, presented from north (segment 1) to south (segment 52); segment IDs are presented below the graphs; positive shoreline change values represent seaward migration and negative values represent landward migration. **Figure 9.** Decadal and 5-yearly changes in volume and shoreline position per segment, presented from north (segment 1) to south (segment 52); segment IDs are presented below the graphs; positive shoreline change values represent seaward migration and negative values represent landward migration.

**Figure 10.** Yearly changes in volume and shoreline position per segment, presented from north (segment 1) to south (segment 52); segment IDs are presented below the graphs; positive shoreline change values represent seaward migration and negative values represent landward migration. **Figure 10.** Yearly changes in volume and shoreline position per segment, presented from north (segment 1) to south (segment 52); segment IDs are presented below the graphs; positive shoreline change values represent seaward migration and negative values represent landward migration.

**Figure 11.** Seasonal changes in volume and shoreline position per segment, presented from north (segment 1) to south (segment 52); segment IDs are presented below the graphs; positive shoreline change values represent seaward migration and negative values represent landward migration. **Figure 11.** Seasonal changes in volume and shoreline position per segment, presented from north (segment 1) to south (segment 52); segment IDs are presented below the graphs; positive shoreline change values represent seaward migration and negative values represent landward migration.

**Figure 12.** Decadal and 5-yearly changes in volume and shoreline position per segment presented according to segment orientation and marked according to the segment type or its location (N-Aguda: segment 40 north of the Aguda breakwater, S-Aguda: segments 43 and 44 south of the Aguda breakwater, N-groyne: segment 13 north of the groyne, S-groyne: segment 15 south of the groyne). **Figure 12.** Decadal and 5-yearly changes in volume and shoreline position per segment presented according to segment orientation and marked according to the segment type or its location (N-Aguda: segment 40 north of the Aguda breakwater, S-Aguda: segments 43 and 44 south of the Aguda breakwater, N-groyne: segment 13 north of the groyne, S-groyne: segment 15 south of the groyne).

**Figure 13.** Yearly changes in volume and shoreline position per segment presented according to the segment orientation and marked according to the segment type or its location (N-Aguda: segment 40 north of the Aguda breakwater, S-Aguda: segments 43 and 44 south of the Aguda breakwater, N-groyne: segment 13 north of the groyne, S-groyne: segment 15 south of the groyne). **Figure 13.** Yearly changes in volume and shoreline position per segment presented according to the segment orientation and marked according to the segment type or its location (N-Aguda: segment 40 north of the Aguda breakwater, S-Aguda: segments 43 and 44 south of the Aguda breakwater, N-groyne: segment 13 north of the groyne, S-groyne: segment 15 south of the groyne).

35

**Figure 14.** Seasonal changes in volume and shoreline position per segment presented according to segment orientation and marked according to the segment type or its location (N-Aguda: segment 40 north of the Aguda breakwater, S-Aguda: segments 43 and 44 south of the Aguda breakwater, N-groyne: segment 13 north of the groyne, S-groyne: segment 15 south of the groyne). **Figure 14.** Seasonal changes in volume and shoreline position per segment presented according to segment orientation and marked according to the segment type or its location (N-Aguda: segment 40 north of the Aguda breakwater, S-Aguda: segments 43 and 44 south of the Aguda breakwater, N-groyne: segment 13 north of the groyne, S-groyne: segment 15 south of the groyne).

### **4. Discussion 4. Discussion**

### *4.1. Overall Morphodynamics 4.1. Overall Morphodynamics*

For the whole study period (November 2009–February 2019), the beach–dune system studied gained 3.6% in volume and 1.1% in area, with the shoreline moving, on average, For the whole study period (November 2009–February 2019), the beach–dune system studied gained 3.6% in volume and 1.1% in area, with the shoreline moving, on average,

1.1 m seaward (Table 1). However, most of this accretion took place during the first year.

volume and shoreline progradation, whereas the second showed a stable volume but a marked shoreline retreat of 3.1 m, despite protection measures (such as the groyne and 1.1 m seaward (Table 1). However, most of this accretion took place during the first year. After that, the beach volume and shoreline position stabilized. Considering the two analysed 5-year periods (2009/2014 and 2014/2019), the first showed a slight increase in volume and shoreline progradation, whereas the second showed a stable volume but a marked shoreline retreat of 3.1 m, despite protection measures (such as the groyne and breakwater) and artificial beach nourishments. The northwest coast of Portugal has been suffering from a scarce sediment input and a high-energy wave climate, which turns this coast into one of the most active in terms of sediment transport fluxes, with the segment just south of the study area between Espinho and Torreira considered one of the most exposed and vulnerable areas to erosion in the country [18]. Erosion hotspots with 2.5 m shoreline retreat have been found in that region at Vagueira beach [25], despite frequent artificial beach nourishments. A retreat of 1–5 m per year was estimated if no mitigation measures were taken [14]. Therefore, although the studied stretch seems to be rather stable in the medium term (decade), there seems to be a trend towards erosion in recent years (though longer time series are needed to confirm whether this trend is persistent). Furthermore, this trend, if confirmed, in combination with the projected sea-level rise [9], suggests that erosion and shoreline retreat is likely to become an even more serious problem in the future.

There are several possible reasons for these differences between early and later periods of the analysed time series, ranging from differences in forcing variables, such as wave and wind intensities and directions, to differences in sediment availability (dependent on river flows, ocean currents, sediment extraction and beach nourishments). Wave and wind conditions during the period between surveys were therefore analysed to look for patterns that may explain the observed morphodynamics. Offshore wave and onshore wind conditions were quite similar for these two 5-year periods (Figure 3), but the strongest and most damaging winter of the time series occurred in 2013/2014. This winter caused widespread and intense erosion on Atlantic coasts, which took months to some years to recover [47,48]. This may have contributed to the sediment loss between May 2010 and June 2014 and to the recovery between June 2014 and May 2018. Nonetheless, the second 5-year period showed a remarkable and concerning average shoreline retreat.

The difference between early and later periods of the data series, representing more accretional and more erosional behaviours, respectively, was also reflected in the marked inter-annual variation found. During the first year (November 2008 to November 2009), the beach–dune system volume increased by more than 4% and the shoreline moved, on average, more than 3 m seaward. Accretion was less intense for the second yearly period analysed (April 2009 to May 2010), and the (approximate) last year (May 2018 to February 2019) showed erosion instead. Looking at the wave and wind conditions during the yearly periods, patterns diverged, reflecting the general medium-term pattern during the first year, with dominant wave directions from the NW and wind directions from NNW, showing more distributed wave directions and a strong southern wind component during the second year and revealing a more frequent wave direction component from the WNW in the third year. Given that the coastal stretch faced roughly WSW, wave incidence in the last year was less oblique for most segments, possibly reducing sediment transport and deposition via the induced N–S littoral drift.

Seasonal analyses showed two very different winter periods (November 2008 to April 2009 and November 2009 to May 2010, respectively), the first presenting overall accretion and shoreline progradation, while the second presented erosion (though with a progressing shoreline). Looking at the corresponding wave and wind climates, the first winter was characterized by NW–WNW waves and winds distributed from N to W and S, whereas the second winter showed dominant W–NW waves and southern winds, suggesting, again, that more oblique wave incidence promoted accretion and less oblique wave incidence promoted erosion.

However, notice that the temporal changes were analysed based on non-regular surveys. Hence, seasonal and annual analyses do not cover the whole time series, as would be desirable. Moreover, the wave data present gaps caused by equipment failure due to ex-

treme weather and wave conditions, meaning that extreme events may be underrepresented in the data. In terms of wind effects, no clear relationship between wind directions and the morpho-sedimentary dynamics could be identified, possibly because the ERA5-Land wind parameters used are model simulations that may not represent local coastal patterns well [45]. Furthermore, morphodynamics results for the higher beach–dune areas, where wind effects are expected to be most noticeable, are likely to be less precise, particularly in the presence of vegetation. Analyses were based on six DEMs, where dune volume reflects sediments, as well as vegetation, and on one DTM, where vegetation was probably filtered out. In vegetated dunes, a DEM will therefore tend to produce higher volumes than a DTM. Furthermore, dune grass and shrub vegetation will increase in volume during the growing season, affecting seasonal budgets, and in some of the dunes, trees were also found (which may cause interannual, 5-yearly or decadal differences due to growth or felling). Detailed analysis of rocky outcrops showed very similar altitudes in the DEM and the DTM, suggesting that they were comparable. However, the fact that the models may have been processed differently and the lower resolution of the DTM have to be considered as potential error sources.

## *4.2. Local Morphodynamics*

Local morphodynamics often diverged from the general pattern of the coastal stretch as a whole, highlighting the importance of local approaches to erosion monitoring. Morphosedimentary dynamics were found to vary temporally, as well as spatially, with patterns depending on beach type and exposure. This was consistent with findings that sediment transport and budget will vary locally along a non-straight coastline due to uneven wave attenuation [27]. Variability between segments was high and the magnitude of change was of a similar order for the different periods analysed, from seasonal to decadal, showing how rapidly coastal morphology can change and recover.

Over the last decade and from north to south, three zones could be distinguished: the rocky segments in the north, which were morphologically stable, as expected; a central zone of sandy beaches with rocky outcrops, displaying variable patterns but a tendency for erosion; and a southern part covering the beaches close to the city of Espinho that tended towards accretion (Figures 7–12). This pattern also occurred during the second winter (2009/2010), the second year (2010), the second 5-year period (2014–2019) and the last year (2018/2019). The first season and first year showed overall accretion in most segments (12–48), and erosion in some rocky segments (4–7) and the most southern segments (51–52), contradicting the pattern of the later periods.

The observed sedimentary dynamics of the rocky shores were unexpected. Detailed analysis showed that differences were found on some sandy patches in the first segments, as well as between rocks, suggesting that there may be artefacts in the DEM due to shadow effects, which are particularly strong in the crevices between rocks (notice that surveys were always done during spring low tide, which may occur at different times of the day, causing more or less shadow) but also near the shoreline, where wave breaking and runup between the rocks and in crevices may have affected DEM precision.

During the first 5-year period, the central part of the study area already showed erosion, but mostly in the upper central segments (16–23). The segments to the south of the Aguda breakwater (43–52) showed accretion, which could be (at least partly) attributed to a protective effect of the breakwater. In the second 5-year period, nearly all of the central segments showed erosion and shoreline retreat. In fact, comparing the two 5-year periods, many segments showed the opposite behaviour, particularly in terms of shoreline position (e.g., 28–31, 34, 40, 44). This suggested a change in the spatial pattern between 2008 and 2019, with a trend towards erosion in the central zone and a trend towards accretion for the southern zone of the studied coastal stretch. Notice, that there is a large groyne just south of the study area at Espinho, which may have contributed to the observed updrift accretion. However, this groyne has been there for decades with more or less success in retaining updrift sediments, suggesting that a change in hydrodynamics and/or sediment

supply through increased artificial beach nourishments [34] supported the accretion in recent years.

The temporal variability observed was likely a result of the different wave and wind patterns, which affected segments differently depending on the local conditions—i.e., geology, exposure, and updrift and downdrift structures. The wave incidence direction and wave impact will depend on the direction of the waves and the segment's orientation. Murray and Asthon [49] evaluated the alongshore sediment flux for different wave angles and concluded that coastline segments with different orientations experience different alongshore sediment fluxes. Their results show that the wave angle that leads to the highest value of sediment transport is not necessarily the most oblique wave in shallow water. However, wave incidence will also depend on local hydrodynamics. The wave data used provide offshore wave heights and directions, but onshore, breaking wave characteristics are determined by wave propagation phenomena, such as shoaling and refraction, which depend on the coastal morphology and setting.

Nonetheless, analyses of the erosion/accretion indicators (sediment budget and shoreline dynamics) in relation to beach orientation (Figures 12–14), showed: (i) the general temporal trends, with periods of more accretion or erosion, that were also seen in the overall and the spatial analyses; (ii) differences between beach types, with rocky beaches being relatively stable, independent of orientation, and sandy beaches showing marked variability; and, (iii) an apparent tendency towards more erosion/shoreline retreat for segments facing more western directions (despite the scatter). The latter was particularly the case for the periods of the last 5 years and the last year, suggesting that shorelines with less obliquely incident waves presented a higher erosion risk. This could be explained by the dominance of cross-shore transport towards the subaerial beach during wave storms, with the sand deposits occurring at or behind the rocky outcrops. The most westward-facing segments tended to have the highest volume losses and shoreline retreats, and erosion/shoreline retreat increased with more western-dominant wave directions, as was the case in the last year. For the 5-year periods, the wave roses seemed rather similar (as could be expected for such a long period), but the WNW component was higher for the second 5-year period than for the first.

In general, more exposed segments (e.g., segments 17, 19, 31; Figure 5) showed frequent sediment loss and shoreline retreat, yet many apparently protected (embayed) segments retreated too (e.g., 18, 20, 23). This may have been due to wave diffraction. The main wave crest orientation was from the NW, inducing a drift current from north to south, which was in some areas inverted due to the presence of obstacles that promoted wave diffraction, causing downdrift erosion. Shoreline retreat was also found downdrift the Canide groyne, which is a typical effect of these structures that are intended to trap sediments in updrift areas. However, even the segments above the groyne did retreat, though less than those downdrift. Segments north (39, 40), as well as south, of the Aguda breakwater (43, 44) moved from accretional during the first years of the survey series to erosional during the last 5 years. This may have been due to changes in wave action, as explained above, but also due to sand beach management operations that removed sediments from north of the breakwater to nourish the beaches of the city Espinho in the south.

There were two flood events in the Douro river during the study period, one in February 2010 and another in March 2018, with peak river flows of more than 6000 m3/s and 4000 m3/s, respectively. These may have affected longshore sediment transport to the sectors south of the river outlet, as the periods comprising these flood events showed consistent erosion in the intermediate segments (16–25) of the study area. However, a detailed hydrodynamic analysis would be needed to confirm this hypothesis.

An interesting issue was that volumetric changes did not always correspond to shoreline dynamics, suggesting sediment transport from the beach to the dunes. There are different modes of beach dynamics—shoreline advance/retreat (translation modes) and beach steepening/flattening (rotation modes)—that can occur simultaneously and are

often linked to sediment redistribution within the beach–dune system [50]. For instance, segments 39 and 40 (located to the north of the Aguda breakwater) presented shoreline regression and stability (−6.1 m and −0.2 m, respectively) but volume stability and accretion (2.1 m3/m and 265.4 m3/m, respectively) for the decade. Simultaneously, a steepening in the profiles was observed, with the beach slope increasing from 6.5 to 7.7◦ and from 3.7 to 6.6◦ , respectively. Notice that the enormous difference in volume in segment 40 was also due to the beach width. Larger segments will have more surface and hence the capacity to present larger sediment budgets.

## **5. Conclusions**

Concluding, despite the limitations of the data time series, the present study demonstrated how local conditions interact with meteo-ocean conditions in shaping local morphosedimentary dynamics, stressing the need for local approaches to monitoring and erosion risk analyses. Dynamics varied markedly in space, with patterns depending on beach type and exposure, as well as on the presence of coastal defence structures that alter local hydrodynamics and, therefore, sediment transport and deposition.

The analysed coastal stretch suffered an average 1.6 m retreat of its shoreline in a decade with a slight increase in its volume (0.6%). However, analyses of shorter periods (annual and five years) revealed greater shoreline retreat and sediment budget values, demonstrating the high temporal dynamics of this coastal stretch and stressing the need for longer monitoring periods (more than 10 years) to assess the main trends of coastal morphodynamics.

Although the studied stretch seemed to be rather stable in the medium term (decade), the northwest coast of Portugal has been suffering from a scarce sediment input and a high-energy wave climate, which turns this coast into one of the most active in terms of sediment transport fluxes. If a projected mean sea-level rise of the order of 1 m by 2100 is confirmed, shoreline retreat will become a major problem in the next few decades.

Coastal management should therefore be based on structural monitoring programs, with surveys at adequate temporal and spatial scales to understand local dynamics and to be able to apply adequate erosion mitigation measures in the right places. Furthermore, given that many coastal systems show sediment deficits, mainly due to anthropogenic interventions in rivers, such as dams, and sediment extraction, sediment budgets should play a central role in the development of coastal defence strategies.

**Author Contributions:** Conceptualization, H.G., J.P. and A.B.; methodology, A.B., J.A.G. and J.P.; software, J.A.G. and A.B.; validation, J.A.G.; formal analysis, J.A.G., A.B. and I.I.; investigation, H.G., J.P., J.A.G., A.B. and I.I.; resources, H.G., J.P. and J.A.G.; data curation, J.A.G.; writing—original draft preparation, A.B.; writing—review and editing, A.B., J.A.G., I.I., H.G., J.P. and L.B.; visualization, A.B.; supervision, L.B.; project administration, H.G., J.P., L.B. and A.B.; funding acquisition, H.G., J.P., L.B. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by the Strategic Funding UIDB/04423/2020 and UIDP/04423/2020 through national funds provided by the FCT—Foundation for Science and Technology and European Regional Development Fund (ERDF). This work was further funded by the projects ATLANTIDA (ref. NORTE-01-0145-FEDER-000040) and Ocean3R (ref. NORTE-01- 0145-FEDER-000064), both of which were supported by the North Portugal Regional Operational Programme (NORTE2020) under the PORTUGAL 2020 Partnership Agreement and through the European Regional Development Fund (ERDF).

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

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** Data were supplied by the European Union MarRISK project: Adaptación costera ante el Cambio Climático: conocer los riesgos y aumentar la resiliencia (0262\_MarRISK\_1\_E) through EP INTERREG V A España-Portugal (POCTEP) program and the project INNOVMAR-Innovation and Sustainability in the Management and Exploitation of Marine Resources (NORTE-01-0145-FEDER- 000035, within Research Line ECOSERVICES), supported by NORTE 2020, under the PORTUGAL 2020 Partnership Agreement, through the ERDF.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article* **Monitoring of Recovery Process at Yeongildae Beach, South Korea, Using a Video System**

**Jung-Eun Oh <sup>1</sup> , Weon-Mu Jeong <sup>1</sup> , Kyong-Ho Ryu <sup>1</sup> , Jin-Young Park <sup>2</sup> and Yeon-S. Chang 1,\***


**Featured Application: The video monitoring system used in this study can be applied for monitoring coastal processes of erosion and recovery and thus can be useful in designing beach protection/prevention plans for damage by storm waves.**

**Abstract:** Once a beach is eroded by storm waves, it is generally recovered under milder wave conditions. To prevent or reduce damage, it is therefore important to understand the characteristics of the site-specific recovery process. Here, we present the results, based on a data set from a video monitoring system and wave measurements, of the recovery process in a pocketed beach located inside a bay where the shoreline retreated harshly (~12 m, on average, of beach width) during Typhoon TAPAH (T1917) in September 2019. It took about 1.5 years for the beach to be recovered to the level before the typhoon. During this period, the erosion and accretion were repeated, with the pattern highly related to the wave power (*Pw*); most of the erosion occurred when *Pw* became greater than 30 kWatt/m, whereas the accretion prevailed when *P<sup>w</sup>* was no greater than 10 kWatt/m. The recovery pattern showed discrepancies between different parts of the beach. The erosion during storm events was most severe in the southern part, whereas the northern shoreline did not significantly change even during TAPAH (T1917). In contrast, the recovery process occurred almost equally at all locations. This discrepancy in the erosion/accretion process was likely due to human intervention, as a shadow zone was formed in the northern end due to the breakwaters, causing disequilibrium in the sediment transport gradient along the shore. The results in this study could be applied in designing the protection plans from severe wave attacks by effectively estimating the size of coastal structures and by correctly arranging the horizontal placement of such interventions or beach nourishment. Although the application of these results should be confined to this specific site, the method using wave energy parameters as criteria can be considered in other areas with similar environments, for future planning of beach protection.

**Keywords:** beach erosion; storm waves; recovery process; video monitoring; pocket beach

## **1. Introduction**

Beaches are important for nearby communities, as they provide areas for recreational activities. They are sources of tourism, and thus well-groomed beaches with a nice coastal environment attract many visitors, even from abroad. Therefore, maintaining these beaches in good condition is necessary to improve their value. Apart from economic reasons, beaches are also environmentally important because they are buffer zones between the ocean and land, protecting on-land facilities by mitigating wave energy. Because of this buffer role, beaches themselves are exposed to attacks of extreme waves under storms such as tropical cyclones. When the wave energy is released in the surf zone due to the breaking of waves, sediments in the beach face and the nearshore seabed are dynamically

**Citation:** Oh, J.-E.; Jeong, W.-M.; Ryu, K.-H.; Park, J.-Y.; Chang, Y.-S. Monitoring of Recovery Process at Yeongildae Beach, South Korea, Using a Video System. *Appl. Sci.* **2021**, *11*, 10195. https://doi.org/10.3390/ app112110195

Academic Editors: Ricardo Castedo, Miguel Llorente Isidro and David Moncoulon

Received: 30 September 2021 Accepted: 28 October 2021 Published: 30 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

transported. Therefore, one of the serious types of damage as a result of storm attacks is beach erosion.

Under energetic wave conditions, nearshore sediments can move offshore rapidly due to wave-induced underwater currents [1], and the shoreline can significantly retreat, even for a short time of hours or days [2,3]. The volume of erosion is directly related to the maximum wave height during the storm, as the near-bed shear stress initiates sediment motion proportional to the square of the wave height, and more sediment moves under higher shear stress once it exceeds a critical value. In addition, the wave period that contributes to the power of storm waves and the wave direction that influences the distribution of wave energy along the shore are important factors in the amount of eroded sediment. In addition to the wave parameters, the duration of a storm [4] or a sequence of consecutive storms [5,6] could play a significant role in the erosional process.

Once a storm event ends, onshore sediment motion prevails, under mild wave conditions, due to wave nonlinearity such as skewness and asymmetry [7], which results in the natural process of shoreline recovery following storm events. Therefore, studies on the recovery process after storms are as important as studies on the erosional process during storms [8]. The time scales of the post-storm nearshore morphological recovery are varying, from days [5,8] to months [9]. However, large storms can cause rapid erosion locally, from which recovery may take many years or even decades if the impacts are sufficiently large [10]. In extreme cases, the damage can be irreversible if the coastal structures that are built to protect the shore are destroyed by the storms [11].

In the recovery process, the interval between multiple storms is important, as it is inversely related to the storm damage. If the interval between two consecutive storms is less than the recovery time of the beach, the next storm could cause greater damage to the beach [12]. For example, storm clusters with small return periods could impose similar erosion impacts to a single storm that had a much longer return period (i.e., much greater storm). If a beach reached equilibrium after the attack of the first storm, additional erosion would occur only when the intensity of the following storms exceeded that of the first one [5]. Therefore, the recovery after storm events is a complex process, not only depending on the post-storm recovery time but also depending on the intensity of the following storms.

The studies on beach recovery after storms have used various approaches. One of the traditional methods is the direct measurement of the elevation and water depth along the lines set perpendicular to the coastline, which is called profiling. The profiling would be useful for long-term monitoring for beach processes with a low initial cost [13]. The profile data could be also used to validate simulation results calculated from numerical models [14]. The model approach is usefully applied for the prediction of beach recovery for either short-term or long-term processes. Short-term prediction models such as XBeach [15] calculate the rapid changes in the seabed and dunes during storms and the post-storm recovery process in the time scale of months [9]. Due to the high computational cost, however, the application of process models is still restricted for year-scale predictions. For such long-term predictions, line models [16] based on the one-line theory [17] or statistical models such as Monte Carlo have been applied [18].

Another tool that is useful for the long-term observation of the beach process is remote sensing. The satellite imagery can be regularly collected over a long period, and it has been used to derive shorelines and identify their spatial patterns [19,20]. However, it is sensitive to cloud cover and has a limited temporal resolution. Similarly, aerial photos are useful to detect geographical changes over decades, although they cannot provide detailed information on detected variation [21]. On the other hand, airborne and in situ LiDAR also provides detailed information on the beach geography, for comparison and validation of other measurements [22,23]. Despite their advantages, the application of aerial photos and LiDAR data is restricted in measuring the beach recovery process continuously due to the high cost. For this reason, Video Monitoring Systems (VMSs) installed at beaches have been successfully applied to observe the recovery process in specific sites [8,9]. The VMS is

not only useful for long-term monitoring [24] but also for short-term observation of the recovery process after storms [5] and for estimation of the wave parameters [25] as well as wave-induced currents [26].

In this study, we also investigated a recovery process after an attack of storm waves, when a severe erosion occurred, by analyzing the shoreline variations in a pocketed beach located inside a bay in the southeast coast of Korea. Our study focused on the conditions that determined the erosional or depositional process at the beach using data sets from a VMS and wave measurements. As previously described, it is a natural process that beaches undergo erosion by attacks of storm waves, and then they are slowly recovered under milder conditions. However, the criteria that decide the shoreline erosion and accretion are still unclear. Although such information obtained at a beach is site-specific and thus cannot be generally applied for other beaches, the analytical methods may still be applied to beaches with different conditions. In addition, if a criterion is provided for a specific beach, then plans can be designed to protect the beach from severe damage that can lead to a loss of its value.

Another focus was examination of the locality in the recovery pattern between different locations along the beach. Previous studies using the VMS generally focused on the time scale of the recovery process [8]. However, it is also important to analyze the discrepancy in the recovery process between different locations, if it exists, because such disequilibrium may lead to the concentration of damage at specific locations if subsequent storms attack the site before fully recovered. The results of the present study are, therefore, useful not only to understand the criteria for erosion/accretion processes at the study site, but also to analyze the causes of disequilibrium during the recovery process. These outcomes are then applied to suggest measures for conservation of the beach, considering the characteristic erosion and recovery pattern, generally or locally.

The paper is organized as follows. The information of the study site and the damage as a result of Typhoon TAPAH (T1917) are described in Section 2.1. The VMS used in this study and its data are introduced in Section 2.2. The general trend of the beach process is described in Section 3.1, and the locality (discrepancies in the process between areas of the beach) is analyzed in Section 3.2. The discussion on the results is provided in Section 4, and the conclusion of the study in Section 5.

## **2. Materials and Methods**

## *2.1. Study Sites*

Yeongildae Beach is in the southeast of the Korean Peninsula (Figure 1a), inside Yeongil Bay, which faces northeast with a mouth width of ~10 km (Figure 1b). Due to its location in the west corner of the bay, the beach has been protected from severe erosion because the wave energy is generally attenuated when reaching the site. Yeongildae Beach is a ~1.7 km long sandy beach where the sediment size ranges from about 0.15 to 0.35 mm. It is a pocketed beach facing ESE (East-Southeast), bounded by two ports—Pohang Port at the southern end and Duho Port at the northern end. Due to the breakwaters that were built for these ports, the sediments were expected to remain within the coastal cell without loss in the longshore direction. However, the beach has been classified as an erosive beach because input sources of sediment to the beach have been lost due to the construction of the ports. In addition, after the breakwaters of Duho Port were constructed, a shadow zone could have been formed in the lee area of the breakwaters, reducing the wave energy by wave diffraction [27]. The shadow zone has caused redistribution of sediments along the beach shore so that the shoreline width has become narrower near the Yeongil Bridge area but thicker at both ends of the beach near the two ports (Figure 1b). The aerial photograph and satellite image in Figure 2 compare the shorelines measured in 1977 (before the construction of the Duho Port breakwaters) and in 2019 (after construction). It shows that the shoreline severely retreated in the middle of the beach, whereas it accredited at both ends in 2019, compared to 1977, indicating the redistribution of sediment to reach an equilibrium by the changed wave energy field after the construction of Duho Port.

*Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 4 of 19

both ends in 2019, compared to 1977, indicating the redistribution of sediment to reach an equilibrium by the changed wave energy field after the construction of Duho Port.

both ends in 2019, compared to 1977, indicating the redistribution of sediment to reach an equilibrium by the changed wave energy field after the construction of Duho Port.

**Figure 1.** (**a**,**b**) Location of Yeongildae Beach, the Republic of Korea, (**c**) map of Yeongildae Beach with the location of the video monitoring system (S1 and S2). The pictures are captured from Google Earth. **Figure 1.** (**a**,**b**) Location of Yeongildae Beach, the Republic of Korea, (**c**) map of Yeongildae Beach with the location of the video monitoring system (S1 and S2). The pictures are captured from Google Earth. **Figure 1.** (**a**,**b**) Location of Yeongildae Beach, the Republic of Korea, (**c**) map of Yeongildae Beach with the location of the video monitoring system (S1 and S2). The pictures are captured from Google Earth.

**Figure 2.** Comparison of shoreline positions of Yeongildae Beach: (**a**) aerial photograph taken in 1977 and (**b**) satellite image measured in 2019. **Figure 2.** Comparison of shoreline positions of Yeongildae Beach: (**a**) aerial photograph taken in 1977 and (**b**) satellite image measured in 2019. **Figure 2.** Comparison of shoreline positions of Yeongildae Beach: (**a**) aerial photograph taken in 1977 and (**b**) satellite image measured in 2019.

The wave data were measured at two locations, M1 and M2, as shown in Figure 3a. At M1, a pressure transducer was moored at the bottom, at an 8.5 m depth, ~1 km away from the shore of the Yeongildae Beach, measuring the wave height (*Hs*) and period (*Tp*). study.

Because the pressure transducer could not provide information on the wave propagation (*Dp*), it was measured at M2, where an Acoustic Wave and Current Profiler (AWAC) was moored at the bottom, at a 21.6 m depth, located in nearly the center of Yeongil Bay, ~7 km away from Yeongildae Beach. The tidal range in the study area was no greater than 0.3 m, and the tide-induced current was no greater than 0.2 m/s. Therefore, the effect by the tide on the beach processes may not be significant, and thus was not considered in this study. Inside Yeongil Bay, the waves mainly propagate from the NE direction, which is similar to the wave propagation outside the bay. Wave conditions in M2 are moderate because the significant wave heights (௦) are generally less than 2.5 m, except for several extreme storm cases. The dominant wave directions in M2 are NE and NNE, which indicates that the wave propagates inside the Yeongil Bay generally in the normal direction to the mouth of the bay.

Figure 3b shows a wave rose diagram that represents the wave conditions at M2.

The wave data were measured at two locations, M1 and M2, as shown in Figure 3a. At M1, a pressure transducer was moored at the bottom, at an 8.5 m depth, ~1 km away from the shore of the Yeongildae Beach, measuring the wave height (௦) and period (). Because the pressure transducer could not provide information on the wave propagation (), it was measured at M2, where an Acoustic Wave and Current Profiler (AWAC) was moored at the bottom, at a 21.6 m depth, located in nearly the center of Yeongil Bay, ~7 km away from Yeongildae Beach. The tidal range in the study area was no greater than 0.3 m, and the tide-induced current was no greater than 0.2 m/s. Therefore, the effect by the tide on the beach processes may not be significant, and thus was not considered in this

*Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 5 of 19

**Figure 3.** (**a**) Map of the two locations of the wave measurements. Wave height and period were measured at M1 using a pressure transducer moored at water depth 8.5 m in the nearshore area of Yeongildae Beach. was measured at M2 using an AWAC moored at water depth 21.6 m in the middle of Yeongil Bay, (**b**) a rose diagram for wave conditions measured at M2 for 2.5 years (from 2018/5/21 to 2020/11/26). **Figure 3.** (**a**) Map of the two locations of the wave measurements. Wave height and period were measured at M1 using a pressure transducer moored at water depth 8.5 m in the nearshore area of Yeongildae Beach. *Dp* was measured at M2 using an AWAC moored at water depth 21.6 m in the middle of Yeongil Bay, (**b**) a rose diagram for wave conditions measured at M2 for 2.5 years (from 21 May 2018 to 26 November 2020).

> Due to the geographical characteristic of Yeongildae Beach being located at the northwest end, inside Yeongil Bay (Figure 1b), the impact of storm waves has not been significant compared to other beaches located outside the bay. In September 2019, however, Typhoon TAPAH (T1917) attacked the study site. Although TAPAH (T1917) was a Category 1 typhoon (max. wind speed 35 m/s), it caused serious damage to many beaches located on the southeast coast of Korea due to the proximity of its path (Figure 4). The damage was also significant at Yeongildae Beach because the beach face was rapidly eroded dur-Figure 3b shows a wave rose diagram that represents the wave conditions at M2. Inside Yeongil Bay, the waves mainly propagate from the NE direction, which is similar to the wave propagation outside the bay. Wave conditions in M2 are moderate because the significant wave heights (*Hs*) are generally less than 2.5 m, except for several extreme storm cases. The dominant wave directions in M2 are NE and NNE, which indicates that the wave propagates inside the Yeongil Bay generally in the normal direction to the mouth of the bay.

> ing the attack of TAPAH (T1917), which provided the motivation of the present study to investigate the characteristic pattern of its recovery process. Due to the geographical characteristic of Yeongildae Beach being located at the northwest end, inside Yeongil Bay (Figure 1b), the impact of storm waves has not been significant compared to other beaches located outside the bay. In September 2019, however, Typhoon TAPAH (T1917) attacked the study site. Although TAPAH (T1917) was a Category 1 typhoon (max. wind speed 35 m/s), it caused serious damage to many beaches located on the southeast coast of Korea due to the proximity of its path (Figure 4). The damage was also significant at Yeongildae Beach because the beach face was rapidly eroded during the attack of TAPAH (T1917), which provided the motivation of the present study to investigate the characteristic pattern of its recovery process.

## *2.2. Video Monitoring System*

At Yeongildae Beach, two VMSs were constructed at S1 and S2, as shown in Figure 1c. The monitoring by the VMS was initiated in November 2018 and in February 2019 at S1 and S2, respectively. The VMSs were built at the top of buildings in both locations to save the cost of constructing new towers to mount the VMSs. At S1, four video cameras were mounted to cover the central and southern part of the beach, and two additional cameras were mounted at S2 to cover the resting northern end of the beach.

**Figure 4.** Track of Typhoon TAPAH (T1917). **Figure 4.** Track of Typhoon TAPAH (T1917).

*2.2. Video Monitoring System*  At Yeongildae Beach, two VMSs were constructed at S1 and S2, as shown in Figure 1c. The monitoring by the VMS was initiated in November 2018 and in February 2019 at S1 and S2, respectively. The VMSs were built at the top of buildings in both locations to save the cost of constructing new towers to mount the VMSs. At S1, four video cameras were mounted to cover the central and southern part of the beach, and two additional cameras were mounted at S2 to cover the resting northern end of the beach. The VMS consists of the cameras, camera controller, data processing, and data transfer systems that send processed data (averaged images) to the main server. Every 30 min from 7:00 a.m. to sunset, each camera took snapshot pictures for 3 min every half-second. In Figure 5a, b, examples of the snapshot images taken at S1 and S2 are shown. The 180 snapshot images collected for 3 min were then averaged to produce an 'averaged image'; in this way, the averaged images were provided twice per hour. The averaged images stored in the controller were transferred to the main server one time per day in the nighttime, when the internet network was not busy. The snapshot images were not saved The VMS consists of the cameras, camera controller, data processing, and data transfer systems that send processed data (averaged images) to the main server. Every 30 min from 7:00 a.m. to sunset, each camera took snapshot pictures for 3 min every half-second. In Figure 5a,b, examples of the snapshot images taken at S1 and S2 are shown. The 180 snapshot images collected for 3 min were then averaged to produce an 'averaged image'; in this way, the averaged images were provided twice per hour. The averaged images stored in the controller were transferred to the main server one time per day in the nighttime, when the internet network was not busy. The snapshot images were not saved in the system unless specified to do so, to save the cost of running the VMS. Once the averaged-image data were transferred to the main server, the locations of the images were corrected in the x-y coordinate system. Because the images were obliquely measured with angles that were different between the cameras, the data needed to be converted into orthogonal images to the ground. After this, the images from each camera were combined into one-image data that covered the whole beach. Figure 5c shows the data in which the orthogonal images are combined into one.

in the system unless specified to do so, to save the cost of running the VMS. Once the averaged-image data were transferred to the main server, the locations of the images were corrected in the x-y coordinate system. Because the images were obliquely measured with angles that were different between the cameras, the data needed to be converted into orthogonal images to the ground. After this, the images from each camera were combined into one-image data that covered the whole beach. Figure 5c shows the data in which the orthogonal images are combined into one. The shoreline positions and beach width were estimated from the image data of the whole beach. For this process, a set of baselines needed to be established, as shown in Figure 5d. The baselines started from the inner line set, along the landward end of the backshore of the beach. Each baseline was set perpendicular to the inner line to be ended in the water so that the position of the shore could be determined manually by comparing the color of the pixels along the baseline. Once the shoreline positions were determined, The shoreline positions and beach width were estimated from the image data of the whole beach. For this process, a set of baselines needed to be established, as shown in Figure 5d. The baselines started from the inner line set, along the landward end of the backshore of the beach. Each baseline was set perpendicular to the inner line to be ended in the water so that the position of the shore could be determined manually by comparing the color of the pixels along the baseline. Once the shoreline positions were determined, the beach width could be calculated by connecting the shoreline positions of a baseline to that of the adjacent baselines. Similarly, the position where a baseline crossed the inner line could be connected to the adjacent positions to form a polygon, as marked with the blue solid lines in Figure 5. The beach width was then obtained by calculating the area inside the polygon. To measure the time variation of the beach width, the inner line and baselines were fixed once determined. In the case of Yeongildae Beach, a total of 34 baselines were established.

the beach width could be calculated by connecting the shoreline positions of a baseline to that of the adjacent baselines. Similarly, the position where a baseline crossed the inner line could be connected to the adjacent positions to form a polygon, as marked with the blue solid lines in Figure 5. The beach width was then obtained by calculating the area inside the polygon. To measure the time variation of the beach width, the inner line and baselines were fixed once determined. In the case of Yeongildae Beach, a total of 34 base-

**Figure 5.** (**a**, **b**) Snapshots of VMS images captured from one of the cameras on S1 and S2, (**c**) orthogonal image in which the pieces of image data captured by the cameras in the VMS system were combined, (**d**) the baseline system to estimate the shoreline positions and beach width. The inner line is set along the landward end of the backshore of the beach, and the baselines are set to start from the inner line and extend seaward perpendicular to the shoreline. A total of 34 baselines are set in Yeongildae Beach, as marked with red solid lines. The shoreline positions are determined by comparing the color of image pixels along each baseline, and the beach width can be calculated from the area of the polygon (marked with blue lines in the figure), which is composed by connecting the crossing points of the baselines to the inner line and shoreline positions horizontally. **Figure 5.** (**a**,**b**) Snapshots of VMS images captured from one of the cameras on S1 and S2, (**c**) orthogonal image in which the pieces of image data captured by the cameras in the VMS system were combined, (**d**) the baseline system to estimate the shoreline positions and beach width. The inner line is set along the landward end of the backshore of the beach, and the baselines are set to start from the inner line and extend seaward perpendicular to the shoreline. A total of 34 baselines are set in Yeongildae Beach, as marked with red solid lines. The shoreline positions are determined by comparing the color of image pixels along each baseline, and the beach width can be calculated from the area of the polygon (marked with blue lines in the figure), which is composed by connecting the crossing points of the baselines to the inner line and shoreline positions horizontally.

### **3. Results 3. Results**

lines were established.

### *3.1. General Pattern of Shoreline Recovery 3.1. General Pattern of Shoreline Recovery*

In Figure 6, time variations of the four parameters measured in the site are compared for about 33 months, since the initiation of the VMS measurement in November 2018. The significant wave height, , and the significant wave period, ௦, were measured at M1 (Figure 3). The time variations of and ௦ in Figure 6a show complex patterns, as In Figure 6, time variations of the four parameters measured in the site are compared for about 33 months, since the initiation of the VMS measurement in November 2018. The significant wave height, *Hm*0, and the significant wave period, *T<sup>s</sup>* , were measured at M1 (Figure 3). The time variations of *Hm*<sup>0</sup> and *T<sup>s</sup>* in Figure 6a show complex patterns, as they highly fluctuate, which may cause difficulties in understanding the impact of these wave parameters on the shoreline change. Therefore, another parameter, wave power (*P<sup>w</sup>* = *ρg* 2 <sup>64</sup>*<sup>π</sup> <sup>H</sup>*<sup>2</sup> *m*0 *Ts*) was calculated (Figure 6b). *P<sup>w</sup>* is an index that directly measures the wave energy flux and can be used as a parameter that combines the effects of the wave height and period. The time variation of *P<sup>w</sup>* is more clearly distinct from those of *Hm*<sup>0</sup> and *T<sup>s</sup>* . Its magnitude is usually less than 10 kWatt/m, with an average value over the

33 months of ~0.9 kWatt/m. However, there were times when *P<sup>w</sup>* sharply increased, exceeding 20 kWatt/m, with its maximum value reaching ~75 kWatt/m. *Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 9 of 19

**Figure 6.** Time variations of (**a**) the significant wave height, , and the significant wave period, ௦, (**b**) the wave power, ௪, (**c**) the amplitude of the near-bed velocity, , and (**d**) the mean beach width, 〈〉, averaged over the 34 baselines. The red rectangles mark the 4 time periods when ௪ exceeded 20 kWatt/m and 〈〉 was decreased, indicating the high correlation between ௪ and 〈〉. The T1 in (**d**) denotes the time just before the beach was severely eroded by Typhoon TAPAH (T1917) in September 2019, T2 denotes the time when the recovery process started in October 2020, and T3 denotes the time when 〈〉 was recovered to the level of T1 in March 2021. **Figure 6.** Time variations of (**a**) the significant wave height, *Hm*0, and the significant wave period, *Ts*, (**b**) the wave power, *Pw*, (**c**) the amplitude of the near-bed velocity, *Um*, and (**d**) the mean beach width, h*y*i, averaged over the 34 baselines. The red rectangles mark the 4 time periods when *P<sup>w</sup>* exceeded 20 kWatt/m and h*y*i was decreased, indicating the high correlation between *P<sup>w</sup>* and h*y*i. The T1 in (**d**) denotes the time just before the beach was severely eroded by Typhoon TAPAH (T1917) in September 2019, T2 denotes the time when the recovery process started in October 2020, and T3 denotes the time when h*y*i was recovered to the level of T1 in March 2021.

The correlation between the beach width and wave power can be more clearly investigated in Figure 7, in which the time variation of ௪ is also compared with that of 〈〉 for the same observation period. To increase the visibility, however, the range of ௪ in the *y*-axis is adjusted so that its upper limit is set to 30 kWatt/m (Figure 7a). In addition, four green rectangles are added in the figure to mark the times when the magnitude of 〈〉 increased (i.e., when the shoreline was advanced). Red/green rectangles are marked over the whole observational period of 33 months so that the erosion/accretion trend can be examined even during the time before Typhoon TAPAH (T1917). Compared to the pattern in the red rectangles, the times of shoreline accretion are also closely related to the wave power in the green rectangles. The maximum magnitude of ௪ was generally higher than 20 kWatt/m in the periods of five red rectangles, whereas it was less than 10 kWatt/m in the periods of four green rectangles. In particular, at four out of five periods of the red rectangles, the maximum magnitude of ௪ exceeded 30 kWatt/m and could reach up to ~75 kWatt/m. During the longest time of the gradual shoreline accretion for ~4 months, which started from T2, the maximum magnitude of ௪ was generally lower than 10 kWatt/m. This indicates that a continuation of low-powered wave conditions is necessary for a recovery process to occur consistently. In contrast, the erosion process occurred in relatively shorter periods, when high-powered waves attacked the site. Specifi-The mean beach width, h*y*i, was estimated by integrating the beach width, *y*, along each baseline, from baseline #1 to #34. Here, the y magnitude of each baseline was obtained by subtracting the average value during the first 2 years from the original VMS beach width (i.e., y is the average removed beach width). As shown in Figure 6d, the beach was significantly eroded when TAPAH (T1917) attacked the site in September 2019. During the period of ~2 days, h*y*i was reduced by ~12 m. In this study, the time of 23 September 2019 was set as T1, as it was the initial time for the erosion to occur due to the typhoon. Since the severe erosion, the beach width hardly recovered, as h*y*i could not reach the original value in T1 (marked with the dashed magenta line in Figure 6d) until 1 March 2021, and this time was set as T3 as an indication of the time of shoreline recovery, although h*y*i touched the original value only for a short time at T3 and decreased again sharply after that. The time variation of h*y*i from T1 and T3 shows an interesting pattern. Once the severe erosion occurred in T1, h*y*i tended to increase slowly (i.e., with much lower speed than the erosion in T1). However, there were times when h*y*i decreased back to the level when the most severe erosion occurred in T1. This accretion and erosion pattern was repeated until 5 October 2020 (this time was set as T2). After T2, h*y*i continued to increase gradually (without severe erosion) until T3, when h*y*i recovered to the level before the erosion in T1. Once it reached the maximum level in T3, h*y*i decreased again, and the gradually increasing process stopped. After this, h*y*i remained level, without showing a clear pattern of erosion/accretion until July 2021, the time of the most recent data.

cally, it is noted that the gradual increase of 〈〉 stopped in T3, when other high-powered waves with ௪ of ~50 kWatt/m (Figure 6b) attacked the site, which indicates that the highpowered waves not only caused the rapid erosion but also could change the accretion The pattern of h*y*i shows a high correlation with that of *Pw*, especially when the erosions occurred from T1 to T3. To increase visibility, four red rectangles are marked to show the time periods when h*y*i and *P<sup>w</sup>* are significantly correlated. In these four

Following the time variation of the total mean beach width 〈〉 in Figure 7b, c dis-

the figure, the solid blue line shows the time variation of averaged beach width for the baselines #2–#5, which is symbolled as 〈〉௦ as it represents the pattern in the southern end area. Similarly, the orange line, 〈〉, represents the variation in the middle of the

trend and bring it to a halt.

periods, the beach width showed a significant decrease (or the increasing trend in h*y*i was halted, as shown in T3), and the maximum magnitude of *P<sup>w</sup>* was observed to exceed 20 kWatt/m, which was significantly higher than observed at other times. In addition to the wave power, the amplitude of the near-bed velocities, *Um*, was estimated using the formula (*U<sup>m</sup>* = 0.0116*Hm*0/*Ts*) by Soulsby (1986) [28]. *U<sup>m</sup>* could be an important factor for sediment motion because cross-shore sediment fluxes are calculated from bed shear stresses that are based on near-bed velocities [29]. As shown in Figure 6c, the time variation pattern of *U<sup>m</sup>* was similar to that of *Pw*; hence, *U<sup>m</sup>* increasing with increasing *P<sup>w</sup>* might trigger the sediment motions. However, the correlation between h*y*i and *U<sup>m</sup>* is less relevant compared to that between h*y*i and *Pw*; that is, *U<sup>m</sup>* peaked on 3 September 2020, not on 23 September 2019, when h*y*i declined considerably. Furthermore, *U<sup>m</sup>* was an indirect parameter estimated from the empirical formula using the measured wave heights and periods from the wave sensor. Therefore, it may not be appropriate to apply *U<sup>m</sup>* to directly measure the erosion/accretion processes in the beach face. In contrast, *P<sup>w</sup>* was directly estimated from the wave measurements and thus can be suggested as a controlling indicator in this study.

The correlation between the beach width and wave power can be more clearly investigated in Figure 7, in which the time variation of *P<sup>w</sup>* is also compared with that of h*y*i for the same observation period. To increase the visibility, however, the range of *P<sup>w</sup>* in the *y*-axis is adjusted so that its upper limit is set to 30 kWatt/m (Figure 7a). In addition, four green rectangles are added in the figure to mark the times when the magnitude of h*y*i increased (i.e., when the shoreline was advanced). Red/green rectangles are marked over the whole observational period of 33 months so that the erosion/accretion trend can be examined even during the time before Typhoon TAPAH (T1917). Compared to the pattern in the red rectangles, the times of shoreline accretion are also closely related to the wave power in the green rectangles. The maximum magnitude of *P<sup>w</sup>* was generally higher than 20 kWatt/m in the periods of five red rectangles, whereas it was less than 10 kWatt/m in the periods of four green rectangles. In particular, at four out of five periods of the red rectangles, the maximum magnitude of *P<sup>w</sup>* exceeded 30 kWatt/m and could reach up to ~75 kWatt/m. During the longest time of the gradual shoreline accretion for ~4 months, which started from T2, the maximum magnitude of *P<sup>w</sup>* was generally lower than 10 kWatt/m. This indicates that a continuation of low-powered wave conditions is necessary for a recovery process to occur consistently. In contrast, the erosion process occurred in relatively shorter periods, when high-powered waves attacked the site. Specifically, it is noted that the gradual increase of h*y*i stopped in T3, when other high-powered waves with *P<sup>w</sup>* of ~50 kWatt/m (Figure 6b) attacked the site, which indicates that the high-powered waves not only caused the rapid erosion but also could change the accretion trend and bring it to a halt.

Following the time variation of the total mean beach width h*y*i in Figure 7b, c displays the time variations of the beach widths in three different parts of the shoreline. In the figure, the solid blue line shows the time variation of averaged beach width for the baselines #2–#5, which is symbolled as h*y*i*<sup>s</sup>* as it represents the pattern in the southern end area. Similarly, the orange line, h*y*i*<sup>m</sup>* , represents the variation in the middle of the beach as averaged for the baselines #15–#17. The yellow line, h*y*i*<sup>n</sup>* , represents the variation in the northern end as averaged for the baselines #30–#33. The variations of each group of beach widths show different changes from time to time. However, all the time variations of each group of beach widths follow the trend of the total mean beach width h*y*i in T1, T2, and T3. That is, h*y*i*<sup>s</sup>* , h*y*i*<sup>m</sup>* , and h*y*i*<sup>n</sup>* started to decline rapidly at T1 and formed troughs around T2, and then they all gradually increased until around T3.

beach as averaged for the baselines #15–#17. The yellow line, 〈〉, represents the variation in the northern end as averaged for the baselines #30–#33. The variations of each group of beach widths show different changes from time to time. However, all the time variations of each group of beach widths follow the trend of the total mean beach width 〈〉 in T1, T2, and T3. That is, 〈〉௦, 〈〉, and 〈〉 started to decline rapidly at T1 and

formed troughs around T2, and then they all gradually increased until around T3.

**Figure 7.** Comparison of time variations between ௪ and 〈〉 as shown in Figure 6. (**a**) the wave power, ௪, the range of the *y*-axis for ௪ is adjusted from 0 to 30 kWatt/m to increase the visibility of the correlation between the two parameters. (**b**) the mean beach width, 〈〉, averaged over the 34 baselines, and (**c**) the mean beach widths of the southern end, the middle, and the northern end of the beach, 〈௦〉, 〈〉, and 〈〉. T1, T2, and T3 are the same time steps marked in Figure 6. The red rectangles mark the period when the maximum magnitude of ௪ became higher than 20 kWatt and 〈〉 decreased. The green rectangles mark the periods when the maximum magnitude of ௪ was no greater than 10 kWatt/m and 〈〉 generally increased. **Figure 7.** Comparison of time variations between *P<sup>w</sup>* and h*y*i as shown in Figure 6. (**a**) the wave power, *Pw*, the range of the *y*-axis for *P<sup>w</sup>* is adjusted from 0 to 30 kWatt/m to increase the visibility of the correlation between the two parameters. (**b**) the mean beach width, h*y*i, averaged over the 34 baselines, and (**c**) the mean beach widths of the southern end, the middle, and the northern end of the beach, h*ys*i, h*ym*i, and h*yn*i. T1, T2, and T3 are the same time steps marked in Figure 6. The red rectangles mark the period when the maximum magnitude of *P<sup>w</sup>* became higher than 20 kWatt and h*y*i decreased. The green rectangles mark the periods when the maximum magnitude of *P<sup>w</sup>* was no greater than 10 kWatt/m and h*y*i generally increased.

### *3.2. Locality in the Erosion/Accretion Process 3.2. Locality in the Erosion/Accretion Process*

The time variation in Figures 6 and 7 shows the general pattern of the erosion and recovery processes of the averaged beach width of the 34 baselines. As described in Section 2.1 with Figures 1c and 2, however, the beach widths show high locality—i.e., a discrepancy in y between different baselines. Due to the construction of Duho Port, a shadow zone was formed in the area near the baselines #29–#34 (Figure 5), where the shoreline was advanced. In contrast, the shoreline at the area near baselines #21–#27 severely retreated even before Typhoon TAPAH (T1917) attacked the site. Similarly, the beach widths in the southern end (near baselines #1–#10) were thicker compared to those in the middle of the beach (near baselines #11–#20). This high locality in the beach width between different areas of the beach might affect the erosion and recovery processes, which is analyzed in this section. The time variation in Figures 6 and 7 shows the general pattern of the erosion and recovery processes of the averaged beach width of the 34 baselines. As described in Section 2.1 with Figures 1c and 2, however, the beach widths show high locality—i.e., a discrepancy in y between different baselines. Due to the construction of Duho Port, a shadow zone was formed in the area near the baselines #29–#34 (Figure 5), where the shoreline was advanced. In contrast, the shoreline at the area near baselines #21–#27 severely retreated even before Typhoon TAPAH (T1917) attacked the site. Similarly, the beach widths in the southern end (near baselines #1–#10) were thicker compared to those in the middle of the beach (near baselines #11–#20). This high locality in the beach width between different areas of the beach might affect the erosion and recovery processes, which is analyzed in this section.

Figure 8 shows the time variation of parameters that are related to the locality in terms of wave power. In Figure 8a, the colored contours represent the beach width () variations of the 34 baselines, as the numbers in the *y*-axis denote the baseline numbers marked in Figure 5. Here, it is noted that the colors in the contours show the variation in along each baseline but cannot be used to compare the magnitude between the baselines because each *y* is the average removed beach width using the first 2-year average. It Figure 8 shows the time variation of parameters that are related to the locality in terms of wave power. In Figure 8a, the colored contours represent the beach width (*y*) variations of the 34 baselines, as the numbers in the *y*-axis denote the baseline numbers marked in Figure 5. Here, it is noted that the colors in the contours show the variation in *y* along each baseline but cannot be used to compare the magnitude between the baselines because each *y* is the average removed beach width using the first 2-year average. It is clear that the *y* values rapidly dropped after TAPAH (T1917). Compared to the northern area (# of baseline >25), however, the *y* values decreased more severely in the southern end (# of baseline <6). The locality was also observed after the typhoon, as *y* values were even reduced in the southern end near baselines < #5 until the recovery process prevailed, starting in November 2020. In contrast, *y* values slightly increased in the northern end (# of baseline >28) by that time. This pattern continued during the recovery process after November

2020, because the *y* values in the northern end became even greater than those before the typhoon since January 2021, whereas they were smaller than those before the typhoon in the southern end. This pattern of locality can be more clearly quantified through h*y*i*<sup>s</sup>* , h*y*i*<sup>m</sup>* , h*y*i*<sup>n</sup>* , and h*y*i in Figure 8c, in which the time variations of the beach widths calculated for three different groups are compared to the mean beach width. In the southern end, the h*y*i*<sup>s</sup>* rapidly decreased ~30 m after TAPAH (T1917) and reached a minimum on 3 September 2020. The h*y*i*<sup>s</sup>* magnitude could not be recovered even during the recovery process, since November 2020, as the maximum h*y*i*<sup>s</sup>* observed in 2021 was still ~10 m smaller than that before the typhoon. In the northern end, h*y*i*<sup>n</sup>* was not significantly reduced by the typhoon and gradually increased after that time as it became greater than the value before the typhoon during the recovery process, since November 2020. In the middle of the beach, h*y*i*<sup>m</sup>* shows the middle course between the two previous cases. It is also noted that the pattern of h*y*i*<sup>m</sup>* is similar to that of h*y*i, indicating the time variation in the middle of the beach can represent the pattern for the whole beach. *Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 12 of 19

**Figure 8.** Time variation of (**a**) beach widths of all 34 baselines, where the *y*-axis represents the baseline numbers, (**b**) wave power of Figure 7a, (**c**) the beach width, 〈〉௦, averaged for baselines #2–#5, representing the southern end area (blue), 〈〉, averaged for baselines #15–#17, representing the middle part of the beach (orange), 〈〉, averaged for baselines #30– #33, representing the northern end area (yellow), and 〈〉, averaged for the all 34 baselines (black), (**d**) ௧ = ሺ/ሻ , the parameter that indicates the locality of erosion/accretion pattern. The colors of the dots in (**d**) represent the range of ௧ magnitude, with reds for positive and blues for negative values. The red and green rectangles in the figure are marked the same as those in Figure 7. (**c**) The grey and purple circles mark the periods when the time variation pattern of the beach width shows locality between the two ends of the beach. **Figure 8.** Time variation of (**a**) beach widths of all 34 baselines, where the *<sup>y</sup>*-axis represents the baseline numbers, (**b**) wavepower of Figure 7a, (**c**) the beach width, <sup>h</sup>*y*i*<sup>s</sup>* , averaged for baselines #2–#5, representing the southern end area (blue), h*y*i*<sup>m</sup>* , averaged for baselines #15–#17, representing the middle part of the beach (orange), h*y*i*<sup>n</sup>* , averaged for baselines #30–#33, representing the northern end area (yellow), and h*y*i, averaged for the all 34 baselines (black), (**d**) *L<sup>t</sup>* = R (*∂y*/*∂t*)*dx*, the parameter that indicates the locality of erosion/accretion pattern. The colors of the dots in (**d**) represent the range of *L<sup>t</sup>* magnitude, with reds for positive and blues for negative values. The red and green rectangles in the figure are marked the same as those in Figure 7. (**c**) The grey and purple circles mark the periods when the time variation pattern of the beach width shows locality between the two ends of the beach.

For a better understanding of the locality in the shoreline variation pattern, a parameter was developed using the formula ௧ = ሺ/ሻ , as plotted in Figure 8d. The dots in the figure show the magnitude of ௧ at each time step, and their colors also indicate the range of ௧ magnitude, with reds for positive and blues for negative values. In the formula of ௧, x is the alongshore direction toward the northern end. Therefore, the negative value of ௧ indicates that / magnitude decreased in the positive x-direction. In other words, the positive/negative values of ௧ indicate the beach width change rate In Figure 8c, the beach width change pattern is compared at two additional time periods, marked with circles. The grey circle covers the period from October 2019 to January 2020, when the beach settled down after TAPAH (T1917). The purple circle covers the recovery period from October 2020 to February 2021. These two periods were chosen because the beach width change pattern shows a discrepancy between the southern and northern ends. During the period of the grey circle, the magnitude of h*y*i*<sup>s</sup>* decreased,

(/) increased/decreased toward the northern end of the beach. For example, the minimum ௧ value (~−4) during the time of TAPAH (T1917) indicates that the erosion rate

might be most severely eroded in the southern end. Figure 9a shows the distribution of ௧ in terms of wave power. Although its probability of occurrence is low (Figure 9b), ௧ usually became negative with increasing magnitude when ௪ was greater than 10 kWatt/m. In particular, the extreme ௧ values (<−3) occurred when TAPAH (T1917) afflicted the beach. This pattern of ௧ distribution indicates that the shoreline was more severely eroded in the southern part of the beach, compared to that in the northern part, when storm waves attacked the site. In Figure 9c, the distribution of wave direction, , is also plotted in terms of ௧. The *x*-axis denotes the wave propagation angles, which increase clockwise from the origin (0°) at the north. Therefore, 45° and 90° denote the NE and E, respectively. The high-powered waves generally approached the beach from ~NE. In the case of Typhoon TAPAH (T1917), the wave propagation direction ranged between whereas that of h*y*i*<sup>n</sup>* increased, showing that the shoreline retreated in the southern end but advanced in the northern end. These results may be interpreted as an indication of a long sediment movement from the southern area toward the northern area. In contrast, both h*y*i*<sup>s</sup>* and h*y*i*<sup>n</sup>* increased during the recovery process, as marked with the purple circle, which indicates that the shoreline was advanced in both areas during the period.

For a better understanding of the locality in the shoreline variation pattern, a parameter was developed using the formula *L<sup>t</sup>* = R (*∂y*/*∂t*)*dx*, as plotted in Figure 8d. The dots in the figure show the magnitude of *L<sup>t</sup>* at each time step, and their colors also indicate the range of *L<sup>t</sup>* magnitude, with reds for positive and blues for negative values. In the formula of *L<sup>t</sup>* , x is the alongshore direction toward the northern end. Therefore, the negative value of *L<sup>t</sup>* indicates that *∂y*/*∂t* magnitude decreased in the positive x-direction. In other words, the positive/negative values of *L<sup>t</sup>* indicate the beach width change rate (*∂y*/*∂t*) increased/decreased toward the northern end of the beach. For example, the minimum *L<sup>t</sup>* value (~−4) during the time of TAPAH (T1917) indicates that the erosion rate was greater in the southern area and decreased in the increasing x-direction so that beach might be most severely eroded in the southern end. Figure 9a shows the distribution of *L<sup>t</sup>* in terms of wave power. Although its probability of occurrence is low (Figure 9b), *L<sup>t</sup>* usually became negative with increasing magnitude when *P<sup>w</sup>* was greater than 10 kWatt/m. In particular, the extreme *L<sup>t</sup>* values (<−3) occurred when TAPAH (T1917) afflicted the beach. This pattern of *L<sup>t</sup>* distribution indicates that the shoreline was more severely eroded in the southern part of the beach, compared to that in the northern part, when storm waves attacked the site. In Figure 9c, the distribution of wave direction, *Dp*, is also plotted in terms of *L<sup>t</sup>* . The *x*-axis denotes the wave propagation angles, which increase clockwise from the origin (0◦ ) at the north. Therefore, 45◦ and 90◦ denote the NE and E, respectively. The high-powered waves generally approached the beach from ~NE. In the case of Typhoon TAPAH (T1917), the wave propagation direction ranged between 50◦ and 60◦ . Considering the shoreline in Figures 1c and 3a, the storm waves during the typhoon directly attacked the southern part, with a slight slope from the normal direction, whereas a shadow zone was likely formed in the northern part of the beach due to the breakwaters of Duho Port. *Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 13 of 19 50° and 60°. Considering the shoreline in Figures 1c and 3a, the storm waves during the typhoon directly attacked the southern part, with a slight slope from the normal direction, whereas a shadow zone was likely formed in the northern part of the beach due to the breakwaters of Duho Port.

**Figure 9.** Distribution of (**a**) wave power, (**b**) probability of occurrence, and (**c**) wave propagation direction in terms of ௧ = ሺ/ሻ . The red circles in (**a**,**c**) mark ௪ and for Typhoon TAPAH (T1917), respectively. **Figure 9.** Distribution of (**a**) wave power, (**b**) probability of occurrence, and (**c**) wave propagation direction in terms of *L<sup>t</sup>* = R (*∂y*/*∂t*)*dx*. The red circles in (**a**,**c**) mark *P<sup>w</sup>* and *D<sup>p</sup>* for Typhoon TAPAH (T1917), respectively.

The discrepancy in the erosion/accretion pattern between the southern and northern parts of the beach can be also observed in the plane view of shoreline changes. In Figure 10, the beach widths of all 34 baselines are compared at three different times, as they correspond to the times of T1, T2, and T3 in Figures 6 and 7. The figure clearly shows that the The discrepancy in the erosion/accretion pattern between the southern and northern parts of the beach can be also observed in the plane view of shoreline changes. In Figure 10, the beach widths of all 34 baselines are compared at three different times, as they correspond to the times of T1, T2, and T3 in Figures 6 and 7. The figure clearly shows that the reduction

reduction in the beach width from T1 to T2 was greatest in the southern end (# of baseline <6), where the maximum reduction distance along baseline #2 was ~40 m. In contrast, the

#29. At the baselines higher than #30, the change was minimal and not observed. During the recovery process from T2 to T3, the beach width increased at most of the baselines. However, the locality is also clearly observed, as its magnitude at T3 was smaller than that at T1 for the baselines #2–#20, indicating the shoreline was not fully recovered to the level before the typhoon attack in the southern part of the beach. In contrast, the beach width at T3 was greater than that at T1 for the baselines higher than #28. This result indicates that the shoreline in this northern area was not eroded significantly during the attacks of storm waves but advanced under milder wave conditions, leading to an accretion of the

In the present study, the wave data were measured in location M1, which led to the assumption that the wave conditions were uniform along the coast of the beach. Therefore, alongshore variation of sediment transport caused by the spatial difference in the wave conditions could not be detected using the present data sets. As described, the locality in the results was induced by the shadowing by breakwaters and the curved shoreline, which might produce alongshore discrepancies in the shoreline evolution pattern, as implied from the locality parameter, ௧, in Figure 8d. The impacts of alongshore variation of wave conditions can be established, in future studies, by employing an array of wave

shoreline in general.

gauges in the longshore direction.

in the beach width from T1 to T2 was greatest in the southern end (# of baseline <6), where the maximum reduction distance along baseline #2 was ~40 m. In contrast, the beach width in the northern end (# of baseline >28) was not significantly changed from T1 to T2, as the maximum change in the beach width was no greater than 5 m along baseline #29. At the baselines higher than #30, the change was minimal and not observed. During the recovery process from T2 to T3, the beach width increased at most of the baselines. However, the locality is also clearly observed, as its magnitude at T3 was smaller than that at T1 for the baselines #2–#20, indicating the shoreline was not fully recovered to the level before the typhoon attack in the southern part of the beach. In contrast, the beach width at T3 was greater than that at T1 for the baselines higher than #28. This result indicates that the shoreline in this northern area was not eroded significantly during the attacks of storm waves but advanced under milder wave conditions, leading to an accretion of the shoreline in general. *Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 14 of 19

**Figure 10.** Beach width variation along the baseline at three different times. Data was calculated by averaging the beach widths over one week, as their centered times were on 23 September 2019 (navy), 5 October 2020 (red), and 1 March 2021 (green). These three times correspond to T1, T2, and T3, as marked in Figures 6 and 7. **Figure 10.** Beach width variation along the baseline at three different times. Data was calculated by averaging the beach widths over one week, as their centered times were on 23 September 2019 (navy), 5 October 2020 (red), and 1 March 2021 (green). These three times correspond to T1, T2, and T3, as marked in Figures 6 and 7.

**4. Discussion**  From the shoreline erosion and accretion patterns shown in Figures 6 and 7, a consistent description can be provided for the recovery process, after the severe erosion by TAPAH (T1917) occurred in T1. The beach width decreased under high wave conditions and increased under low conditions, which corresponds to the general understanding of the shoreline evolution as previously described in Section 1. However, the criteria that determine the range of accretion and erosion are still unclear. The results in the present study can provide these criteria, although their application should be only confined to this study site. In the present study, the wave data were measured in location M1, which led to the assumption that the wave conditions were uniform along the coast of the beach. Therefore, alongshore variation of sediment transport caused by the spatial difference in the wave conditions could not be detected using the present data sets. As described, the locality in the results was induced by the shadowing by breakwaters and the curved shoreline, which might produce alongshore discrepancies in the shoreline evolution pattern, as implied from the locality parameter, *L<sup>t</sup>* , in Figure 8d. The impacts of alongshore variation of wave conditions can be established, in future studies, by employing an array of wave gauges in the longshore direction.

### Once the severe erosion by TAPAH (T1917) occurred in T1, the recovery process **4. Discussion**

mainly started in T2, about 1 year after the erosion, and continued until T3, about 5 months after T2. The factors that affected short-term (days) processes can be further examined by posing a question about the recovery process—why didn't the recovery process start earlier rather than starting ~1 year, until T2, after the typhoon? The answer may be simple: there were earlier times in which the recovery actually occurred under mild wave conditions. However, it could not continue because another set of high waves attacked the site, restraining the shoreline accretion. Some of the high waves were powerful enough to reverse the accretion trend into erosion. Therefore, the recovery process could not continue further, while the accretion and erosion processes were repeated until T2. During the pe-From the shoreline erosion and accretion patterns shown in Figures 6 and 7, a consistent description can be provided for the recovery process, after the severe erosion by TAPAH (T1917) occurred in T1. The beach width decreased under high wave conditions and increased under low conditions, which corresponds to the general understanding of the shoreline evolution as previously described in Section 1. However, the criteria that determine the range of accretion and erosion are still unclear. The results in the present study can provide these criteria, although their application should be only confined to this study site.

riod of ~5 months from T2 to T3, the wave power remained low so that the recovery process could continue without severe disturbances. The next question on the process is which wave condition determines the trend of erosion or accretion? Based on the results of this study, a rough estimation can be pro-Once the severe erosion by TAPAH (T1917) occurred in T1, the recovery process mainly started in T2, about 1 year after the erosion, and continued until T3, about 5 months after T2. The factors that affected short-term (days) processes can be further examined by posing a question about the recovery process—why didn't the recovery process start

vided: the shoreline was advanced when the wave power was no greater than 10 kWatt/m. In contrast, the erosion likely occurred when the wave power became greater than 20

kWatt/m. During the total observational period of 33 months, the maximum wave power reached up to 75 kWatt/m, whereas the average wave power during the period was ~0.9 kWatt/m. When the wave power was too low in this site, the sediments did not move, and neither erosion nor accretion could occur. However, the criteria for this 'no sediment mo-

It should be noted that there are critical cases in which the criteria may not be applicable. For example, once the beach width reached the minimum value at T1, its magnitude did not fall below this value during the rest of the observation period. This indicates that the erosional status reached a critical limit with this minimum width, and thus further erosion was prevented unless more extreme events occurred to beak the equilibrium

tion' condition could not be determined based on the results of this study.

earlier rather than starting ~1 year, until T2, after the typhoon? The answer may be simple: there were earlier times in which the recovery actually occurred under mild wave conditions. However, it could not continue because another set of high waves attacked the site, restraining the shoreline accretion. Some of the high waves were powerful enough to reverse the accretion trend into erosion. Therefore, the recovery process could not continue further, while the accretion and erosion processes were repeated until T2. During the period of ~5 months from T2 to T3, the wave power remained low so that the recovery process could continue without severe disturbances.

The next question on the process is which wave condition determines the trend of erosion or accretion? Based on the results of this study, a rough estimation can be provided: the shoreline was advanced when the wave power was no greater than 10 kWatt/m. In contrast, the erosion likely occurred when the wave power became greater than 20 kWatt/m. In fact, most rapid erosions occurred when the wave power exceeded 30 kWatt/m. During the total observational period of 33 months, the maximum wave power reached up to 75 kWatt/m, whereas the average wave power during the period was ~0.9 kWatt/m. When the wave power was too low in this site, the sediments did not move, and neither erosion nor accretion could occur. However, the criteria for this 'no sediment motion' condition could not be determined based on the results of this study.

It should be noted that there are critical cases in which the criteria may not be applicable. For example, once the beach width reached the minimum value at T1, its magnitude did not fall below this value during the rest of the observation period. This indicates that the erosional status reached a critical limit with this minimum width, and thus further erosion was prevented unless more extreme events occurred to beak the equilibrium again. On the other hand, there was a period of ~3 months from mid-February to mid-May 2019 (between the first green rectangle and the first red rectangle in Figure 7) when the wave power was generally no greater than 10 kWatt/m, but the beach width did not increase. This is likely because the beach width was thicker than the other period, and thus additional accretion of the shoreline might be disturbed, i.e., the beach width already reached saturation.

The results from locality analysis on the erosion/accretion process showed that the southern part of the beach was more severely eroded, especially when high-powered waves afflicted the site, whereas the shoreline in the northern part of the beach was relatively well protected. During the recovery process after October 2020, however, the imbalance between the two ends was reduced, and the shoreline was advanced almost equally at all parts—southern, middle, and northern—of the beach. The reason for this is still unclear, but the northern end could be protected from erosion due to the shadow zone formed by the breakwaters of Duho Port, in which the wave energy was attenuated, especially when the high-powered storm waves attacked the site.

In the southern end, the breakwater of Pohang Port does not form a shadow zone, considering that most of the waves approached from the NE or NNE, as shown in Figure 3b. Figure 9c also shows that the propagating direction of waves during Typhoon TAPAH (T1917) or other storm waves ranged from 50◦ to 60◦ , which confirms that the damage by storm waves could be focused on the southern part of the beach. In contrast, the recovery process occurred under milder wave conditions, by transporting the sediments equally at all parts of the beach. However, there was a possibility of alongshore sediment movement from the southern part to the middle and northern areas. As shown in Figure 8c, h*y*i*<sup>s</sup>* decreased but h*y*i*<sup>m</sup>* and h*y*i*<sup>n</sup>* increased during the period just after TAPAH (T1917), as marked with the grey circle. This locality of different shoreline evolution processes implies that the sediments could be transported in the longshore direction from the south to the north.

The results of this study can be usefully applied in designing the protection/prevention plans against damage to the beach and nearby coastal areas, not only from the disastrous storm events but also from the processes that may break the equilibrium in the beach status. For example, the suggested criteria that could determine the erosion/accretion conditions

could be useful in planning coastal structures. The size of these structures can be effectively estimated if information on the wave conditions that affect the erosion or accretion in this specific site is known in advance. If the size of such a structure is too large, it may be cost-ineffective. In addition, the beach recovery process could be disturbed if the structure reduces too much wave energy; in this case, the minimum power required for accretion cannot be met. Similarly, the structure would not correctly function if its size is too small to reduce the wave power below the criteria for erosion. The information on the locality in the shoreline evolution process could be also useful in planning hard and soft coastal structures. If the longshore transport from the southern to the northern part is quantitatively understood (although it has not been done in this study), structures such as groins can be designed to reduce the loss of sediment in specific areas. In addition, the results of the locality analysis in this study would be useful in planning beach nourishment—in determining where to put the sand to maximize the effect of the beach fill by considering the imbalance in the beach process. For example, beach nourishment and placing some coastal structures have been planned in the study area, although they have not been finalized. The results of this study are, therefore, expected to provide useful information in determining such plans, to maximize the effect of the structures and nourishment. *Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 16 of 19 reach equilibrium under lower wave energy compared to beaches on the open coast. Figure 11 compares the wave power between Yeongildae Beach (8.5 m depth) and Hwajin Beach (32.0 m depth), located ~20 km north of Yeongildae Beach. Because Hwajin Beach is located on the open coast outside Yeongil Bay, its wave force is much greater than that observed at Yeongildae Beach despite, the difference in water depth. As a result, the sediments in this site could respond more easily at lower wave energy levels, which might lead to a different pattern of erosion and recovery process. For example, the sand size in the study site ranged from 0.15 to 0.35 mm, as it is categorized as fine

It should be noted that the results from this study are confined only to this specific site. This small beach is characterized as a pocketed beach, where the sediments cannot be added from or lost to other coastal cells. In addition, the beach is located inside a bay, so the waves are usually attenuated when reaching the study site such that the shoreline can reach equilibrium under lower wave energy compared to beaches on the open coast. Figure 11 compares the wave power between Yeongildae Beach (8.5 m depth) and Hwajin Beach (32.0 m depth), located ~20 km north of Yeongildae Beach. Because Hwajin Beach is located on the open coast outside Yeongil Bay, its wave force is much greater than that observed at Yeongildae Beach despite, the difference in water depth. to medium sand, which was smaller than the sand size measured in the two beaches located in the open coast outside the bay, where it was ~0.5 mm at both beaches (categorized as medium to coarse sand) [30,31]. The criteria for erosion/accretion suggested in this study can hardly be determined in other beaches on the open coast, where nearshore crescentic sandbars are actively developed. Because the horns of the sandbars are coupled with the shoreline pattern [32], the erosion/accretion of the shoreline has a strong locality rather than showing the general trend as observed in this study site, where no crescentic sandbars developed. For this reason, the information on the locality from these study results should be also confined to this specific site.

**Figure 11.** Comparison of the observed wave power between the two beaches of (**a**) Yeongildae Beach (this study site) and (**b**) Hwajin Beach, located in the open coast outside Yeongil Bay. Although the distance between the two beaches is only ~20 km, the wave power in Hwajin Beach is much greater than that in Yeongildae Beach. **Figure 11.** Comparison of the observed wave power between the two beaches of (**a**) Yeongildae Beach (this study site) and (**b**) Hwajin Beach, located in the open coast outside Yeongil Bay. Although the distance between the two beaches is only ~20 km, the wave power in Hwajin Beach is much greater than that in Yeongildae Beach.

**5. Conclusions**  In this study, we employed a data set of VMS and wave parameters measured over 33 months to investigate the erosional and depositional processes in Yeongildae Beach, located inside Yeongil Bay on the southeastern coast of the Korean Peninsula. The beach As a result, the sediments in this site could respond more easily at lower wave energy levels, which might lead to a different pattern of erosion and recovery process. For example, the sand size in the study site ranged from 0.15 to 0.35 mm, as it is categorized as fine to medium sand, which was smaller than the sand size measured in the two beaches located

was severely eroded when Typhoon TAPAH (T1917) hit the site in late September 2019,

beach width to be recovered to the level before the attack of the typhoon. The study then analyzed the recovery process, during which shoreline erosion and accretion repeated,

The study site is a pocketed beach, as the breakwaters of two ports were built at both ends of the beach, blocking the input or loss of sediments in the longshore direction crossing the lateral boundary of the beach. The beach is also located at the west corner, inside

during which the beach width retreated ~12 m on average; it took about 1.5 years for the 57

corresponding to the wave conditions.

in the open coast outside the bay, where it was ~0.5 mm at both beaches (categorized as medium to coarse sand) [30,31]. The criteria for erosion/accretion suggested in this study can hardly be determined in other beaches on the open coast, where nearshore crescentic sandbars are actively developed. Because the horns of the sandbars are coupled with the shoreline pattern [32], the erosion/accretion of the shoreline has a strong locality rather than showing the general trend as observed in this study site, where no crescentic sandbars developed. For this reason, the information on the locality from these study results should be also confined to this specific site.

## **5. Conclusions**

In this study, we employed a data set of VMS and wave parameters measured over 33 months to investigate the erosional and depositional processes in Yeongildae Beach, located inside Yeongil Bay on the southeastern coast of the Korean Peninsula. The beach was severely eroded when Typhoon TAPAH (T1917) hit the site in late September 2019, during which the beach width retreated ~12 m on average; it took about 1.5 years for the beach width to be recovered to the level before the attack of the typhoon. The study then analyzed the recovery process, during which shoreline erosion and accretion repeated, corresponding to the wave conditions.

The study site is a pocketed beach, as the breakwaters of two ports were built at both ends of the beach, blocking the input or loss of sediments in the longshore direction crossing the lateral boundary of the beach. The beach is also located at the west corner, inside the bay, and the distance from the beach to the mouths of the bay is 8–15 km. Due to its location on the beach, the wave energy was lower and the sediment size was finer, compared to those observed on the open coast outside the bay, because the waves were attenuated when reaching the site.

The results of the analysis corresponded to the general understanding of beach processes—that beach widths decreased/increased under high/low wave energy conditions. In this study site, however, the pattern of shoreline evolution was found to be highly correlated with *Pw*, the wave power that combined the impacts of wave height and period, rather than the wave propagating direction, which was observed to be similar for high-powered waves. In particular, the beach width generally increased when *P<sup>w</sup>* was no greater than 10 kWatt/m and tended to decrease under the condition of *P<sup>w</sup>* greater than 20 kWatt/m. However, most erosional events occurred when *P<sup>w</sup>* exceeded 30 kWatt/m. Especially after TAPAH (T1917), it took a long time (~1.5 years) for the beach width to be fully recovered to the previous level because high-powered waves occasionally disturbed the recovery trend. Therefore, the full recovery of the beach width could be reached over the last five months in the recovery period, during which *P<sup>w</sup>* was kept continuously lower than 10 kWatt/m.

The locality of beach recovery process was also analyzed by examining the discrepancy in the erosion/accretion pattern between different locations within the beach. It was observed that the erosion by TAPAH (T1917) was most severe in the southern part of the beach, whereas the change in the beach width was minimal in the northern part during the same period. This pattern of locality was similarly observed under other storm conditions. Considering the wave propagation direction was similar (~NE) for the storm waves, it is likely that a shadow zone was formed in the northern end due to the breakwaters of Duho Port, such that the sediments in the lee area of the breakwaters were protected because the wave energy was attenuated. On the contrary, the breakwater of Pohang Port located in the southern end did not form a shallow zone, and thus the storm waves directly attacked the southern part, causing erosion. In the phase of recovery under milder wave conditions, the shadow zone did not play an important role, and the shoreline was advanced almost equally at all parts of the beach. However, it was also observed that there was a period just after TAPAH (T1917) when erosion occurred in the southern part while accretion occurred in the northern part, indicating a possible alongshore sediment movement from the south to the north.

The results from this study can be usefully applied to design the protection and prevention plans of the beach from storm damage, considering the general and local erosion and recovery patterns. For example, the criteria that determined the trends of erosion and accretion based on wave power can be used to effectively estimate the size of coastal structures, such as detached breakwaters for shore protection. In addition, the results of locality analysis can be useful to plan the horizontal arrangement of such hard structures or beach nourishment, because they can provide information on where to place the structure or sand to maximize their ability to stabilize the beach. It is noted that the application of these results should be confined only to this specific site and may not be directly applicable to other beaches on the open coasts. For example, the suggested criteria are valid where the wave energy is relatively lower than that on the open coasts. The conditions in this specific beach might result in characteristic shoreline responses. In the case of other beaches located on the open coasts, where various nearshore sandbars develop, the shoreline could be coupled with sandbar positions, and the erosion/accretion pattern would be more complex than that observed in this study site. Regardless of these limitations, the method in this study can be usefully applied in areas with a similar environment for designing prevention plans for disasters due to storms.

**Author Contributions:** Conceptualization, J.-E.O., K.-H.R., W.-M.J. and Y.-S.C.; methodology, J.-E.O., K.-H.R. and J.-Y.P.; formal analysis, J.-E.O.; investigation, J.-E.O., W.-M.J. and Y.-S.C.; resources, J.-E.O. and J.-Y.P.; data curation, J.-E.O.; writing—original draft preparation, J.-E.O. and Y.-S.C.; writing review and editing, J.-E.O. and Y.-S.C.; visualization, J.-E.O. and Y.-S.C.; project administration, K.-H.R.; funding acquisition, W.-M.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Korea Institute of Ocean Science and Technology (KIOST), grant number PE99932, and the Ministry of Oceans and Fisheries, grant number PG52320.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article* **A Multi-Criteria Decision Analysis (MCDA) Approach for Landslide Susceptibility Mapping of a Part of Darjeeling District in North-East Himalaya, India**

**Abhik Saha <sup>1</sup> , Vasanta Govind Kumar Villuri 1,\* , Ashutosh Bhardwaj <sup>2</sup> and Satish Kumar <sup>1</sup>**


**Abstract:** Landslides are the nation's hidden disaster, significantly increasing economic loss and social disruption. Unfortunately, limited information is available about the depth and extent of landslides. Therefore, in order to identify landslide-prone zones in advance, a well-planned landslide susceptibility mapping (LSM) approach is needed. The present study evaluates the efficacy of an MCDA-based model (analytical hierarchy process (AHP)) and determines the most accurate approach for detecting landslide-prone zones in one part of Darjeeling, India. LSM is prepared using remote sensing thematic layers such as slope, rainfall earthquake, lineament density, drainage density, geology, geomorphology, aspect, land use and land cover (LULC), and soil. The result obtained is classified into four classes, i.e., very high (11.68%), high (26.18%), moderate (48.87%), and low (13.27%) landslide susceptibility. It is observed that an entire 37.86% of the area is in a high to very high susceptibility zone. The efficiency of the LSM was validated with the help of the receiver operating characteristics (ROC) curve, which demonstrate an accuracy of 96.8%, and the success rate curve showed an accuracy of 81.3%, both of which are very satisfactory results. Thus, the proposed framework will help natural disaster experts to reduce land vulnerability, as well as aid in future development.

**Keywords:** landslide; multi-criteria decision analysis (MCDA); analytical hierarchy process (AHP); receiver operating characteristics; area under the curve (AUC)

## **1. Introduction**

Disastrous natural hazards such as landslides have brought enormous casualties and economic losses in the past. In hilly areas, the essential requirement for ensuring people's safety is the identification of high-risk landslide-prone zones where development should not be conducted, and where soil stabilization works should be conducted. Nearly 12.60% of landslides have occurred in the Indian Himalayan region, and all areas can be detected as being landslide-prone. In developing countries such as India, various terrain hills are experiencing increasing population density and rapid infrastructural development, which may increase the vulnerability to potential landslides and thus socio-economic losses. The study area is situated in the Himalayan region, where landslides are a common occurrence due to the area's topography and geology. The region encompasses several major cities, including Darjeeling, Ghum, Rangbull, and Sonada, which attract a large number of tourists daily. The population in the area is dense, and significant land-use changes have taken place over the years, further exacerbating the risk of landslides. By conducting landslide susceptibility mapping in these areas, it is possible to identify highrisk zones and inform land-use planning decisions to mitigate the impact of landslides on the local population and infrastructure. Quantitative and qualitative approaches are

**Citation:** Saha, A.; Villuri, V.G.K.; Bhardwaj, A.; Kumar, S. A Multi-Criteria Decision Analysis (MCDA) Approach for Landslide Susceptibility Mapping of a Part of Darjeeling District in North-East Himalaya, India. *Appl. Sci.* **2023**, *13*, 5062. https://doi.org/10.3390/ app13085062

Academic Editors: Miguel Llorente Isidro, Ricardo Castedo and David Moncoulon

Received: 30 March 2023 Revised: 11 April 2023 Accepted: 12 April 2023 Published: 18 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

used to classify landslide susceptibility mapping (LSM) [1–4]. To prevent these losses and improve the efficiency of LSM, new methods and different hybrid models are introduced every day by researchers [5–10]. Machine learning approaches have also been introduced in this field in recent years [11–14]. Different ensemble methods and hybrid models are also used to assess and identify LSM [8,15]. The overall incidences of landslides are increasing, and the factors causing the landslides are the growth of urbanization, and development in landslide-prone areas.

Changing climate patterns, increased local rainfall, etc., are the reasons behind constant changes in landslide-prone zones [16,17]. Most losses would be avoided by identifying the problematic areas, land deforestation, etc. Hence, recognizing the existing and potentially unstable slopes is most significant task. Due to mass changes in urbanization and growing awareness of its socio-economic impact in the hilly region, significant attention is focused on the study of landslides [18]. Hence, LSM of any terrain helps to classify areas into different classes concerning the degree of potential hazards. Recently, various GIS tools have helped to identify conditioning factors that may influence an area's vulnerability to landslides. Rainfall and earthquakes significantly trigger landslides due to their external and temporal aspects [19]. For effective disaster management and future improvement, developers and engineers use landslide susceptibility maps [12,20,21]. The identification of risk areas, coupled with a proper landslide assessment plan, can help disaster management planners and developers to minimize accidents and losses. Many parts of the Himalayas are not physically accessible but geospatially can be accessible everywhere, so assessing landslides using the geospatial technique is easy [22,23]. Geospatial and geotechnical modeling has been used by different researchers in the last few years for landslide susceptibility and hazard assessment using different parameters [24–26]. Land-use change is a major contributor to landslides, causing significant environmental and socioeconomic impacts [27]. In order to achieve sustainability goals and reduce landslide risks, understanding the spatiotemporal evolution and influencing factors of land-use change is essential [28–30]. Land-use change can create slope instability, increasing the susceptibility of an area to landslides. Deforestation, urbanization, and mining are examples of land-use changes that can alter soil properties and destabilize slopes [31]. The impacts of land-use change on landslide susceptibility are complex and influenced by factors such as topography, geology, climate, land-use policies, and population growth [26]. A comprehensive approach is necessary to achieve sustainability goals, which includes considering the impacts of land-use change on landslide susceptibility. Sustainable land-use planning can play a critical role in reducing landslide risks by preserving natural ecosystems, restoring degraded land, and promoting low-impact development practices [32,33]. Accurate landslide susceptibility mapping is crucial for identifying at-risk areas and making informed land-use planning decisions. With advances in remote sensing and geographic information systems (GIS) technology, sophisticated landslide susceptibility models have been developed to capture the complex interactions between land-use change and landslide susceptibility [23,34,35]. Understanding the spatiotemporal evolution and influencing factors of land-use change is the key to achieving sustainability goals and mitigating landslide risks. Sustainable landuse planning and effective landslide susceptibility mapping can contribute to promoting a low-carbon, resilient future. The LSM has multiple uses, identifying unstable areas in advance, disallowing new construction in hazard-prone areas, relief operations, etc. An analytical hierarchy process (AHP)-based method has been attempted in the study area to prepare LSM due to its simplifying method of assigning ranks and weights [21,29,36]. Different natural hazards such as floods, gully erosion, earthquakes, and liquification can also be classified through AHP methodology [37–39]. It is a method to derive ratio scales from paired comparisons from the principal eigenvectors [29].

The study of landslide susceptibility mapping in Darjeeling is due to a combination of various factors in the area. A high-resolution DEM is used for the delineation of various thematic layers. With a high population density, significant land-use changes, and challenging topography, the study region provides an interesting case study for analyzing the

*2.1. Study Area* 

*2.2. Datasets Used* 

complex interactions between topography, land use, and population growth contributing to landslide risk in the Himalayan region. This study has the potential to provide valuable insights for informing land-use planning decisions and reducing the impact of landslides on local communities. The current study focuses on identifying the possible landslide areas for stabilizing the citizens' safety and the area's future development. To aid in the development of more accurate and reliable landslide susceptibility models that can be applied to other regions with similar challenges, the AHP method is adopted to recognize the landslide prone areas. The published result can be utilized by West Bengal tourism for future development and controlled habitation in hazard-prone areas. aid in the development of more accurate and reliable landslide susceptibility models that can be applied to other regions with similar challenges, the AHP method is adopted to recognize the landslide prone areas. The published result can be utilized by West Bengal tourism for future development and controlled habitation in hazard-prone areas. **2. Materials and Methods** 

challenging topography, the study region provides an interesting case study for analyzing the complex interactions between topography, land use, and population growth contributing to landslide risk in the Himalayan region. This study has the potential to provide

landslide areas for stabilizing the citizens' safety and the area's future development. To

### **2. Materials and Methods** The study region is located in the Darjeeling district of West Bengal, India, in the

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 3 of 23

### *2.1. Study Area* Eastern Himalayan mountain environment (Figure 1). The area is situated in latitudes

The study region is located in the Darjeeling district of West Bengal, India, in the Eastern Himalayan mountain environment (Figure 1). The area is situated in latitudes 88◦171<sup>0</sup> E to 88◦343<sup>0</sup> E, 26◦933<sup>0</sup> N to 27◦082<sup>0</sup> N and encompasses around 256 Km<sup>2</sup> of the region. The maximum elevation in this region is 2600 m from mean sea level (MSL). The most dreaded and widespread landslide event in the Darjeeling Himalayas occurred in 1968, triggering many shallow and deep-seated landslides [40]. It is believed that many prominent and large landslides in the Darjeeling hills are still active today [41]. Thus, the region requires a systematic study for LSM. The landslide's study area induces inherent factors such as highly dissected slopes, barren land, etc. Tea plantation is the primary land use operation in this area. The study region had a population density of 586 per sq·Km, while the Darjeeling district's population was 1,846,823 as per the 2011 census. 88°171′ E to 88°343′ E, 26°933′ N to 27°082′ N and encompasses around 256 Km2 of the region. The maximum elevation in this region is 2600 m from mean sea level (MSL). The most dreaded and widespread landslide event in the Darjeeling Himalayas occurred in 1968, triggering many shallow and deep-seated landslides [40]. It is believed that many prominent and large landslides in the Darjeeling hills are still active today [41]. Thus, the region requires a systematic study for LSM. The landslide's study area induces inherent factors such as highly dissected slopes, barren land, etc. Tea plantation is the primary land use operation in this area. The study region had a population density of 586 per sq·Km, while the Darjeeling district's population was 1,846,823 as per the 2011 census.

**Figure 1.** The study area shows a part of the Darjeeling district. **Figure 1.** The study area shows a part of the Darjeeling district.

pared using various remotely sensed data from multiple sources. The GIS software (ArcGIS 10.4) was used for preparing layers such as slope, aspect, drainage density, and

## *2.2. Datasets Used*

Thematic layers significant to the causative factors presented in Table 1 were prepared using various remotely sensed data from multiple sources. The GIS software (ArcGIS 10.4) was used for preparing layers such as slope, aspect, drainage density, and lineament density from ALOS PALSAR DEM. The soil and geomorphology maps of West Bengal were collected from the National Bureau of Soil Survey and Land Use Planning (NBSS and LUP), Government of India. A geology map was compiled from the Geological Survey of India (GSI). The Catrosat-2 series (MX) dataset generated a land use and land cover (LULC) map. The drainage density and lineament density maps were made from the ALOS PALSAR DEM data. Data on rainfall were gathered from the Current Research Unit (CRU). The inverse distance weighting (IDW) method for interpolation was used to prepare a rainfall map. An earthquake map was prepared from the data procured from the National Centre for Seismology, New Delhi, India.


**Table 1.** Data used for generating landslide conditioning factor.

## *2.3. Methodology*

The basic methodology adopted for preparing LSM is shown in Figure 2. The initial stage of a susceptibility assessment involves gathering all available information and data about the study area. This stage is crucial as it lays the foundation for further analysis of the relationship between landslide occurrence and conditioning parameters. However, assessing the cause–effect relationships can be challenging since landslides are typically caused by multiple factors. To create a landslide susceptibility map for the Darjeeling region of West Bengal, ten inputs were chosen based on the most significant factors influencing landslides, including slope, aspect, earthquake activity, drainage density, rainfall, land use/land cover (LULC), lineament density, geology, geomorphology, and soil type [12,29]. Our study expanded upon the number of factor classes considered compared to prior research. We reviewed the literature by other authors and sought to assign non-homogeneous weights to each factor class, in contrast to the uniform weighting used in previous studies. The various layer classes are rated in the range between 0 and 9, where higher ratings represent more decisive landslide influence. After preparing the final LSM, the study area was differentiated into different susceptibility zones, such as low, high, moderate, and very high.

**Figure 2.** Graphical representation of overall methodology of landslide susceptibility mapping. **Figure 2.** Graphical representation of overall methodology of landslide susceptibility mapping.

**Formatted:** Not Highlight

## *2.4. Landslide Inventory*

*2.4. Landslide Inventory*  The position and extent of the landslide must be appropriately identified when creating the maps of landslide susceptibility. Landslide inventory is a critical component of all types of landslide zoning classification, including hazard, susceptibility, and risk zoning [42]. It has to deal with the time, location, sort, amount, and travel distance of land sliding in a specific area. Different methods are available for identifying landslides. Field observations, satellite photographs, book studies for details on previous landslides, and aerial photography are a few of them [43]. The landslide inventory map was produced using a combination of visual analysis of aerial imagery, field data, and satellite photography [14]. The whole data on landslides in India were made available by the Geological Survey of India (GSI). Resampled public data from the GSI were used to create the inventory map for the study area. One hundred and fourteen landslides have been found by this study using the available data (Figure 3). Most of the landslides recorded in this area have cliffs, highways, and banks of rivers. Erosion, steep terrain, and a lack of flora are The position and extent of the landslide must be appropriately identified when creating the maps of landslide susceptibility. Landslide inventory is a critical component of all types of landslide zoning classification, including hazard, susceptibility, and risk zoning [42]. It has to deal with the time, location, sort, amount, and travel distance of land sliding in a specific area. Different methods are available for identifying landslides. Field observations, satellite photographs, book studies for details on previous landslides, and aerial photography are a few of them [43]. The landslide inventory map was produced using a combination of visual analysis of aerial imagery, field data, and satellite photography [14]. The whole data on landslides in India were made available by the Geological Survey of India (GSI). Resampled public data from the GSI were used to create the inventory map for the study area. One hundred and fourteen landslides have been found by this study using the available data (Figure 3). Most of the landslides recorded in this area have cliffs, highways, and banks of rivers. Erosion, steep terrain, and a lack of flora are often the causes of landslides along cliffsides.

often the causes of landslides along cliffsides.

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 6 of 23

**Figure 3.** Map of landslide inventory of study area. **Figure 3.** Map of landslide inventory of study area.

*2.5. Landslide Conditioning Factor (LCF) 2.5. Landslide Conditioning Factor (LCF)*

### 2.5.1. Slope 2.5.1. Slope

The slope is one of the triggering factors for landslides in any hilly region. Steeper slopes are more hazard-prone to instability compared with lesser slopes. Shear stress in soil and other uncemented materials gradually rises as the slope angle increases [44,45]. In a steep slope, the gravitational force works more than in a moderate slope, though the shear stress is the same [31]. Depending on the slope angle, it is classified into five different classes with the help of a natural break classification technique: <15°, 14.6°–22.45°, 22.46°– 29.74°, 29.75°–37.88°, >37.89°. (Figure 4a). Landslides usually happen in areas with high slopes. Hence, areas with high and extremely high slopes have a higher likelihood of experiencing landslides. On the other hand, landslide incidence is more or less stable in areas with a moderate to low slope angle [46]. The slope is one of the triggering factors for landslides in any hilly region. Steeper slopes are more hazard-prone to instability compared with lesser slopes. Shear stress in soil and other uncemented materials gradually rises as the slope angle increases [44,45]. In a steep slope, the gravitational force works more than in a moderate slope, though the shear stress is the same [31]. Depending on the slope angle, it is classified into five different classes with the help of a natural break classification technique: <15◦ , 14.6–22.45◦ , 22.46–29.74◦ , 29.75–37.88◦ , >37.89◦ . (Figure 4a). Landslides usually happen in areas with high slopes. Hence, areas with high and extremely high slopes have a higher likelihood of experiencing landslides. On the other hand, landslide incidence is more or less stable in areas with a moderate to low slope angle [46].

**Formatted:** Highlight

**Figure 4.** *Cont*.

**Figure 4.** *Cont*.

**Figure 4.** Landslide conditioning factor: (**a**) slope, (**b**) rainfall, (**c**) earthquake, (**d**) lineament density, (**e**) drainage density, (**f**) geology, (**g**) geomorphology, (**h**) aspect, (**i**) LULC, (**j**) soil. **Figure 4.** Landslide conditioning factor: (**a**) slope, (**b**) rainfall, (**c**) earthquake, (**d**) lineament density, (**e**) drainage density, (**f**) geology, (**g**) geomorphology, (**h**) aspect, (**i**) LULC, (**j**) soil.

### 2.5.2. Rainfall 2.5.2. Rainfall

As rainfall is one of the triggering factors for landslides, May to July is the most vulnerable period [47,48]. The inverse distance weighting (IDW) method for interpolation is used to prepare a rainfall map (Figure 4b). For the rainfall map, the last ten years of data from the Darjeeling and Kalimpong areas were collected from the Climate Research Unit (CRU). The range of rainfall data are 2005–2143 mm, and are further divided into five classes: 2005–2038 mm, 2039–2062 mm, 2063–2088 mm, 2089–2114 mm, and 2115–2143 As rainfall is one of the triggering factors for landslides, May to July is the most vulnerable period [47,48]. The inverse distance weighting (IDW) method for interpolation is used to prepare a rainfall map (Figure 4b). For the rainfall map, the last ten years of data from the Darjeeling and Kalimpong areas were collected from the Climate Research Unit (CRU). The range of rainfall data are 2005–2143 mm, and are further divided into five classes: 2005–2038 mm, 2039–2062 mm, 2063–2088 mm, 2089–2114 mm, and 2115–2143 mm.

### mm. 2.5.3. Earthquake

2.5.3. Earthquake A landslide due to an earthquake can be triggered either by increased shear stress due to acceleration horizontally, or decreased strength of materials [49]. When any hilly area is seismically activated, there is potential for an earthquake to occur [49]. The Eastern Himalayas region is very active seismically. A map of earthquakes has been created using point data from the National Centre for Seismology in New Delhi, which spans more than 200 years. The earthquake map (Figure 4c) has been created using the Inverse Distance A landslide due to an earthquake can be triggered either by increased shear stress due to acceleration horizontally, or decreased strength of materials [49]. When any hilly area is seismically activated, there is potential for an earthquake to occur [49]. The Eastern Himalayas region is very active seismically. A map of earthquakes has been created using point data from the National Centre for Seismology in New Delhi, which spans more than 200 years. The earthquake map (Figure 4c) has been created using the Inverse Distance Weighted tool in a GIS environment, and is divided into three classes.

#### Weighted tool in a GIS environment, and is divided into three classes. 2.5.4. Lineament Density

2.5.4. Lineament Density A study of landslide density helps us to understand the causative elements of landslides [50]. From the literature, it is understood that in the Himalayan region, landslide circumstances are generally very close to the local geological lineament [50–53]. Landslides are further vulnerable in the joint, fractured, and faulted zones [28]. The thematic map lineament density (Figure 4d) was generated from ALOS PALSAR DEM images after A study of landslide density helps us to understand the causative elements of landslides [50]. From the literature, it is understood that in the Himalayan region, landslide circumstances are generally very close to the local geological lineament [50–53]. Landslides are further vulnerable in the joint, fractured, and faulted zones [28]. The thematic map lineament density (Figure 4d) was generated from ALOS PALSAR DEM images after interpretation in a GIS environment. Rectilinear inclinations of morphological features, linear stream courses, structural alignments, and tonal contrast aid in interpreting the lineaments.

interpretation in a GIS environment. Rectilinear inclinations of morphological features, linear stream courses, structural alignments, and tonal contrast aid in interpreting the Despite discovering huge lineaments, no substantial thrusts or faults have been observed in the study area [12].

## 2.5.5. Drainage Density

Drainage density (DD) is a triggering factor for landslides, where drainage controls the landslide, and its densities correspond to the nature of geotechnical features as well as soil properties [40]. River Tista is the biggest river in the study area. Almost all streamlines initiated from high ridges flow down by making valleys and meet the Tista river. Infiltration is inversely related to drainage density [54]. Equation (1) is used to generate drainage density:

$$\text{DD} = \frac{\text{L}\_{\text{d}}}{\text{S}\_{\text{d}}} \tag{1}$$

where L<sup>d</sup> is the measurement term of the drainage system and S<sup>d</sup> represents drainage basin size. The Euclidean distance approach generates drainage density maps in GIS and is classified into three classes (Figure 4e).

## 2.5.6. Geology

The structure and content of different geology determine the strength of the rock [14]. Compared to soft rocks, the more substantial rocks provide excellent protection from the main thrusts and are less prone to landslides [55,56]. Geological rock properties represent texture, color, grain size, or properties. The geology map (Figure 4f) was prepared using the geology of the study area published by GSI.

## 2.5.7. Geomorphology

Geomorphology plays a significant role in the command of landslides. Inflated altitude regions are more vulnerable to landslides than lower altitude regions. The geomorphology map (Figure 4g) was extracted from the published map from the BHUKOSH portal (https: //www.bhukosh.gsi.gov.in (accessed on 31 May 2022)) in the GIS environment.

## 2.5.8. Aspect

The slope aspect represents the orientation of the slope angle [57]. The south-east direction is more vulnerable than the rest of the directions [56]. East- and south-facing slopes are more likely to have landslides, according to the distribution of landslides [58]. Generally, north-facing slopes have more vegetation density than south-facing slopes [40]. Due to sun rotation, in the afternoon, west-facing slopes experience the hottest time of the day [59]. The slope aspect map (Figure 4h) is classified into 1. flat, 2. north, 3. north-east, 4. east, 5. south-east, 6. south, 7. south-west, 8. west, 9. north-west, and this map was generated in the GIS platform. The slope facing south and east makes up a significant component of the study area. The west-facing slope covers the lower part of the area, while a slope angle pointing north comprises a moderate to low slope zone.

## 2.5.9. Land Use and Land Cover

A key element contributing to the incidence of landslides is land use and land cover. Barren regions are more prone to landslides than lush trees. Deforestation is another essential factor that causes landslides. Areas covered with vegetation have seen that the big woody trees with long root systems have helped to improve slope stability in the area. Vegetation density is inversely proportional to landslides [55]. The LULC map (Figure 4i) was extracted using CARTOSAT-2 (MX) series imagery. The LULC map is organized into many categories, including rural and urban areas, sparse forests, tea plantations, agricultural land, barren land, and waterbodies. The southern portion of the area is where you will find agricultural land and developed regions with excellent road access. The bulk of the basin is covered in forest. In the southern zone, human intervention is more common than in the northern region, owing to accessibility issues.

## 2.5.10. Soil

The occurrence of landslides impacts topsoil cover on a slope. Regarding soil texture, rocky and sandy loam have higher landslide potential than silt loam, fine sandy loam, gravelly silt loam, and loam [22,60]. The soil map (Figure 4j) for the study area was created from a local soil map made by NBSS and LUP data on a scale of 1:500,000. Three different textural groups: fine loamy, coarse loamy, and loamy soil skeletal, are found in this region.

## *2.6. Method*

## 2.6.1. Analytical Hierarchy Process (AHP)

The multi-criteria decision analysis (MCDA) method, AHP, is a systematic way to arrange and evaluate difficult mathematical choices. Since then, it has undergone substantial research and development. The factors affecting landslide susceptibility must be added while creating the susceptibility map. The layers have been weighted for this based on their significance. In AHP, the pair-wise comparison matrix method is implemented for the MCDM structure. Qualitative (subjective) and quantitative (objective) decision-making analysis is performed using the AHP method. The number of columns and rows in the comparison matrix is equal, where value 1 is placed on the diagonal of the matrix, and one side of the diagonal stores the scores. Each layer must rate against the other to construct a pair-wise confusion matrix. The rating value ranges from 1 to 9 (Table 2). Each pair-wise comparison matrix value represents the importance of two factors. It was thought that if attribute A is given a score of 9, meaning it is more important than attribute B, then B must receive a score of 1/9, meaning it is less important than A.

The consistency ratio (CR) and consistency index (CI) need to be calculated to validate the pair-wise comparison matrix [36,61,62]. We need to know the value of the 1970 Saatyprojected random consistency index (RI) (Table 3) to compute the CR [59]. If the value of CR < 0.10, we will accept the pair-wise comparison matrix. The following formula used for calculating CI and CR [63]:

$$\text{Consistency Index} (\text{CI}) = \frac{\lambda\_{\text{max}} - \mathbf{n}}{\mathbf{n} - 1} \tag{2}$$

**Table 2.** The scale of weightage for AHP given by Saaty [63].


**Table 3.** RI value for calculating CR by Saaty [63].


The principle eigenvalue is λmax and n is the number of factors.

$$\text{Consistency Ratio} (\text{CR}) = \frac{\text{Consistency index} (\text{CI})}{\text{Random Consumption} \,\text{Index} (\text{RI})} \tag{3}$$

2.6.2. Multi-Collinearity Analysis of LFC

In the least squares regression analysis, the variance inflation factor (VIF) is used to assess the level of multi-collinearity. The exponent shows the multicollinearity-based increases in the coefficient [64]. The amount of multi-collinearity may be evaluated using the VIF's value. An observational rule states that multi-collinearity is high if the VIF value exceeds 5. The second approach to investigating multi-collinearity uses the tolerance (T) margin of error. Tolerance is a comparatively universal kind of multiple correlation coefficient [65]. Fully multi-collinear variables have no margin of error since they are entirely predictable from all independent variables. If a variable's tolerance value is one, it has zero correlation with any of the independent variables [33,65]. The criteria VIF and T revealed that the study was multi-collinear. It is shown below.

$$\text{Tolerance} = 1 - \mathbb{R}\_{\text{j}}^{2} \tag{4}$$

$$\text{VIF} = \frac{1}{\text{Tolerance}}\tag{5}$$

Here, R<sup>j</sup> 2 is the coefficient of determination (R-squared) of the model of the descriptive variable j as the response variable, and the other explanatory variables as the independent variable.

## 2.6.3. Model Accuracy Evaluation Method

A range of statistical indicators may be used to evaluate the efficacy of landslide hazard zonation models. To assess the performance of the prediction model, many validation model evaluation approaches were employed in this study, including sensitivity, receiver operating characteristics (ROC), specificity, area under curve (AUC), and accuracy. The effectiveness of landslide susceptibility prediction has recently been extensively studied using the ROC\_AUC technique [26,66]. The inputs used to create the ROC curve were true positive, which refers to a landslide that was correctly anticipated on the axis-"X", and false positive, which refers to a landslide that was not correctly predicted on the axis-"Y". The models were statistically compared using AUC, a measure of the inclusive effectiveness of the models. ROC-AUC's accuracy is assessed as poor when the value is between 0.5 and 0.6, moderate when the value is between 0.6 and 0.7, high when the value is between 0.7 and 0.8, and exceptional when the value is between 0.8 and 0.9 [67]. Consequently, AUC values may be used to evaluate the accuracy of a prediction model.

The receiver operating characteristics (ROC) curve was used to assess the map and the model's precision [38,68–70]. The ROC curve was shown concerning the axis-"X" and axis-"Y", where the axis-"X" is the false positive rate (1-specificity) and the "Y" axis is the true positive rate (sensitivity) [25,31,70,71].

$$\text{X} - \text{axis} = \frac{\text{TN}}{\text{TN} + \text{FP}} = 1 - \text{Specificity} \tag{6}$$

$$\text{Y} - \text{axis} = \frac{\text{TP}}{\text{TP} + \text{FN}} = \text{Sensitivity} \tag{7}$$

where TP = true positive, FP = false positive, TN = true negative, FN = false negative. A ROC curve was used for the assessment of the hazard map.

## **3. Result**

## *3.1. Multi-Collinearity Analysis of LCF*

Table 4 displays the findings of the multi-collinearity connection between the ten independent components employed in this investigation. The lineament density has the highest value (VIF = 1.423), according to the findings of multi-collinearity between these variables. In contrast, the land use and land cover component variance factor has the lowest value (VIF = 1.029). In addition, all variables have VIF values less than five, and all factors have high T values. Since each independent variable uniquely impacts the dependent variable, there is no multi-collinearity among the variables used in this study. Hence, with less collinearity, every layer is used to compute the final LSM map.


**Table 4.** Multi-collinearity analysis of landslide conditioning factors.

## *3.2. Assessment with AHP for Generating Landslide Susceptibility Mapping*

The LSM map was prepared by integrating all ten layers: drainage density, rainfall, lineament density, slope, geology, LULC, geomorphology, soil, aspect, and earthquake. Each sub-factor of LCF is assigned weights for generating LSM [72]. Table 5 represents the weights of each sub-factor of LCFs. LSM is created by the AHP method in the GIS environment with ten selected LCFs by the formula

$$\text{LSM} = \sum\_{i=1}^{n} \mathbf{w}\_{i} \ast \mathbf{l}\_{i} \tag{8}$$

where l are the individual layers, w are their corresponding weights, and n is the number of layers (Table 6).


**Table 5.** LSM weights and rating system of different landslide conditioning factors.


**Table 5.** *Cont.*

**Table 6.** The final pair-wise comparison matrix of landslide conditioning factors.


Where SL = Slope, Rainfall = RF, Earthquake = EQ, LD = Lineament Density, DD = Drainage Density, GEO = Geology, GEM = Geomorphology, AS = Aspect, LULC = land use and land cover, SO = Soil.

From the outcome of the AHP model, the results were utilized to generate the final LSM and classified into low susceptibility, moderate susceptibility, high susceptibility, and very high susceptibility zones (Figure 5). Susceptibility mapping for the study area shows that the very high susceptibility zone covers 14.40% (19.312 km<sup>2</sup> ), the high susceptibility zone covers 32.20% (43.180 km<sup>2</sup> ), and the moderate susceptibility zone covers 36.77% (49.321 km<sup>2</sup> ). The low susceptibility zone covers 16.63% (22.303 km<sup>2</sup> ) of the total area. So, we observed that 83.37% of the total area comes under the moderate to very high zone. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 15 of 23

**Figure 5.** Final landslide susceptibility map of the part of the Darjeeling district. **Figure 5.** Final landslide susceptibility map of the part of the Darjeeling district.

#### *3.3. Validation of LSM and Comparison with Field Data 3.3. Validation of LSM and Comparison with Field Data*

AHP prepares an LSM of the study area. The AUC curve is drawn using a validating dataset. The dataset is divided into training (70%) and testing (30%) datasets. The testing dataset is used for the validation procedure. The AUC value obtained for this study region is 96.8%. An acceptable model has an AUC value from 0.70 to 0.8, a good hazard model has an AUC from 0.8 to 0.9, and an outstanding model has an AUC greater than 0.9. Our findings show that our model had an AUC of 96.8% (Figure 6), suggesting an extradentary model for susceptibility and functioning [57]. To evaluate the accuracy of a predictive model, it is important to compare the model's AHP prepares an LSM of the study area. The AUC curve is drawn using a validating dataset. The dataset is divided into training (70%) and testing (30%) datasets. The testing dataset is used for the validation procedure. The AUC value obtained for this study region is 96.8%. An acceptable model has an AUC value from 0.70 to 0.8, a good hazard model has an AUC from 0.8 to 0.9, and an outstanding model has an AUC greater than 0.9. Our findings show that our model had an AUC of 96.8% (Figure 6), suggesting an extradentary model for susceptibility and functioning [57].

predictions to real-world data. In the case of a landslide susceptibility map developed using the analytic hierarchy process (AHP) model, one method of validation is to use a success rate curve. The success rate curve is constructed by plotting the cumulative per-

area in decreasing the landslide susceptibility index (LSI) values [21,73]. The area under the success rate curve can be used to assess the accuracy of the prediction [53]. To validate

of 81.30% (Figure 7). The results show good accuracy of the model.

**Figure 6.** Overall model's accuracy is represented in AUC-ROC curve. **Figure 6.** Overall model's accuracy is represented in AUC-ROC curve.

To evaluate the accuracy of a predictive model, it is important to compare the model's predictions to real-world data. In the case of a landslide susceptibility map developed using the analytic hierarchy process (AHP) model, one method of validation is to use a success rate curve. The success rate curve is constructed by plotting the cumulative percentage of observed landslide occurrence against the cumulative percentage of the study area in decreasing the landslide susceptibility index (LSI) values [21,73]. The area under the success rate curve can be used to assess the accuracy of the prediction [53]. To validate the AHP model's landslide susceptibility map, 30% of the data were randomly selected and used for model validation. The results of the success rate curve analysis for the AHP model showed an area under the curve of 0.813, corresponding to a prediction accuracy of 81.30% (Figure 7). The results show good accuracy of the model. **Figure 6.** Overall model's accuracy is represented in AUC-ROC curve.

the AHP model's landslide susceptibility map, 30% of the data were randomly selected and used for model validation. The results of the success rate curve analysis for the AHP model showed an area under the curve of 0.813, corresponding to a prediction accuracy

**Figure 7.** Accuracy of the model represented in a success rate curve. and development must be performed cautiously. One of these dangers is the possibility **Figure 7.** Accuracy of the model represented in a success rate curve.

Since hilly and mountainous areas are more vulnerable to natural dangers, planning

**4. Discussion** 

## **4. Discussion**

Since hilly and mountainous areas are more vulnerable to natural dangers, planning and development must be performed cautiously. One of these dangers is the possibility of a landslide. Landslides are likely to occur under similar conditions as those observed in previous studies [74,75]. Landslide susceptibility mappings are crucial in these areas because they provide decision-makers and planners with a better understanding of the process, allowing them to take the first line of defensive action [76]. It is challenging to create an accurate LSM that can be used to designate areas prone to landslides [77]. As a consequence, several solutions are being developed daily all around the globe to handle these challenges of accuracy and reliability [78]. Due to the exhaustive research of LSM, new methodologies have been created. Our research employs MCDA algorithms to develop an accurate model of landslide hazard zonation. LCFs help to initiate landslides. It is critical to choose the appropriate LCFs when developing an accurate landslide susceptibility model. LSM models provide strong predictive capability with less inaccuracy as a consequence. There are numerous LCFs, and they vary depending on the characteristics locally and globally. These landslide indicators are linked to the land use, climatic, geological, and geomorphological variables that regulate landslides. The area's edge is dominated by the convex slope, the middle section by the straight slope, and the bottom portion by the concave slope. There are no standards for selecting LCFs based on the variety of indicators and the characteristics of the region [79]. Much effort has been spent in choosing the most appropriate and stressful factors. The multi-collinearity test helps us to find the correlation between LCFs that may impact the overall accuracy of the models. Ten (10) LCFs were chosen as independent layers in the present research to assess the susceptibility of the area to landslides. The VIF was used to evaluate the LCFs' multi-collinearity. Our findings show that the layers have no multi-collinearity between them. As a result, all variables were included in this model. According to the study's results, one of the reliable approaches for LSM is AHP, which describes the weighted-overlay analysis technique with a multicriteria decision approach. To classify the area into different susceptibility zones, ten (10) thematic layers were used, including rainfall, slope, earthquake, aspect, lithology, drainage density, lineament density, LULC, soil, and geomorphology. Figure 5 depicts a statistically generated landslide susceptibility map interpreted using the landslide hazard index (LHI) value. The model's min and max LHI values were 1393 and 5391, respectively, with a mean of 2606.88 and a standard deviation (s.d.) of 584.93 (Figure 8). Since this histogram showed that the values were unevenly distributed, the natural-break categorization approach was employed for zonation mapping [80]. As a result, four landslide hazard zones were identified and mapped: low susceptibility, moderate susceptibility, high susceptibility, and extremely high susceptibility (Figure 5). The analytical area % of the extremely high susceptibility zone for AHP is 11.68. (Figure 9). According to the hazard area, AHP falls within the high to extremely high zone with 37.86 percent (Figure 10). As a result, we may assume that the whole study area is in the high-risk zone. The AHP methodology is widely used in landslide susceptibility mapping. However, it has some limitations in the context of the selected area. One limitation of using AHP methodology for landslide susceptibility mapping in Darjeeling is that it relies heavily on expert opinion, which may be subjective and may vary from person to person. The quality and availability of data can also impact the accuracy of the AHP model, especially in areas where data are limited. Another limitation of AHP methodology for landslide susceptibility mapping in Darjeeling is that it requires extensive data inputs, which can be challenging to obtain and verify in remote and mountainous regions. The AHP method assumes that the weights assigned to each criterion remain constant over time, which may not be the case in dynamic and rapidly changing environments such as Darjeeling. The AHP method does not consider the temporal and spatial variability of landslide triggering factors, which can vary significantly depending on the location and time of year. The AHP method does not account for the complex interactions and feedback mechanisms between factors, which can lead to underestimating or overestimating the risk of landslides in certain areas. Moreover, the trustworthiness of

the LSM map is dependent on the findings obtained from the ROC-AUC and success rate curve. The GSI's released data were used to create a landslide inventory map. A total of 114 landslides were recorded, with 79 (70%) serving as training data and 35 (30%) serving as testing data. The AHP model generated maps with AUC values of 96.8 percent. The success rate curve was also plotted for determining the accuracy of the model, which was found to be 81.30 percent. Since time is an important factor in hazard research, this finding is useful in an emergency. One may conclude that the models' accuracy is comparable. We recommend using the AHP model in landslide investigations because it can generate excellent and reliable landslide hazard maps for risk reduction and management planning. serving as training data and 35 (30%) serving as testing data. The AHP model generated maps with AUC values of 96.8 percent. The success rate curve was also plotted for determining the accuracy of the model, which was found to be 81.30 percent. Since time is an important factor in hazard research, this finding is useful in an emergency. One may conclude that the models' accuracy is comparable. We recommend using the AHP model in landslide investigations because it can generate excellent and reliable landslide hazard maps for risk reduction and management planning. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 18 of 23 serving as training data and 35 (30%) serving as testing data. The AHP model generated maps with AUC values of 96.8 percent. The success rate curve was also plotted for determining the accuracy of the model, which was found to be 81.30 percent. Since time is an important factor in hazard research, this finding is useful in an emergency. One may conclude that the models' accuracy is comparable. We recommend using the AHP model in landslide investigations because it can generate excellent and reliable landslide hazard

maps for risk reduction and management planning.

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 18 of 23

**Figure 8.** Histogram analysis of the LSM model. **Figure 8.** Histogram analysis of the LSM model. **Figure 8.** Histogram analysis of the LSM model.

**Figure 9.** The graph shows the area percentage of different landslide susceptibility classes. **Figure 9.** The graph shows the area percentage of different landslide susceptibility classes.

**Figure 10.** The map shows the high to very high susceptibility zone of the study area. **Figure 10.** The map shows the high to very high susceptibility zone of the study area.

### **5. Conclusions 5. Conclusions**

One of the deadliest natural dangers in the mountainous terrain is a landslide. Due to man-made and natural events, the highly damaged mountainous landscape is vulnerable to landslides (such as earthquakes, climate change, and human intervention, respectively). Because of these threats, communities in West Bengal's northern Himalayan highlands live in a never-ending nightmare, leading to socio-economic losses. To protect property and livelihood, interim and long-term solutions to reduce the risk of landslides in this region are required. To aid in infrastructural development and socio-economic development planning, we must identify and map the most vulnerable areas. LSM may be a critical tool for assessing risk management in rugged terrain. For many years, the application of MCDA models for landslide hazard assessment has produced amazingly efficient and exact results. The primary goal of this study is to evaluate the efficacy of MCDA models AHP and to determine the most accurate and helpful approach for detecting landslideprone places in the research region. The result is based on AHP, classified into four landslide susceptibility classes. The classes are very high, high, moderate, and low landslide susceptibility. The results show that 26.18% of the total area falls under high landslide susceptibility zones, and 11.68% falls into very high susceptibility areas. So, the entire 37.86% area is in a high to very high susceptibility zone. Our findings (based on LSM and AUC) show that the usual strategy is successful. The AHP MCDA algorithm achieved 96.8 percent for the AUC-ROC curve and 81.3 percent for the success rate curve. From the result, we can identify that most areas which fall under the high to very high susceptibility classes are in the eastern, south-eastern, and southern parts of the area. The pattern of landslide sites reveals that landslides often occur on heavily dissected slopes. Less debrisfilled vegetation areas can increase the risk of landslides. The area's combined landslip One of the deadliest natural dangers in the mountainous terrain is a landslide. Due to man-made and natural events, the highly damaged mountainous landscape is vulnerable to landslides (such as earthquakes, climate change, and human intervention, respectively). Because of these threats, communities in West Bengal's northern Himalayan highlands live in a never-ending nightmare, leading to socio-economic losses. To protect property and livelihood, interim and long-term solutions to reduce the risk of landslides in this region are required. To aid in infrastructural development and socio-economic development planning, we must identify and map the most vulnerable areas. LSM may be a critical tool for assessing risk management in rugged terrain. For many years, the application of MCDA models for landslide hazard assessment has produced amazingly efficient and exact results. The primary goal of this study is to evaluate the efficacy of MCDA models AHP and to determine the most accurate and helpful approach for detecting landslide-prone places in the research region. The result is based on AHP, classified into four landslide susceptibility classes. The classes are very high, high, moderate, and low landslide susceptibility. The results show that 26.18% of the total area falls under high landslide susceptibility zones, and 11.68% falls into very high susceptibility areas. So, the entire 37.86% area is in a high to very high susceptibility zone. Our findings (based on LSM and AUC) show that the usual strategy is successful. The AHP MCDA algorithm achieved 96.8 percent for the AUC-ROC curve and 81.3 percent for the success rate curve. From the result, we can identify that most areas which fall under the high to very high susceptibility classes are in the eastern, south-eastern, and southern parts of the area. The pattern of landslide sites reveals that landslides often occur on heavily dissected slopes. Less debris-filled vegetation areas can increase the risk of landslides. The area's combined landslip categorization is

debris slide. Hence, it can be concluded that the most likely place for a landslide to occur in this study area is a severely dissected hill with debris and low to moderate vegetation, usually sparse forest, barren land, or a tea plantation area. This research successfully identified the landslide-prone locations, and the best location for planners for development in the hilly region of Darjeeling, West Bengal, India. Finally, the LSM established in this research may be used by decision-makers, land-use planners, and other governmental and non-governmental entities as an efficient tool for optimizing resource management, infrastructure development, and human activity in the studied area.

**Author Contributions:** Conceptualization, A.S.; methodology, A.S., V.G.K.V. and A.B.; software, A.S.; validation, A.S., V.G.K.V. and A.B.; formal analysis, A.S.; investigation, A.S., V.G.K.V., A.B. and S.K.; resources, A.S., V.G.K.V. and A.B.; writing—original draft preparation, A.S.; writing—review and editing, A.S., V.G.K.V., A.B. and S.K.; visualization, A.S., V.G.K.V., A.B. and S.K.; supervision, V.G.K.V. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are thankful to the Department of Mining Engineering, Indian Institute of Technology (Indian School of Mines), Dhanbad, for providing an adequate environment for the research work. Authors thanks Lakshya Tripathi for the support during the study.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Revalidation Technique on Landslide Susceptibility Modelling: An Approach to Local Level Disaster Risk Management in Kuala Lumpur, Malaysia**

**Elanni Affandi <sup>1</sup> , Tham Fatt Ng 1,\* , Joy J. Pereira <sup>2</sup> , Ferdaus Ahmad <sup>3</sup> and Vanessa J. Banks <sup>4</sup>**


**Abstract:** Landslide susceptibility modelling in tropical climates is hindered by incomplete inventory due to rapid development and natural processes that obliterate field evidence, making validation a challenge. Susceptibility modelling was conducted in Kuala Lumpur, Malaysia using a new spatial partitioning technique for cross-validation. This involved a series of two alternating east-west linear zones, where the first zone served as the training dataset and the second zone was the test dataset, and vice versa. The results show that the susceptibility models have good compatibility with the selected landslide conditioning factors and high predictive accuracy. The model with the highest area under curve (AUC) values (SRC = 0.92, PRC = 0.90) was submitted to the City Council of Kuala Lumpur for land use planning and development control. Rainfall-induced landslides are prominent within the study area, especially during the monsoon period. An extreme rainfall event in December 2021 that triggered 122 landslides provided an opportunity to conduct retrospective validation of the model; the high predictive capability (AUC of PRC = 0.93) was reaffirmed. The findings proved that retrospective validation is vital for landslide susceptibility modelling, especially where the inventory is not of the best quality. This is to encourage wider usage and acceptance among end users, especially decision-makers in cities, to support disaster risk management in a changing climate.

**Keywords:** landslide susceptibility; validation; predictive capability; disaster risk; tropical climate; Malaysia

## **1. Introduction**

Landslides triggered by rainfall have resulted in the highest number of fatalities in Asia, and whilst its attribution to climate change is limited by deficient records, landslide occurrence is expected to increase in many regions [1,2]. Development in steep hills and unstable slopes, in combination with seasonally dry periods followed by excessive rainfall are contributing factors for frequent landslides in Southeast Asia [3]. The number of people in cities and urban settlements that will be exposed to landslides is expected to increase in Southeast Asia due to increased monsoon rainfall as a result of climate change [1,4]. The expected increase of landslides calls for reliable demarcation of areas susceptible to this hazard for supporting disaster risk management in cities.

Quantitative techniques predominate in landslide susceptibility modelling. The techniques could be deterministic (analytically driven), heuristic (knowledge-driven), or data-driven statistics, encompassing bivariate, multivariate, and machine learning approaches [5,6]. The use of statistical approaches in actual practice is dependent on a reliable inventory that reflects the actual factors that cause a slope to fail [7,8]. Statistical approaches

**Citation:** Affandi, E.; Ng, T.F.; Pereira, J.J.; Ahmad, F.; Banks, V.J. Revalidation Technique on Landslide Susceptibility Modelling: An Approach to Local Level Disaster Risk Management in Kuala Lumpur, Malaysia. *Appl. Sci.* **2023**, *13*, 768. https://doi.org/10.3390/ app13020768

Academic Editors: Ricardo Castedo, Miguel Llorente Isidro and David Moncoulon

Received: 2 December 2022 Revised: 21 December 2022 Accepted: 28 December 2022 Published: 5 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

are reliable for shallow landslides when drawing on inventories that integrate field investigation, high resolution digital elevation models (DEM) and images, combined with expert inputs for more nuanced representation of local conditions [9–12].

Validation is critical for susceptibility modelling to ensure the quality and reliability of outputs for wider usage among the decision makers [13–15]. Model validation generally involves two aspects, evaluating compatibility and assessing predictive capability [5,16,17]. The area under curve (AUC) is computed from the receiver operating characteristics (ROC) curve to report results. The evaluation of compatibility (i.e., agreement between the outputs and observed data) is from success rate curves (SRC) while the assessment of predictive capability is through prediction rate curves (PRC). Partitioning of the inventory into test and training datasets is the basis for validation [15,18].

Partitioning of the inventory is handled using three types of techniques. These are random, spatial, and temporal techniques [13,14,19,20]. Temporal validation uses dataset of landslides that occurred at different time periods to assess the predicted model capability. Random partitioning involves separation of the landslide inventory into two subsets, where one subset is arbitrarily removed to assess predictability. Spatial partitioning is a comparison of two mutually exclusive sub-regions of the area under investigation, or to another area with similar setting. This technique has limitations in lithologically heterogeneous areas, where some parameters that influence landslides may be absent in one of the sub-regions, leading to low predictive capability [15,21]. Retrospective validation has recently been introduced to assess the actual accuracy of prediction from susceptibility modelling [18]. The application of this approach is expected to increase the reliability, transparency, and acceptance of susceptibility modelling.

Malaysia has tropical climate with mean daily temperatures of between 26 ◦C to 28 ◦C and annual average rainfall ranging from 2000 mm to 4000 mm, with two distinct monsoon seasons [22]. The hot and humid conditions have resulted in deep weathered soil profiles that are prone to landslides. Rainfall-induced shallow landslides are common in cities primarily due to land use change and other human activity, and several have been observed to be reactivated failures [23–25]. In some cases, failures have been linked to inadequate engineering design during construction and poor slope maintenance [23]. Landslide mapping is a challenge due to vegetation growth, active erosion, and rapid development that obliterates evidence of landslide occurrence. This is compounded by inadequate public records, which is a prevailing issue for many developing countries due to limited resources [26].

Landslide susceptibility modelling has been conducted sporadically in Malaysia over the past decade. Statistical approaches are prevalent in areas where inventories are available such as Cameron Highlands, Penang, Putrajaya, and eastern Selangor [27–36]. Generally, the inventories are limited, and the best available data are utilised. More recently, a bivariate statistics approach combined with expert inputs was used to improve the susceptibility model [9]. In many cases, model validation is not mentioned [28–30,32,33], while in others only the success rates are reported [9,34]. Investigations that report on predictive capability have relied on random partitioning and spatial correlation techniques [27,31,35,36]. All the investigations have been conducted for academic purposes and retrospective validation has not been reported.

Kuala Lumpur, the capital city of Malaysia, suffers from numerous slope failures especially after a severe rainfall episode, where residential infrastructure and business complexes are most affected [23]. Many of the incidents occur in cut and fill slopes, road embankments, highways, and hill slope development projects [28,37,38]. The situation is expected to worsen with climate change and continuous expansion of the city. This paper highlights a new spatial partitioning technique for landslide susceptibility modelling that was conducted in Kuala Lumpur, Malaysia. The technique is introduced for assessing the predictive capability of landslides, to overcome the limitations associated with lithological heterogeneity. The ensuing landslide susceptibility map was then formally submitted to the City Council of Kuala Lumpur (DBKL) to support decision-making. An extreme

rainfall event with a recorded precipitation of more than 250 mm occurred in December 2021 exceeding the monthly rainfall in over a day. The unusual phenomenon was caused by the northeast monsoon flow factors and a low-pressure weather system which caused continuous heavy rain over the west coast of peninsula regions besides the usually affected central east coast and northern peninsula regions. The event caused massive flooding of 100-year return period and triggered 122 landslides in the city. This provided an opportunity to conduct retrospective validation to confirm the robustness of the partitioning technique and increase confidence in the landslide susceptibility map.

## **2. Materials and Methods**

## *2.1. Study Area*

Kuala Lumpur covers an area of 243 km<sup>2</sup> and is the most urbanised and densely populated territory in Malaysia. The city is located in west-central Peninsular Malaysia, underlain by a variety of rock types, geological structures and geomorphological features [39]. The Dinding Schist present in northeastern Kuala Lumpur consists of meta-volcanics and quartz-mica schist (Figure 1). It is overlain by the Hawthornden Schist comprising mainly of graphitic quartz-mica schist. The schistosity of these two rock units is trending approximately north-south. The northern and east-central area is underlain by the Kuala Lumpur Limestone. The limestone is fully covered by alluvium, fill material, and mine tailing with thickness up to 66 m. Kenny Hill Formation is mainly found in the southern and central parts of the city. It comprises interbedded phyllite and quartzite, and the strata are mainly trending north–south. All the above metasedimentary sequences were intruded by granite, which is present mainly in the western and southeastern Kuala Lumpur. Soils developed over the schists, Kenny Hill Formation, and granite average 13 m, 9 m, and 15 m thick, respectively. The ground elevation obtained from the digital terrain model (DTM) of Kuala Lumpur ranges from 9 m to 320 m; about 90% of the city is below 100 m. Flat alluvial plains occur along the Kelang River and its tributaries. The alluvium has been largely mined for tin in the past. The alluvial plain is flanked by hills at the east and west. The hills in the Kenny Hill Formation areas have low relief (~50 m) and are generally elongated in the north–south direction. The hills in the granite and schist areas have higher relief (~100 m) and are also elongated in the north–south direction.

## *2.2. Input Data*

The landslide inventory consists of 650 landslide points sourced from Department of Mineral and Geoscience Malaysia (JMG). It is the sole official landslide repository that was established in 2014 through field mapping and includes only rainfall-triggered landslides. Although it contains only partial information on location and type of slope failure (manmade or natural), the original landslide points were refined through quality assessment to ensure the evidence of landslides location and its accuracy. The points were overlain on the DTM, orthorectified aerial photographs and satellite imageries from Google Earth, and repositioned at the crown of the landslides.

The DTM was derived from LiDAR data from DBKL. The 2014 LiDAR data have a resolution of 1 m, and they were resampled to a pixel size of 5 m for the study. Aerial photographs were acquired from the Department of Survey and Mapping Malaysia (JUPEM) and DBKL. Topographic maps used were at the scales of 1:50,000, 1:25,000 and 1:10,000 published by JUPEM. Bedrock geology maps obtained from JMG were updated to incorporate information on surface geology and geological structures. The additional data were obtained from boreholes and field visits as well as interpretation from aerial photographs and topographic maps.

**Figure 1.** Geological map of the Kuala Lumpur area. The city of Kuala Lumpur is located within the Selangor state, in central-west Peninsular Malaysia. **Figure 1.** Geological map of the Kuala Lumpur area. The city of Kuala Lumpur is located within the Selangor state, in central-west Peninsular Malaysia.

## *2.3. Landslide Conditioning Factors*

*2.2. Input Data*  The landslide inventory consists of 650 landslide points sourced from Department of Mineral and Geoscience Malaysia (JMG). It is the sole official landslide repository that was established in 2014 through field mapping and includes only rainfall-triggered landslides. Although it contains only partial information on location and type of slope failure (manmade or natural), the original landslide points were refined through quality assessment to ensure the evidence of landslides location and its accuracy. The points were overlain A total of fourteen intrinsic landslide conditioning factors were prepared from the geological maps, DTM, topographic maps, and satellite imageries. All the geoprocessing was done using the software ArcGIS 10.5. Bedrock geology and surface geology were extracted from the geological data. Geomorphological factors such as slope gradient, slope aspect, and slope curvature were derived directly from the DTM. The standard deviation of elevation calculated from a 10 m by 10 m moving window was used to derive the surface roughness [40]. The topographical position index (TPI) [41,42] was calculated as the difference between the elevation at a point and the mean elevation within a 250 m radius circular moving window.

on the DTM, orthorectified aerial photographs and satellite imageries from Google Earth, and repositioned at the crown of the landslides. The DTM was derived from LiDAR data from DBKL. The 2014 LiDAR data have a resolution of 1 m, and they were resampled to a pixel size of 5 m for the study. Aerial photographs were acquired from the Department of Survey and Mapping Malaysia (JU-PEM) and DBKL. Topographic maps used were at the scales of 1:50,000, 1:25,000 and 1:10,000 published by JUPEM. Bedrock geology maps obtained from JMG were updated The topographic wetness index (TWI) is a contributing factor that affects the topography and used to quantify the topographic control on hydrological processes, expressed by Equation (1), whereas the stream power index (SPI) measures the erosive power of water flow as defined by Equation (2). The value α is calculated by multiplying the flow accumulation with the cell area. The distance to stream, distance to road and distance to lineament maps were created by making multiple ring buffer around the streams, roads and lineament, respectively. Landsat 8 imagery was used to produce the normalised difference vegetation index (NDVI) map. The value can be calculated using Equation (3).

$$\text{TWI} = \ln(\frac{a}{\tan \beta}) \tag{1}$$

geological maps, DTM, topographic maps, and satellite imageries. All the geoprocessing was done using the software ArcGIS 10.5. Bedrock geology and surface geology were extracted from the geological data. Geomorphological factors such as slope gradient, slope aspect, and slope curvature were derived directly from the DTM. The standard deviation of elevation calculated from a 10 m by 10 m moving window was used to derive the

photographs and topographic maps.

*2.3. Landslide Conditioning Factors* 

$$\text{SPI} = \ln(\alpha.\tan\beta) \tag{2}$$

where *α* is the local upslope contributing area, and *β* is the slope [43].

$$\text{NDVI} = \frac{NIR - Red}{NIR + Red} \tag{3}$$

where *NIR* is the near-infrared band, and *Red* is the red band.

An expert consultation process was carried out to select the best landslide conditioning factors for the area [9]. The conditioning factors went through an iterative selection and elimination procedure by experts with experience working on landslides in tropical terrain. The evaluation was based on redundancy and relevance of the factors to the local geological and terrain conditions, as well as the initial results of bivariate statistical analysis. Only seven of the original conditioning factors were selected as shown in Table 1. They are surface geology, slope gradient, elevation, distance to lineament, distance to road, surface roughness, and TPI (Figures 2 and 3).

Most of the ground is covered by surface geology and shows a strong influence towards landslide initiation compared to bedrock geology, caused by the deep weathering profile which generates thick residual soil over most area [44]. The residual soil of Kenny Hill formation holds the highest weightage of 1.00 with 68.7% of total landslide, due to the interlayering of metasedimentary rocks of variable strength and properties while the granite represents 127 (20%) landslides with a weightage value of −0.18. Slope gradient is usually regarded as the primary causative factor for the onset of slope failure associated to shear strength of the material on slope [45]. The analysis of the parameter indicates that the landslide density and weightage values increase with increasing slope gradient. Elevation of the slope is not a conditioning factor by itself but affects the overall surface of the terrain and topographic features which controls the vegetation distribution [31]. The highest weightage is represented by the 100–150 m class with a value of 0.55 and has 10% of the total landslides while 50–100 m class represents 68.8% of the total landslides with a lower weightage of 0.49 due to the larger area. The correlation between distance to lineament parameter and landslides is explained by weakened earth material by deformation along faults. The structural stability of the surrounding area is reduced by induced regional perturbations in the fracturing occurrence and enhanced weathering in the rocks [46]. The distance to lineament factor and landslide shows direct correlation with the highest weightage (0.93) is within 0–250 m proximity. Slope which has closer proximity to roads especially near cut and fill slope could affect the stability by increase of stress in its base and accumulation of water from nearby slope [47,48]. The highest weightage of 0.18 produced from the 25–50 m class interval representing 25% of total landslides while the majority of landslides (47%) occurs within 0–25 m distance from a roadway yields a low weightage of 0.02 due to larger class area. Surface roughness generally correlates to excavation, surface erosion and the density of vegetation cover and could signify the type of land cover of that particular area and has been classified into high, moderate and low values [46]. Half of the landslides (53%) are categorised as having moderate surface roughness with a weightage of 0.07, however the highest weightage of 0.29 represented by slope with high roughness due to smaller areal class. TPI gives an indication to where the landslide point is located with reference to the topographical position or sometimes expressed as landscape position [49,50]. The highest percentage and weightage of landslides is within the middle slope with a value of 61.17% and 1.13, respectively, with the upper slope class having a weightage of −0.01 with 15.08% of total landslides.


**Table 1.** Landslide conditioning factors and their definition as used in the bivariate statistical analysis.

**Figure 2.** Seven landslide conditioning factors were selected according to the steps highlighted in the flow chart (**a**). The factors are surface geology (**b**), slope gradient (**c**), and elevation (**d**). **Figure 2.** Seven landslide conditioning factors were selected according to the steps highlighted in the flow chart (**a**). The factors are surface geology (**b**), slope gradient (**c**), and elevation (**d**).

**Figure 3.** Map of the landslide conditioning factors on distance to lineaments (**a**), distance to roads (**b**), topographic position index (**c**), and surface roughness (**d**). **Figure 3.** Map of the landslide conditioning factors on distance to lineaments (**a**), distance to roads (**b**), topographic position index (**c**), and surface roughness (**d**).

#### *2.4. Bivariate Statistical Analysis 2.4. Bivariate Statistical Analysis*

Landslide susceptibility modelling was conducted using the bivariate statistical approach that relies on the correlation between landslide population and the landslide conditioning factors. Each conditioning factor is subdivided into several classes, and a weight is calculated for each class from its statistical spatial relationship with landslide distribu-Landslide susceptibility modelling was conducted using the bivariate statistical approach that relies on the correlation between landslide population and the landslide conditioning factors. Each conditioning factor is subdivided into several classes, and a weight is calculated for each class from its statistical spatial relationship with landslide distribution [5,12,51,52].

tion [5,12,51,52]. The method proposed by [53] involves calculation of the natural logarithm of landslide density of each class within each factor divided by the overall landslide density in The method proposed by [53] involves calculation of the natural logarithm of landslide density of each class within each factor divided by the overall landslide density in the area expressed as the Equation (4).

$$W\_i = \ln\left(\frac{Density}{Density}\right) = \ln\left(\frac{\frac{Npix(S\_i)}{Npix(N\_i)}}{\frac{\sum Npix(S\_i)}{\sum Npix(N\_i)}}\right) \tag{4}$$

∑i() ∑i() where, *Wi* = the weight given to a certain factor class. *Densclass* = the landslide density where, *W<sup>i</sup>* = the weight given to a certain factor class. *Densclass* = the landslide density within the factor class. *Densmap* = the landslide density within the entire map.

within the factor class. *Densmap* = the landslide density within the entire map. i() =

*Np*i*x*(*Si*) = the number of pixels with landslide occurrence in a certain factor class. *Np*i*x*(*Ni*) = total number of pixels in a certain factor class.

The weight indicates the correlation of the class factor with the landslide occurrence and refers to landslide density of each factor class. Negative weights indicate that the landslide density is lower than average, while a positive value signifies that it is higher. A zero weight is obtained when there is an absence of landslide occurrences within a parameter class [53]. Since the analysis is based on statistical calculation of the pixels in a map, raster maps are required for the processing. Vector maps in the form of polygon shapefiles were converted into raster files for further processing.

The raster factor maps were reclassified by assigning the weight to all the factor classes. The landslide susceptibility is a summation of the reclassified factor maps. The susceptibility classes can be defined based on the landslide density [54,55]. The classification was carried out manually based on the percentage of landslide points within each class. The classes are very high (>50%), high (20–50%), moderate (10–20%), low (2–10%) and very low (<2%) susceptibility, respectively.

## *2.5. Model Evaluation and Validation*

The study employed the ROC method for model evaluation and validation. The AUC is used as the metric to evaluate the quality of the model. The susceptibility map was divided into 25 equal area classes with decreasing susceptibility and the percentage of landslide in each of class is calculated. The cumulative landslide percentage was plotted against the cumulative area to derive the ROC curve.

The landslide inventory was partitioned into two groups, i.e., training dataset and test dataset. The SRC is acquired when ROC curve of the same susceptibility map is plotted against the training set. The PRC is generated from the ROC curve derived from the susceptibility map generated from the training set against landslide points in the test set. The performance of the model is represented by a value range under the curve of 0.5 to 1.0, where 0.5 represents a test with accuracy no better than chance while value close to 1.0 which suggests an ideal model with a perfect fit [19].

The whole 2014 landslide inventory was used to generate a landslide susceptibility map and the SRC is plotted to assess the model's compatibility. Spatial cross-validation was done on the same inventory. The study area was partitioned into alternating linear zones (zone A and zone B) measuring 1 km in width as to represent the smallest geological unit in Kuala Lumpur known as schist. Another reason for the selection of the dimension is to ensure consistency since the topographic map of Malaysia has a 1 km grid. The zones are orientated east–west, perpendicular to the north-south geological and geomorphological trend. This is to ensure the inclusivity of the factors involved in both sub-regions. The validation was done twice. First zone A was used as the training set and zone B as the test set, and in second validation, the role of zones A and B is reversed. The SRC and PRC were plotted for both cases.

## *2.6. Retrospective Validation*

In December 2021, an extreme rainfall event occurred in Kuala Lumpur and its surrounding areas causing massive flood events. The event also triggered widespread landslides. Fieldwork was carried out in January 2022 to record the landslides. The research team mapped 122 landslides more precisely than the older events in the original inventory. The position of the landslides was recorded directly in the field using the application Avenza Maps, where the DTM used as the base map. Retrospective validation was conducted using a quantitative approach. The susceptibility maps and the landslide points were compared by calculating the AUC values of PRC to evaluate the performance of the landslide susceptibility maps that was generated from the original inventory.

## **3. Results**

## *3.1. Landslide Inventories*

The 2014 landslide inventory comprising 650 events is represented mainly by small rotational landslides that occurred in the Kuala Lumpur granite and rotational to translational failures observed in the Kenny Hill formation. The mechanism of the landslides is recognised as similar as they often occur together and caused sliding of material, but data for each location are not available in the inventory. The landslides within the area are rainfall-induced which occurred on both natural and man-made slope. About 70% of the landslide are less than 10 m in length, 27% between 10 m and 20 m, and only 3% are more than 20 m. The landslides occurred mainly in areas under the Kenny Hill Formation (70%) and granite (20%). Landslides are common in the hilly areas on slope with gradient of 15–25◦ (35%) and 25–35◦ (31%), in the west-central, south, and northeast Kuala Lumpur. The distribution of the 2021 landslides is similar, but there is a noticeable increase of landslides in granite (31%). Spatially partitioning the 2014 inventory for model validation revealed that 58% of landslides occurred in the first subset (zone A) while the remaining 42% were in the second subset (zone B).

## *3.2. Landslide Susceptibility Models*

Three landslide susceptibility models were produced from bivariate statistical analysis using 7 landslide conditioning factors (Figures 2 and 3). The weights calculated for these three models are listed in Table S1. In the first model (model 1), all the landslide points from the 2014 inventory were used and there was no spatial partitioning of the dataset (Figure 4a). The AUC of the SRC of this model is 0.91 (Figure 5a). The areas classified as very high and high susceptibility cover 7% and 9% of the city, respectively. These two susceptibility classes are concentrated in the northeastern, west-central, and southern parts of Kuala Lumpur; much of this area is not developed. About 7% of the area is classified as moderately susceptible. The moderately susceptible areas are distributed mainly in the vicinity of the two earlier classes. The low susceptibility areas (20%) are scattered around the earlier classes and within the very low susceptibility class. Covering 57% of the area, the very low susceptibility class is the largest. It is distributed mainly in the northern and eastern Kuala Lumpur, in relatively flat areas underlain by alluvium and mine tailings.

The next two models were used for spatial cross-validation, where the study area is partitioned into a series of two alternating linear zones (zones A and B) trending eastwest. In the initial partitioned model (model A), the susceptibility map is generated using landslide points in zone A as the training dataset and landslide points in zone B as the test dataset (Figure 4b). The role of landslide points in zones A and B are reversed in the next model (model B) (Figure 4c). The distribution of the very high to moderate susceptibility classes of these two models are similar to the unpartitioned model (model 1), which obtained using the entire dataset of the inventory (Figure 4a–c). Compared to model 1, the low susceptibility areas of both models A and B are higher (27%), and the very low susceptibility areas are lower (48%).

## *3.3. Model Evaluation and Validation*

The spatial cross-validation produced results that are similar. The AUC of the SRC for model A and model B is 0.92 and 0.89, respectively (Figure 5a). The high AUC values of these two models and the model 1 indicate good compatibility of the models derived from the 7 selected landslide conditioning factors. It also indicates a high level of inclusivity where all the parameter classes are well represented in the three models. Model 1 and model A produced similar high SRC values, which indicates excellent model fitness. Although the 0.01 difference may not be statistically significant, a possible explanation for this finding is the higher distribution of landslides within zone A subset (58%) generated a model with the closest SRC value to that of model 1.

**Figure 4.** (**a**) Landslide susceptibility map of model 1, where all the landslide points in the 2014 inventory were used to generate the model. (**b**) Landslide susceptibility map of model A, where landslides in zone A were used as training dataset and landslides in zone B as test dataset. (**c**) Landslide susceptibility map of model B where landslides in zone B were used as training dataset and landslides in zone A as test dataset. (**d**) Landslides occurred in 2021 that were used for retrospective validation are plotted on model A. **Figure 4.** (**a**) Landslide susceptibility map of model 1, where all the landslide points in the 2014 inventory were used to generate the model. (**b**) Landslide susceptibility map of model A, where landslides in zone A were used as training dataset and landslides in zone B as test dataset. (**c**) Landslide susceptibility map of model B where landslides in zone B were used as training dataset and landslides in zone A as test dataset. (**d**) Landslides occurred in 2021 that were used for retrospective validation are plotted on model A.

The PRC for model A and model B also produced AUC values of 0.90 and 0.89, respectively, indicating that both models have similarly high predictive capability (Figure 5a). The similar AUC values show that the spatial partitioning method used is suitable for model validation. The AUC of model A is marginally better than model B. In 2020, this model was selected to be formally submitted to DBKL, the local council of Kuala Lumpur, to support decision-making on land use planning and development control.

**Figure 5.** (**a**) ROC curves of spatial model evaluation and validation with associated AUC values for model 1, model A, and model B. (**b**) ROC curves of retrospective validation with associated AUC values for the same models. **Figure 5.** (**a**) ROC curves of spatial model evaluation and validation with associated AUC values for model 1, model A, and model B. (**b**) ROC curves of retrospective validation with associated AUC values for the same models.

*3.3. Model Evaluation and Validation*  The spatial cross-validation produced results that are similar. The AUC of the SRC for model A and model B is 0.92 and 0.89, respectively (Figure 5a). The high AUC values of these two models and the model 1 indicate good compatibility of the models derived from the 7 selected landslide conditioning factors. It also indicates a high level of inclusivity where all the parameter classes are well represented in the three models. Model 1 and model A produced similar high SRC values, which indicates excellent model fitness. Although the 0.01 difference may not be statistically significant, a possible explanation for The retrospective validation of all the three models using the 122 landslide events that occurred in 2021 showed a slight increase in confidence level (Figure 4d). The AUC of the PRC for model 1, model A and model B is 0.93, 0.93, and 0.92, respectively (Figure 5b). The values proved to be marginally higher than the spatial partitioning method. Comparing the PRC to that of model 1, the same AUC values could be attributed to the small sample size of the 2021 landslide events and similar landslide distribution in both models. The revalidation using the retrospective method reaffirms the high predictive ability of all three landslide susceptibility models. This indicates that the product that was submitted to the city council is reliable for decision-making purposes.

### this finding is the higher distribution of landslides within zone A subset (58%) generated **4. Discussion**

a model with the closest SRC value to that of model 1. The PRC for model A and model B also produced AUC values of 0.90 and 0.89, respectively, indicating that both models have similarly high predictive capability (Figure 5a). The similar AUC values show that the spatial partitioning method used is suitable for model validation. The AUC of model A is marginally better than model B. In 2020, this model was selected to be formally submitted to DBKL, the local council of Kuala Lumpur, to support decision-making on land use planning and development control. The retrospective validation of all the three models using the 122 landslide events that occurred in 2021 showed a slight increase in confidence level (Figure 4d). The AUC of the PRC for model 1, model A and model B is 0.93, 0.93, and 0.92, respectively (Figure 5b). The values proved to be marginally higher than the spatial partitioning method. Comparing the PRC to that of model 1, the same AUC values could be attributed to the small sample size of the 2021 landslide events and similar landslide distribution in both models. The revalidation using the retrospective method reaffirms the high predictive ability of all three landslide susceptibility models. This indicates that the product that was submitted The sole official landslide inventory established in 2014 contained only 650 landslide points that could be verified despite the numerous slope failures have been reported in Kuala Lumpur for decades, especially after a severe rainfall. This is primarily because much of the evidence in the field was obliterated. In comparison, after one unusually heavy rainfall incident, 122 events were recorded by the research team. This indicates that the number of historical landslides is most likely higher and has gone unrecorded. Furthermore, information in the original 2014 landslide inventory was limited to event locations, while more complete information was obtained from the fresh landslide scars in 2021; indicating the need for urgency in field data collection. Though the 2014 inventory was manually improved to correctly position and interpret the type and size of landslides using DTM and orthorectified aerial photographs; much time and expertise on local conditions was required. An automated approach has been advocated for the identification of landslides using high resolution DEM [7]. However, the accuracy of this method in cities is yet to be reported. In order to obtain complete and better-quality datasets, landslide inventories in tropical cities need to be routinely updated after heavy rainfalls.

to the city council is reliable for decision-making purposes. **4. Discussion**  The sole official landslide inventory established in 2014 contained only 650 landslide points that could be verified despite the numerous slope failures have been reported in Kuala Lumpur for decades, especially after a severe rainfall. This is primarily because Landslide inventories with information on landslide masses reportedly produced more accurate susceptibility maps compared to the use of landslide mass centre point [7]. However, it has also been reported that the scarp has higher predictive accuracy than the mass centre as it closely represents the pre-failure conditions [56]. The 2014 inventory utilised the position of the landslide scarp centre as it was a more prominent feature that was more readily identified. In most cases, the body and toe of the landslide is not clear due to erosion and vegetation. In this study, the inventory sourced from JMG was not of the best

much of the evidence in the field was obliterated. In comparison, after one unusually heavy rainfall incident, 122 events were recorded by the research team. This indicates that the number of historical landslides is most likely higher and has gone unrecorded. quality with respect to data on the type and extent of the landslides. Notwithstanding, a consistent sampling technique using information on landslide scarp resulted in sufficiently accurate susceptibility models.

The landslide conditioning factors used for modelling were restricted to 7, although a total of 14 factors were initially produced. This smaller selection was considered appropriate through expert consultation and resulted in a susceptibility model that had high prediction accuracy. In this investigation, the high resolution (1 m) Light Detection and Ranging (LiDAR) derived DTM was resampled to 5 m, to reduce the noise, artifacts, and processing time. This was also considered sufficient as the topographic data (DTM) was of relatively high resolution. It has been demonstrated that finer resolution of topographic data leads to more accurate and precise susceptibility models. This study supports recent findings that as long as suitable input data and techniques are selected based on the data quality and purpose, most landslide susceptibility models result in sound overall prediction accuracy regardless of the approach [7].

Retrospective validation conducted in this study indicates that the spatial partitioning technique along alternating linear zones trending east–west that was used to cross validate the susceptibility models is able to resolve the issue of heterogeneity caused by north–south trending geological and geomorphological features. All the landslide conditioning factor classes were represented in both zones A and B (Supplementary, Table S1). The difference in the area for 70% of the factor classes between the two zones is less than 10% (Supplementary Table S2). It appears that spatial partitioning techniques that are context specific, taking into account local geological and geomorphological features, result in high confidence susceptibility maps. Nevertheless, additional analysis is required using linear zones trending in other directions to confirm this outcome. Other random partitioning and spatial correlation techniques that have been previously reported in the country should also be investigated [27,31,35,36].

Notwithstanding, the landslide susceptibility map that was submitted to DBKL is reliable for decision-making as indicated by the high predictive value (AUC = 0.93) from retrospective validation. The map indicates that that 16% of the Kuala Lumpur area has very high to high susceptibility. Fortunately, much of this comprises rugged terrain that is currently not developed. The population density map of district level produced using the 2020 census data indicates the highest population per km<sup>2</sup> is distributed at the northwest region found in Wangsa Maju followed by the second highest class represented by Setiawangsa, Batu and Seputeh in southwest region (Figure 6) [57]. With reference to the landslide susceptibility model, we recognise that part of the hilly areas in Wangsa Maju, Setiawangsa and Seputeh located in high to very high susceptibility class, requiring urgent effort in hazard mitigation, disaster preparedness, and stringent actions for development. In the current scenario, targeted disaster risk management strategies could be formulated by DBKL to manage the zones that have already been developed and the vulnerable areas. This could include investing in monitoring and early warning measures and conducting preparedness programs which aim to build the resilience of communities. The undeveloped areas could be either be gazetted for protection or stringent development, where developers are required to conduct detailed geological and geotechnical investigation, with provision of appropriate disaster mitigation measures.

Resource constraints are a contributing factor to the lack of effort to maintain and regularly update the inventory in Kuala Lumpur. A complete and frequently updated inventory would facilitate further investigation including the development of hazard maps to determine the probability of landslide events. To evaluate the accuracy of bivariate statistical model in comparison with other established methods, future work using a semiquantitative approach, namely the analytic hierarchy process (AHP), is recommended to further improve the findings. In addition, the determination of rainfall threshold values that trigger landslides could be used to establish early warning for disaster risk management [58,59]. Further investigation is also required in Kuala Lumpur to understand the significant increase of landslides in areas underlain by granites in 2021, whether this is

due to unevenness in rainfall distribution or inherent characteristics of the terrain. This is critical in light of the worsening situation expected with climate change and continuous expansion of the city. Innovative partnerships could be considered, encompassing academia, the government agencies, and the local council, with involvement of the insurance and banking sectors, to maintain and conduct further investigation of landslide risks in Kuala Lumpur. Such public–private partnerships could be scaled up to cover more cities in the country, if successful. formulated by DBKL to manage the zones that have already been developed and the vulnerable areas. This could include investing in monitoring and early warning measures and conducting preparedness programs which aim to build the resilience of communities. The undeveloped areas could be either be gazetted for protection or stringent development, where developers are required to conduct detailed geological and geotechnical investigation, with provision of appropriate disaster mitigation measures.

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 14 of 18

**Figure 6.** Population density map of Kuala Lumpur at district level based on 2020 census data. **Figure 6.** Population density map of Kuala Lumpur at district level based on 2020 census data.

### Resource constraints are a contributing factor to the lack of effort to maintain and **5. Conclusions**

regularly update the inventory in Kuala Lumpur. A complete and frequently updated inventory would facilitate further investigation including the development of hazard maps to determine the probability of landslide events. To evaluate the accuracy of bivariate statistical model in comparison with other established methods, future work using a semiquantitative approach, namely the analytic hierarchy process (AHP), is recommended to further improve the findings. In addition, the determination of rainfall threshold values that trigger landslides could be used to establish early warning for disaster risk management [58,59]. Further investigation is also required in Kuala Lumpur to understand the significant increase of landslides in areas underlain by granites in 2021, whether this is due to unevenness in rainfall distribution or inherent characteristics of the terrain. This is critical in light of the worsening situation expected with climate change and continuous expansion of the city. Innovative partnerships could be considered, encompassing academia, the government agencies, and the local council, with involvement of the insurance and banking sectors, to maintain and conduct further investigation of landslide risks in The projected intensification of monsoon rainfall over Southeast Asia and continuous growth of the city is expected to increase landslide incidents in Kuala Lumpur, Malaysia. This calls for reliable demarcation of areas within the city that are susceptible to landslides to support disaster risk management. Kuala Lumpur is underlain by heterogeneous lithology and north–south trending geological structures and geomorphological features, which limits the effectiveness of conventional spatial partitioning techniques for validation in landslide susceptibility modelling. This limitation was overcome by employing a new spatial partitioning technique using a series of two alternating east–west linear zones, where the first zone served as the training dataset and the second zone was the test dataset, and vice versa. The portioning resulted in an inclusive distribution of landslide conditioning factor classes; the difference in the area of factor classes between the two zones is mainly less than 10%. The ROC curves show that the susceptibility models have good compatibility with the selected landslide conditioning factors and have high predictive accuracy. Model A with the highest AUC values (SRC = 0.92, PRC = 0.90) was submitted to DBKL for

land use planning and development control. Retrospective revalidation reaffirmed the high predictive ability (AUC of PRC = 0.93) of the landslide susceptibility model. The findings indicates that the product is reliable for decision-making purposes, supporting local level hazard and risk mitigation efforts and better ground prediction in future city development The authors suggest future work to be done using other established methods where comparison of the findings can be done, as data acquisition of landslide occurrences and precipitation data within the study area improves.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/app13020768/s1, Table S1: (a) The weights of the 14 landslide conditioning factors. Only the first seven factors were used in model 1. (b) The weights for model A where all the 377 landslide points in zone A were used as the training dataset. (c) The weights for model B where all the 274 landslide points in zone B were used as the training dataset. Table S2: Comparison of the area extent of the landslide conditioning factor classes in zone A and zone B.

**Author Contributions:** Funding acquisition, J.J.P.; resources, J.J.P.; conceptualisation, T.F.N., E.A., V.J.B. and F.A.; data collection (2014 landslide inventory), F.A.; field investigation, E.A. and T.F.N.; data analysis, E.A.; writing, E.A., T.F.N. and J.J.P.; editing, T.F.N., J.J.P. and V.J.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to thank the Newton-Ungku Omar Fund (Application ID: 59348- 455144, UM Project No. IF002–2017, UKM Project No. XX-2017-002) for supporting this research.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The datasets generated during and analysed are available from the corresponding author on reasonable request.

**Acknowledgments:** The authors would like to thank the Editor for the opportunity to be part of this special issue. The assistance of project members from multiple institutions is gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interests.

## **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **GIS-Multicriteria Analysis Using AHP to Evaluate the Landslide Risk in Road Lifelines**

**Giovanni Leonardi 1,\* , Rocco Palamara <sup>1</sup> , Francesco Manti <sup>2</sup> and Antonio Tufano <sup>3</sup>**


**Abstract:** The present paper proposes a new methodology to characterize the landslide susceptibility of the Reggio Calabria metropolitan area. For this purpose, various factors were used, such as land use, slope, rainfall, elevation, lithology, distance from roads and rivers, and thanks to the use of GIS devices and the AHP method, the landslide risk was defined for the whole territory. The values obtained were classified into four categories: low, moderate, high, and very high. They were then exported into the GIS environment to produce a landslide susceptibility map. The study carried out demonstrates the fragility of the Calabrian territory. From the results obtained, in fact, 66% of the metropolitan territory of Reggio Calabria appears to have a medium–high landslide risk.

**Keywords:** landslide risk; natural phenomena; Reggio Calabria; GIS; AHP method

## **1. Introduction**

Landslides are natural phenomena as well as man-made disasters that frequently lead to loss of human life and property, as well as causing serious damage to natural resources throughout the world. These occurs because of different associated natural hazards or anthropogenic activities. Natural phenomena include meteorological changes, such as intense rainfall, prolonged rainfall, or snowmelt and again earthquakes or volcanic eruptions. Human disturbances instead include land use alteration, deforestation, carrying out excavations, changes in the slope profile, and infrastructure constructions [1].

The effects associated with the occurrence of landslides include injuries, loss of human life, significant damage to natural resources and infrastructure, such as roads, bridges, and communication lines. To limit this, the realization of the landslide susceptibility map is fundamental to identify possible preventive measures and the planning of evacuation strategies to avoid significant damage and victims. Moreover, landslide susceptibility, hazard and risk maps are of great help to planners for selecting suitable areas to implement development schemes in any area [2,3].

Despite tremendous progress in science and technology, landslides considerably affect the socio-economic conditions of all regions of the globe. The susceptibility to landslides can be defined as the probability of the spatial occurrence of a landslide event based on the relationships between the occurrence distribution and a set of predisposing factors in a given area, such as geo-environmental thematic variables [4]. Landslides cannot be completely prevented but an accurate prediction of the landslide-susceptible areas is of particular importance since this can lead to saving natural resources as well as human lives.

Landslides occur and are controlled by one or more conditioning factors, such as topography, slope, failure mechanism, intensity of rainfall, land cover, geological formation, strength of rocks, and many more [5–7].

**Citation:** Leonardi, G.; Palamara, R.; Manti, F.; Tufano, A. GIS-Multicriteria Analysis Using AHP to Evaluate the Landslide Risk in Road Lifelines. *Appl. Sci.* **2022**, *12*, 4707. https://doi.org/10.3390/app 12094707

Academic Editors: Miguel Llorente Isidro, Ricardo Castedo and David Moncoulon

Received: 15 March 2022 Accepted: 3 May 2022 Published: 7 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The landslide susceptibility assessment can be tricky because of the difficult evaluation of both the spatial and temporal distribution of past events for large areas, mainly due to the limitations and gaps of historical records and geographic information [8].

Many different types of landslide risk assessment zonation techniques have been developed over the last decades. These methods can typically be labeled as being either qualitative or quantitative procedures, also called direct and indirect methods [9–11].

The qualitative approaches are inventory-based and knowledge-driven methods such as distribution analysis or inventory [12], geomorphic mapping [13], and map integration models [14]. They yield direct and countable results but are limited by data availability. By contrast, the qualitative approach is less complex, but the subjective estimation of experts renders its degree of confidence and the accuracy of its results questionable [15].

Direct mapping has the purpose of determining the degree of landslide susceptibility based on a qualitative approach and provides an estimate of "where" landslides are expected.

In this approach, the experience, knowledge of the ground boundary conditions, and the evaluation of pre-existing terrain maps are taken into consideration. Precisely for these reasons, direct methods present results that are subjective since they are linked to the experience and knowledge of the personnel and present many difficulties in transcribing the phenomena due to the occurrence of particular environmental changes in certain areas [16]. This methodology usually uses factors characteristic of the locations where landslides have occurred in the past, together with an assessment of the areas with a potential to experience land sliding in the future, but with no assessment of the landslide's occurrence frequency.

The reliability of landslide-susceptibility maps mainly depends on the scale of work, the quality of the data and the analysis methodology.

On the other hand, indirect mapping is based on quantitative models, such as statistical and deterministic ones, which aim to find numerical correlations between the various factors that favor the occurrence of landslides and their distribution on the earth's surface and based on information obtained from the interrelation between landslide conditioning factors and the landslide distribution. For large areas, statistical methods are the most widely used [17]. The statistical methods include bivariate statistical analysis [18], logistic regression [19–21], multivariate regression [22,23]; multivariate adaptive regression spline [24], statistical index [25,26], analytical hierarchy process analytical hierarchy process [27,28], weight of evidence [29,30], and evidential belief function [31,32].

In recent times, quantitative approaches have been widely used for landslide susceptibility and hazard evaluation. However, it is mostly very difficult to include temporal probability in this analysis of large areas, due to the heterogeneity of subsoil conditions, the numerous conditioning factors, and the absence of a complete historical record of the events that can determine the occurrence of landslides (for example precipitation, earthquakes, and landslides themselves).

Furthermore, these models, however, show a high sensitivity of the results to the quality and accuracy of the thematic data and the specificity of some sensitive factors that could end up not being taken into consideration.

Landslide susceptibility and hazard assessments often use multi-criteria decision analysis (MCDA) techniques since in most cases, the types and format of data that are available are qualitative and quantitative, thereby requiring a semi-quantitative method that incorporates both types of data [33].

The semi-quantitative landslide assessment methods can be considered an effective expert's tool for weighting and ranking the chosen factors, which represent the main causes for landslide susceptibility of the study area, in an objectively optimal and simple way.

Some qualitative methods become semi-quantitative by incorporating ranking and weighting [34,35], as is the case of the analytic hierarchy process (AHP), a multi-objective and multi-criteria decision-making methodology which has been widely applied for the solution of decision problems [36].

This method is based on the analytical hierarchy of involved factors and the comparison between the various pairs of them in order to enable the assignment of a relevant ratio for each factor. Thanks to this method, it is therefore possible to estimate the weight of each factor considered, through the linear correlation of each factor with respect to the others. The ability to correlate different factors has made this method a valuable tool for many researchers in compiling landslide susceptibility maps, which are obtained by correlating and comparing a large number of factors [37].

The management of a large number of factors to be correlated and estimated for the definition of landslide susceptibility is therefore carried out, in most cases, through the use of the geographical information system (GIS) [35].

The increasing popularity of GIS and its ease of use over the last decades has led many studies to use indirect susceptibility mapping approaches.

GIS, which uses data-integration techniques, is a very suitable tool for landslide susceptibility mapping. In fact, with the increasing availability of high-resolution spatial data sets, GIS, remote sensing, and computers with large and fast processing capacity, it has becoming possible to partially automate the evaluation of landslide hazard and susceptibility mapping process and thus minimize fieldwork.

The reliability of these maps depends mostly on the applied methodology as well as on the available data used for the hazard risk estimation. Moreover, GIS is an excellent and useful tool for mapping the susceptibility of an area prone to landslide manifestation.

This study focused on producing a landslide susceptibility map of the Reggio Calabria metropolitan area (Italy) by combining GIS techniques and the AHP method.

The main objectives of the present work can be synthesized to create a landslide susceptibility map for the province of Reggio Calabria, using the GIS software by means of the weighted combination of various factors, such as the slope, lithology, elevation, rainfall, land use, distance of the road and river. From the obtained map, it is possible to evaluate the landslide susceptibility in the areas that are particularly relevant because of the connections between the internal urban areas and the main towns and services located along the coasts and to analyze the degree of landslide risk of the connecting infrastructures (lifelines). This method could be used in other areas of the Calabria Region for the realization of the landslide susceptibility map in order to be able to use the maps obtained to plan the safety of the slopes or, in particular, identify the main road infrastructures more at landslides risk, comparing the results obtained with the various inventories of landslides that occurred in the past.

## **2. Study Area**

The Reggio Calabria metropolitan city, in the Calabria region in the south of Italy, was selected as the study area because the territory is subject to phenomena of hydrogeological and seismic risk; the present study area has suffered from every landslide type.

The territory of Calabria is geologically young and often subject to natural modifications, so hydro-geological disaster (landslides and floods) is one of the risk factors to which Calabria is exposed. This is due, among other conditions, to the physical conformation of the region and the climatic conditions.

Landslides are the most common consequence of the soil instability due to the seismic activity and the hydrogeological problems in Reggio Calabria metropolitan city. Extraordinary conditions, such as earthquakes, and more ordinary meteorological conditions, such as rainfalls, can both induce landslides, rock fall, or debris flow.

The considered study area with 97 cities and towns and with a total population of 530,000 inhabitants is the area of Calabria with the highest population density equal to 165.57 residents/km<sup>2</sup> .

The Reggio Calabria metropolitan city covers an area of 3183 km<sup>2</sup> , of which 1685 km<sup>2</sup> (52.95%) is represented by hilly terrain, 1.275 km<sup>2</sup> (40.07%) is mountainous and the remaining 223 km<sup>2</sup> (6.97%) are represented by land lowland.

In regions where urban residential areas coincide with mountainous terrains, like in this area, the risk is higher for people, and the economic costs include relocating communities and repairing physical structures.

The Autorità di Bacino of Calabria Region, considering the various mentioned problems, has made available several thematic maps for a complete hydrogeological balance (PAI—Piano di Assetto Idrogeologico, 2001) [38]. These maps, updated in 2016, are accessible through Quantum—GIS, and they describe the landslide hazard and risk areas in the whole region, but they do not give any information about the landslide risk assessment in the strategical infrastructures identified as lifelines [2]. Furthermore, these assessments were carried out without taking into consideration the main connecting and emergency road infrastructures, which, however, play a role particularly relevant during an emergency to allow rapid and efficient access, assistance, and rescue, guaranteeing evacuation and, more in general, maintaining access for all emergency services [2,3].

Therefore, starting from the data provided by the PAI, containing the landslides that have occurred in the past together with an assessment of the areas with a potential to experience land-sliding in the future, and with the integration of factors of particular importance, such as distance from roads and from the rivers, the objective of the present work is to create a landslide susceptibility map for the metropolitan area of Reggio Calabria by means of the weighted combination of various indices using the AHP method.

## **3. Materials and Method**

The first step in every susceptibility assessment consists of collecting all available information and data for the study area, and this stage may be the most important part.

Evaluating the relationship between the landslide occurrence and the conditioning parameters becomes very important for the landslide susceptibility mapping; further analysis of the cause–effect relationships is not always simple, as a landslide is seldom linked to a single cause. So, to produce the landslide susceptibility map in the metropolitan area of Reggio Calabria, a total of 7 inputs were selected for the model, considering the main influencing factors of landslides: the slope, the lithology, the elevation, the rainfall, the land use, the distance of the road, and the river.

In the study carried out, the number of factor classes was increased compared to that of previously published studies, analyzing what was done in the literature by other authors and wanting to obtain a non-homogeneous weight of each class of factor, unlike what was obtained in the previous studies.

## *3.1. Slope Factor*

The slope and the aspect of the slopes play an important role in the occurrence of landslides because they represent the result of the combined influence of many agents. Slope is an important factor in the analysis of landslide; in fact, as the slope increases, the probability of the occurrence of landslide increases because as the slope angle increases, the shear stress of the soil increases.

A digital elevation model (DEM) was utilized using 3D extension of ArcGIS, and the slope angle was extracted from it, using contours with 5 m intervals digitized from topographic sheets and saved as a line layer.

Six slope categories were used as factor classes for the analysis of landslides risk as shown in Table 1.


**Table 1.** Weight of the individual classes of landslide hazard factors.

## *3.2. Lithology Factor*

Lithology is a further important parameter with regard to landslide manifestation. It is widely recognized that geology greatly influences the occurrence of landslides because lithological and structural variations often lead to a difference in the strength and permeability of rocks and soils.

For the study area, three factor classes were determined by grouping the different geological formations according to their geological engineering behavior and according to their physical and mechanical characteristics.

Thus, lithology includes three classes as follows: sedimentary deposits, clays and sandstones and rocks.

## *3.3. Elevation Factor*

The altitude does not contribute directly to landslide manifestation, but in relation to the other parameters, such as precipitation, the altitude contributes to landslide manifestation and influences the whole system.

Elevation is useful to classify the local relief and to locate points of maximum and minimum heights within terrains. The elevation is considered the fourth important parameter in the classification of landslide risk. The grid maps of the altitude with cell size 5 × 5 m

were produced from the DEM. The separation of the altitude into 5 classes was as follows: <100 m, 101–300 m, 301–500 m, 250–500 m, and >1000 m. The highest density corresponds to the class with an elevation range of 500–1000 m, representing hilly areas.

## *3.4. Rainfall Factor*

As it is well known, precipitation is among the most usual triggering factors for landslide manifestation. For the necessities of this study, the precipitation map was produced, using the annual average precipitation data of the main meteorological stations well distributed in the study area.

This map was separated into 5 classes as shown in Table 1.

## *3.5. Land Use Factor*

The data for the land use in an area are a parameter that seriously affects the slope failures, as slope stability is very sensitive to changes in vegetation. For the necessities of this study, the land use, which reflects the vegetation covering, was classified into 4 categories as follows: artificial and waterproof surfaces, arable agricultural areas with sparse vegetation, agroforestry areas with scattered trees and shrubs, wooded areas, or areas with important vegetation.

## *3.6. Distance of Road and River Factor*

As it is obvious, the artificial and natural parts of the slopes around a road are more sensitive in landslide manifestation. In addition, some lifeline roads have to guarantee effectiveness and efficiency in the immediate aftermath of a natural disaster; in some cases, they have to be maintained also during the event to allow rapid and efficient access and assistance and rescue, guaranteeing evacuation and, more in general, maintaining access for all emergency services. Therefore, the road network was chosen as a principal parameter that consists of the road network of the province of Reggio Calabria, composed of 1850 km of roads.

Buffer zones were created around road lifelines at distances of 100, 200, 300, 400 and 500 m.

Similar to the road network, the hydrographic network was digitized and saved as line layers in the GIS database, using the topographic sheets as a data source, and buffer zones were created around the bed of the rivers and the streams of the area, at distances of 50, 100, 200, and 250 m. These distances of river are measured from the river's bed boundaries from both sides.

Each factor was characterized in classes whose weight (**Wi,L**) was determined on the basis of the portion of the territory occupied by the class of the factor in the entire metropolitan area of Reggio Calabria (Table 1).

This method was used to identify the amount of territory that the factor class occupies to then be put into relation with the weights of the same calculated with other methods (**Wi,PAI, Wi,AHP**) and to be able to determine the degree of landslide susceptibility.

The analysis was carried out by means of GIS devices, which allow to easily analyze a considerable amount of data and to convert the map pixels into data sets.

The liability and accuracy of the collected data also influence the success of the applied methodology. For this reason, raster images were used with an accuracy of a pixel equal to 5 m. The images obtained are shown in Figure 1.

(**a**)

**Figure 1.** *Cont*.

(**b**)

**Figure 1.** *Cont*.

(**c**)

**Figure 1.** *Cont*.

(**d**)

**Figure 1.** *Cont*.

(**e**)

**Figure 1.** *Cont*.

(**f**)

**Figure 1.** *Cont*.

(**g**)

**Figure 1.** *Cont*.

(**h**)

**Figure 1.** Topographical parameter maps: (**a**) slope degree; (**b**) land use; (**c**) elevation; (**d**) rainfall; (**e**) geology; (**f**) river distance; (**g**) road distance, (**h**) PAI landslide risk areas. **Figure 1.** Topographical parameter maps: (**a**) slope degree; (**b**) land use; (**c**) elevation; (**d**) rainfall; (**e**) geology; (**f**) river distance; (**g**) road distance, (**h**) PAI landslide risk areas.

So, the weight of each class was determined by the examination of the value attained by each factor class in the landslide areas identified by the PAI, with respect to the past landslide events (**Wi,PAI**) that occurred in the territory of the metropolitan area (Table 2). This parameter allows to determine the importance of each individual class based on historical data relating to landslides that occurred in the past. So, the weight of each class was determined by the examination of the value attained by each factor class in the landslide areas identified by the PAI, with respect to the past landslide events (**Wi,PAI**) that occurred in the territory of the metropolitan area (Table 2). This parameter allows to determine the importance of each individual class based on historical data relating to landslides that occurred in the past.


**Table 2.** Weight of classes of factor in the landslide areas identified by the PAI.

Later using the AHP method, we were able to obtain a relative significance of the relevant factors after the pairwise comparison matrix was constructed.

The pairwise comparison method was developed by Saaty in the context of the analytical hierarchy process. The AHP approach allows assessing the relative weight of multiple elements in an initiative manner [39].

After defining the weight of each factor, with respect to the past landslide events identified by the PAI, its relative weight was calculated with the AHP method. The rating score of relative significance was set from 1 to 9, indicating less important to much more important factors.

The pairwise comparison matrix is shown in Table 3 using a 7 × 7 matrix, where diagonal elements are equal to 1.

The row describes the importance of factor, and the values of each row are compared with each column to define the relative importance to obtain a rating score. For example, the slope is significantly more important than the elevation and therefore assigned the value 8. Conversely, the importance of the elevation with respect to the slope is equal to 1/8.


**Table 3.** Pairwise comparison matrix.

After evaluating and confirmed the consistency of the created eigenvector matrix for the AHP method, the weight for each factor using this method was obtained (Table 4).

**Table 4.** Weight of normalized landslide hazard factors through AHP method.


The landslides index for each individual pixel was finally obtained by the weighted sum of the weight determined on the basis of the portion of territory occupied by the factor class (**Wi,L**), the weight reached by each factor class in the landslide areas identified by the PAI (**Wi,PAI**), and the weight of the same factor obtained from the AHP method (**Wi,AHP**).

$$\mathbf{LI} = \sum\_{\mathbf{i}=1}^{n} \mathbf{W\_{\bar{i},L}} + \mathbf{W\_{\bar{i},PAI}} + \mathbf{W\_{\bar{i},AHP}}$$

The values thus obtained for the assessment of susceptibility were classified into five ranges; in particular, values below 0.17 represent the portions of territory with low risk, values between 0.1701 and 0.25 represent the portions of territory with moderate risk, values between 0.2501 and 0.3 represent the portions of territory with high risk, and values greater than 0.301 represent the portions of the territory with very high risk.

The above-mentioned results highlight that 61% of the whole territory of the metropolitan area is affected by a high landslide susceptibility value. In particular, the surface at low risk of landslide is equal to 7% of the entire metropolitan area, at moderate risk is equal to 32%, at high risk is equal to 36%, and at very high risk is equal to 25%.

The values obtained for the assessment of susceptibility were classified into five ranges in such a way that, having identified the lower and higher values in each range considered, the number of pixels with respect to the total was equalized between the various classes. This ensures that each class range has approximately the same number of values in each class and that the change between intervals is fairly consistent [40].

From the obtained map (Figure 2), it is possible to notice that the highly and very highly susceptible zones are the hills, i.e., the areas between the coasts and the central part of the province. These areas are particularly relevant because of the connections within the various internal urban areas and the main towns and services located along the coasts. Further, the landslide susceptibility evaluation procedure adopted resulted to be in agreement with the work performed by the Calabria Region Basin Authority: 80% of the highly and very highly susceptible areas coincide with the landslide areas identified by the PAI maps.

**Figure 2.** Landslides susceptibility map. **Figure 2.** Landslides susceptibility map.

### **4. Results 4. Results**

the PAI maps.

The study carried out demonstrates the fragility of the Calabrian territory. From the results obtained, in fact, 38% of the metropolitan territory of Reggio Calabria appears to have a medium landslide risk and 28%, a high landslide risk. The results obtained from the model are congruent with the landslide risk areas identified by the PAI and are more reliable than previous studies [5]. In particular, the model produces an overestimation of only 8% of the areas at risk of medium–high landslide, compared to the areas surveyed by the PAI, thus also managing to determine the risk of landslides in those areas not classified by the PAI due to the absence of exposed elements. These areas, however, are still of fundamental importance, as they are crossed by connecting infrastructures. The study carried out demonstrates the fragility of the Calabrian territory. From the results obtained, in fact, 38% of the metropolitan territory of Reggio Calabria appears to have a medium landslide risk and 28%, a high landslide risk. The results obtained from the model are congruent with the landslide risk areas identified by the PAI and are more reliable than previous studies [5]. In particular, the model produces an overestimation of only 8% of the areas at risk of medium–high landslide, compared to the areas surveyed by the PAI, thus also managing to determine the risk of landslides in those areas not classified by the PAI due to the absence of exposed elements. These areas, however, are still of fundamental importance, as they are crossed by connecting infrastructures.

the various internal urban areas and the main towns and services located along the coasts. Further, the landslide susceptibility evaluation procedure adopted resulted to be in agreement with the work performed by the Calabria Region Basin Authority: 80% of the highly and very highly susceptible areas coincide with the landslide areas identified by

This method could be used in other areas of the Calabria Region for the realization of the landslide susceptibility map in order to be able to use the maps obtained to plan the This method could be used in other areas of the Calabria Region for the realization of the landslide susceptibility map in order to be able to use the maps obtained to plan the safety of the slopes or, in particular, identify the main road infrastructures on which to carry out work risk mitigation (lifelines). The results of this preliminary analysis can give important information on the relative criticality of the different road sectors, thereby allowing attention and economic budgets to be shifted toward the most critical aspects, where structural and non-structural mitigation measures could be implemented. In such a

sense, they supply useful information regarding the intervention priorities and, above all, a scale of inspection surveys and analyses of a structural nature which can find the effective level of susceptibility of the infrastructure, as well as the most opportune mitigation actions for the reduction of risk levels. This latter aspect will be the subject of future studies.

**Author Contributions:** Data curation, F.M.; Methodology, G.L., R.P., F.M. and A.T.; Writing—original draft, G.L., R.P. and A.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article* **Landslide Susceptibility Assessment Considering Landslide Volume: A Case Study of Yangou Watershed on the Loess Plateau (China)**

**Hang Gao <sup>1</sup> and Xia Zhang 2,\***


**Abstract:** Because of the special geological conditions on the Loess Plateau, Landslide erosion is not only the main goal of prevention and control of geological disasters, but also an important erosion mode of soil and water loss in the basin. Thus, landslide susceptibility assessment before only considering landslide frequency is not far enough for a decision-maker. The study aims to consider both frequency and scale of landslides for a better landslide susceptibility evaluation. Taking the Yangou small watershed as an example, this study used a VR model, RIRA method, and the GIS method to comprehensively consider frequency and scale to analyze landslide susceptibility of the small watershed. Based on the detailed analysis of the existing literature, slope, elevation, NDVI, land-use, lithology, amount distant to road, amount distant to river, profile curvature, and rainfall as landslide are selected as the conditioning factors (CFs) of the landslide, to draw the sensitivity map. The map of landslide susceptibility was classified into five zones: very low, low, medium, high, and very high, and the cover areas occupy 6.90, 12.81, 12.83, 9.42, and 5.87 km<sup>2</sup> , respectively. A total of 60% of the landslide occurred in the zones of high and very high susceptibility, accounting for 87% of the total volume in the study area. The very high susceptibility is the area with a larger relief and along the river and road. The findings will help decision makers to formulate scientific comprehensive policies that take into account disaster prevention and soil conservation measures in specific regions.

**Keywords:** landslide susceptibility; volume ratio; reserve increase-rate-analysis method; GIS; Loess Plateau

## **1. Introduction**

Landslides, a mass failure on steep slopes, facilitated by gravity are an important process controlling the sedimentary structures and growth patterns of the steep slope [1], and major natural disasters that often cause human and economic losses due to natural forces and human actions [2]. There are mainly three manifestations such as falling, sliding, or flowing [3]. China is covered by a vast area that is considered to be one of the most landslide-prone regions in the world, while landslide phenomena cause an estimated 700 to 1000 deaths every year and more than RMB 10 billion annually in infrastructure and property damage [4]. Landslides, especially on the Loess Plateau, have been widespread due to their topography and geographical location. In order to prevent and mitigate the damage, it is necessary to need comprehensive information regarding landslides, such as occurrence, spatial distribution, and susceptibility. Landslides also deliver huge amounts of sediment to rivers in the mountainous and hilly watersheds [5,6]. In particular, most landslides occur on slopes, whether man-made or natural, where loose material is always transported to the toe and downhill, which may have an out-of-field effect on sediment

**Citation:** Gao, H.; Zhang, X. Landslide Susceptibility Assessment Considering Landslide Volume: A Case Study of Yangou Watershed on the Loess Plateau (China). *Appl. Sci.* **2022**, *12*, 4381. https://doi.org/ 10.3390/app12094381

Academic Editors: Miguel Llorente Isidro, Ricardo Castedo and David Moncoulon

Received: 24 March 2022 Accepted: 24 April 2022 Published: 26 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

transport through river flows to the downhill [7]. On the Loess Plateau, the phenomenon is particularly serious, which is one of the reasons for the poor soil erosion in the region due to unmanageable, heavily containing sediment flow [8]. It can be seen that landslides have a great effect on water and soil loss in the basin. Older land use policies may not always reflect optimal planning for land use prone to landslides [9]. A reasonable zoning prediction of regional landslides can therefore assist planners and engineers in making decisions on the use of such lands [6].

Landslide susceptibility is recognized as the propensity of an area to generate landslides [10]. Generally, landslide occurrence in each region is a function of various factors. Each of the factors and function has a different influence. In order to assess susceptibility from gravity, it is, therefore, necessary to identify and analyze the factors leading to gravity [11]. The occurrence of landslides is controlled both by a series of predisposing factors, e.g., geological, geomorphological, climatic and hydro-geological, and triggering factors, including seismicity, heavy rainfall, human activities, and freeze-thaw [12]. Although it is still difficult to predict a landslide event in space and time, a region may be divided into sub-regions with homogenous properties, which are classified and ranked according to the potential hazard degrees of group movement considering the prerequisite factors explained above. A landslide susceptibility map is considered as an effective solution, and it is also an important task for decision making on regional planning and protection.

A variety of methods can be employed to develop the landslide susceptibility assessment. There are two main types: a "knowledge-driven model" that uses relevant domain knowledge for analysis and a "data-driven model" that uses intelligent methods to automatically find patterns from data. The former methods include fuzzy logic [13], fuzzy comprehensive evaluation [14], and the analytic hierarchy process [15]. The latter methods include logistic regression analysis [16], decision tree [17], random forest [18], artificial neural network [19], and support vector machine [20]. In general, the knowledge-driven model is more subjective while the data-driven model is relatively objective and can accurately reflect the correlation between landslide susceptibility and its basic environmental factors. However, the data-driven model is complex in the model training and testing process and requires a large number of training samples. Both types of models can get a more reasonable result, however; they mainly provide information on the probability of landslide occurrence, but no information on the size of a potential landslide once occurring. In other words, the landslide susceptibility map applied by those methods only reflected the possibility of landslide occurrence, while it cannot express the probability of a landslide. Research has shown that there is a significant difference in spatial variation between the susceptibility and erosion risk of a mass movement in the watershed [21]. Thus, the results only consider landslide frequency as incomplete for the decision-makers to implement land planning or other policy in the region. Therefore, it will also be meaningful to understand the sediment yield and prevent soil erosion in the catchment when both considering the scale and the frequency of the landslide.

If the volume of the landslide is considered in the assessment of landslide susceptibility, it is necessary to choose a suitable method. Whatever the relationship between landslide and the landslide-related factors, it can be inferred from the relationship between landslide and the landslide-related factors that had occurred. Thus, volume ratio (VR), which reflected the relationship between the volume of the landslide and each landslide conditioning factors, is used in the study, and the revised increase-rate-analysis method (RIRA) method is also a suitable way for quantitative analysis of the relationship between landslide scale and various factors. The RIRA method is an improvement on the increaserate-analysis method (IRA), which is used to evaluate the sensitivity of the landslide to a single variable [1]. On this basis, Lucas [22] has proposed a revised increase-rate-analysis method (RIRA) based on the importance ranking of influencing factors and sensitivity analysis for a large number of data-influencing factors. The method has been successfully applied in the analysis of the influencing factors of the silt dam conditions in the Culiacan basin in northwestern Mexico. Recently, as a basic analysis tool for landslide hazards,

ArcGIS can effectively perform spatial data management and manipulation for the analysis and has been used in many studies [16]. It is also an indispensable tool for landslide susceptibility mapping. a landslide susceptibility map for the Yangou small watershed on the Loess Plateau of China; (c) assessing the landslide susceptibility and providing suggestions for landslide prevention and soil erosion of the small catchment.

basin in northwestern Mexico. Recently, as a basic analysis tool for landslide hazards, ArcGIS can effectively perform spatial data management and manipulation for the anal‐ ysis and has been used in many studies [16]. It is also an indispensable tool for landslide

The objective of this study included: (a) applying the VR model, the RIRA method, and GIS in the spatial event prediction based on considering landslide scale; (b) creating

*Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 3 of 17

The objective of this study included: (a) applying the VR model, the RIRA method, and GIS in the spatial event prediction based on considering landslide scale; (b) creating a landslide susceptibility map for the Yangou small watershed on the Loess Plateau of China; (c) assessing the landslide susceptibility and providing suggestions for landslide prevention and soil erosion of the small catchment. **2. Study Area** The study focuses on the Yangou watershed, a small catchment located in Yanan city in the north Shanxi province of China (Figure 1). It is enclosed within the latitudes 36°28′

### **2. Study Area** N and 36°34′ N and longitudes 109°27′ E and 109°33′ E, covering an area of about 48 km2.

**3. Methods**

susceptibility mapping.

The study focuses on the Yangou watershed, a small catchment located in Yanan city in the north Shanxi province of China (Figure 1). It is enclosed within the latitudes 36◦280 N and 36◦34<sup>0</sup> N and longitudes 109◦27<sup>0</sup> E and 109◦33<sup>0</sup> E, covering an area of about 48 km<sup>2</sup> . The catchment is typical of the gullied-hill zone and is located in the central area of the Loess Plateau, which is a sub-watershed in the Yanhe basin. The river is a secondary branch of the Yan River and flows from southeast to northwest. The watershed exhibits complex topographic variations, with a gully density of 4.8 km km−<sup>2</sup> and an elevation ranging from 988 m to 1404 m. The terrain gradient of the basin is mostly composed of steep hill slopes. The slope analysis showed that the slope distribution in the study area ranged from 0 ◦ to 51◦ . The climate belongs to the transitional zone from warm semi-dry to semi-wet. The average annual precipitation is 550 mm, and 70% of the annual precipitation occurs from June to September. Generally, the southern section had higher values than the north. The geotechnical medium in the area can be divided into rock mass and soil mass, mainly sandstone and loess. The saturated bulk density (γ), Poisson's ratio, cohesion (C), and friction angle (ϕ) of sandstone are 2.68 g/cm<sup>3</sup> , 0.22, 4.23 MPa, and 41.6◦ , respectively. The moisture content (ω), dry bulk density (γ), cohesion (C), and friction angle (ϕ) of loess are 15.6%, 1.51 g/cm<sup>3</sup> , 39.48 MPa, and 26.95◦ . The catchment is typical of the gullied‐hill zone and is located in the central area of the Loess Plateau, which is a sub‐watershed in the Yanhe basin. The river is a secondary branch of the Yan River and flows from southeast to northwest. The watershed exhibits complex topographic variations, with a gully density of 4.8 km km−<sup>2</sup> and an elevation ranging from 988 m to 1404 m. The terrain gradient of the basin is mostly composed of steep hill slopes. The slope analysis showed that the slope distribution in the study area ranged from 0° to 51°. The climate belongs to the transitional zone from warm semi‐dry to semi‐wet. The average annual precipitation is 550 mm, and 70% of the annual precipi‐ tation occurs from June to September. Generally, the southern section had higher values than the north. The geotechnical medium in the area can be divided into rock mass and soil mass, mainly sandstone and loess. The saturated bulk density (γ), Poisson's ratio, co‐ hesion (C), and friction angle (φ) of sandstone are 2.68 g/cm3, 0.22, 4.23 MPa, and 41.6°, respectively. The moisture content (ω), dry bulk density (γ), cohesion (C), and friction angle (φ) of loess are 15.6%, 1.51 g/cm3, 39.48 MPa, and 26.95°.

**Figure 1.** Study area in China. (**a**) Shaanxi province in China; (**b**) Yan'an city in the Shaanxi province; (**c**) Yangou watershed in the Yan'an city. **Figure 1.** Study area in China. (**a**) Shaanxi province in China; (**b**) Yan'an city in the Shaanxi province; (**c**) Yangou watershed in the Yan'an city.

The landslide susceptibility maps were prepared using VR and RIRA models to as‐ sess the landslide susceptibility in the Yangou catchment. Figure 2 showed the flow chart

## **3. Methods**

The landslide susceptibility maps were prepared using VR and RIRA models to assess the landslide susceptibility in the Yangou catchment. Figure 2 showed the flow chart of the application method, including the data collection, gathering of the landslide conditioning factors, determination of relationship between each condition factors' subclasses, and landslide inventory map using VR, application of the RIRA method in the landslide sensibility, landslide susceptibility assessment of numerical data integration, and landslide susceptibility mapping. of the application method, including the data collection, gathering of the landslide condi‐ tioning factors, determination of relationship between each condition factors' subclasses, and landslide inventory map using VR, application of the RIRA method in the landslide sensibility, landslide susceptibility assessment of numerical data integration, and land‐ slide susceptibility mapping.

*Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 4 of 17

**Figure 2.** Flowchart of the used methodology. **Figure 2.** Flowchart of the used methodology.

The landslide inventory database, which included 27 locations, volumes, and types were provided by the Xi'an Geological Survey Center of the China Geological Survey. The landslide is divided into four sections: ≤1 × 104, 1 × 104–10 × 104, 10 × 104–100 × 104, and ≥100 × 104 m3, corresponding to tiny, small, medium, and large landslides. The amount of the landslide and its percentage in each interval were counted. In this study, nine param‐ eters were chosen according to the literature review and the general environment of the study area. These conditioning factors are land‐use, slope, profile curvature, altitude/ele‐ vation, distance to roads, distance to rivers, lithology, normalized difference vegetation index (NDVI), and rainfall. The altitude/elevation, slope, and profile curvature were ex‐ tracted from the Digital Elevation Model (DEM) (cell size: 25 m × 25 m), which is converted from a topographic map with a scale of 1:100,000. They were also collected from the "Xi'an Geological Survey Center of China Geological Survey". The NDVI was extracted from the The landslide inventory database, which included 27 locations, volumes, and types were provided by the Xi'an Geological Survey Center of the China Geological Survey. The landslide is divided into four sections: <sup>≤</sup><sup>1</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> , 1 <sup>×</sup> <sup>10</sup>4–10 <sup>×</sup> <sup>10</sup><sup>4</sup> , <sup>10</sup> <sup>×</sup> <sup>10</sup>4–100 <sup>×</sup> <sup>10</sup><sup>4</sup> , and <sup>≥</sup><sup>100</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , corresponding to tiny, small, medium, and large landslides. The amount of the landslide and its percentage in each interval were counted. In this study, nine parameters were chosen according to the literature review and the general environment of the study area. These conditioning factors are land-use, slope, profile curvature, altitude/elevation, distance to roads, distance to rivers, lithology, normalized difference vegetation index (NDVI), and rainfall. The altitude/elevation, slope, and profile curvature were extracted from the Digital Elevation Model (DEM) (cell size: 25 m × 25 m), which is converted from a topographic map with a scale of 1:100,000. They were also collected from the "Xi'an Geological Survey Center of China Geological Survey". The NDVI was extracted from the Landsat 8/OLI data set, which is provided by the Geospatial Data Cloud site, Computer Network Information Center [23]. The computational formula of NDVI is

$$\text{NDVI} = (NIR - R) / (NIR + R) \tag{1}$$

NDVI / *NIR R NIR R* (1) where *NIR* is the reflection value of the near‐infrared band, and *R* is the reflection value where *NIR* is the reflection value of the near-infrared band, and *R* is the reflection value of the red band. The landuse data were provided by the Resource and Environment Science and Data Center [24]. It was decided that the study area was mainly covered by eight

Landsat 8/OLI data set, which is provided by the Geospatial Data Cloud site, Computer

of the red band. The landuse data were provided by the Resource and Environment Sci‐

land, and grassland. The lithology of the Yangou Watershed was delineated using the China Geological Map at a 1:250,000‐scale. The distance to rivers was extracted using hy‐ drological analysis of DEM and neighborhood statistic in the ArcGIS. The distance to roads was developed using China Road Network data at a 1:100,000‐scale. The rainfall land-use types, namely water, wetland, farmland, resident land, bareland, forest, shrubland, and grassland. The lithology of the Yangou Watershed was delineated using the China Geological Map at a 1:250,000-scale. The distance to rivers was extracted using hydrological analysis of DEM and neighborhood statistic in the ArcGIS. The distance to roads was developed using China Road Network data at a 1:100,000-scale. The rainfall data were from the China meteorological background data set of a 500 m × 500 m pixel size, which is provided by the Resource and Environment Data Cloud Platform [25]. Then, these nine conditioning factors were selected in the present study and were standardized to the same size of 25 m × 25 m for further analyses. The study area covers 355 columns and 447'rows. The spatial analysis tool in the ArcGIS software was used to count the corresponding values of each factor at the landslide location. Also, the occurrence frequency ratio of the different classes was calculated in the land use and geology as the corresponding values of the landslide.

To show the relationship between each factor's subclasses and landslide volumes, the VR approach is proposed. It is a variant of the probabilistic method that is based on the observed relationships between the volume of the landslide and each landslide conditioning factor. The volume ratios for the classes or types of each conditioning factor were calculated by dividing the landslide volume by the typical value ratio of the study area. Each factor's ratio value was calculated using Equation (2):

$$VR\_{jk} = \frac{V\_{jk}^\*}{V\_{jk}} / \frac{T^\*}{T} \tag{2}$$

*VRjk* = volume ratio of each factor; *Vjk*\* = volume of observed landslide of class "*k*" of factor "*j*"; *Vjk* = volume of class "*k*" of factor "*j*"; *T* = typical value of the map, such as the total area of the map, the average distance to roads or rivers in each class, the average elevation of each class, and the average rainfall of each class. *T*\* = typical value of the observed landslide of the map, such as the total area of each class, the average distance of the landslide to roads or rivers, the average elevation of the landslide, and the average rainfall of the landslide. A value of 1 is the neutral value, and the values higher than 1 show a high positive correlation with a certain factor of the landslide [26]. A higher VR value shows that a higher probability of the landslide volume occurs.

The RIRA is a method of sensitivity analysis used by calculating the relationship between dependent variables and independent variables. RIRA can assess the relative contribution of each variable by placing greater emphasis on known variables contributing to landslides, and which are sorted by contribution [22]. In addition, each variable has a specific value in the method, which ensures the objectivity of the data calculation. To eliminate the dimensional impact between the independent variables, all data of the factors were standardized using the z-scores method before calculation. The dependent variables (such as the amount of landslide) and the independent variables (such as profile curvature, slope, altitude, NDVI, land-use, geology, distance to rivers, distance to roads, and rainfall) were listed separately. The variables were sorted in order from small to large, and the sensitivity coefficients of each influencing factor were calculated. In order to quantify the impact of landslide-related factors on the landslide, the calculation steps of RIRA were as follows: first, the landslide increase rate (%) was calculated using Equation (3):

$$RP\_i = \frac{2(P\_i - P\_{i-1})}{P\_i + P\_{i-1}} \tag{3}$$

where *P<sup>i</sup>* is the amount of the *i*-th landslide (104 m<sup>3</sup> ). Then, the increase rate (%) of each impact factor on the landslide was calculated as follows:

$$R t\_i = \frac{2(t\_i - t\_{i-1})}{t\_i + t\_{i-1}} \tag{4}$$

where *<sup>t</sup><sup>i</sup>* and *<sup>t</sup>i*−*<sup>1</sup>* are the values of the explanatory variable that are sequentially connected. *Rt<sup>i</sup>* is the increase rate calculated for two successively independent variables ordered by the values. Therefore, the Absolute Sensitivity Parameter (s), towards an independent variable, *t*, was calculated by its mean growth rate:

$$s\_{aj} = \left| \overline{(R P\_{i,j} / R t\_{i,j})}\_{i=1,N} \right| \tag{5}$$

Being *j* the *j*-th independent variable, *N* the total number of landslide data, namely slope, aspect, altitude, NDVI, land use, geology, distance to rivers and distance to roads. A factor with a larger absolute value of the s represents the higher sensitivity of the amount of the landslide to the related factor. In the study, the revised increase-rate-analysis (RIRA) method is used to assess the relationship between the factors and landslide.

At last, the landslide susceptibility of each subclass was calculated by Equation (6):

$$S\_{ak} = s\_{a\dot{j}} \times V R\_{\dot{j}k} \tag{6}$$

*Sak* is the susceptibility of each *k-th* class. Then, each class susceptibility map can be obtained using the *Sak* value in the ArcGIS software. At last, all thematic layers were combined using the *Sak* value to derive the landslide susceptibility map of the Yangou watershed. Calculated susceptibility values were classified into areas of very low, low, medium, high, and very high susceptibility using the natural break method available within the ArcGIS software.

## **4. Results and Discussion**

## *4.1. Landslide Inventory*

Landslides in the study area primarily lie along the northwest and southeast direction. There are two types of landslides in the region: slides and falls. About 15 observed landslides belong to slides and 12 falls. The amount of erosion is 704 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> and 14.4 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , accounting for 98% and 2% of the total landslide, respectively. The slide is the main type of landslide in the region, which had an important contribution to the frequency of landslide and the total amount of slope movements in the region. The triggering mechanism of falls is different from that of slides. When the water in the gully scours and transports the lower accumulation, and the gravity moment generated by the upper steep wall under the action of gravity is greater than the tensile force moment of the soil, the upper soil layer loses its balance to fall [27]. There were many small falls in the study area. Most of the fall-bodies fell and broke, which were difficult to preserve for a long time. Most of them have been transported to form sand production, which is difficult to investigate [28].

The volume of the smallest landslide is approximately 0.2 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , the largest is around 270 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , and the average is estimated to be 26.6 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> . Figure 3 showed the amount and frequency of landslides in each interval, which were tiny, small, medium, and large landslides. The total amount of erosion caused by large landslides of <sup>≥</sup><sup>100</sup> <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> was the largest, which was 503.6 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , accounting for 70.2% of the total erosion and followed by the erosion of the medium landslide, which was 174.7 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , accounting for 24.4% of total erosion. The small landslide was 37.1 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , accounting for 5.2% of total erosion, and the tiny landslide had the least landslide, at 2.7 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> , accounting for 0.2% of total erosion. In the case of frequency of landslides, there were 19 small-scale erosions with a volume of fewer than 10 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> in the study area, accounting for 70% of the total frequency of erosion, that is, the occurrence of tiny-scale and small-scale landslides were higher in the area. While the amount of medium-sized and large-scale collapse accounted for 94.6% of the total, their frequencies of occurrence were only 5 and 3 times, accounting for 19% and 11% of the total frequency, respectively. It can be seen that the number of tiny–small landslides in the study area was more frequent. The frequency of medium and large-scale landslides was little, but they had a major contribution to the total amount of landslides. It can also be seen from the statistics that the contribution of frequency

and the volume of landslides are not the same, and it is more meaningful to consider both effects in the study. The cumulative amount of landslide erosion in the study area reached 718.1 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> . A large amount of landslide erosion has an important impact on sediment yield and sediment transport in the basin. The landslide-triggered erosion enters the channel and then flows into the river system by runoff, directly or indirectly transporting a large amount of sediment to the channel, which is one of the main processes of soil erosion and river sediment sources. Large slides also cause road damage and soil erosion. The volume of the falls is not large but their migration speed is fast. Because the occurrence to stop movement does not exceed 1 min, people are often unable to hide and suffer the risk of casualties; the harm is no less than the landslide. Between 1980 and 2015, there were 53 landslides over the Loess Plateau, killing and missing 717 people [29]. *Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 7 of 17 enters the channel and then flows into the river system by runoff, directly or indirectly transporting a large amount of sediment to the channel, which is one of the main processes of soil erosion and river sediment sources. Large slides also cause road damage and soil erosion. The volume of the falls is not large but their migration speed is fast. Because the occurrence to stop movement does not exceed 1 min, people are often unable to hide and suffer the risk of casualties; the harm is no less than the landslide. Between 1980 and 2015, there were 53 landslides over the Loess Plateau, killing and missing 717 people [29].

**Figure 3.** Amount and frequency of landslides in the Yangou watershed. **Figure 3.** Amount and frequency of landslides in the Yangou watershed.

### *4.2. Application of the VR Model 4.2. Application of the VR Model*

Slope angle (°)

In order to study, in detail, the relationship between each factor and the scale of the landslide, each factor is divided into several classes. The results of the spatial relationship between landslide volumes and the classes of landslide conditioning factors using the VR model are shown in Table 1. In order to study, in detail, the relationship between each factor and the scale of the landslide, each factor is divided into several classes. The results of the spatial relationship between landslide volumes and the classes of landslide conditioning factors using the VR model are shown in Table 1.

In the case of slope angle, the results showed that the VR values also increased as the slope increased, indicating that the slope angle has a good positive correlation with the **Table 1.** Spatial relationship between landslide conditioning factors and landslide by volume ratio model.


0–15 0.7294 0.8629 0.8453 15–25 0.2427 0.1272 1.9089 25–35 0.0259 0.0094 2.7538 >35 0.0019 0.0006 3.5264


**Table 1.** *Cont.*

In the case of slope angle, the results showed that the VR values also increased as the slope increased, indicating that the slope angle has a good positive correlation with the probability of the landslide volume. The slope is a major factor in landslide occurrence because it relates to drainage, fault line, and road networks [2]. In the Yangou watershed, the 0◦–15◦ slope area constituted 20.34% of the total area, the 15◦–25◦ area accounted for 35.49%, the 25◦–35◦ area accounted for 36.57%, and slope areas of more than 35◦ were 7.61% of the total area. The findings revealed a wide range of slope changes across the whole research area. According to the survey, 78% of the landslides occurred on slopes greater than 15◦ in the study area. Generally, the mass movement shifts from sliding to slumping with an increasing angle [30,31]. Thus, landslides are created in a hilly region or mountainous areas with steep slopes. The slope affects the scale and intensity of surface material flow and energy conversion. For example, the slope affects the area and amount of rainfall on the slope, thereby affecting the magnitude of runoff, infiltration, and runoff kinetic energy [32].Thus, the greater the slope, the greater the gravity erosion [20,33].

There is no correlation between landslide volume and elevation rise in terms of elevation. There are many reasons for the influence of elevation on the development of landslide disasters, such as that different elevation ranges have different climatic characteristics and local depressions are due to slope differences, whether there is a sliding surface and an intensity of human activities in different elevation ranges [34]. Thus, elevation has no direct impact on the distribution of gravity erosion disaster points in the region. The 1200–1300 m area had the highest VR value of 0.5755, followed by 1100–1200 m (0.1896), and <1000 m (0.1450). The 1000–1100 m area has the lowest value of 0.959. In the Yangou watershed, 96.3% of the landslides occurred on the altitude ranging from 1000 m to 1300 m. This showed that landslides usually occurred at intermediate elevations. The results can also be

found in the other places since slopes tend be covered by a layer of thin colluvium that is prone to landslides [35,36].

The NDVI is a measure of surface reflectance to give a quantitative estimate of biomass and vegetation growth [37]. The value of NDVI in the study area ranges from −0.27 to 0.59 in the summer. Generally, the larger the NDVI value is, the higher the vegetation coverage is. In the case of NDVI, there is no landslide in the classes of no cover and high cover. In the other classes, the VR value reflected the probability of landslide volume increased with the increase of the plant cover. In other words, the vegetation didn't play important roles in the decrease of the landslide amount. There are positive and negative effects of vegetation cover on landslides. Vegetation roots can reduce the creep of the soil layer by increasing the shear strength of soil, which is negatively correlated with the amount of landslide erosion. On the other hand, in the field investigation of small watersheds, it was found that with the increase of vegetation coverage, the runoff on the slope increased, thereby increasing the gully erosion of the watershed [38]. Plant root splitting may also promote the occurrence of landslide erosion.

Land-use is another factor used to consider the natural and man-made environmental impacts on the land surface. The main land use types in the study area were forestland and cropland, accounting for 40% and 31% of the total study area, respectively. Secondly, the grassland accounted for 23% of the total study area. Other classes accounted for a combined 6% of the total study area. The bareland had the highest VR value of 7.4512, followed by resident land (3.5780) and grassland (1.7180). The class of forest and farmland had the VR values of 0.9734 and 0.0498, respectively. In the shrubland, wetland, and water, the VR values are zero. The results indicated surface attachment plays an important role in slope stability compared with bareland, and in our case, the cultivated areas are less prone to landslides when compared to other land cover types. This may be due to positive impacts of land management with better water and soil control measures on slope stability [39].

Different geological formations have different compositions and structure, which contribute to the strength of the material. Stronger rocks are more resistant to the driving force of landslides and less prone to damage than weaker rocks [40]. The main lithology units of the study area consisted of Jurassic sandstone and Quaternary loess. Vertical joints, structural joints, collapse joints, and weathering joints of loess are generally developed, which are potential geological factors leading to geological disasters such as loess slides and falls. The deformation and failure caused by loess collapsibility provide a channel for precipitation collection and rapid infiltration, and often lead to geological disasters such as landslides. The vertical joints and weathered joints of sandstone in the slope zone are developed, and the development and expansion of joint cracks form the dangerous rock mass of slope is seen, which often leads to fall disasters [25]. The joints in all the study area are very developed. From the VR values, it is seen that the high landslide volume occurred in the lithological units of sandstone with the VR value of 1.4640, followed by the lithological unit of loess with the VR value of 0.8152. The Jurassic sandstone had a higher value because it is mainly in the valley area; as the drainage area, the hydrological conditions are more complicated.

The stability of the soil mass might be harmed by road cutting caused by slope excavations. As a result, one of the influencing elements for landslides is closeness to a road. In the case of distance to roads, distance to roads of 0–50 m had the highest VR values, while the VR value in the distance to roads of 300–500 m had saltation. The result may be that the population density is high in the study area [41], and thus the road net is also density. In the range of 300–500 m, large collapses had occurred. However, the overall trend is that the VR value decreased as the distance increased, indicating that the distance to roads has a high influence on landslide occurrence and scale. Actually, the amount of landslide occurrence within 150 m of the road accounted for more than 78% in the study area.

There are reverse correlations between landslide volume and parameters based on the VR value of the distance to river. In the study area, 89% of the landslide occurred within 150 m of streams. At the same time, the hydrographic axes of the third and fourth order streams are considered to be the main factor of the occurrence of landslides [28]. This is because the scouring effect of river flow will erode the slope. At the same time, the seepage effect of water flow makes the slope soil saturated. A large amount of free water is accumulated in collapsible loess, and the shear strength of soil is reduced, which has a negative impact on the slope stability. In addition, human activities are mostly distributed in the valley, which reduces the slope stability and is also the cause of frequent disasters and serious landslide erosion on both sides of the river.

Profile curvature measures the rate of change of slope, which controls erosion and deposition by affecting the acceleration and deceleration of the flow across the surface. Convex surface is positive and concave surface is negative [42]. For the profile curvature, 44% of the landslide occurred on the concave slopes, while 56% occurred on the convex slopes, and the VR value for convex slope (1.0505) is greater than that of the concave slope (0.9567), indicating that the convex slope has a little more effect on the volume of landslides. The concave slope easily collects rainwater, and rainwater infiltrates along the slope, which reduces the shear strength of soil and causes the instability and deformation of soil. At the same time, the potential energy of convex slope soil is closer to the limit value of mechanical equilibrium, and is also vulnerable to stronger weathering, which may also lead to the occurrence of slope instability. Stochastic statistics of 300 unstable slopes in China have shown that the instability of convex slopes is better than that of concave slopes [43]. Besides, in the experiment, it is verified the larger volume of landslides easily occurred in the convex slope [44]. The failure characteristics of loess soil slope are creep, and soil creep tends to produce convex terrain [45]. Therefore, the convex slope is dominant in the steep loess slope, and it is also prone to landslides [46].

The precipitation intensity is an important factor to affect the occurrence of landslides in the semi-arid region, although the interval between minimum and maximum rainfall is small [40]. In the case of rainfall, the medium rainfall area had the highest VR value of 0.3799, followed by more rainfall area (0.2416), less rainfall area (0.2052), and heavy rainfall area (0.1755), which indicated that there is no cear relationship between the amount of this factor and landslide volume. On the Loess Plateau, long-term and strong early effective rainfall both triggered clusters of shallow landslides in the previous study [47].

## *4.3. Application of the RIRA Method*

The sensitivity coefficient of different related factors on landslides is quite different. The results of the sensitivity of the slope, elevation, NDVI, land-use, lithology, distance to road, distance to river, profile curvature, and rainfall on the volume of landslides performed by RIRA are reported in Table 2. The larger the sensitivity of an independent variable, the more prominently this variable influences the volumes of landslide. Profile curvature was found to be the most important influential factor, as it is the most sensitive parameter (0.3616), followed by the distance to river (0.2567) and distance to road (0.0908). Rainfall was the fourth in the ranking of the relative sensitivity, and the value is 0.0628. The relative sensitivity of elevation, slope, and NDVI are comparable to each other, although these variables are less important than the three variables motioned above (0.0523, 0.0345, and 0.0203, respectively). Finally, lithology and land use have relatively little relative sensitivity.

**Table 2.** Sensitivity analysis of landslides and related factors in the Yangou watershed (China).


In summary, the independent variables, profile curvature, distance to river, distance to road, rainfall, elevation, and slope, have an important impact on the triggering and scale of landslides. Especially, both the relative sensitivity of two variables, profile curvature, and the distance to river, are very high. Convex and concave slope profiles are common in natural hill slopes. In the Yangou watershed, it is the most import factor. That means

landslides in areas with large undulating terrain are subject to key monitoring. Therefore, the influence of slope morphology should be taken into account when humans are carrying out activities. In the case of the relationship between landslides and distance to streams, the VR value rises as the distance to a stream decreases. The river has a strong erosive effect on the exposed area of the bedrock, which affects the stability of the valley slope. Especially in smaller valleys, both the down erosion and side erosion of flowing water exist, and the valley slopes on both sides are steep, which are in the erosion of flowing water and are high-risk areas of slides and falls. In addition, the most direct impact of long-term surface water accumulation is the rise of the groundwater level, which forms groundwater in a large range and affects slope stability through groundwater. The groundwater on both sides of the basin passes can aggravate the slope instability by softening and eroding the earth and soil, reducing the strength of the rock and soil and simultaneously generating dynamic and hydrostatic pressures and pore water pressures, exerting floating force on the rock and soil and increasing the weight of the rock and soil [25]. Thus, the distance to river has the second sensitivity of all factors. In the study area, most of the collapse of the loess hillside slope, or the initiation of loess landslides, after having been excavated during the construction of an expressway, was induced by excavation and rainfall [48]. The sensitivity of distance to road on the scale of the landslide is of third importance. In the construction of rural road projects, rough construction, such as cut-off or slope cutting and filling, is blindly carried out without the stability analysis of the slopes on both sides of the valley, which makes the slope steeper, loses support, forms a steep slope, and causes hidden dangers for the occurrence of landslide disasters [49]. Precipitation is the main recharge source for soil moisture in the hill and gully region of the Loess Plateau [50]. According to the survey, 40% of the catastrophic landslides are trigged by rainfall on the Loess Plateau of China [26]. It is not easy to infiltrate after precipitation, however, once the effective strength of the loess is reduced by rainfall, the landslide is triggered. There is a close correspondence between landslide occurrence and the pattern of precipitation, and 84.6% of landslides have happened in the monsoon in the Baota district since 1985 [51]. Therefore, rainfall is the fourth important factor in the study area. Elevation data reflected the location of the landslide, which means the volume of landslides has higher relative sensitivity on the variable of elevation in the study area. The slope is an important performance of the slope geometry, which determines the state and distribution of stresses in the slope mass and controls the stability and mode of the instability [51]. However, the slope has relatively little impact on the occurrence and scale of landslides in the study area. The increased vegetation coverage has a certain inhibitory effect on the amount of landslides, but considering the comprehensive effect of other factors, the influence of plants is not obvious. In addition, landslides with different scales occurred in different lithology. Thus, the influence of lithology on landslides is also relatively rare in the study area. The impact of environmental and human factors on the surface, such as construction activities, road cutting, and natural resources exploitation, may lead to landslides [34]. However, during the 9th Five-Year Plan period (1996–2000), Yangou watershed was selected as a demonstration region for integrated research for the development of ecological agriculture in the Loess Plateau, and the land use structure changed rapidly [41]. Nowadays, the land cover tends to be stable, so the absolute sensitivity of landslides to land use is lowest.

## *4.4. Results of the Landslide Susceptibility Mapping*

Figures 4 and 5 showed the different space layers of conditioning factors and the landslide susceptibility map constructed by the VR model and RIRA method, respectively. In Figure 4, higher susceptibility values were in the layers with the convex slope and concave slope, as well as the distance of less than 50 m to the river, the distance of less than 50 m to the road, the slope of 25–35◦ , and the slope higher than 35◦ ; the values were all more than 10. Among the nine related factors examined in the study, elevation, distance to river, lithology, and rainfall are hard to change, and the other factors, slope, profile curve, distance to road, land-use, and NDVI, can be modified to reduce the landslide

susceptibility. Slope gradients can be reduced by converting slope land into terraces [21]. In terms of profile curve and distance to road, human activity and road construction should be carried out in a planned way. Especially, a preliminary reconnaissance survey should be implemented on the topography and geology of the area. In terms of land use, many landslides were observed on the bareland. Thus, the reduction of ground exposure may reduce landslide susceptibility. For the NDVI, it can be seen that soil bioengineering is not omnipotent and did not apply to all situations, e.g., deep-seated, rock, or extremely steep slopes [52]. Therefore, the comprehensive effect of the conditioning factors should be considered in the process of prevention and control of landslides. In the study area, reconstruction of the eco-environment had already been conducted since 1997. In the past decade, grasses or shrubs are the main types of plantings on the top of hills or upper part of hill slopes, while trees were mainly selected in the hilly slope with a moderate gradient, and most of the slope of cropland was converted to woodland or grassland [53]. Nevertheless, a quantitative assessment of the effectiveness of those comprehensive measures in inhabiting landslides is still lacking. *Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 13 of 17

**Figure 4.** Susceptibility value of spatial layers influencing landslide occurrence in the study area: (**a**) profile curvature, (**b**) distance to river, (**c**) distance to road, (**d**) rainfall, (**e**) elevation, (**f**) slope, (**g**) NDVI, (**h**) lithology, and (**i**) land use. **Figure 4.** Susceptibility value of spatial layers influencing landslide occurrence in the study area: (**a**) profile curvature, (**b**) distance to river, (**c**) distance to road, (**d**) rainfall, (**e**) elevation, (**f**) slope, (**g**) NDVI, (**h**) lithology, and (**i**) land use.

**Figure 5.** Landslide susceptibility map by the VR model and RIRA method: (**a**) landslide suscepti‐ bility classes in the Yangou watershed, (**b**) distribution of susceptibility over the landslide at sus‐ ceptibility map (%), and (**c**) histogram of landslide susceptibility index value. **Figure 5.** Landslide susceptibility map by the VR model and RIRA method: (**a**) landslide susceptibility classes in the Yangou watershed, (**b**) distribution of susceptibility over the landslide at susceptibility map (%), and (**c**) histogram of landslide susceptibility index value.

**5. Conclusions** To understand the role of landslide disasters in small watersheds more comprehen‐ sively, a landslide susceptibility assessment that considers the amount of landslides is necessary. The applied analysis of the VR model, RIRA method, and ArcGIS within the framework of the present study achieved more comprehension for disaster prevention and land planning with respect to the location and volume of the past landslides. The relationship among landslide occurrence, landslide volume, and nine landslide condition‐ ing factors such as slope, elevation, NDVI, land use, lithology, distance to road, distance to river, profile curvature, and rainfall are evaluated using the above models. Our analysis demonstrated that the sub‐classes in the same conditioning factors had different VR val‐ ues, as indicated that sub‐regions with homogenous properties contributed differently to the landslide scales. In addition, the results showed that the various parameters of land‐ slides had different sensitive values to the landslide volume and occurrence. Among those factors, profile curvature and distance to river were the two most highly sensitive, fol‐ lowed by distance to road, rainfall, elevation, slope, and NDVI. The lithology and land‐ use were the two lowest sensitive factors. The sensitive values are 0.3616, 0.2567, 0.0908, 0.0628, 0.0523, 0.0345, 0.0203, 0.0099, and 0.0097, respectively. Also, the landslide suscep‐ tibility map considering landslide volume classified the study area into five zones, with susceptibility degrees of very high, high, medium, low, and very low, which showed that about 60% of the landslides occurred in a high and very high susceptibility area, account‐ As shown in Figure 5c, a very low susceptibility zone occupies 6.90 km<sup>2</sup> . Similarly, low, medium, high, and very high susceptibility zones occupy 12.81, 12.83, 9.42, and 5.87 km<sup>2</sup> , respectively. As depicted in Figure 5a, the susceptibility levels decreased with distance from rivers and road. Very high susceptibility was in the region with a larger relief and along the river and road, which was also in line with the site investigation. In the whole province, landslides were mostly distributed on both sides of the roads and rivers. For example, collapses often take place on both sides of the Yan and the Fenchuan Rivers, especially along the roads and railways, [51]. This is due to the fact that Yan'an City and urban population centers are located in major river valleys and that intensive human engineering activities are often responsible for many landslides and related serious damages [51]. Therefore, in the road construction and other engineering works on the Loess Plateau, the engineers and researchers should be committed to the prevention and control of loess landslides [48], and, also, the road and the river located in the gullies. As a result, the down-cutting of channels during rainstorms, as well as the high content of soil moisture adjacent to a channel, make it more prone to landslides [21]. Moreover, in the regional susceptibility assessment of the Baota district, the specific area was extended along the Nanchuan River, and its tributaries to Shitougou village got a high susceptibility [28]. Figure 5b revealed that no landslide is located within the very low susceptibility areas; 15% of the landslide are situated within the low susceptibility areas; 26% within the moderate suitability areas; 30% within the high susceptibility areas; 30% within the very high susceptibility areas. The result of the landslide assessment highlighted that the majority of the occurred landslides

ing for 87% of the total volume of landslides. Landslides such as slides and falls are im‐

(60%), are situated within the zones of high and very high susceptibility, while the part of landslide volume was accounting for 87% of the total landslides. Landslide susceptibility is predicting the probability of landslides' occurrence regardless of their magnitude which is, however, a critical piece of information in mass movement erosion control. However, the previous study revealed small and medium mass movements often have much higher transport rate of debris than large ones [21]. Both the volume and frequency of landslides were considered in our study; the landslide mitigation and control areas can be further prioritized in the watershed from the map. The high and very high susceptibility area must receive higher priority in the watershed for large landslides, while the low susceptibility risk area can be excluded in the observation. Moreover, the area of low and medium susceptibility can be paid more attention to, as well as the soil erosion caused by landslides. The results of a zoning plan can also reduce the workload of field investigations and move the focus of the actual work in watershed management.

## **5. Conclusions**

To understand the role of landslide disasters in small watersheds more comprehensively, a landslide susceptibility assessment that considers the amount of landslides is necessary. The applied analysis of the VR model, RIRA method, and ArcGIS within the framework of the present study achieved more comprehension for disaster prevention and land planning with respect to the location and volume of the past landslides. The relationship among landslide occurrence, landslide volume, and nine landslide conditioning factors such as slope, elevation, NDVI, land use, lithology, distance to road, distance to river, profile curvature, and rainfall are evaluated using the above models. Our analysis demonstrated that the sub-classes in the same conditioning factors had different VR values, as indicated that sub-regions with homogenous properties contributed differently to the landslide scales. In addition, the results showed that the various parameters of landslides had different sensitive values to the landslide volume and occurrence. Among those factors, profile curvature and distance to river were the two most highly sensitive, followed by distance to road, rainfall, elevation, slope, and NDVI. The lithology and land-use were the two lowest sensitive factors. The sensitive values are 0.3616, 0.2567, 0.0908, 0.0628, 0.0523, 0.0345, 0.0203, 0.0099, and 0.0097, respectively. Also, the landslide susceptibility map considering landslide volume classified the study area into five zones, with susceptibility degrees of very high, high, medium, low, and very low, which showed that about 60% of the landslides occurred in a high and very high susceptibility area, accounting for 87% of the total volume of landslides. Landslides such as slides and falls are important contents of geological disaster prevention and control, and this phenomenon is also the main form of soil erosion. Small landslide erosion is more likely to cause sediment yield. Thus, the spatial distribution of the existing landslide susceptibility can be analyzed in order to identify the components that are located in large landslide scale areas and small landslide scale areas, so as to set the priorities to prevent and control landslide disasters and soil erosion. Engineers, decision-makers, and environmental managers may implement the analysis that we implemented in the present study during new or existing planning projects and produce maps that will make possible the adoption of policies and strategies aimed at multi-hazard mitigation.

**Author Contributions:** Writing—original draft preparation, H.G.; writing—review and editing, X.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded and supported by the National Natural Science Foundation of China (Grant No: 52009103, 42177346) and the Innovation Capability Support Program of Shaanxi (2019TD-040).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We are grateful to the anonymous reviewers for their valuable suggestions in improving the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Yuming Wu <sup>1</sup> , Aohua Tian 1,2,3 and Hengxing Lan 1,\***


**Abstract:** Numerical simulation is one of the methods to assess landslide movement processes, which is beneficial for engineering design and urban planning. With the development of computer technology, GIS has gradually become the mainstream platform for landslide simulation due to data availability and algorithm integrability. However, the dynamic processes of landslides are complicated, which makes integration difficult on GIS platforms. Some assumptions are applied to simplify these dynamic processes and solve this problem. Generally, there are two main types of numerical models on GIS platforms: models based on the Eulerian description and models based on the Lagrangian description. Case studies show that Eulerian models are suitable for flow-like movement, and Lagrangian models are suitable for discrete rigid bodies movement. Different models face different problems: the Eulerian-based models show numerical diffusion and oscillation, and the Lagrangian-based model needs to consider complicated shear and collision processes. In addition, the 3-D model can describe more details in the *z*-direction, while the depth-averaged model can obtain a reasonable range of motion, depth, and speed quickly. From the view of numerical simulation, inappropriate models, assumptions, and numerical schemes will produce errors. The landslide type refers to several forms of mass wasting associated with a wide range of ground movements, which guides establishing dynamic models and numerical schemes on GIS platforms and helps us obtain results accurately.

**Keywords:** landslide dynamics; landslide classification; dynamic models; depth-averaged model; GIS

## **1. Introduction**

Landslides cause a large number of casualties and property losses every year. Quantitative risk analysis (QRA) is essential for reducing property damage and loss of life. Numerical simulation is one of the methods for assessing the landslide movement process and potential hazard levels [1,2], which is part of quantitative risk analysis [3]. Many numerical models are used to obtain the landslide runout distance, thus providing a basis for urban planning [4–7].

From dynamics, all landslides follow Newton's second law [2] and, thus, mass conservation, momentum conservation, and energy conservation are the governing equations in numerical simulations. In terms of force, gravity is the primary driving force [8], and friction is the major resistance force during motion. However, different shear processes and collision processes make landslide dynamics complex.

There are two model types: discrete models and continuum models. The discrete models describe the interactions and shocks well, but the continuum models describe the flow well. The Lagrangian functions describe the translation and rotation of each discrete element for the discrete models. Eulerian functions are the best to calculate the field for the continuum models. With the development of computation, meshless methods such as

**Citation:** Wu, Y.; Tian, A.; Lan, H. Comparisons of Dynamic Landslide Models on GIS Platforms. *Appl. Sci.* **2022**, *12*, 3093. https://doi.org/ 10.3390/app12063093

Academic Editors: Ricardo Castedo, Miguel Llorente Isidro and David Moncoulon

Received: 18 February 2022 Accepted: 16 March 2022 Published: 17 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

SPH are proposed to solve problems dominated by complex boundary dynamics for the continuum models. In the SPH method, finding neighboring nodes is difficult. Models with different accuracies are used to describe the processes [7]. These models include Rockfall Analyst [9], Rocfall [10], and Rockyfor3D [11] for rockfalls, and Massflow [12], Flo-2D [13], r.avaflow [14], and Titan2D [15] for mass flow.

In some simulations, different software programs yield different trajectories based on methods, assumptions, and numerical schemes [16]. Models may face various challenges. For example, fluid mechanics models may encounter numerical diffusions and oscillations [17]; rigid body models need to consider complicated shear and collision processes [18]. Therefore, a suitable model is a prerequisite for obtaining accurate results, and the physical insights provided by models can help us understand the landslide process.

GIS provides some functions that allow users to build numerical models [19]. It makes programming convenient and easy. There are two types of landslide simulations on GIS platforms: add-in programs and stand-alone programs. Add-in programs use the GUI and functions on GIS platforms to build a dynamic landslide model. This combination is easy to implement, but the GIS platform may limit some capacities. Stand-alone programs use open libraries such as GDAL. They are computationally efficient but require a lot of development time. The landslide dynamics programs on GIS platforms need to be simplified to accommodate the GIS file structure. Therefore, appropriate assumptions and suitable models are the keys to simulation on GIS platforms.

Landslide classification uses simple words to describe the characteristics of landslide phenomena, which may guide us to select a suitable model for simulation on GIS platforms. Several authors, including Heim [20], Zaruba and Mencl [21], and Sharpe [22], proposed different landslide classification systems. In most of these classifications, the name of the landslide is a combination of the material and motion type [23–26]. The most wellknown classification system is Varnes' classification [27]. In this classification, material types include rock, debris, and earth, and motion types include falling, toppling, sliding, spreading flow, and complex. Subsequently, Hungr [23] updated the Varnes' classification to assess landslide phenomena, and the material types were more detailed than those in the original classification. Landslide types are often indicators of the critical movement process, and a suitable numerical model must consider the appropriate variables.

However, the relationship between the landslide movement process and the physical mechanism of landslides is not very clear. It results in difficulty in obtaining accurate runout zones on GIS platforms. In this paper, we compared numerical models on GIS platforms and established relationships between landslide types and the selection of numerical models on GIS platforms.

## **2. Models**

## *2.1. Descriptions*

Landslides contain two main motion processes: rigid body motion and flow-like motion. The rigid body motion assumes that the block does not deform or change shape [28]; the flow-like motion assumes that the material is continuous mass rather than discrete particles. The driving force is gravity, and the resistance force is mainly friction during the movement (Figure 1).

In dynamics, there are five major categories in the Varnes' movement type: falling, toppling, sliding, spreading, and flowing. The movement type is related to the critical processes that occur during motion. All the processes follow Newton's laws of motion [29,30], which are given by

$$\frac{dmv}{dt} = f\tag{1}$$

where *m* is the mass, *v* is the velocity, and *f* is the resultant force, which depends on the movement type. The left-hand side (LHS) of Equation (1) represents the motion process, and the right-hand side (RHS) is the resultant force.

**Figure 1.** Two main motion processes: (**A**) the rigid body motion, and (**B**) the flow-like motion.

In dynamics, there are five major categories in the Varnes' movement type: falling,

toppling, sliding, spreading, and flowing. The movement type is related to the critical

*Appl. Sci.* **2022**, *12*, x FOR PEER REVIEW 3 of 16

**Figure 1.** Two main motion processes: (**A**) the rigid body motion, and (**B**) the flow-like motion. **Figure 1.** Two main motion processes: (**A**) the rigid body motion, and (**B**) the flow-like motion. There are two common descriptions of motion in numerical models: the Eulerian description and the Lagrangian description [31] (Figure 2). In the Eulerian description, a

In dynamics, there are five major categories in the Varnes' movement type: falling, toppling, sliding, spreading, and flowing. The movement type is related to the critical processes that occur during motion. All the processes follow Newton's laws of motion [29,30], which are given by There are two common descriptions of motion in numerical models: the Eulerian description and the Lagrangian description [31] (Figure 2). In the Eulerian description, a model is not concerned about the location or velocity of any particular particle or columns, and the focus is on the grid. In the Lagrangian description, individual particles or columns are marked, and their positions and velocities are described as a function of time. model is not concerned about the location or velocity of any particular particle or columns, and the focus is on the grid. In the Lagrangian description, individual particles or columns are marked, and their positions and velocities are described as a function of time.
