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